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.


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,6,7] Therefore, having these discrete velocity values, a particle moves from one lattice site to a unique neighboring site per constant time interval (say ∆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-toone advection may also have non-trivial implications in simulations of flows at finite Knudsen numbers.[7,23,24] 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.[8,9,10]).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.[11] 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 Rieman-nian manifold can be theoretically described.[12,13] 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, [14,15,16]).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.[12,13] 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 spatiotemporal 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 Appendices A, B and C, a basic theoretical description of the continuum Boltzmann kinetic theory in curvilinear coordinates as in the literature, [13] some fundamental properties, [11] as well as a derivation of the Navier-Stokes hydrodynamic equation on a general curvilinear coordinate system.

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.

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 x be any spatial point in 3D Euclidean space.We choose a curvilinear coordinate system, so that x = x(q) is uniquely defined by q.Here q ≡ (q 1 , q 2 , q 3 ) 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 x now be only defined on sites (i.e., vertices) of the curvilinear mesh.We choose the coordinate system {q} that is a one-to-one mapping between x and q.In other words, for any site x on the curvilinear mesh, there is a unique value q associated with it, that is x = x(q) is uniquely defined.In addition, we construct the coordinate values of q on the mesh as follows: For any site x(q), its nearest neighbor site along the ith (i = 1, 2, 3) coordinate curve in the positive or negative direction is a spatial point x ±i .Due to the unique mapping, we have x ±i = x(q ±i ), where q ±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 ±i and q only differ in their ith coordinate component values by a constant d 0 .That is, q ±i = (q 1 ±i , q 2 ±i , q 3 ±i ), and q j ±i − q j = ±d 0 δ j i (i, j = 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 {q} actually form a simple uniform 3D cubic Cartesian lattice with the lattice spacing unity.One may interpret this "Cartesian" lattice {q} as a result of deformation (bending, twisting and stretching/compressing) of the original curvilinear mesh {x} 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 {x} on the mesh are specified, and hence the distance from any one vertex to another on a such mesh is also fully determined.Let us define the distance vector from x(q) to one of its neighbors x(q ±i ) (i = 1, 2, 3) as For instance, D ±1 (q) ≡ x(q 1 ± 1, q 2 , q 3 ) − x(q 1 , q 2 , q 3 ).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, D ±i (q) 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, D i (q) = −D −i (q).For example, according to the definition of (1), D −1 (q) = x(q 1 − 1, q 2 , q 3 ) − x(q 1 , q 2 , q 3 ) = −D 1 (q 1 − 1, q 2 , q 3 ) = −D 1 (q) = −(x(q 1 + 1, q 2 , q 3 ) − x(q 1 , q 2 , q 3 )) and the inequality only becomes an equality everywhere if the curvilinear mesh is a uniform lattice (so that |D i | = const, independent of spatial coordinate value q).Following the basic differential geometry concept, we now construct the basis tangent vectors at x(q) along each of the three coordinate directions, where the scalar factor β i (q) is and β i (q) = 1 when the two distance vectors D i (q) and −D −i (q) are parallel.With such a construction, the parity symmetry is achieved, so that Notice, in (3) ∆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 g i (q) (i = 1, 2, 3) of a curvilinear coordinate system are not orthonormal in general.That is, g i (q) • g j (q) = δ ij .Hence, we should further define the corresponding metric tensor based on these basis tangent vectors, as well as the volume J of the cell centered at x(q), and we can always choose a proper handedness so that J(q) > 0. Clearly J(q) is a constant in space for a uniform Bravais lattice.One can easily verify that with det[g ij (q)] being the determinant of the metric tensor [g ij (q)].Furthermore, we can define the co-tangent basis vectors g i (q) (i = 1, 2, 3) accordingly, Indeed, the basis tangent vectors and the co-tangent vectors are orthonormal to each other, g i (q) • g j (q) = δ j i , i, j = 1, 2, 3 where δ j i is the Kronecker delta function.Similarly, we can define the inverse metric tensor, Clearly the inverse metric tensor is the inverse of the metric tensor, [ 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, 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 ∆x = 1 and ∆t = 1.The constant number c i α is either a positive or negative integer or zero, and it is the ith component value of the three dimensional coordinate array c α ≡ (c 1 α , c 2 α , c 3 α ).For example, in the so called D3Q19 with the Greek index α running from 0 to 19, [5] c α ∈ {(0, 0, 0), (±1, 0, 0), (0, ±1, 0), (0, 0, ±1), (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1)} 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.[17,7,6,18] These are, when there exists a proper set of constant weights {w α ; α = 1, . . ., b}, the set of lattice component vectors admit moment isotropy up to the 6th order, namely where the constant T 0 depends on the choice of a set of lattice vectors and δ ij is the Kronecker delta function.For example, T 0 = 1/3 for the so called D3Q19 lattice (albeit it only satisfies (10) up to the 4th order).Notice, one should not confuse the three dimensional array c α with the lattice vector e α (q) in ( 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, It is easily seen that Θ i j (q + c α , q) 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.

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, [19,20,1,21,5] where N α (q, t) is the number of particles belonging to the discrete direction c α in the cell q at time t.We have assumed in (12) a unity time increment (i.e., ∆t = 1) without loss of generality.Ω α (q, t) in ( 12) is the collision term that satisfies local mass and momentum conservations, The particle density distribution function f α (q, t) is related to N α (q, t), via where J(q) is the volume of cell centered at q, as defined previously.The fundamental fluid quantities are given by the standard hydrodynamic moments, where ρ(q, t) and u(q, t) are fluid density and velocity at the location x(q) and time t.Using the relationship in (9), the velocity moment above can also be rewritten as and the velocity component value in the curvilinear coordinate system is given by, ρ(q, t)U i (q, t) Or for simplicity of notation, we define a three-dimensional fluid velocity array U (q, t) ≡ (U 1 (q, t), U 2 (q, t), U 3 (q, t)).Therefore, ( 17) can be equivalently expressed as ρ(q, t)U (q, t) It is immediately seen that ( 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 (13) can also be written as, Often in LBM the collision term takes on a linearized form, [1,18] namely where M αβ and f eq β (q, t) represent a collision matrix and the equilibrium distribution function, respectively.In particular, the so called Bhatnagar-Gross-Krook (BGK) form corresponds to with τ being the collision relaxation time.[22,20,21,5] In order to recover the correct Navier-Stokes hydrodynamics, besides ( 13) and ( 19), the collision matrix needs to satisfy an additional condition, [18,23,24,25,26,27] α Obviously the BGK form trivially satisfies such an additional property.
The extra term δN α (q, t) in ( 12) represents the change of particle distribution due to an effective inertial force associated with the curvature and nonuniformity of a general curvilinear mesh.Obviously δN α (q, t) 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 where N ′ α (q, t) is the post-collide distribution at (q, t) that is equal to the right side of eqn.(12).Due to curvilinearity, though the amount of particles advected from cell q is exactly equal to what is arrived at cell q + c α (as defined in (22)), the momentum is in fact changed along the way.Indeed, in general e α (q + c α )N α (q + c α , t + 1) = e α (q)N ′ α (q, t) 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, where the right side of the unequal sign in (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 ( 23) and ( 22), we see that the net momentum change via advection from all the neighboring cells into cell q is given by, Similarly, we can realize that the net momentum change via advection out of cell q to all its neighboring cells is given by Subsequently, if we impose the constraints on δN α (q, t) below, 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 χ(q, t) = [χ I (q, t) + χ o (q, t)]/2.More specifically, the first constraint in (26) means that no mass source is introduced by δN α (q, t).On the other hand, the second constraint in (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 Using the geometric quantities defined in the previous subsection, it can be immediately shown that where the geometric function Θ i j (q + c α , q) is defined in (11).From ( 27) and (11), one can immediately see that F i (q, t) vanishes if the curvilinear mesh is a regular uniform Cartesian lattice, as is obviously the case for a basic LBM.The second constraint in (26) can also be expressed in coordinate component form as 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, with A specific form of δN α (q, t) can be chosen below With a simple algebra, one can verify that it satisfies the moment constraints of ( 26), ( 28) and (29).Notice, due to the appearance of N ′ α (q, t) in ( 27), the overall collision process for determining N ′ α (q, t) 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, with ρ(q, t)a i (q, t) ≡ F i (q, t).It is straightforward to show that these fundamental conditions are met by the following equilibrium distribution form, The equilibrium distribution form above is analogous to that of a low Mach number expansion of the Maxwell-Boltzmann distribution expressed in curvilinear coordinates.[12] It reduces to the standard LBM equilibrium distribution if the curvilinear mesh is a uniform Cartesian lattice so that g ij = δ ij .[18,7,17] With all the quantities and dynamic properties defined above, we can show that the lattice Boltzmann equation ( 12) computed on the (non-Euclidean) uniform Cartesian lattice {q} obeys the Navier-Stokes hydrodynamics in curvilinear coordinates.Consequently, mapping of the resulting fluid values, ρ(q, t) and U i (q, t), onto the curvilinear mesh is given by a simple transformation below, ρ(x(q), t) = ρ(q, t) , u(x(q), t) = U i (q, t)g i (q) (34)

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 Equation ( 12) below, Expanding it in both time and space up to the second order leads to Now we introduce multiple scales in time and space based on the conventional Chapman-Enskog expansion procedure, [28,19] for they introduce no source of mass.Equating the same powers of ǫ, eqn.(36) leads to the following two equations and where the collision term has taken the linearized form in (20).Eqn.(38) can be directly inverted, Therefore, eqn.(39) can also be written as, β + δN (1)   α (41) Taking the mass moment of eqn.( 38) and (41), and using mass conservation properties of M αβ and δN (n) above, we obtain for the leading order Based on the equilibrium moment definitions (37), eqn.( 42) is equivalent to For the next order, we have Combining ( 42) and (44), we get Substitute the definitions in (37), eqn.(45) becomes where is to be further discussed later in this section.Define the fluid velocity as and since J = J(q) is not a function of time, eqn.( 46) is in fact the standard mass continuity equation expressed in general curvilinear coordinates, A is the divergence of a vector A in a general curvilinear coordinate system, eqn.( 49) is simply the familiar mass continuity equation in the coordinate free representation, Taking the momentum moment of eqn.(38), we get, Define (see also in (32)) together with (47), eqn.(51) becomes Similarly, taking the momentum moment of eqn.(41), and use relation (40), we can derive Using the collision matrix property (21) and its inverse, namely then eqn.(56) becomes Substitute JΠ ij, (1) of the above into (54), we get After some simple cancellations and rearrangements, it becomes Using definitions in (32), N eq α = Jf eq α , and (47), eqn.(60) can be equivalently expressed as where in the above we have defined Jρa (1),i ≡ α c i α δN (1) α .Now let us examine the forcing terms Jρa (0),i and Jρa (1),i .Recall from the momentum conservation condition, the overall forcing term in eqn.( 12) is given by ( 28), that is with JF i given by (27), where Θ i j (q + c α , q) is defined in (11).Taking the continuum limit, we have for the leading order of Θ i j (q + c α , q) Consequently, when substitute ( 63) and ( 64) in (62), we obtain, to the leading order where is the Christoffel symbol as defined in differential geometry.[11] Also in the above, the relationship ( 22) is used, namely Furthermore, since we have from (65) the following Comparing terms of the same order in ǫ, we recognize that and Replacing N (1) α with that in (40), with the use of property (57) and after some straightforward algebra, eqn.(71) becomes α c i α δN (1) From the definition (52) and (69), we see that Similarly, let us also define where Hence, from (72), we have Taking (73) to eqn.(53), we have for the leading order the equation below On the other hand, taking (74) and (76) to eqn.(61), we get in the viscous order Combining ( 77) and (78) together and using the definition in (48), one obtains formally the Cauchy's transport equation, namely where Π jk = Π jk,eq + Π jk,neq .Equation (79) can be further expressed in terms of hydrodynamic quantities.From the moment properties of the equilibrium distribution function in (32), we have which has exactly the same form as that from the continuum kinetic theory.On the other hand, from (74) we have where, from (32), Π ij,eq is given by (80) and Q ijk,eq is given below which has exactly the same form as that from the continuum kinetic theory.However, the term τ α c i α c j α δN (0) α needs to be further evaluated in terms of the hydrodynamic quantities.
From (29) and (30), we have with Taking the long wave length limit, (84) in the continuum limit becomes Therefore, )J(q)[Γ i kl (q)Q jkl,eq (q, t) + Γ j kl (q)Q ikl,eq (q, t)] (86) So we finally obtain from (81) that with Π ij,eq and Q ijk,eq given by ( 80) and (82).Substituting eqns.(80) and (87) into eqn.(79),comparing with eqn.(122) in Appendix B, together with the definitions of (80) and (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 (τ − 1 2 ) 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 Appendices A and B).Notice, same as the derivation in the continuum kinetic theory, in (87) the 0th order time derivative term ∂ t0 Π ij,eq is to be replaced by the leading (Euler) order equations on the right hand sides of eqns.(77) and (43).Without repeating the rest of algebra, here we present the final form -the Navier-Stokes equation in a general curvilinear coordinate system, where the pressure p = ρT 0 , µ = (τ − 1 2 )ρT 0 , and with the standard covariant derivative of the velocity component defined as Using the standard definitions of differential operators in differential geometry, [11] we can recognize that eqn.( 88) is indeed the familiar Navier-Stokes equation in a generic coordinate-free operator representation,

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 [12] to discrete space and time.Unlike some previous works, [14,15,16] here we use a volumetric representation so that conservation laws are exactly ensured.[9] 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, [14,15,16] 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,23,24] 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.
where the integral operator above is defined as dv ≡ dv 1 dv 2 dv 3 .Substituting (95) into (96), we get Define a particle density function then we have the hydrodynamic moments specified below dvN (q, v, t) = J(q) dvf (q, v, t) = J(q)ρ(q, t) dvv i N (q, v, t) = J(q) dvv i f (q, v, t) = J(q)ρ(q, t)u i (q, t) (100) Taking the moment integral of the Boltzmann equation (98), and using the collision properties of (97), we obtain the two continuity equations, corresponding to mass and momentum conservations respectively The term dv ∂ ∂v i (v j v k Γ i jk N ) = 0 via integration by parts, and using definition in (100) we get the mass continuity equation as follows or in the more familiar form since the volume J is not dependent on time t.
Integrate by parts, and use ∂v i ∂v j = δ i j , Hence, the second equation in (101) becomes where the momentum flux tensor Π ij = Π ij (q, t) is defined by Eqn.( 104) is known as the Cauchy's transport equation.We can separate the momentum flux tensor into two parts associated, respectively, to the equilibrium and the non-equilibrium parts of the distributions f (q, v, t) = f eq (q, v, t)+f neq (q, v, t), so that Π ij (q, t) = Π ij,eq (q, t)+Π ij,neq (q, t).The equilibrium distribution is given by the Maxwell-Boltzmann form, where U ≡ v − u and θ is the temperature.In terms of curvilinear coordinates, Hence, we can rewrite (106) in terms of curvilinear coordinates below and the normalization factor Here, the inverse metric tensor [g ij ] is defined such that g ik g kj = δ i j in differential geometry.det[g ij ] is the determinant of [g ij ].Using a few basic properties of Gaussian integral, we immediately obtain from (105) and (107) that ρ = dvf eq , ρu i = dvv i f eq , Π ij,eq = dvv i v j f eq = g ij ρθ + ρu i u j (109) Therefore, at the Euler order, in which the momentum flux tensor only includes the equilibrium contribution, eqn.(104) reduces to However, since the underlying space is "flat" (Euclidean), the metric tensor obeys the following property Substituting (112) into eqn.(111),we arrive at a more standard form of the Euler equation, with the pressure defined by an ideal gas equation of state, p = ρθ.

B Derivation of the Navier-Stokes Equations in General Coordinates
To derive the Navier-Stokes hydrodynamics up to the viscous order, we use the Chapman-Enskog expansion procedure,[28] Here ǫ (<< 1) denotes a small number.Thus, the Boltzmann equation (98) leads to the following two equations, and where, for simplicity, in the above we have used the BGK collision operator form Ω = −(N − N eq )/τ .[22] Taking the mass and momentum moments over (114), and use the properties in (109) as well as conservation of mass and momentum by the collision in (97), we immediately obtain the leading (Euler) hydrodynamics as given by ( 103) and (113), with an ideal gas equation of state, p ≡ ρθ.
Eqn.(114) can be inverted to give, With the definition of particle density distribution function (99), N (q, v, t) = J(q)f (q, v, t), and J depends on q only, hence eqn.( 117) is equivalent to With (116), it is easily checked via straightforward algebra that N (1) (and f (1) ) gives vanishing mass and momentum moments, namely dvN (1) (q, v, t) = 0, On the other hand, taking the momentum flux moment, we have from (117) the following, The last term in (119) above can be further simplified below, Substituting (120) into (119), we obtain Or equivalently, we have where the equilibrium heat flux Q ijk,eq is defined as The integration (123) is straightforward to perform with f eq (given by ( 107)), and it yields Therefore, together with (109), eqn.(122) becomes Next, let us take into account the properties of the metric tensor in the flat space, these are Then we can further simply (125) to + ∂ ∂q k (Jρu i u j u k ) + Γ i kl Jρθg jk u l + Γ i kl Jρu j u k u l + Γ j kl Jρθg ki u l + Γ j kl Jρu i u k u l } (127) For the purpose of deriving hydrodynamics for isothermal fluid, we set θ = const, in addition we substitute the Euler order relations from (116) for the terms with ∂ t0 , then after some straightforward algebra, eqn.(127) becomes Π ij,neq = −τ ρθ{g ik ( ∂u j ∂q k + Γ j kl u l ) + g jk ( Using the standard definition of covariant derivative in differential geometry, together with the rate of strain definition as the symmetric form of velocity derivative, C Incompressibility in Curved Space with Underlying Euclidean Metrics The size of a volume element in a 3 2 single particle phase-space is denoted as J p dqdv, where dq ≡ dq 1 dq 2 dq 3 and dv ≡ dv 1 dv 2 dv 3 .Clearly, for a curvilinear coordinate system, J p = J 2 = g, with g ≡ det[g ij ] being the determinant of the metric tensor, and J = J(q) is the Jacobian of the curvilinear coordinates.The transport of a particle phase-space density function φ = φ(q, v, t), without collision, follows a continuity equation in phase space, ∂ t (J p φ) + ∂ ∂q i (J p φv i ) + ∂ ∂v i (J p φ vi ) = 0 (133) Since a particle moves along a straight path (geodesics) for the underlying Euclidean space according to the 1st law of Newton, we have (see Appendix A) We can show that the motion of particles in phase-space is incompressible.Define the velocity field in phase-space as, Ẋ ≡ ( q, v) = (v, v) (135) then, the divergence of the velocity field in phase-space is given by Substitute (134) in eqn.(136), we get But since both J p (= J 2 ) and Γ i jk are functions of q only, the second term in (137) becomes where the last equality in the above is due to the symmetry of the Christoffel symbol, Γ i jk = Γ i kj .Furthermore, using the fundamental property then (138) becomes, ∂J ∂q j v j = ∂J p ∂q j v j = ∂ ∂q j (J p v j ) (139) the last equality is because ∂ ∂q j v j = 0. Plug the result of (139) into (137), hence we have proved the incompressibility property of the velocity field in the phase-space.That is, With the incompressibility property of the phase space velocity field, the continuity equation (133) takes on a form of the Vlasov equation, It is worth to point out, the incompressibility property of the phase-space velocity field is an intrinsic property of the particle motion via Newtonian mechanics in Euclidean space.Hence this is true in any coordinate systems whether Cartesian or curvilinear.Nonetheless, it is useful to directly show that such a property is preserved in a curvilinear coordinate system.
2 (g jk u i | k + g ik u j | k ) then eqn.(128) further simplies to Π ij,neq = −2µS ij(129)where µ ≡ τ ρθ is the dynamic viscosity.Taking mass and momentum moments of eqn.(115), and the vanishing righthand side, we obtain Combining the first order hydrodynamics of (130) with that of the Euler order (116), we finally arrive at the full Navier-Stokes hydrodynamics in a curvilinear coordinate system (Jρu i u j ) + Γ i jk ρu j u k = −g ij ∂p ∂q j with an ideal gas equation of state, p ≡ ρθ.Eqns.(131) are in fact the mass continuity and the Navier-Stokes equation in coordinate-free operator forms