Volumetric Lattice Boltzmann Models in General Curvilinear Coordinates: Theoretical Formulation

A theoretical formulation of lattice Boltzmann models on a general curvilinear coordinate system is presented. It is based on a volumetric representation so that mass and momentum are exactly conserved as in the conventional lattice Boltzmann on a Cartesian lattice. In contrast to some previously existing approaches for arbitrary meshes involving interpolation approximations among multiple neighboring cells, the current formulation preserves the fundamental one-to-one advection feature of a standard lattice Boltzmann method on a uniform Cartesian lattice. The new approach is built on the concept that a particle is moving along a curved path. A discrete space-time inertial force is derived so that the momentum conservation is exactly ensured for the underlying Euclidean space. We theoretically show that the new scheme recovers the Navier-Stokes equation in general curvilinear coordinates in the hydrodynamic limit, along with the correct mass continuity equation.

1 Introduction

Lattice Boltzmann Methods (LBM) have been developed as an advantageous method for computational fluid dynamics during past few decades [1, 2]. The underlying dynamics of a lattice Boltzmann system resides on the kinetic theory, in that it involves motion of many particles according to the Boltzmann equation [3, 4]. There are two fundamental dynamical processes in a basic Boltzmann kinetic system, namely advection and collision. In the advection process, a particle moves from one location to another according to its microscopic velocity, while in the collision process particles interact with each other obeying conservation laws and relax to equilibrium. In a standard LBM model, particle velocity takes on a discrete set of constant values, and the latter form exact links from one lattice site to its neighboring sites on a simple Bravais lattice corresponding to a three dimensional (3D) uniform cubical Cartesian mesh [5–7]. Therefore, having these discrete velocity values, a particle moves from one lattice site to a unique neighboring site per constant time interval (say $\text{Δ}t=1$). This exact one-to-one advection of a particle on the lattice effectively realizes a Lagrangian advection characteristics. That is, a particle location during the advection process is exactly predicted at any later time from its location at a given time. Such a Lagrangian advection characteristics is the reason for one of the key advantages of LBM, besides its simplicity, the advection process produces an extremely low numerical dissipation (due to no smearing of locations of a single particle). Indeed, if such a one-to-one advection is relaxed, numerical dissipation has been seen to significantly increase even with the usage of much higher order interpolation schemes. In addition, the Lagrangian one-to-one advection may also have non-trivial implications in simulations of flows at finite Knudsen numbers [7–9]. In contrast, unlike particle velocity, a fluid velocity varies in both time and space, hence an one-to-one advection process is impossible for a Navier-Stokes equation based numerical scheme on a spatially fixed mesh.

However, a uniform and cubical Cartesian mesh poses fundamental limitations. First of all, often in realistic fluid study one deals with a solid geometry with curved surfaces, and obviously a Cartesian mesh does not smoothly conform with such a geometric shape. Secondly, a flow usually has small scale structures concentrated in certain spatial locations and directions. For instance, in the turbulent boundary layer, the flow scale in the direction normal to the wall is much smaller than in the tangential direction or in the fluid region outside of the boundary layer. Consequently, the requirement on spatial resolution is significantly higher in the normal direction inside a boundary layer. A cubic Cartesian mesh does not provide the flexibility with different spatial resolutions in different directions. Therefore, it is of fundamental importance to extend the standard LBM for a general non-cubic Cartesian mesh.

There have been various attempts in the past to extend LBM on arbitrary meshes, including curvilinear mesh (cf [10–12]). The essential idea is to relax the aforementioned exact one-to-one advection between a pair of lattice sites while keeping the advection of particles according to the original constant discrete velocity values. As a consequence, on an arbitrary mesh, a particle after advecting from its original mesh site does not in general land on a single neighboring mesh site. Thus its location needs to be distributed onto a set of neighboring mesh sites via interpolation. As mentioned above, such an effective “one-to-many” advection process destroys the preciseness of a Lagrangian advection characteristics thus resulting in a significantly increased numerical dissipation.

It is of fundamental interest to formulate an LBM on a general mesh that preserves the exact one-to-one feature of the particle advection process. A possible way to accomplish this is to construct the micro-dynamic process on a non-Euclidean space represented by a general curvilinear coordinate system based on Riemann geometry [13]. In such a non-Euclidean space, a constant particle velocity corresponds to a curved and spatially varying path in the underlying Euclidean space. Indeed, the continuum Boltzmann kinetic theory on a Riemannian manifold can be theoretically described [14, 15]. The key difference between a constant particle velocity on a Euclidean and a non-Euclidean space is that the latter is accompanied with an inertial force, due to the 1st law of Newton in Euclidean space. According to this concept, it is entirely conceivable to formulate LBM models on a cubic Cartesian lattice in non-Euclidean space that corresponds to a general curvilinear mesh in the Euclidean space. The first such attempts were made by Mendosa et al. in a series of papers (see [16–18]). The significance of their work is that it established this fundamental approach to LBM on curvilinear meshes via a Riemann geometric framework.

As Mendosa et al., we formulate LBM on a general curvilinear coordinate system via the concept of Riemannian geometry. There are several critical differences between our approach presented here and that of Mendosa et al. First of all, we adopt a volumetric formulation so that the mass and momentum conservations are exactly ensured. The resulting continuity equation of mass is automatically shown to have the correct form in curvilinear coordinates, and there is no need to introduce any artificial mass source term to correct for any artifacts in the resulting hydrodynamics. Secondly, as in the continuum kinetic theory on a manifold, the only external source term in the extended LBM is associated with an effective inertial force due to curvilinearity. It contributes no mass source. Moreover, such a force term is constructed fully self-consistently within the discrete lattice formulation. This force form in discrete space and time recovers asymptotically the one in the continuum kinetic theory in the hydrodynamic limit. It does not rely on any outside analytical forms borrowed from the continuum kinetic theory [14, 15]. Thirdly, the inertial force enforces the exact momentum conservation for the underlying Euclidean space in the discrete space-time LBM model. Lastly, we demonstrate that the force term must be constructed properly so that it adds momentum in the system at a proper discrete time moment in order to produce the correct resulting Navier-Stokes hydrodynamics to the viscous order. Indeed, the correct Navier-Stokes equation is fully recovered in the hydrodynamic limit with the new formulation without introducing any extra correction terms.

To summarize, our goal is to theoretically demonstrate how discrete kinetic theories can be formulated in curvilinear coordinates in a way that mass and momentum are exactly conserved, and also the correct hydrodynamics is recovered in the microscopic limit. We hope the reader may find such discrete non-equilibrium statistical systems interesting, similar to discrete systems found in other branches of physics, e.g., spin models on a lattice in the equilibrium statistical theory, lattice field theories, etc. In addition to fundamental interest, non-equilibrium discrete systems like LBM may be useful for applications in numerical analysis and for practical use. Computational aspects such as scheme stability, numerical dissipation, performance, and boundary conditions treatment, as well as numerical verification of theoretical results, are outside of the scope of this study. We also restrict ourselves to one of the simplest situations, namely the low speed isothermal flow, which admits a closed self-contained analytical treatment which explicates the main features of volumetric formulation in curvilinear coordinates. Possible extensions and applications of this work are briefly discussed in Section 4 below.

The description of non-equilibrium dynamics and transport at the spatio-temporal scales not treatable by hydrodynamics can often benefit from using kinetic theory methods. Here we show how to construct lattice kinetic theories in a general coordinate system that possess exact conservation laws at both kinetic and continuum scales, as well as correct macroscopic limit behavior. In addition to the fundamental interest, this can be useful for describing interactions at interfaces between different fluids or materials with complex geometry.

The subsequent sections are organized as follows. We describe in Section 2 the formulation of the new LBM on a general curvilinear mesh. In Section 3, we provide a detailed theoretical derivation and show that the new LBM indeed produces the correct Navier-Stokes equation in general curvilinear coordinates. In order for the paper to be more self-contained, we provide, in Supplementary Appendices SA1–SA3 basic theoretical description of the continuum Boltzmann kinetic theory in curvilinear coordinates as in the literature [15], some fundamental properties [13], as well as a derivation of the Navier-Stokes hydrodynamic equation on a general curvilinear coordinate system.

2 Formulation of Lattice Boltzmann in Curvilinear Coordinates

This section is divided into two subsections. First, we construct the geometric quantities that are necessary for defining a general curvilinear coordinate system. Then we present the volumetric lattice Boltzmann formulation on a general curvilinear mesh.

2.1 Construction of Geometric Quantities on a Curvilinear Mesh

We describe how to construct a coordinate system for formulating kinetic theory on a curvilinear mesh. Let $\mathbf{x}$ be any spatial point in 3D Euclidean space. We choose a curvilinear coordinate system, so that $\mathbf{x}=\mathbf{x}\left(q\right)$ is uniquely defined by q. Here $q\equiv \left(q^1,q^2,q^3\right)$ are the coordinate values along three non-collinear congruences of parameterized basis curves. The above forms the basic starting point for a general curvilinear mesh. For when a 3D curvilinear mesh with its entire spatial layout is given, then all its vertex locations in space are known and specified. Let $\mathbf{x}$ now be only defined on sites (i.e., vertices) of the curvilinear mesh. We choose the coordinate system $\left\{q\right\}$ that is a one-to-one mapping between $\mathbf{x}$ and q. In other words, for any site $\mathbf{x}$ on the curvilinear mesh, there is a unique value q associated with it, that is $\mathbf{x}=\mathbf{x}\left(q\right)$ is uniquely defined. In addition, we construct the coordinate values of q on the mesh as follows: For any site $\mathbf{x}\left(q\right)$, its nearest neighbor site along the ith ($i=\mathrm{1,2,3}$) coordinate curve in the positive or negative direction is a spatial point ${\mathbf{x}}_{\pm i}$. Due to the unique mapping, we have ${\mathbf{x}}_{\pm i}=\mathbf{x}\left(q_{\pm i}\right)$, where $q_{\pm i}$ is a unique coordinate value for the neighboring site. It is entirely possible to choose the spatial variation of q in such a way that $q_{\pm i}$ and q only differ in their ith coordinate component values by a constant $d_0$. That is, $q_{\pm i}=\left(q_{\pm i}^1,q_{\pm i}^2,q_{\pm i}^3\right)$, and $q_{\pm i}^j-q^j=\pm d_0{\delta }_i^j$ ($i,j=\mathrm{1,2,3}$). Here we choose the constant $d_0$ to be unity without loss of generality.

It is easy to see that, under such a construction, the coordinate values $\left\{q\right\}$ actually form a simple uniform 3D cubic Cartesian lattice with the lattice spacing unity. One may interpret this “Cartesian” lattice $\left\{q\right\}$ as a result of deformation (bending, twisting and stretching/compressing) of the original curvilinear mesh $\left\{\mathbf{x}\right\}$ in the Euclidean space. Thus, its topological structure is the same as the original curvilinear mesh, but the Cartesian lattice is on a non-Euclidean space as a result of distortion of the original Euclidean space.

When a curvilinear mesh is provided, spatial locations of all the vertices $\left\{\mathbf{x}\right\}$ on the mesh are specified, and hence the distance from any one vertex to another on a such mesh is also fully determined. Let’s define the distance vector from $\mathbf{x}\left(q\right)$ to one of its neighbors $\mathbf{x}\left(q_{\pm i}\right)$ ($i=\mathrm{1,2,3}$) as

${\mathbf{D}}_{\pm i}\left(q\right)\equiv \mathbf{x}\left(q_{\pm i}\right)-\mathbf{x}\left(q\right);i=\mathrm{1,2,3}(1)$

For instance, ${\mathbf{D}}_{\pm 1}\left(q\right)\equiv \mathbf{x}\left(q^1\pm 1,q^2,q^3\right)-\mathbf{x}\left(q^1,q^2,q^3\right)$. Notice that, due to spatial non-uniformity of a general curvilinear mesh, the spatial distance from one mesh site to its nearest neighbor site is in general changing from location to location. In other words, ${\mathbf{D}}_{\pm i}\left(q\right)$ is a function of q. Furthermore, one can realize that the distance value in the positive direction along the ith coordinate curve is in general not equal to that in the negative direction. Explicitly, in terms of the distance vectors, ${\mathbf{D}}_i\left(q\right)\ne -{\mathbf{D}}_{-i}\left(q\right)$. For example, according to the definition of Eq. 1,

$\begin{array}{c}{\mathbf{D}}_{-1}\left(q\right)=\mathbf{x}\left(q^1-1,q^2,q^3\right)-\mathbf{x}\left(q^1,q^2,q^3\right)=-{\mathbf{D}}_1\left(q^1-1,q^2,q^3\right)\\ \ne -{\mathbf{D}}_1\left(q\right)=-(\mathbf{x}(q^1+1,q^2,q^3)-\mathbf{x}(q^1,q^2,q^3))\end{array}(2)$

and the inequality only becomes an equality everywhere if the curvilinear mesh is a uniform lattice (so that $\left|{\mathbf{D}}_i\right|=const$, independent of spatial coordinate value q).

Following the basic differential geometry concept, we now construct the basis tangent vectors at $\mathbf{x}\left(q\right)$ along each of the three coordinate directions,

${\mathbf{g}}_i\left(q\right)\equiv {\beta }_i\left(q\right)\left[{\mathbf{D}}_i\left(q\right)-{\mathbf{D}}_{-i}\left(q\right)\right]/2\text{Δ}x;i=\mathrm{1,2,3}(3)$

where the scalar factor ${\beta }_i\left(q\right)$ is

${\beta }_i\left(q\right)\equiv \frac{\left[\left|{\mathbf{D}}_i\left(q\right)\right|+\left|{\mathbf{D}}_{-i}\left(q\right)\right|\right]}{\left|{\mathbf{D}}_i\left(q\right)-{\mathbf{D}}_{-i}\left(q\right)\right|}$

and ${\beta }_i\left(q\right)=1$ when the two distance vectors ${\mathbf{D}}_i\left(q\right)$ and $-{\mathbf{D}}_{-i}\left(q\right)$ are parallel. With such a construction, the parity symmetry is achieved, so that

${\mathbf{g}}_i\left(q\right)=-{\mathbf{g}}_{-i}\left(q\right)$

Notice, in Eq. 3$\text{Δ}x$ is a length scale corresponding to a representative vertex spacing of the mesh. It is well known that, unlike a Cartesian coordinate system in Euclidean space, the basis tangent vectors ${\mathbf{g}}_i\left(q\right)$ ($i=\mathrm{1,2,3}$) of a curvilinear coordinate system are not orthonormal in general. That is, ${\mathbf{g}}_i\left(q\right)\cdot {\mathbf{g}}_j\left(q\right)\ne {\delta }_{ij}$. Hence, we should further define the corresponding metric tensor based on these basis tangent vectors,

$g_{ij}\left(q\right)\equiv {\mathbf{g}}_i\left(q\right)\cdot {\mathbf{g}}_j\left(q\right),i,j=\mathrm{1,2,3}(4)$

as well as the volume J of the cell centered at $\mathbf{x}\left(q\right)$,

$J\left(q\right)\equiv \left({\mathbf{g}}_1\left(q\right)\times {\mathbf{g}}_2\left(q\right)\right)\cdot {\mathbf{g}}_3\left(q\right)(5)$

and we can always choose a proper handedness so that $J\left(q\right)>0$. Clearly $J\left(q\right)$ is a constant in space for a uniform Bravais lattice. One can easily verify that

$g\left(q\right)\equiv det\left[g_{ij}\left(q\right)\right]=J^2\left(q\right)(6)$

with $det\left[g_{ij}\left(q\right)\right]$ being the determinant of the metric tensor $\left[g_{ij}\left(q\right)\right]$. Furthermore, we can define the co-tangent basis vectors ${\mathbf{g}}^i\left(q\right)$ ($i=\mathrm{1,2,3}$) accordingly,

$\begin{array}{l}{\mathbf{g}}^1\left(q\right)\equiv {\mathbf{g}}_2\left(q\right)\times {\mathbf{g}}_3\left(q\right)/J\left(q\right)\\ {\mathbf{g}}^2\left(q\right)\equiv {\mathbf{g}}_3\left(q\right)\times {\mathbf{g}}_1\left(q\right)/J\left(q\right)\\ {\mathbf{g}}^3\left(q\right)\equiv {\mathbf{g}}_1\left(q\right)\times {\mathbf{g}}_2\left(q\right)/J\left(q\right)\end{array}(7)$

Indeed, the basis tangent vectors and the co-tangent vectors are orthonormal to each other,

${\mathbf{g}}_i\left(q\right)\cdot {\mathbf{g}}^j\left(q\right)={\delta }_i^j,i,j=\mathrm{1,2,3}$

where ${\delta }_i^j$ is the Kronecker delta function. Similarly, we can define the inverse metric tensor,

$g^{ij}\left(q\right)\equiv {\mathbf{g}}^i\left(q\right)\cdot {\mathbf{g}}^j\left(q\right),i,j=\mathrm{1,2,3}(8)$

Clearly the inverse metric tensor is the inverse of the metric tensor, $\left[g^{ij}\left(q\right)\right]={\left[g_{ij}\left(q\right)\right]}^{-1}$, or

${\displaystyle \sum }_{k=1}^3{g}_{ik}\left(q\right)g^{kj}\left(q\right)={\delta }_i^j$

and

$det\left[g^{ij}\left(q\right)\right]=1/det\left[g_{ij}\left(q\right)\right]$

Having these basic geometric quantities defined above, we can now introduce the lattice Boltzmann velocity vectors on a general curvilinear mesh, similar to the ones on a standard Cartesian lattice,

${\mathbf{e}}_{\alpha }\left(q\right)\equiv c_{\alpha }^i{\mathbf{g}}_i\left(q\right)\frac{\text{Δ}x}{\text{Δ}t}(9)$

Here and thereafter, the summation convention is used for repeated Roman indices. For convenience, in subsequent derivations we adopt the ‘lattice units’ convention so that $\text{Δ}x$ = 1 and $\text{Δ}t=1$. The constant number $c_{\alpha }^i$ is either a positive or negative integer or zero, and it is the ith component value of the three dimensional coordinate array $c_{\alpha }\equiv \left(c_{\alpha }^1,c_{\alpha }^2,c_{\alpha }^3\right)$. For example, in the so called D3Q19 with the Greek index α running from 0 to 19,⁵

$c_{\alpha }\in \left\{\left(\mathrm{0,0,0}\right),\left(\pm \mathrm{1,0,0}\right),\left(0,\pm \mathrm{1,0}\right),\left(\mathrm{0,0},\pm 1\right),\left(\pm 1,\pm \mathrm{1,0}\right),\left(\pm \mathrm{1,0},\pm 1\right),\left(0,\pm 1,\pm 1\right)\right\}$

As shown in the subsequent sections, a set of necessary moment isotropy and normalization conditions must be satisfied in order to recover the correct full Navier-Stokes hydrodynamics [6, 7, 19, 20]. These are, when there exists a proper set of constant weights $\left\{w_{\alpha };\alpha =1,\dots ,b\right\}$, the set of lattice component vectors admit moment isotropy up to the 6th order, namely

$\begin{array}{l}{\displaystyle \sum _{\alpha }{w}_{\alpha }=1}\\ {\displaystyle \sum _{\alpha }{w}_{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j}=T_0{\delta }^{ij}\equiv T_0{\text{Δ}}^{\left(2\right),ij}\\ {\displaystyle \sum _{\alpha }{w}_{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j{c}_{\alpha }^k{c}_{\alpha }^l}=T_0^2\left[{\delta }^{ij}{\delta }^{kl}+{\delta }^{ik}{\delta }^{jl}+{\delta }^{il}{\delta }^{jk}\right]\equiv T_0^2{\text{Δ}}^{\left(4\right),ijkl}\\ {\displaystyle \sum _{\alpha }{w}_{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j{c}_{\alpha }^k{c}_{\alpha }^l{c}_{\alpha }^m{c}_{\alpha }^n}\\ =T_0^3\left[\begin{array}{l}{\delta }^{ij}{\text{Δ}}^{\left(4\right),klmn}+{\delta }^{ik}{\text{Δ}}^{\left(4\right),jlmn}+{\delta }^{il}{\text{Δ}}^{\left(4\right),jkmn}\\ +{\delta }^{im}{\text{Δ}}^{\left(4\right),jkln}+{\delta }^{in}{\text{Δ}}^{\left(4\right),jklm}\end{array}\right]\equiv T_0^3{\text{Δ}}^{\left(6\right),ijklmn}\end{array}(10)$

where the constant $T_0$ depends on the choice of a set of lattice vectors and ${\delta }^{ij}$ is the Kronecker delta function. For example, $T_0=1/3$ for the so called D3Q19 lattice (albeit it only satisfies Eq. 10 up to the 4th order). Notice, one should not confuse the three dimensional array $c_{\alpha }$ with the lattice vector ${\mathbf{e}}_{\alpha }\left(q\right)$ in Eq. 9. The former is simply an array of three constant (integer) numbers, while the latter is a vector (with a well defined direction and magnitude) at q on the curvilinear mesh.

Lastly, we define a set of specific geometric quantities below that are essential for the extended LBM model,

$\begin{array}{l}{\text{Θ}}_j^i\left(q+c_{\alpha },q\right)\equiv \left[{\mathbf{g}}_j\left(q+c_{\alpha }\right)-{\mathbf{g}}_j\left(q\right)\right]\cdot {\mathbf{g}}^i\left(q\right)\\ i,j=\mathrm{1,2,3};\alpha =\mathrm{0,1},\dots ,b\end{array}(11)$

It is easily seen that ${\text{Θ}}_j^i\left(q+c_{\alpha },q\right)$ vanish if the mesh is a uniform Cartesian lattice.

Once a curvilinear mesh is specified, all the geometric quantities above are fully determined and can thus be pre-computed before starting a dynamic LBM simulation.

2.2 Volumetric Lattice Boltzmann Model on a Curvilinear Mesh

Although the basic theoretical framework of the work is more general, for simplicity of describing the basic concept, we present the formulation for the so called isothermal LBM in this section. Similar to the standard lattice Boltzmann equation (LBE), we write the evolution of particle distribution below [1, 5, 21–23]

$N_{\alpha }\left(q+c_{\alpha },t+1\right)=N_{\alpha }\left(q,t\right)+{\text{Ω}}_{\alpha }\left(q,t\right)+\delta N_{\alpha }\left(q,t\right)(12)$

where $N_{\alpha }\left(q,t\right)$ is the number of particles belonging to the discrete direction $c_{\alpha }$ in the cell q at time t. We have assumed in Eq. 12 a unity time increment (i.e., $\text{Δ}t=1$) without loss of generality. ${\text{Ω}}_{\alpha }\left(q,t\right)$ in Eq. 12 is the collision term that satisfies local mass and momentum conservations,

$\begin{array}{l}{\displaystyle \sum _{\alpha }{\text{Ω}}_{\alpha }\left(q,t\right)}=0\\ {\displaystyle \sum _{\alpha }{\mathbf{e}}_{\alpha }\left(q\right){\text{Ω}}_{\alpha }\left(q,t\right)=0}\end{array}(13)$

The particle density distribution function $f_{\alpha }\left(q,t\right)$ is related to $N_{\alpha }\left(q,t\right)$, via

$J\left(q\right)f_{\alpha }\left(q,t\right)=N_{\alpha }\left(q,t\right)(14)$

where $J\left(q\right)$ is the volume of cell centered at q, as defined previously. The fundamental fluid quantities are given by the standard hydrodynamic moments,

$\begin{array}{l}\rho \left(q,t\right)={\displaystyle \sum _{\alpha }{f}_{\alpha }\left(q,t\right)}\\ \rho \left(q,t\right)\mathbf{u}\left(q,t\right)={\displaystyle \sum _{\alpha }{\mathbf{e}}_{\alpha }\left(q\right)f_{\alpha }\left(q,t\right)}\end{array}(15)$

where $\rho \left(q,t\right)$ and $\mathbf{u}\left(q,t\right)$ are fluid density and velocity at the location $\mathbf{x}\left(q\right)$ and time t. Using the relationship in Eq. 9, the velocity moment above can also be rewritten as

$\begin{array}{c}\rho \left(q,t\right)\mathbf{u}\left(q,t\right)={\displaystyle \sum _{\alpha }{c}_{\alpha }^i{\mathbf{g}}_i\left(q\right)f_{\alpha }\left(q,t\right)}\\ =\rho \left(q,t\right)U^i\left(q,t\right){\mathbf{g}}_i\left(q\right)\end{array}(16)$

and the velocity component value in the curvilinear coordinate system is given by,

$\rho \left(q,t\right)U^i\left(q,t\right)={\displaystyle \sum }_{\alpha }{c}_{\alpha }^i{f}_{\alpha }\left(q,t\right)(17)$

Or for simplicity of notation, we define a three-dimensional fluid velocity array $U\left(q,t\right)\equiv \left(U^1\left(q,t\right),U^2\left(q,t\right),U^3\left(q,t\right)\right)$. Therefore, Eq. 17 can be equivalently expressed as

$\rho \left(q,t\right)U\left(q,t\right)={\displaystyle \sum _{\alpha }{c}_{\alpha }{f}_{\alpha }\left(q,t\right)}(18)$

It is immediately seen that Eq. 18 has the same form for the fluid velocity as that in the standard Cartesian lattice based LBM. Similarly, the momentum conservation of the collision term in Eq. 13 can also be written as,

${\displaystyle \sum _{\alpha }{c}_{\alpha }{\text{Ω}}_{\alpha }\left(q,t\right)}=0(19)$

Often in LBM the collision term takes on a linearized form [1, 20], namely

${\text{Ω}}_{\alpha }\left(q,t\right)={\displaystyle \sum _{\beta }J\left(q\right)M_{\alpha \beta }\left[f_{\beta }\left(q,t\right)-f_{\beta }^{eq}\left(q,t\right)\right]}(20)$

where $M_{\alpha \beta }$ and $f_{\beta }^{eq}\left(q,t\right)$ represent a collision matrix and the equilibrium distribution function, respectively. In particular, the so called Bhatnagar-Gross-Krook (BGK) form corresponds to

$M_{\alpha \beta }=-\frac{1}{\tau }{\delta }_{\alpha \beta }$

with τ being the collision relaxation time [5, 22–24]. In order to recover the correct Navier-Stokes hydrodynamics, besides Eqs 13,19 the collision matrix needs to satisfy an additional condition [8, 9, 20, 25–27],

${\displaystyle \sum _{\alpha }{c}_{\alpha }{c}_{\alpha }{M}_{\alpha \beta }}=-\frac{1}{\tau }{c}_{\beta }{c}_{\beta }(21)$

Obviously the BGK form trivially satisfies such an additional property.

The extra term $\delta N_{\alpha }\left(q,t\right)$ in Eq. 12 represents the change of particle distribution due to an effective inertial force associated with the curvature and non-uniformity of a general curvilinear mesh. Obviously $\delta N_{\alpha }\left(q,t\right)$ vanishes in the standard LBM on a Cartesian lattice. We explain and construct its explicit form below.

We define the advection process as an exact one-to-one hop from one site to another as in the standard LBM, namely

$N_{\alpha }\left(q+c_{\alpha },t+1\right)=N_{\alpha }^{\prime }\left(q,t\right)(22)$

where $N_{\alpha }^{\prime }\left(q,t\right)$ is the post-collide distribution at $\left(q,t\right)$ that is equal to the right side of Eq. 12. Due to curvilinearity, though the amount of particles advected from cell q is exactly equal to what is arrived at cell $q+c_{\alpha }$ (as defined in Eq. 22), the momentum is in fact changed along the way. Indeed, in general

${\mathbf{e}}_{\alpha }\left(q+c_{\alpha }\right)N_{\alpha }\left(q+c_{\alpha },t+1\right)\ne {\mathbf{e}}_{\alpha }\left(q\right)N_{\alpha }^{\prime }\left(q,t\right)$

In the above, the left side of the inequality sign is the momentum value at the end of advection process while the right side is the value at the beginning of the process. The inequality is there because the path of particles is curved (as well as stretched or compressed), so that its velocity at the end of the advection is changed from its original value. This is fundamentally different from that on a uniform Cartesian lattice, in that the particles have a constant velocity throughout the advection process. Consequently, we have the following inequalities in the overall momentum values,

${\displaystyle \sum _{\alpha }{\mathbf{e}}_{\alpha }\left(q\right)N_{\alpha }\left(q,t\right)}\ne {\displaystyle \sum _{\alpha }{\mathbf{e}}_{\alpha }\left(q-e_{\alpha }\right)N_{\alpha }^{\prime }\left(q-c_{\alpha },t-1\right)}(23)$

where the right side of the unequal sign in Eq. 23 represents the total amount of momentum advected out of all the neighboring cells, while the left side is the total momentum arrived at cell q after the advection along the curved paths. Thus from Eqs 22,23, we see that the net momentum change via advection from all the neighboring cells into cell q is given by,

$J\left(q\right){\mathit{\chi }}^I\left(q,t\right)=-{\displaystyle \sum }_{\alpha }[{\mathbf{e}}_{\alpha }\left(q\right)-{\mathbf{e}}_{\alpha }\left(q-c_{\alpha }\right)]N_{\alpha }\left(q,t\right)(24)$

Similarly, we can realize that the net momentum change via advection out of cell q to all its neighboring cells is given by

$J\left(q\right){\mathit{\chi }}^o\left(q,t\right)=-{\displaystyle \sum _{\alpha }[{\mathbf{e}}_{\alpha }\left(q+c_{\alpha }\right)-{\mathbf{e}}_{\alpha }\left(q\right)]N_{\alpha }^{\prime }\left(q,t\right)}(25)$

Subsequently, if we impose the constraints on $\delta N_{\alpha }\left(q,t\right)$ below,

$\begin{array}{l}{\displaystyle \sum }_{\alpha }\delta N_{\alpha }\left(q,t\right)=0\\ {\displaystyle \sum }_{\alpha }{\mathbf{e}}_{\alpha }\left(q\right)\delta N_{\alpha }\left(q,t\right)=J\left(q\right)\mathit{\chi }\left(q,t\right)\end{array}(26)$

then we preserve an exact mass conservation as well as recover the exact momentum conservation in discrete space-time for the underlying Euclidean space. Here $\mathit{\chi }\left(q,t\right)=\left[{\mathit{\chi }}^I\left(q,t\right)+{\mathit{\chi }}^o\left(q,t\right)\right]/2$. More specifically, the first constraint in Eq. 26 means that no mass source is introduced by $\delta N_{\alpha }\left(q,t\right)$. On the other hand, the second constraint in Eq. 26 introduces an “inertial force” that equals exactly to the amount needed for achieving the momentum conservation in the underlying Euclidean space at any lattice site q and time t. The mechanism is in analogy with the continuum kinetic theoretic description in a curved space. Writing in the coordinate component form, we have

$F^i\left(q,t\right)=\mathit{\chi }\left(q,t\right)\cdot {\mathbf{g}}^i\left(q\right)$

Using the geometric quantities defined in the previous subsection, it can be immediately shown that

$J\left(q\right)F^i\left(q,t\right)=-\frac{1}{2}{\displaystyle \sum _{\alpha }{c}_{\alpha }^j\left\{{\text{Θ}}_j^i\left(q+c_{\alpha },q\right)N_{\alpha }^{\prime }\left(q,t\right)-{\text{Θ}}_j^i\left(q-c_{\alpha },q\right)N_{\alpha }\left(q,t\right)\right\}}(27)$

where the geometric function ${\text{Θ}}_j^i\left(q+c_{\alpha },q\right)$ is defined in Eq. 11. From Eqs 11, 27, one can immediately see that $F^i\left(q,t\right)$ vanishes if the curvilinear mesh is a regular uniform Cartesian lattice, as is obviously the case for a basic LBM. The second constraint in Eq. 26 can also be expressed in coordinate component form as

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }\left(q,t\right)}=J\left(q\right)F^i\left(q,t\right)(28)$

As to be demonstrated in the next section, in order to recover the full viscous Navier-Stokes equation, an additional constraint on the momentum flux also needs to be imposed below,

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }\left(q,t\right)}=J\left(q\right)\left[\delta {\text{Π}}^{ij}\left(q,t\right)+\delta {\text{Π}}^{ji}\left(q,t\right)\right](29)$

with

$\delta {\text{Π}}^{ij}\left(q,t\right)\equiv -\frac{1}{2}\left(1-\frac{1}{2\tau }\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^k\left[{\text{Θ}}_k^j\left(q+c_{\alpha },q\right)-{\text{Θ}}_k^j\left(q-c_{\alpha },q\right)\right]f_{\alpha }^{eq}\left(q,t\right)}(30)$

A specific form of $\delta N_{\alpha }\left(q,t\right)$ can be chosen below

$\delta N_{\alpha }\left(q,t\right)=w_{\alpha }J\left(q\right)[\frac{c_{\alpha }^j{F}^j\left(q,t\right)}{T_0}+\left(\frac{c_{\alpha }^j{c}_{\alpha }^k}{T_0}-{\delta }^{jk}\right)\frac{\delta {\text{Π}}^{jk}\left(q,t\right)}{T_0}](31)$

With a simple algebra, one can verify that it satisfies the moment constraints of Eqs 26, 28, 29. Notice, due to the appearance of $N_{\alpha }^{\prime }\left(q,t\right)$ in Eq. 27, the overall collision process for determining $N_{\alpha }^{\prime }\left(q,t\right)$ defines an implicit relationship, so that various ways to handle it need to be explored depending on types flows of interest.

The final part for complete specifying the extended LBM on a curvilinear mesh is the form of the equilibrium distribution function, which needs to be defined appropriately in order to recover the correct Euler equation as well as the Navier-Stokes equation in curvilinear coordinates in the hydrodynamic limit. In particular, the following fundamental conditions on hydrodynamic moments must be realized,

$\begin{array}{l}{\displaystyle \sum _{\alpha }{f}_{\alpha }^{eq}}=\rho ,\\ {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{f}_{\alpha }^{eq}}=\rho U^i,\\ {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j{f}_{\alpha }^{eq}}\equiv {\text{Π}}^{ij,eq}=g^{ij}\rho T_0+\rho {\tilde{U}}^i{\tilde{U}}^j,\\ {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j{c}_{\alpha }^k{f}_{\alpha }^{eq}}\equiv Q^{ijk,eq}=\left[g^{ij}{\tilde{U}}^k+g^{jk}{\tilde{U}}^i+g^{ki}{\tilde{U}}^j\right]\rho T_0+\rho {\tilde{U}}^i{\tilde{U}}^j{\tilde{U}}^k\end{array}(32)$

In the above,

${\tilde{U}}^i\left(q,t\right)=U^i\left(q,t\right)+\frac{1}{2}{a}^i\left(q,t\right)$

with $\rho \left(q,t\right)a^i\left(q,t\right)\equiv F^i\left(q,t\right)$. It is straightforward to show that these fundamental conditions are met by the following equilibrium distribution form,

$f_{\alpha }^{eq}=\rho w_{\alpha }\begin{array}{c}\{1+\frac{c_{\alpha }^i{U}^i}{T_0}+\frac{1}{2T_0}\left(\frac{c_{\alpha }^i{c}_{\alpha }^j}{T_0}-{\delta }^{ij}\right)\left[\left(g^{ij}-{\delta }^{ij}\right)T_0+{\tilde{U}}^i{\tilde{U}}^j\right]+\frac{1}{6T_0^3}\left(c_{\alpha }^i{c}_{\alpha }^j{c}_{\alpha }^k-T_0\left(c_{\alpha }^i{\delta }^{jk}+c_{\alpha }^j{\delta }^{ki}+c_{\alpha }^k{\delta }^{ij}\right)\right)\\ \left[T_0\left[\left(g^{ij}{\tilde{U}}^k-{\delta }^{ij}{U}^k\right)+\left(g^{jk}{\tilde{U}}^i-{\delta }^{jk}{U}^i\right)+\left(g^{ki}{\tilde{U}}^j-{\delta }^{ki}{U}^j\right)\right]+{\tilde{U}}^i{\tilde{U}}^j{\tilde{U}}^k\right]\}\end{array}(33)$

The equilibrium distribution form above is analogous to that of a low Mach number expansion of the Maxwell-Boltzmann distribution expressed in curvilinear coordinates [14]. It reduces to the standard LBM equilibrium distribution if the curvilinear mesh is a uniform Cartesian lattice so that $g^{ij}={\delta }^{ij}$ [7, 19, 20].

With all the quantities and dynamic properties defined above, we can show that the lattice Boltzmann Eq. 12 computed on the (non-Euclidean) uniform Cartesian lattice $\left\{q\right\}$ obeys the Navier-Stokes hydrodynamics in curvilinear coordinates. Consequently, mapping of the resulting fluid values, $\rho \left(q,t\right)$ and $U^i\left(q,t\right)$, onto the curvilinear mesh is given by a simple transformation below,

$\begin{array}{l}\rho \left(\mathbf{x}\left(q\right),t\right)=\rho \left(q,t\right),\\ \mathbf{u}\left(\mathbf{x}\left(q\right),t\right)=U^i\left(q,t\right){\mathbf{g}}_i\left(q\right)\end{array}(34)$

3 Derivation of the Navier-Stokes Hydrodynamics in Curvilinear Coordinates

In this section, we provide a detailed derivation and show that the extended LBM presented above indeed produces the correct Navier-Stokes equations in general curvilinear coordinates. We rewrite the lattice Boltzmann Eq. 12) below,

$N_{\alpha }\left(q+c_{\alpha },t+1\right)=N_{\alpha }\left(q,t\right)+{\text{Ω}}_{\alpha }\left(q,t\right)+\delta N_{\alpha }\left(q,t\right)(35)$

Expanding it in both time and space up to the second order leads to

$\left[{\partial }_t+c_{\alpha }^i\frac{\partial }{\partial q^i}+\frac{1}{2}{\left({\partial }_t+c_{\alpha }^i\frac{\partial }{\partial q^i}\right)}^2\right]N_{\alpha }={\text{Ω}}_{\alpha }+\delta N_{\alpha }(36)$

Now we introduce multiple scales in time and space based on the conventional Chapman-Enskog expansion procedure, [21, 28]

${\partial }_t=\mathit{ϵ}{\partial }_{t_0}+{\mathit{ϵ}}^2{\partial }_{t_1};\frac{\partial }{\partial q^i}=\mathit{ϵ}\frac{\partial }{\partial q^i}$

and

$N_{\alpha }=N_{\alpha }^{eq}+\mathit{ϵ}{N}_{\alpha }^{\left(1\right)}+{\mathit{ϵ}}^2{N}_{\alpha }^{\left(2\right)}+\cdots$

Here ϵ ($\ll 1$) denotes a small number. Conservation laws require,

$\begin{array}{l}{\displaystyle \sum _{\alpha }{N}_{\alpha }^{eq}}={\displaystyle \sum _{\alpha }{N}_{\alpha }}=J\rho ,{\displaystyle \sum _{\alpha }{c}_{\alpha }^i{N}_{\alpha }^{eq}}={\displaystyle \sum _{\alpha }{c}_{\alpha }^i{N}_{\alpha }}=J\rho U^i\\ {\displaystyle \sum _{\alpha }{N}_{\alpha }^{\left(n\right)}}={\displaystyle \sum _{\alpha }{c}_{\alpha }^i{N}_{\alpha }^{\left(n\right)}}=0;n>0\end{array}(37)$

Likewise, $\delta N_{\alpha }=\mathit{ϵ}\delta N_{\alpha }^{\left(0\right)}+{\mathit{ϵ}}^2\delta N_{\alpha }^{\left(1\right)}+\cdots$, and

${\displaystyle \sum _{\alpha }\delta N_{\alpha }^{\left(n\right)}}=0;n=\mathrm{0,1},\dots$

for they introduce no source of mass. Equating the same powers of ϵ, Eq. 36 leads to the following two equations

$\left({\partial }_{t_0}+c_{\alpha }^i\frac{\partial }{\partial q^i}\right)N_{\alpha }^{eq}={\displaystyle \sum _{\beta }{M}_{\alpha \beta }{N}_{\beta }^{\left(1\right)}+\delta N_{\alpha }^{\left(0\right)}}(38)$

and

$\begin{array}{c}{\partial }_{t_1}{N}_{\alpha }^{eq}+\frac{1}{2}{\left({\partial }_{t_0}+c_{\alpha }^i\frac{\partial }{\partial q^i}\right)}^2{N}_{\alpha }^{eq}+\left({\partial }_{t_0}+c_{\alpha }^i\frac{\partial }{\partial q^i}\right)N_{\alpha }^{\left(1\right)}\\ ={\displaystyle \sum _{\beta }{M}_{\alpha \beta }{N}_{\beta }^{\left(2\right)}+\delta N_{\alpha }^{\left(1\right)}}\end{array}(39)$

where the collision term has taken the linearized form in Eq. 20. Eq. 38 can be directly inverted,

$N_{\alpha }^{\left(1\right)}={\displaystyle \sum _{\beta }{M}_{\alpha \beta }^{-1}\left[\left({\partial }_{t_0}+c_{\beta }^i\frac{\partial }{\partial q^i}\right)N_{\beta }^{eq}-\delta N_{\beta }^{\left(0\right)}\right]}(40)$

Therefore, Eq. 39 can also be written as,

$\begin{array}{l}{\partial }_{t_1}{N}_{\alpha }^{eq}+\frac{1}{2}\left({\partial }_{t_0}+c_{\alpha }^i\frac{\partial }{\partial q^i}\right)\left[{\displaystyle \sum _{\beta }{M}_{\alpha \beta }{N}_{\beta }^{\left(1\right)}+\delta N_{\alpha }^{\left(0\right)}}\right]\\ +\left({\partial }_{t_0}+c_{\alpha }^i\frac{\partial }{\partial q^i}\right)N_{\alpha }^{\left(1\right)}\\ ={\displaystyle \sum _{\beta }{M}_{\alpha \beta }{N}_{\beta }^{\left(2\right)}+\delta N_{\alpha }^{\left(1\right)}}\end{array}(41)$

Taking the mass moment of Eqs 38, 41, and using mass conservation properties of $M_{\alpha \beta }$ and $\delta N^{\left(n\right)}$ above, we obtain for the leading order

${\displaystyle \sum _{\alpha }\left({\partial }_{t_0}+c_{\alpha }^i\frac{\partial }{\partial q^i}\right)N_{\alpha }^{eq}}=0(42)$

Based on the equilibrium moment definitions Eqs 37, 42 is equivalent to

${\partial }_{t_0}\left(J\rho \right)+\frac{\partial }{\partial q^i}\left(J\rho U^i\right)=0(43)$

For the next order, we have

${\partial }_{t_1}{\displaystyle \sum _{\alpha }{N}_{\alpha }^{eq}}+\frac{1}{2}\frac{\partial }{\partial q^i}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(0\right)}}=0(44)$

Combining Eq. 43 and Eq. 44, we get

${\partial }_t{\displaystyle \sum }_{\alpha }{N}_{\alpha }^{eq}+\frac{\partial }{\partial q^i}{\displaystyle \sum }_{\alpha }{c}_{\alpha }^i{N}_{\alpha }^{eq}+\frac{1}{2}\frac{\partial }{\partial q^i}{\displaystyle \sum }_{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(0\right)}=0(45)$

Substitute the definitions in Eqs 37, 45 becomes

${\partial }_t\left(J\rho \right)+\frac{\partial }{\partial q^i}\left[J\rho \left(U^i+\frac{a^{\left(0\right),i}}{2}\right)\right]=0(46)$

where

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(0\right)}}\equiv J\rho a^{\left(0\right),i}(47)$

is to be further discussed later in this section. Define the fluid velocity as

${\tilde{U}}^i\equiv U^i+\frac{1}{2}{a}^{\left(0\right),i}(48)$

and since $J=J\left(q\right)$ is not a function of time, Eq. 46 is in fact the standard mass continuity equation expressed in general curvilinear coordinates,

${\partial }_t\rho +\frac{1}{J}\frac{\partial }{\partial q^i}\left(J\rho {\tilde{U}}^i\right)=0(49)$

Recognizing that $\frac{1}{J}\frac{\partial }{\partial q^i}\left(J{A}^i\right)=\nabla \cdot \mathbf{A}$ is the divergence of a vector $\mathbf{A}$ in a general curvilinear coordinate system, Eq. 49 is simply the familiar mass continuity equation in the coordinate free representation,

${\partial }_t\rho +\nabla \cdot \left(\rho \tilde{u}\right)=0(50)$

Taking the momentum moment of Eq. 38, we get,

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i\left({\partial }_{t_0}+c_{\alpha }^j\frac{\partial }{\partial q^j}\right)N_{\alpha }^{eq}}={\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(0\right)}}(51)$

Define (see also in Eq. 32)

$J{\text{Π}}^{ij,eq}\equiv {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j{N}_{\alpha }^{eq}}(52)$

together with Eq. 47, Eq. 51 becomes

${\partial }_{t_0}\left(J\rho U^i\right)+\frac{\partial }{\partial q^j}\left(J{\text{Π}}^{ij,eq}\right)=J\rho a^{\left(0\right),i}(53)$

Similarly, taking the momentum moment of Eq. 41,

$\begin{array}{l}{\partial }_{t_1}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i{N}_{\alpha }^{eq}}+{\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\frac{\partial }{\partial q^j}}{\displaystyle \sum _{\alpha }\left[{\delta }_{\alpha \beta }+\frac{1}{2}{M}_{\alpha \beta }\right]N_{\beta }^{\left(1\right)}}\\ +\frac{1}{2}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i\left({\partial }_{t_0}+c_{\alpha }^j\frac{\partial }{\partial q^j}\right)\delta N_{\alpha }^{\left(0\right)}}\\ ={\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(1\right)}}\end{array}(54)$

Let

$J{\text{Π}}^{ij,\left(1\right)}\equiv {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j}{\displaystyle \sum _{\alpha }\left({\delta }_{\alpha \beta }+\frac{1}{2}{M}_{\alpha \beta }\right)N_{\beta }^{\left(1\right)}}(55)$

and use relation Eq. 40, we can derive

$\begin{array}{l}J{\text{Π}}^{ij,\left(1\right)}={\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j}{\displaystyle \sum _{\beta }\left({\delta }_{\alpha \beta }+\frac{1}{2}{M}_{\alpha \beta }\right)}{\displaystyle \sum _{\gamma }{M}_{\beta \gamma }^{-1}\left[\left({\partial }_{t_0}+c_{\gamma }^k\frac{\partial }{\partial q^k}\right)N_{\gamma }^{eq}-\delta N_{\gamma }^{\left(0\right)}\right]}\\ ={\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j}{\displaystyle \sum _{\beta }\left(M_{\alpha \beta }^{-1}+\frac{1}{2}{\delta }_{\alpha \beta }\right)\left[\left({\partial }_{t_0}+c_{\beta }^k\frac{\partial }{\partial q^k}\right)N_{\beta }^{eq}-\delta N_{\beta }^{\left(0\right)}\right]}\end{array}(56)$

Using the collision matrix property Eq. 21 and its inverse, namely

${\displaystyle \sum _{\alpha }{c}_{\alpha }{c}_{\alpha }{M}_{\alpha \beta }}=-\frac{1}{\tau }{c}_{\beta }{c}_{\beta };{\displaystyle \sum _{\alpha }{c}_{\alpha }{c}_{\alpha }{M}_{\alpha \beta }^{-1}}=-\tau c_{\beta }{c}_{\beta }(57)$

then Eq. 56 becomes

$J{\text{Π}}^{ij,\left(1\right)}=-\left(\tau -\frac{1}{2}\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\left[\left({\partial }_{t_0}+c_{\alpha }^k\frac{\partial }{\partial q^k}\right)N_{\alpha }^{eq}-\delta N_{\alpha }^{\left(0\right)}\right]}(58)$

Substitute $J{\text{Π}}^{ij,\left(1\right)}$ of the above into Eq. 54, we get

$\begin{array}{l}{\partial }_{t_1}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i{N}_{\alpha }^{eq}}-\frac{\partial }{\partial q^j}\left\{\left(\tau -\frac{1}{2}\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\left[\left({\partial }_{t_0}+c_{\alpha }^k\frac{\partial }{\partial q^k}\right)N_{\alpha }^{eq}-\delta N_{\alpha }^{\left(0\right)}\right]}\right\}\\ +\frac{1}{2}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i\left({\partial }_{t_0}+c_{\alpha }^j\frac{\partial }{\partial q^j}\right)\delta N_{\alpha }^{\left(0\right)}}\\ ={\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(1\right)}}\end{array}(59)$

After some simple cancellations and rearrangements, it becomes

$\begin{array}{l}{\partial }_{t_1}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i{N}_{\alpha }^{eq}}+\frac{1}{2}{\partial }_{t_0}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(0\right)}}\\ -\frac{\partial }{\partial q^j}\left[\left(\tau -\frac{1}{2}\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\left({\partial }_{t_0}+c_{\alpha }^k\frac{\partial }{\partial q^k}\right)N_{\alpha }^{eq}}\right]\\ +\frac{\partial }{\partial q^j}\left[\tau {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }^{\left(0\right)}}\right]\\ ={\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(1\right)}}\end{array}(60)$

Using definitions in Eq. 32, $N_{\alpha }^{eq}=J{f}_{\alpha }^{eq}$, and Eq. 47, Eq. 60 can be equivalently expressed as

$\begin{array}{l}{\partial }_{t_1}\left(J\rho U^i\right)+\frac{1}{2}{\partial }_{t_0}\left(J\rho a^{\left(0\right),i}\right)\\ -\frac{\partial }{\partial q^j}\left[\left(\tau -\frac{1}{2}\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\left({\partial }_{t_0}+c_{\alpha }^k\frac{\partial }{\partial q^k}\right)N_{\alpha }^{eq}}\right]\\ +\frac{\partial }{\partial q^j}\left[\tau {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }^{\left(0\right)}}\right]\\ =J\rho a^{\left(1\right),i}\end{array}(61)$

where in the above we have defined $J\rho a^{\left(1\right),i}\equiv {\displaystyle \sum }_{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(1\right)}$.

Now let us examine the forcing terms $J\rho a^{\left(0\right),i}$ and $J\rho a^{\left(1\right),i}$. Recall from the momentum conservation condition, the overall forcing term in Eq. 12 is given by Eq. 28, that is

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }\left(q,t\right)}=J\left(q\right)F^i\left(q,t\right)$

with $J{F}^i$ given by Eq. 27,

$J\left(q\right)F^i\left(q,t\right)=-\frac{1}{2}{\displaystyle \sum _{\alpha }{c}_{\alpha }^j\left\{{\text{Θ}}_j^i\left(q+c_{\alpha },q\right)N_{\alpha }^{\prime }\left(q,t\right)-{\text{Θ}}_j^i\left(q-c_{\alpha },q\right)N_{\alpha }\left(q,t\right)\right\}}(62)$

where ${\text{Θ}}_j^i\left(q+c_{\alpha },q\right)$ is defined in Eq. 11. Taking the continuum limit, we have for the leading order of ${\text{Θ}}_j^i\left(q+c_{\alpha },q\right)$

$\begin{array}{c}{\text{Θ}}_j^i\left(q+c_{\alpha },q\right)\equiv \left[{\mathbf{g}}_j\left(q+c_{\alpha }\right)-{\mathbf{g}}_j\left(q\right)\right]\cdot {\mathbf{g}}^i\left(q\right)\\ \approx c_{\alpha }^k\frac{\partial {\mathbf{g}}_j\left(q\right)}{\partial q^k}\cdot {\mathbf{g}}^i\left(q\right)\end{array}(63)$

Similarly,

$\begin{array}{c}{\text{Θ}}_j^i\left(q-c_{\alpha },q\right)\equiv \left[{\mathbf{g}}_j\left(q-c_{\alpha }\right)-{\mathbf{g}}_j\left(q\right)\right]\cdot {\mathbf{g}}^i\left(q\right)\\ \approx -c_{\alpha }^k\frac{\partial {\mathbf{g}}_j\left(q\right)}{\partial q^k}\cdot {\mathbf{g}}^i\left(q\right)\end{array}(64)$

Consequently, when substitute Eqs 63, 64 in Eq. 62, we obtain, to the leading order

$J\left(q\right)F^i\left(q,t\right)=-\frac{1}{2}{\text{Γ}}_{jk}^i\left(q\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k\left[N_{\alpha }^{\prime }\left(q,t\right)+N_{\alpha }\left(q,t\right)\right]}(65)$

where

${\text{Γ}}_{jk}^i\left(q\right)\equiv \frac{\partial {\mathbf{g}}_j\left(q\right)}{\partial q^k}\cdot {\mathbf{g}}^i\left(q\right)(66)$

is the Christoffel symbol as defined in differential geometry [13]. Also in the above, the relationship Eq. 22 is used, namely

$N_{\alpha }\left(q+c_{\alpha },t+1\right)=N_{\alpha }^{\prime }\left(q,t\right)(67)$

Furthermore, since

$N_{\alpha }\left(q+c_{\alpha },t+1\right)\approx N_{\alpha }\left(q,t\right)+\left({\partial }_{t_0}+c_{\alpha }^k\frac{\partial }{\partial q^k}\right)N_{\alpha }\left(q,t\right)$

we have from Eq. 65 the following

$\begin{array}{c}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }\left(q,t\right)}=J\left(q\right)F^i\left(q,t\right)\\ =-{\text{Γ}}_{jk}^i\left(q\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k{N}_{\alpha }\left(q,t\right)}\\ -\frac{1}{2}{\text{Γ}}_{jk}^i\left(q\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k\left({\partial }_{t_0}+c_{\alpha }^l\frac{\partial }{\partial q^l}\right)N_{\alpha }\left(q,t\right)}\end{array}(68)$

Comparing terms of the same order in ϵ, we recognize that

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(0\right)}}=J\rho a^{\left(0\right),i}=-{\text{Γ}}_{jk}^i{\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k{N}_{\alpha }^{eq}}(69)$

and

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(1\right)}}=J\rho a^{\left(1\right),i}=-{\text{Γ}}_{jk}^i{\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k{N}_{\alpha }^{\left(1\right)}}(70)$

$-\frac{1}{2}{\text{Γ}}_{jk}^i{\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k\left({\partial }_{t_0}+c_{\alpha }^l\frac{\partial }{\partial q^l}\right)N_{\alpha }^{eq}}(71)$

Replacing $N_{\alpha }^{\left(1\right)}$ with that in Eq. 40, with the use of property Eq. 57 and after some straightforward algebra, Eq. 71 becomes

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i\delta N_{\alpha }^{\left(1\right)}}=J\rho a^{\left(1\right),i}={\text{Γ}}_{jk}^i\left\{\left(\tau -\frac{1}{2}\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k\left({\partial }_{t_0}+c_{\alpha }^l\frac{\partial }{\partial q^l}\right)N_{\alpha }^{eq}-\tau {\displaystyle \sum }_{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k\delta N_{\alpha }^{\left(0\right)}}\right\}(72)$

From the definition Eqs 52, 69, we see that

$J\rho a^{\left(0\right),i}=-{\text{Γ}}_{jk}^{i}J{\text{Π}}^{jk,eq}(73)$

Similarly, let us also define

$\begin{array}{c}J{\text{Π}}^{ij,neq}\equiv -\left(\tau -\frac{1}{2}\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\left({\partial }_{t_0}+c_{\alpha }^k\frac{\partial }{\partial q^k}\right)N_{\alpha }^{eq}}\\ +\tau {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }^{\left(0\right)}}\\ =-\left(\tau -\frac{1}{2}\right)\left[{\partial }_{t_0}\left(J{\text{Π}}^{ij,eq}\right)+\frac{\partial }{\partial q^k}\left(J{Q}^{ijk,eq}\right)\right]\\ +\tau {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }^{\left(0\right)}}\end{array}(74)$

where

$J{Q}^{ijk,eq}\equiv {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j{c}_{\alpha }^k{N}_{\alpha }^{eq}}(75)$

Hence, from Eq. 72, we have

$J\rho a^{\left(1\right),i}=-{\text{Γ}}_{jk}^{i}J{\text{Π}}^{jk,neq}(76)$

Taking Eq. 73 to Eq. 53, we have for the leading order the equation below

${\partial }_{t_0}\left(J\rho U^i\right)+\frac{\partial }{\partial q^j}\left(J{\text{Π}}^{ij,eq}\right)+{\text{Γ}}_{jk}^{i}J{\text{Π}}^{jk,eq}=0(77)$

On the other hand, taking Eqs 74, 76 to Eq. 61, we get in the viscous order

${\partial }_{t_1}\left(J\rho U^i\right)+\frac{1}{2}{\partial }_{t_0}\left(J\rho a^{\left(0\right),i}\right)+\frac{\partial }{\partial q^j}\left(J{\text{Π}}^{ij,neq}\right)+{\text{Γ}}_{jk}^{i}J{\text{Π}}^{jk,neq}=0(78)$

Combining Eqs. 77, 78 together and using the definition in Eq. 48, one obtains formally the Cauchy’s transport equation, namely

${\partial }_t\left(\rho {\tilde{U}}^i\right)+\frac{1}{J}\frac{\partial }{\partial q^j}\left(J{\text{Π}}^{ij}\right)+{\text{Γ}}_{jk}^i{\text{Π}}^{jk}=0(79)$

where ${\text{Π}}^{jk}={\text{Π}}^{jk,eq}+{\text{Π}}^{jk,neq}$.

Eq. 79 can be further expressed in terms of hydrodynamic quantities. From the moment properties of the equilibrium distribution function in Eq. 32, we have

${\text{Π}}^{ij,eq}=g^{ij}\rho T_0+\rho {\tilde{U}}^i{\tilde{U}}^j(80)$

which has exactly the same form as that from the continuum kinetic theory. On the other hand, from Eq. 74 we have

$\begin{array}{c}{\text{Π}}^{ij,neq}=-\left(\tau -\frac{1}{2}\right)\left[{\partial }_{t_0}{\text{Π}}^{ij,eq}+\frac{1}{J}\frac{\partial }{\partial q^k}\left(J{Q}^{ijk,eq}\right)\right]\\ +\frac{\tau }{J}{\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }^{\left(0\right)}}\end{array}(81)$

where, from Eq. 32, ${\text{Π}}^{ij,eq}$ is given by Eq. 80 and $Q^{ijk,eq}$ is given below

$Q^{ijk,eq}=\left[g^{ij}{\tilde{U}}^k+g^{jk}{\tilde{U}}^i+g^{ki}{\tilde{U}}^j\right]\rho T_0+\rho {\tilde{U}}^i{\tilde{U}}^j{\tilde{U}}^k(82)$

which has exactly the same form as that from the continuum kinetic theory. However, the term $\tau {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }^{\left(0\right)}}$ needs to be further evaluated in terms of the hydrodynamic quantities.

From Eqs 29, 30, we have

${\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }}\left(q,t\right)=J\left(q\right)\left[\delta {\text{Π}}^{ij}\left(q,t\right)+\delta {\text{Π}}^{ji}\left(q,t\right)\right](83)$

with

$\delta {\text{Π}}^{ij}\left(q,t\right)\equiv -\frac{1}{2}\left(1-\frac{1}{2\tau }\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k\left[{\text{Θ}}_k^i\left(q+c_{\alpha },q\right)-{\text{Θ}}_k^i\left(q-c_{\alpha },q\right)\right]f_{\alpha }^{eq}\left(q,t\right)}(84)$

Taking the long wave length limit, Eq. 84 in the continuum limit becomes

$\begin{array}{c}\delta {\text{Π}}^{ij}\left(q,t\right)\approx -\left(1-\frac{1}{2\tau }\right){\text{Γ}}_{kl}^i\left(q\right){\displaystyle \sum _{\alpha }{c}_{\alpha }^j{c}_{\alpha }^k{c}_{\alpha }^l{f}_{\alpha }^{eq}\left(q,t\right)}\\ =-\left(1-\frac{1}{2\tau }\right){\text{Γ}}_{kl}^i\left(q\right)Q^{jkl,eq}\left(q,t\right)\end{array}(85)$

Therefore,

$\tau {\displaystyle \sum _{\alpha }{c}_{\alpha }^i{c}_{\alpha }^j\delta N_{\alpha }^{\left(0\right)}}=-\left(\tau -\frac{1}{2}\right)J\left(q\right)[{\text{Γ}}_{kl}^i\left(q\right)Q^{jkl,eq}\left(q,t\right)+{\text{Γ}}_{kl}^j\left(q\right)Q^{ikl,eq}\left(q,t\right)](86)$

So we finally obtain from Eq. 81 that

$\begin{array}{c}{\text{Π}}^{ij,neq}=-(\tau -\frac{1}{2})\{{\partial }_{t_0}{\text{Π}}^{ij,eq}+\frac{1}{J}\frac{\partial }{\partial q^k}\left(J{Q}^{ijk,eq}\right)\\ +[{\text{Γ}}_{kl}^i\left(q\right)Q^{jkl,eq}+{\text{Γ}}_{kl}^j\left(q\right)Q^{ikl,eq}]\}\end{array}(87)$

with ${\text{Π}}^{ij,eq}$ and $Q^{ijk,eq}$ given by Eqs 80, 82.

Substituting Eqs 80, 87 into Eq. 79, comparing with eqn (B.9) in Supplementary Appendix SA2, together with the definitions of Eqs 80, 82 (having the same forms as that from the continuum kinetic theory), we finally arrive at exactly the same form as derived out of a continuum kinetic theory in curvilinear coordinates. The only difference between the two is that the collision time τ is replaced by $\left(\tau -1/2\right)$ here. The rest of derivation for obtaining the Navier-Stokes hydrodynamics is a straightforward algebra and is no different from that of the continuum kinetic theory in curvilinear coordinates (provided in Supplementary Appendices SA1, SA2). Notice, same as the derivation in the continuum kinetic theory, in Eq. 87 the 0th order time derivative term ${\partial }_{t_0}{\text{Π}}^{ij,eq}$ is to be replaced by the leading (Euler) order equations on the right hand sides of Eqs 43, 77. Without repeating the rest of algebra, here we present the final form - the Navier-Stokes equation in a general curvilinear coordinate system,

$\begin{array}{c}{\partial }_t\left(\rho {\tilde{U}}^i\right)+\frac{1}{J}\frac{\partial }{\partial q^j}\left(J\rho {\tilde{U}}^i{\tilde{U}}^j\right)+{\text{Γ}}_{jk}^i\rho {\tilde{U}}^j{\tilde{U}}^k\\ =-g^{ij}\frac{\partial p}{\partial q^j}+\frac{1}{J}\frac{\partial }{\partial q^j}\left(2J\mu S^{ij}\right)+2{\text{Γ}}_{jk}^i\mu S^{jk}\end{array}(88)$

where the pressure $p=\rho T_0$, $\mu =\left(\tau -1/2\right)\rho T_0$, and

$S^{ij}=\frac{1}{2}\left[{\tilde{U}}^i{|}_k{g}^{kj}+{\tilde{U}}^j{|}_k{g}^{ki}\right]$

with the standard covariant derivative of the velocity component defined as

${\tilde{U}}^i{|}_k\equiv \frac{\partial {\tilde{U}}^i}{\partial q^k}+{\text{Γ}}_{kl}^i{\tilde{U}}^l$

Using the standard definitions of differential operators in differential geometry [13], we can recognize that Eq. 88 is indeed the familiar Navier-Stokes equation in a generic coordinate-free operator representation,

${\partial }_t\left(\rho \tilde{u}\right)+\nabla \cdot \left(\rho \tilde{u}\tilde{u}\right)=-\nabla p+\nabla \cdot \left(2\mu \mathbf{S}\right)(89)$

4 Discussion

In this paper, we present a theoretical formulation of lattice Boltzmann models in a general curvilinear coordinate system. The formulation is an application of the Riemannian geometry for kinetic theory [14] to discrete space and time. Unlike some previous works [16–18], here we use a volumetric representation so that conservation laws are exactly ensured [11]. Furthermore, in the current formulation, we find that the main and the only additional source term in the extended LBM model is corresponding to the inertial force to ensure the exact momentum conservation in the underlying Euclidean space. This is the same as in the continuum kinetic theory. On the other hand, this forcing term needs to be applied at an appropriate discrete time in order to realize the correct viscous fluid effect associated with non-equilibrium physics. The equilibrium distribution function also needs to be properly modified that is directly analogous to the Maxwell-Boltzmann distribution on a curved space. Unlike the previous formulations [16–18], there are no other terms or treatments added to cancel any discrete artifacts. Through a detailed analysis, we have shown that the current LBM formulation recovers the correct Navier-Stokes behavior in the hydrodynamic limit, as long as a discrete lattice velocity set satisfying a sufficient order of isotropy is used. Extensive numerical validations of this extended LBM for various flows on various curvilinear meshes are to be presented in future publications. The main benefit of this kind of theoretical formulation in a general curvilinear coordinate system is its preservation of the key LBM one-to-one Lagrangian nature of particle advection. Not only this is desirable in algorithmic simplicity, as a standard LBM on a Cartesian lattice, it also has non-trivial implications for flows at finite Knudsen number [7–9]. Although the specific LBM model constructed here is for an isothermal fluid, the fundamental framework is directly extendable to more general fluid flow situations. Possible extensions in the future may include transport of scalars, complex fluids, finite Knudsen flows, as well as higher speed flows with energy and temperature dynamics in curved space. The latter is essential for study of highly compressible flows and flows with substantial temperature variations. Another interesting possible extension of the current theoretical formulation in the future is for a time varying coordinate system. This is useful particularly for studying of fluid flows around a dynamically deforming solid object.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Conflict of Interest

Author HC was employed by the company Dassault Systemes.

Acknowledgments

I am grateful to Raoyang Zhang, Pradeep Gopalakrishnan, Ilya Staroselsky and Alexei Chekhlov for insightful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2021.691582/full#supplementary-material

References

1. Benzi, R, Succi, S, and Vergassola, M. The Lattice Boltzmann Equation: Theory and Applications. Phys Rep (1992). 222:145–97. doi:10.1016/0370-1573(92)90090-m

CrossRef Full Text | Google Scholar

2. Chen, S, and Doolen, GD. Lattice Boltzmann Method for Fluid Flows. Annu Rev Fluid Mech (1998). 30:329–64. doi:10.1146/annurev.fluid.30.1.329

CrossRef Full Text | Google Scholar

3. Huang, K. Statistical Mechanics. 2nd ed. India: John Wiley (1987).

4. Cercignani, C. The Boltzmann Equation and its Applications. New York: Springer-Verlag (1975).

5. Qian, YH, D'Humières, D, and Lallemand, P. Lattice BGK Models for Navier-Stokes Equation. Europhys Lett (1992). 17:479–84. doi:10.1209/0295-5075/17/6/001

CrossRef Full Text | Google Scholar

6. Chen, H, Goldhirsch, I, and Orszag, SA. Discrete Rotational Symmetry, Moment Isotropy, and Higher Order Lattice Boltzmann Models. J Sci Comput (2008). 34:87–112. doi:10.1007/s10915-007-9159-3

CrossRef Full Text | Google Scholar

7. Shan, X, Yuan, X-F, and Chen, H. Kinetic Theory Representation of Hydrodynamics: a Way beyond the Navier-Stokes Equation. J Fluid Mech (2006). 550:413. doi:10.1017/s0022112005008153

CrossRef Full Text | Google Scholar

8. Chen, H, Zhang, R, Staroselsky, I, and Jhon, M. Recovery of Full Rotational Invariance in Lattice Boltzmann Formulations for High Knudsen Number Flows. Physica A: Stat Mech its Appl (2006). 362:125–31. doi:10.1016/j.physa.2005.09.008

CrossRef Full Text | Google Scholar

9. Zhang, R, Chen, H, and Shan, X. Efficient Kinetic Method for Fluid Simulation beyond the Navier-Stokes Equation. Phys Rew. E (2006). 74:046703. doi:10.1103/physreve.74.046703

PubMed Abstract | CrossRef Full Text | Google Scholar

10. He, X, Luo, L-S, and Dembo, M. Some Progress in Lattice Boltzmann Method Part I Nonuniform Mesh Grids. J Comput Phys (1996). 129:357–63. doi:10.1006/jcph.1996.0255

CrossRef Full Text | Google Scholar

11. Chen, H. Volumetric Formulation of the Lattice Boltzmann Method for Fluid Dynamics: Basic Concept. Phys Rev E (1998). 58:3955–63. doi:10.1103/physreve.58.3955

CrossRef Full Text | Google Scholar

12. Barraza, JAR, and Dieterding, R. Towards a Generalised Lattice-Boltzmann Method for Aerodynamic Simulations. J Comput Phys (2019). 45:101182. doi:10.1016/j.jocs.2020.101182

Google Scholar

13. Aris, R. Vectors, Tensors, and the Basic Equations of Fluid Mechanics. New York: Dover (1962).

14. Vedenyapin, V, Sinitsyn, A, and Dulov, E. Kinetic Boltzmann, Vlasov and Related Equations. London: Elsevier (2011).

15. Love, PJ, and Cianci, D. From the Boltzmann Equation to Fluid Mechanics on a Manifold. Phil Trans R Soc A (2011). 369:2362–70. doi:10.1098/rsta.2011.0097

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Mendoza, M, Debus, J-D, Succi, S, and Herrmann, HJ. Lattice Kinetic Scheme for Generalized Coordinates and Curved Spaces. Int J Mod Phys C (2014). 25:1441001. doi:10.1142/s0129183114410010

CrossRef Full Text | Google Scholar

17. Debus, J-D, Mendoza, M, and Herrmann, HJ. Dean Instability in Double-Curved Channels. Phys Rev E (2014). 90:053308. doi:10.1103/physreve.90.053308

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Debus, J-D, Mendoza, M, Succi, S, and Herrmann, HJ. Poiseuille Flow in Curved Spaces. Phys Rev E (2016). 93:043316. doi:10.1103/physreve.93.043316

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Chen, H, and Shan, X, Fundamental Conditions for N-Th-Order Accurate Lattice Boltzmann Models, Physica D 237, issues: 14-17, 2003–8. (2007). doi:10.1016/j.physd.2007.11.010

Google Scholar

20. Chen, H, Teixeira, C, and Molvig, K. Digital Physics Approach to Computational Fluid Dynamics: Some Basic Theoretical Features. Intl J Mod Phys C (1997). 9(8):675. doi:10.1142/s0129183197000576

CrossRef Full Text | Google Scholar

21. Frisch, U, D’Humieres, D, Hasslacher, B, Lallemand, P, Pomeau, Y, and Rivet, JP. Lattice Gas Hydrodynamics in Two and Three Dimensions. Complex Syst (1987). 1:649–707.

Google Scholar

22. Chen, S, Chen, H, Martnez, D, and Matthaeus, W. Lattice Boltzmann Model for Simulation of Magnetohydrodynamics. Phys Rev Lett (1991). 67:3776–9. doi:10.1103/physrevlett.67.3776

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Chen, H, Chen, S, and Matthaeus, W. Recovery of the Navier-Stokes Equations Using a Lattice-Gas Boltzmann Method. Phys.Rev A (1992). 45:5339–42. doi:10.1103/physreva.45.r5339

CrossRef Full Text | Google Scholar

24. Bhatnagar, PL, Gross, EP, and Krook, M. A Model for Collision Processes in Gases fsI Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys Rev (1954). 94(3):511–25. doi:10.1103/physrev.94.511

CrossRef Full Text | Google Scholar

25. Latt, J, and Chopard, B. Lattice Boltzmann Method with Regularized Pre-collision Distribution Functions. Math Comput Simulat (2006). 72(2-6):165. doi:10.1016/j.matcom.2006.05.017

CrossRef Full Text | Google Scholar

26. Chen, H, Gopalakrishnan, P, and Zhang, R. Recovery of Galilean Invariance in thermal Lattice Boltzmann Models for Arbitrary Prandtl Number. Int J Mod Phys C (2014). 25(10):1450046. doi:10.1142/s0129183114500466

CrossRef Full Text | Google Scholar

27. Chen, H, Zhang, R, and Gopalakrishnan, P. Filtered Lattice Boltzmann Collision Formulation Enforcing Isotropy and Galilean Invariance. Phys Scr (2020). 95:034003. doi:10.1088/1402-4896/ab4b4d

CrossRef Full Text | Google Scholar

28. Chapman, S, and Cowling, T. The Mathematical Theory of Non-uniform Gases. 3rd ed. Cambridge: Cambridge Univ. Press (1990).