Three-Dimensional Central Moment Lattice Boltzmann Method on a Cuboid Lattice for Anisotropic and Inhomogeneous Flows

We present a new 3D lattice Boltzmann (LB) algorithm based on central moments for the D3Q27 lattice using a cuboid grid, which is parameterized by two grid aspect ratios that are related to the ratios of the particle speeds with respect to that along a reference coordinate direction. The use of the cuboid lattice grid enables the method to compute flows having different characteristic length scales in different directions more efficiently. It is constructed to simulate the Navier-Stokes equations consistently via introducing counteracting corrections to the second order moment equilibria obtained via a Chapman-Enskog analysis that eliminate the errors associated with the grid anisotropy and the non-Galilean invariant terms. The implementation is shown to be compact and modular, with an interpretation based on special matrices, admitting ready extension of the standard algorithm for the cubic lattice to the cuboid lattice via appropriate scaling of moments based on grid aspect ratios before and after collision step and equilibria corrections. The resulting formulation is general in that the same grid corrections developed for the D3Q27 lattice for recovering the correct viscous stress tensor is applicable for other lattice subsets, and a variety of collision models, including those based on the relaxation of raw moments, central moments and cumulants, as well as their special case involving the distribution functions. The cuboid central moment LBM is validated against a variety of benchmark flows, and when used in lieu of the corresponding raw moment formulation for simulating shear flows, we show that it results in significant improvements in numerical stability. Finally, we demonstrate that our cuboid LB approach is efficient in simulating anisotropic shear flow problems with significant savings in computational cost and memory storage when compared to that based on the cubic lattice.


Introduction
The lattice Boltzmann (LB) methods, which arise as minimally discretized numerical schemes of the Boltzmann transport equation -a cornerstone formulation in kinetic theory, has attracted much attention in the recent decades [1,2,3,4,5]. As a mesoscopic approach, it has enriched the variety of computational fluid dynamics (CFD) techniques that are being developed and has been applied to a wide range of fluid flows successfully [6,7,8,9]. It involves the evolution of the distribution of the particle populations by a collision step, which is then followed by their lock-step advection along discrete directions, referred to as the streaming step. The former is often modeled by a relaxation of the distributions (e.g., [10]) or their moments (e.g., [11]) to equilibria. The discrete particle velocity directions are referred to as the lattice, which are usually designed to obey the associated physical symmetry and isotropy of the fluid flow being simulated.
The lattice is usually referred to by the notation DdQq, where d represents the number of spatial dimensions and q denotes the number of discrete particle directions. The macroscopic fields such as the fluid velocity are obtained from the lower order moments of the distribution functions, while the higher order kinetic moments can be constructed to evolve so as to facilitate the numerical robustness of the method. The linear advection and locally nonlinear collision of the LB method, along with its ability to naturally represent the physics of complex fluids and flows based on kinetic models, are among the important assets of this approach.
Anisotropic and inhomogeneous fluid motions involve relatively large spatial variations in the characteristic features of the fluid flow in one or more directions when compared to that in the other directions. These include shear flows, flows around boundary layers and in compacted geometrically disordered media, where the characteristic length scales or spatial gradients in the flow are dominant in certain directions relative to others. Such multidimensional flows governed by the Navier-Stokes equations can be more efficiently computed using methods which use grids that naturally conform with the direction-specific non-homogeneities inherent to the problem of interest. The LB methods, on the other hand, are usually constructed using uniform lattice grids, such as using a square lattice in two-dimensions (2D) and a cubic lattice in three-dimensions (3D) in order to satisfy their symmetry and isotropy requirements. One approach to overcome this issue is to employ nonuniform grids with the LB method augmented using interpolations (e.g., [12,13]), which however introduce considerable additional numerical dissipation compromising the accuracy of the approach [14]. Alternatively, conventional discretizations, such as finite difference, finite volume or finite element approaches (e.g., [15,16]) 2 could be used which however entail additional complexity and overhead to the LB schemes. It is thus desirable to use the standard LB discretization along the particle characteristics that preserves the lock-step or perfect streaming with relatively low attendant numerical dissipation while allowing the use of different particle speeds along different directions. Such types of LB schemes are associated with the use of rectangular lattice grids in 2D and cuboid lattice grids in 3D.
Starting from the initial work of Koelman [17] much focus has been given to the construction of the LB methods on rectangular grids during the last two decades using different collision models with some necessary modifications to the algorithm to satisfy the Navier-Stokes (NS) equations. In the simpler rectangular LB versions, the single-relaxation-time (SRT) model has been modified by including additional particle velocities and whose equilibria were constructed via solving a quadrature problem [18], or by extending the equilibrium distribution functions with additional corrections to restore the isotropy effects [19], or by including some counteracting forcing terms [20]. On the other hand, by exploiting the additional degrees of freedom existing in the multiple-relaxation-time (MRT) collision model based on raw moments, rectangular LB schemes were constructed in the following four different ways: (i) by coupling between various relaxation parameters and the grid aspect ratio via a linear stability analysis [21], (ii) by keeping the matrix of the transformation between the moment space and the velocity space of the distribution functions independent of the grid aspect ratio [22], (iii) by using an additional adjustable parameter that determines the relative orientation in the energy-normal stress subspace via an inverse design analysis based on the Chapman-Enskog expansion [23], and (iv) by extending the equilibrium moments to include the stress components for restoring the isotropy of the recovered macroscopic equations [19]. It should be pointed out that Zong et al [23] demonstrated that the emergent continuum limit equations of some of the earlier rectangular LB formulations [17,21,22] were not fully consistent with the NS equations. Moreover, although the latter rectangular LB schemes [23,19] yielded physically correct formulations, they require cumbersome implementations involving the specification of many free parameters and complicated expressions for the transport coefficients and mapping matrices dependent on the grid aspect ratio, and, like the other consistent schemes [18,20], are subject to major stability issues at lower grid aspect ratios when simulating flows at relatively large characteristic flow velocities or low viscosities.
The rationale for the limitations of the prior rectangular LB schemes were clarified in our recent work [24,25,26]. These include their choice of the orthogonal moment basis, construction of the discrete equilibria involving only the lower order velocity terms and without correcting for the non-Galilean invariant cubic velocity errors arising from aliasing effects, and the use of collision models based on raw moments. Moreover, it has been pointed out in an earlier work for standard lattices (square or cubic) [27,28] that the use of central moments and avoiding the orthogonalization of the moment basis leads to significant stability improvements. The construction of the collision step in the LB formulations using central moments based on the peculiar velocity [29,30] naturally maintains the Galilean invariance of the moments independently supported by the lattice and its advantages have been demonstrated for a variety of fluid dynamical problems 3 (see e.g., [31,32,33,34,35,36]). Based on these considerations, we proposed a 2D central moment rectangular LB scheme recently and demonstrated its superior numerical features for simulating flows at higher Reynolds numbers using relatively small grid aspect ratio when compared to the other existing LB methods based on the rectangular lattice [24,25,26].
Since anisotropic and inhomogeneous flows in situations of practical interest are often 3D in nature, it would be beneficial to develop LB algorithms on cuboid lattices. Their exist few prior studies in this regard.
For example, Hegele et al [18] presented a SRT-LB scheme using a D3Q23 lattice that involves the use of two additional particle velocities in each of the two Cartesian directions embedded to the D3Q19 lattice and constructed the equilibrium distribution functions by solving a quadrature problem to ensure the desired isotropy. It evaluated the accuracy of the resulting approach for simulating cylindrical waves by using grid sizes with only small deviations from the cubic lattice. Later, Jiang and Zhang [37] presented a cuboid LB scheme using a D3Q19 lattice for simulations of porous media. However, due to severe numerical stability restrictions, it was able to use grid aspect ratios with minor variations from unity. More recently, Wang et al [38] presented a raw moment-based MRT-LB scheme on the D3Q19 cuboid lattice which utilized an orthogonal moment basis. The equilibria and the transport coefficients were constructed via an inverse design analysis to satisfy the Navier-Stokes equations and involved the specification of many free parameters entailing a cumbersome implementation of the method. The resulting numerical approach was validated for some benchmark flow problems at moderate grid aspect ratios. Here, it should be pointed out that the choice of the equilibria plays a crucial role in maintaining the physical isotropy of the fluid flow equations [39].
This can be naturally constructed by matching with the continuous Maxwell distribution as done in central moment LB formulations rather than involving complicated fitting of parameters. Moreover, as mentioned above, the use of central moments is expected to provide significant improvements to the state-of-the-art in LBM based on non-cubic grids. However, a 3D central moment LBM based on a cuboid lattice does not currently exist in the literature, which can be constructed via an extension and significant modification of our recent work on the rectangular lattice [26] and is the main objective of this present paper.
In this study, we will present a new 3D cuboid central moment LB method on a D3Q27 lattice and is referred to as the 3DCCM-LBM in what follows, which is parameterized by two grid aspect ratios representing the ratios of the characteristic particle speeds along the two directions with respect to those in the remaining direction. The D3Q27 lattice is chosen since it provides greater possible isotropy and accuracy when compared to its D3Q15 and D3Q19 lattice subsets [40]; however, the equilibria corrections needed for the D3Q27 cuboid lattice to satisfy the NS equations (see below for more details) are applicable for these lattice subsets as well.
In contrast to the prior approaches, our cuboid LB scheme will be constructed by using a natural moment basis that avoids orthogonalization and the discrete central moment equilibria are specified by matching with those corresponding to the continuous Maxwell distribution function. Moreover, the prior LB algorithms on stretched lattice grids based on the various collision models discussed above, including our recent 2D rectangular central moment LB scheme [26], prescribe the transformations between the distribution functions 4 and raw moments before and after collision by using a moment basis that separates the trace of the second order moments (related to bulk viscosity) from its other components (related to shear viscosity). The use of this strategy in the context of the 3D formulations based on the cuboid lattice would result quite complicated transformation matrices dependent on the grid aspect ratios. This is obviated in our present approach by naturally separating the bulk and shear viscosity effects associated with the viscous stress only within the collision step and not for such transformations, and the effect of the grid aspect ratios are introduced via much simpler pre-and post-collision diagonal scaling matrices in such mappings. We will show the resulting compact approach can be naturally interpreted based on special matrices. The truncation errors resulting from the grid anisotropy associated with the cuboid lattice as well as those due to the aliasing effects will be eliminated by deriving the necessary correction terms from a Chapman-Enskog analysis, which will be augmented to the second order moment equilibria. The resulting 3DCCM-LBM is modular in construction in that the existing 3D central moment (and also its special case involving raw moment) based algorithms developed for the cubic lattices can be readily extended to cuboid lattices by using the corrections to the moment equilibria derived in this work and introducing the pre-and post-collision diagonal scaling matrices based on the grid aspect ratios. Our cuboid LB approach will be first validated against the analytical solutions and/or numerical results for some standard benchmark flow problems. The advantages of this scheme in efficiently simulating an example inhomogeneous and anisotropic flow problem with significant savings in memory storage and computational effort will then be demonstrated. Moreover, we will also show numerical stability improvements in the use of our 3D central moment scheme when compared to that based on raw moments for computation of shear flows at relatively large flow velocities and/or low viscosities.
The organization of this paper is as follows. Section 2 presents a Chapman-Enskog multiscale analysis of a 3D LB equation based on the non-orthogonal moment basis using a D3Q27 cuboid lattice. Since the viscous stress tensor in the NS equations is based on the second order non-equilibrium raw moments that are the same as the corresponding second order non-equilibrium central moments, it suffices to present our analysis based on the simpler raw moment formulation; and as part of this, we will identify the correction terms involving velocity gradients and grid aspect ratios for eliminating the grid anisotropy and non-Galilean invariant cubic velocity related truncation errors and show consistency with the 3D NS equations. Formulas for computing the local strain rate tensor based on non-equilibrium moments and parameterized by the grid aspect ratios will also be derived in this section. In Sec. 3, we will discuss the complete details of the implementation aspects of our 3DCCM-LBM. Additional algorithmic details in this regard are given in various appendices (see Appendix A -Appendix D). Formulas for the momentum-augmented bounce-back scheme for the D3Q27 cuboid lattice that are dependent on the grid aspect ratios for simulating shear flows due to moving boundaries are derived and summarized in Sec. 4. Section 5 presents a numerical validation study of the 3DCCM-LBM using different grid aspect ratios for a variety of benchmark flow problems. This is followed by a demonstration of the computational effectiveness of using cuboid grids (via 3DCCM-LBM) in lieu of utilizing cubic grids for simulating inhomogeneous and anisotropic flows in Sec. 6 and the numerical stability improvements achieved with using central moments rather than raw moments on the D3Q27 cuboid lattice at different grid aspect ratios in Sec. 7. Finally, the main conclusions of this work are given in Sec. 8.

Chapman-Enskog Analysis on a D3Q27 Cuboid Lattice: Isotropy Corrections, Macroscopic
Flow Equations, and Local Formulas for the Strain Rate Tensor 2.1. Cuboid lattice parameters, moment basis, and definitions of central moments and raw moments The cuboid grid based on the three dimensional, twenty seven velocities (D3Q27) lattice used in deriving the formulation in what follows is shown in Fig. 1. Figure 1: Three-dimensional, twenty seven particle velocities (D3Q27) cuboid lattice. The components of the particle velocities e = (e x , e y , e z ) in the y and z directions are related to the grid aspect ratios r and s, respectively, which is defined below. In other words, these two free parameters represent the ratios of the characteristic particle speeds along the y and z directions with respect to those in the x direction, and formalize the flexibility accorded by the cuboid lattice. Thus, if ∆x, ∆y and ∆z are the space steps in the x, y and z directions, respectively, for evolving over a time step ∆t that would then fix the particle speeds in the respective directions (i.e., c x = ∆x/∆t, c y = ∆y/∆t and c z = ∆z/∆t), we can then define the grid aspect ratios as r = ∆y/∆x, and s = ∆z/∆x. As is standard in LB formulations, we work with the usual lattice units for simplicity in the following, i.e., the reference space step in the x direction and the time step are taken to be unity: ∆x = 1 and ∆t = 1. Thus, ∆y = r and ∆z = s. We now clarify the notations used in Fig. 1, where the magnitude of each of the particle velocity is the distance between the head and tail of an arrow representing that direction; in lattice units, with a unit time step, the Cartesian components of the length of each lattice direction or the streaming distance are bounded by (1, r, s). Now, since in Fig. 1, for every particle velocity direction, its opposite counterpart is also depicted, the total distance for such a pair encompasses a length with components (2, 2r, 2s) as shown. Based on these considerations, we can then list the Cartesian components of the particle velocities for the D3Q27 cuboid lattice as follows: |e y = (0, 0, 0, r, −r, 0, 0, r, r, −r, −r, 0, 0, 0, 0, r, −r, r, −r, r, r, −r, −r, r, r, −r, −r) † , 6 where |· denotes a column vector based on the standard 'ket' notation and † refers to taking the transpose of any array. We will also need the following 27-dimensional vector with unit elements in what follows: (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 We then define a linearly independent set of non-orthogonal basis vectors for the D3Q27 cuboid lattice using a combination of the monomials of the type |e m x e n y e p z for integer exponents m,n and p [30] as follows: where [30] |T 0 = |1 , As discussed in the introduction section, we will not orthogonalize the above set of basis vectors further to maintain simplicity and robustness and retain them in their natural forms in the following derivation based on the 3D cuboid lattice analogous to our rectangular central moment LB formulation [26]. It should be noted that while the above basis vectors segregate the trace of the second order components |e 2 x + e 2 y + e 2 z from the other second order components in order to enable the independent specification of the bulk viscosity from shear viscosity, as mentioned earlier, in the algorithmic implementation of the central moment approach discussed in the next section (see Sec. 3), these operations will be confined only within the collision step and not for performing any mappings between distribution functions and moments. These considerations are crucial in avoiding cumbersome transformations and corrections for eliminating any truncation errors in order to develop an efficient algorithm on a cuboid lattice.
Then, we list the vectors of the distribution functions f , their equilibria f eq and the source terms S accounting for the effect of any body force with components F = (F x , F y , F z ) experienced by the motion of the fluid with density ρ and velocity u = (u x , u y , u z ) for the D3Q27 cuboid lattice as f = (f 0 , f 1 , f 2 , . . . , f 26 ) † , f eq = (f eq 0 , f eq 1 , f eq 2 , . . . , f eq 26 ) † , S = (S 0 , S 1 , S 2 , . . . , S 26 ) † . In anticipation of the developments in Sec. 3, we first define the central moments of the distribution functions f α and their equilibria f eq , as well as the source term S α of order (m + n + p) using the weights as the peculiar velocity components, i.e., the components of 7 the particle velocity shifted by those of the fluid velocity, as follows: Also, we define the bare or raw moments of the above three quantities of order (m + n + p) by using just the particle velocity components as the weights, which are given by Note that here and in the following any quantity with a prime notation ( ) is used to distinguish it as a raw moment (as opposed to a central moment). For convenience, we now enumerate these raw moments as , n eq = k eq 000 , k eq 100 , k eq 010 , k eq 001 , k eq 110 , k eq 101 , k eq 011 , (k eq 200 − k eq 020 ), (k eq 200 − k eq 002 ), (k eq 200 + k eq 020 + k eq 002 ), (k eq 120 + k eq 102 ), (k eq 210 + k eq 012 ), (k eq 201 + k eq 021 ), (k eq 120 − k eq 102 ), (k eq 210 − k eq 012 ), (k eq 201 − k eq 021 ), k eq 111 , (k eq 220 + k eq 202 + k eq 022 ), (k eq 220 + k eq 202 − k eq 022 ), (k eq 220 − k eq 202 ), k eq 211 , k eq 121 , k eq 112 , k eq 122 , k eq 212 , k eq 221 , k eq Here, the distribution functions and raw moments are related via T, i.e., the basis vectors (Eq. (4)), as n = Tf , n eq = Tf eq , Ψ = TS,

The lattice Boltzmann equation
It should be noted that the use of a cuboid lattice would result in an anisotropic viscous tensor given in terms of the grid aspect ratios r and s via the second order non-equilibrium moments, and such grid-related anisotropy effects need be eliminated via carefully designed correction terms. In constructing such counteracting corrections, since the second order non-equilibrium raw components are the same as the corresponding central moments by definition, it is enough to perform an analysis on the following simpler lattice Boltzmann equation (LBE) on the raw moments, i.e., the MRT-LBE [26]: where Λ = diag ω 0 , ω 1 , ω 2 , . . . , ω 26 is the diagonal relaxation time matrix and I is a 27 × 27 identity matrix.
Here, the right side of Eq. (8) prescribes the relaxation of different moments n j to their equilibria n eq j at a rate given by the parameter ω j , which is augmented by the effect of the source term via Ψ j (j = 0, 1, 2, . . . , 26) and the results are then mapped back into the velocity space via the inverse mapping T −1 . This is followed by their perfect advection of the distribution functions f α to the nearest node along the particle characteristics e α (α = 0, 1, 2, . . . , 26) as implied by the left side Eq. (8). Then, the hydrodynamic fields are computed via appropriate moments as
where c s is the speed of sound, which is an adjustable parameter of the collision model and will be related to the transport coefficients via a Chapman-Enskog analysis later. It should be noted, however, that the analysis in what follows that provides corrections for the second order moments, and the algorithm devised in the next section (Sec. 3) for the cuboid lattice are still applicable for other collision models. These include those based on the factorization property [41,30] and cumulant collision model [28], which differ from the Maxwellianbased central moment model in the evolution of the higher order moments. Then, the corresponding discrete raw moment equilibria are obtained from Eq. (10) via binomial transforms, which read as follows: k eq 000 = ρ, k eq 100 = ρu x , k eq 010 = ρu y , k eq 001 = ρu z , k eq 110 = ρu x u y , k eq 101 = ρu x u z , k eq 011 = ρu y u z , In addition, the central moments of the source terms for recovering the NS equations in 3D are given by [30] σ 000 = 0, σ 100 = F x , σ 010 = F y , σ 001 = F z , σ mnp = 0 if (m + n + p) ≥ 2, and the corresponding raw moments follow from their binomial transform as [30] σ 000 = 0, 2.4. Chapman-Enskog Analysis: Identification of truncation errors due to grid anisotropy and non-Galilean invariance from aliasing effects on the D3Q27 cuboid lattice We will now perform an analysis based on the Chapman-Enskog (C-E) multiscale expansion [42], written for the LBE in the matrix form by d'Humières [43], to identify the truncation errors due to anisotropic effects arising from using the cuboid lattice grid and non-Galilean invariant (GI) cubic velocity errors resulting from the aliasing effects associated with the D3Q27 lattice. Corrections to eliminate these errors will then be derived in what follows. The analysis for a 3D LBE formulation given above follows the approach presented in Premnath and Banjeree [44] and Hajabdollahi and Premnath [45], which was adopted for constructing a rectangular LB algorithm recently [26]. In this version of the C-E analysis, we expand the moments about their equilibria by including the non-equilibrium effects as higher order perturbations and the time derivative as a multiple scale expansion as follows: where = ∆t is the perturbation parameter that serves the purpose of bookkeeping and isolating terms of the different orders. Moreover, we apply a multivariate Taylor where E i = T (e i I)T −1 and e i = |e i with i ∈ (x, y, z). We then re-express Eqs. (15b)and (15c) in the long form, which respectively read as O( ) moment system: O( 2 ) moment system: (2) . (17) Substituting the moment equilibria n (0) listed in Eq. (11) into the O( ) Eq. (16), we now write explicitly all the equations up to the second order moment components which are relevant to determining the fluid dynamical behavior as follows: ) are influenced by the grid aspect ratios r and s. Hence, it follows that such underlined terms will impact the viscous stress tensor and hence the hydrodynamics. Thus, in order to obtain the complete picture, we show in the following the evolution of the conserved moments at the O( 2 ) level, i.e., the leading four moment components of Eq. (17), which manifest the effect of the non-equilibrium parts on the evolution of density and momentum components at the slower time scale t 1 : Evidently, the hydrodynamical behavior given in Eqs. (19b)-(19d) is influenced by the grid aspect ratios via the non-equilibrium moments n (1) 7 , n (1) 8 and n (1) 9 , which needs to eliminated. Before accomplishing this, we need the complete expressions of the second order non-equilibrium moments in order to isolate the truncation errors from terms that correspond to physics. Hence, we rewrite Eqs. (18e)-(18j) after segregating the terms associated with the grid aspect ratios r and s (see underlined terms) from those associated with the standard cubic lattice (terms within the brackets {· · · }) as follows: It should be noted here that diagonal parts of the non-equilibrium moments n (1) contain contributions from non-GI cubic velocity terms due to aliasing effects associated with the D3Q27 lattice (e.g., α f α e 3 αi = α f α e αi , where i ∈ {x, y, z}) which appear together with the physical terms within those enclosed within the brackets {· · · }) of Eqs. (20d)-(20f) in addition to the grid-dependent anisotropic errors.
The grid-anisotropy related error terms will be denoted by E js , while the non-GI cubic velocity truncation errors will be referred to as E jg for j = 7, 8 and 9. The structures of the non-equilibrium moments can then be obtained from Eqs. (20a)-(20f) after substituting for the time derivatives of the conserved moments by means of terms related to the spatial derivatives using Eqs. (18a)-(18d) and simplifying the resulting equations by retaining terms of O(u 3 i ) (following Hajabdollahi and Premnath [45]). The off-diagonal second order nonequilibrium moments n (1) Finally, the diagonal second order non-equilibrium moments n 002 and n (1) 002 are given as follows: where the truncation errors due to grid-anisotropy E 7s , E 8s and E 9s and non-GI aliasing effects E 7g , E 8g and E 9g can be written as 2.5. Elimination of errors due to lattice stretching and non-GI cubic velocity terms via extended moment equilibria In order to eliminate the truncation errors identified above in Eqs. (23a)-(23f) associated with Eqs. (22a)-(22c), we now extend the moment equilibria n eq to an effective moment equilibria n eq,eff as n eq,eff = n eq + ∆tn eq(1) , where n eq(1) are the required correction terms. Since the truncation errors appear in the diagonal parts of the second order non-equilibrium moments only, it suffices to consider non-zero correction terms only for those corresponding components identified by the indices 7, 8 and 9 in the effective equilibria. Hence, we write The expressions in Eq. (25) are inspired by an inspection of Eqs. (23a)-(23f) and involve the diagonal components of the velocity gradient tensor ∂ x u x ∂ y u y and ∂ z u z and the density gradients ∂ x ρ, ∂ y ρ and ∂ z ρ and their unknown coefficients θ 7x , θ 7y , θ 8x , θ 8z , θ 9x , θ 9y , and θ 9z and λ 7x , λ 7y , λ 8x , λ 8z , λ 9x , λ 9y , and λ 9z . The formulas for these coefficients will now be established by performing a modified C-E analysis that accounts for the effective moment equilibria introduced in Eq. (24) into the expansion. That is, n = n eq,eff + n (1) + 2 n (2) + · · · = n eq + n eq(1) + n (1) + 2 n (2) + · · · , Then, repeating the steps performed in Sec. 2.4, the respective moment equations of O( k ) for k = 0, 1 and 2, respectively, in Eqs. (15a)-(15c) are replaced with the following: where E jg and E js (j = 7, 8 and 9) are truncation error terms associated with grid anisotropy and non-GI due to aliasing, which are given in Eqs. (23a)-(23f) and n eq(1) 7 , n eq(1) 8 and n eq(1) 9 are the corresponding yet to be determined corrections.
It now remains to obtain the desired expressions for the correction terms. In this regard, we combine the moment equations (Eq. (26b)), in particular, those for the hydrodynamic moments (mass and momentum components) evolving at time scale t 0 with times the corresponding equations (Eq. (26c)) based on the time scale t 1 variations. Then, using ∂ t = ∂ t0 + ∂ t1 , we obtain the effective evolution equations for all the moments, including the macroscopic fields related to the mass and momentum components. In particular, the momentum equations will then contain the truncation error terms and the counteracting correction terms along with those associated with the physics. We can then isolate the terms related to the truncation errors and corrections from those related to the desired NS equations and set the combined effect of the former to be zero. This would then lead to constraint relations between the errors and the required corrections as discussed in detail in Hajabdollahi and Premnath [45]. Thus, if we define the following vector containing the truncation errors as where Then, it follows from Eqs. (26b),(26c) and Eqs. (27a)-(27f), the required constraint equation between the moment equilibria corrections vector n eq(1) identified in Eq. (24) and the vector of truncation errors Ξ is given by Specifically, it reduces to n eq(1) j In Eq. ((31)), we do not assume the summation convention of repeated indices. Now, evaluating Eq. (31) by applying Eq. (25) for j = 7 and using Eqs. (23a) and (23b), we get Comparing, we finally obtain the following required expressions for the coefficients θ 7x , θ 7y , λ 7x and λ 7y : Repeating these steps by applying Eqs. (31) and (25) together for j = 8 and j = 9 and invoking Eqs. (23c) and (23d), and Eqs. (23e) and (23f), respectively, i.e., we arrive at the required formulas for the coefficients θ 8x , θ 8z , λ 8x , λ 8z , and θ 9x , θ 9y , θ 9z , λ 9x , λ 9y , and λ 9z as follows: and When the extended moment equilibria (Eq. (24)) using the corrections (Eq. 25)) with these above coefficients are used in the LBE formulation based on the D3Q27 cuboid lattice, it recovers the 3D NS equations given where p = c 2 s ρ represents the pressure field, and the bulk viscosity ξ and shear viscosity ν are related to the relaxation parameters of the second order moments as Before concluding this section, the following important remarks related to our above derivation are in order: • The expressions for the moment equilibria corrections and the transport coefficients involve just the minimally required set of free parameters, viz., the grid aspect ratios r and s and the speed of sound c s , and are considerably simpler than those used in a recent work [38] that used a set of orthogonalized moment basis on the D3Q19 cuboid lattice and many additional parameters in the specification of the equilibria. Moreover, unlike the previous study [38], the equilibria used in our approach given in Eq. (11) involve higher order velocity terms obtained from matching with the moments of the continuous Maxwell distribution thereby naturally eliminating the non-GI truncation errors. This aspect along with the use of non-orthogonal moment basis in the present derivation that avoids a spurious coupling of the evolution of the higher order moments with those of the lower moments existing in approaches that use orthogonal moment basis [27] result in a robust LBE formulation on a cuboid lattice.
• The results derived here for D3Q27 cuboid lattice for the counteracting corrections and the transport coefficients are also equally valid for the D3Q15 and D3Q19 cuboid lattice sets when they use subsets of the same non-orthogonal moment basis.
• The formulas derived in this section for the equilibria corrections necessary to eliminate the grid anisotropy and non-GI error are general and applicable for a wide variety of collision models. For example, if k mnp , k mnp and c mnp correspond to raw moments, central moments and cumulants, respectively, of order (m + n + p), owing to the definitions of these quantities, it follows that their second order non-equilibrium components are identical, i.e., k mnp for (m + n + p) = 2. From this it follows that the formulas obtained above can be directly incorporated into either raw moment-, central moment-or cumulant-based 3D LB algorithms extended to using a cuboid lattice. Moreover, we point out that our derivation can be readily used to construct even an SRT-based LB scheme using a cuboid lattice by setting all the relaxation parameters equal to one another and then constructing the equilibrium distribution functions including the necessary corrections in the velocity space via using f eq,eff = T −1 n eq,eff = T −1 (n eq + ∆tn eq (1) ). However, as demonstrated in our recent work on the rectangular lattice [26], the SRT-and raw moment-based LB schemes are less stable compared to those based on central moments; hence the present investigation will focus on further developing the derivation presented here into an efficient and robust 3D LB algorithm on the D3Q27 cuboid lattice and their numerical study in the following sections. Such an implementation for a cumulant collision kernel on a cuboid lattice is a subject for a future study.
• When r = s = 1 and c 2 s = 1/3, all the correction terms shown above become equal to zero, and the previous formulations applicable for the cubic lattice sets are recovered.
• Finally, when the interest in simulating relatively very low Mach number flows, the non-GI cubic velocity errors can become relatively small and can thus be neglected. Under such conditions, the formulas derived in the above for the corrections lead to further simplifications. In particular, the 18 coefficients determined above reduce to the following: and, thus, we do not require the computation of the density gradients; moreover, the local expressions for the velocity gradients derived in the next section will become more compact and can be readily re-derived when required.

Strain rate tensor components based on non-equilibrium moments for the D3Q27 cuboid lattice
The moment equilibria corrections in Eq. (25) depend on the diagonal parts of the strain rate tensor ∂ x u x , ∂ y u y and ∂ z u z . These components, along with the off-diagonal components, i.e., 1 2 (∂ x u y + ∂ y u x ), 1 2 (∂ x u z + ∂ z u x ) and 1 2 (∂ y u z + ∂ z u y ) can be obtained locally from the second-order non-equilibrium moments as follows. From Eq. (27d) and Eq. (31) for j = 7 and using Eqs. (23a) and (23b), and after rearranging, we Finally, using Eq. (27f) and Eq. (31) together for j = 9 along with Eqs. (23e) and (23f), we obtain the following: We can use the above three equations to solve for the diagonal parts of the strain rate tensor as follows.
First, we define the following variables: Here, the density gradients ∂ x ρ, ∂ y ρ and ∂ z ρ are based on a (isotropic) second order finite-difference scheme (see the work by Kumar [46] and Leclaire [47] for details). The Eqs.
9 , respectively, which can be obtained from either using raw moments or central moments as follows: 9 = (k 200 + k 020 + k 002 ) − k eq 200 + k eq 020 + k eq 002 = (k 200 + k 020 + k 002 ) − (k eq 200 + k eq 020 + k eq 002 ) , where the equilibrium central moments k eq 200 , k eq 020 and k eq 002 are given in Eq. (10), while the corresponding raw moments are listed in Eq. (11). For example, Based on these considerations, we can rewrite the right and left sides Eqs. (38)- (40), respectively, by means of the following variables: 9 + e 9ρ = (k 200 + k 020 + k 002 ) − 3ρc 2 20 and, Then, Eqs. (38)- (40) reduce to the following system of equations which can be readily solved as follows, thus providing the local expressions for the diagonal components of the strain rate tensor on a cuboid lattice: For completeness, we note that the off-diagonal components of the strain rate tensor can be obtained from Eqs. (27a)-(27c) as Before concluding this section, we provide a guideline in choosing the speed of sound c s for the cuboid lattice-based LB formulations. The optimal value for the cubic lattice c 2 s = 1/3. However, in general cases for the cuboid lattice, the particle speeds in the y and z directions are, respectively, r and s times the particle speed in the x direction. With different particle speeds in different coordinate directions, based on physical considerations, the speed of sound is then chosen as c 2 s = (1/3)min(r 2 , s 2 ). This choice maintains the physically consistent isotropy requirements, reduces the number of spurious terms to be eliminated via the corrections derived earlier, and recovers the optimal value for the cubic lattice.
3. 3D cuboid central moment LBM (3DCCM-LBM) using a non-orthogonal moment basis on the D3Q27 lattice Before constructing a 3D central moment LBM on the cuboid lattice grid using the derivation presented in the previous section, we first note a challenge in directly using the moment basis T given in Eqs. (3) and (4) and the resulting definition of the moments in Eq. (7), and then present its simple resolution that is more suitable for devising an efficient algorithm. In particular, the moment basis T (Eqs. (3) and (4)) involves linear combinations of moments of their second order diagonal components that are written to separate the trace of the diagonal elements from the others in order to allow an independent specification of the bulk and shear viscosities and are also accordingly parameterized by the grid aspect ratios r and s. This was necessary in showing consistency of our approach to the 3D NS equations and in the derivation of the required corrections for using the cuboid lattice. Now, since the collision step to be performed in terms of the relaxations of either the raw moments or the central moments and whose effect need to be transformed back into the velocity space in terms of the distribution functions via T −1 (see e.g., Eq. (8) for the raw moment representation of the LBE), a direct inversion of T would lead to cumbersome expressions involving the grid aspect ratios r and s. However, we will now show that by using a re-defined moment basis, the segregation of the evolution of the moments for independent adjustments of bulk and shear viscosities can be still be achieved by confining it only within the collision step for the relaxation process and not for the mappings, which would have the same overall effect as the original moment basis. Such equivalent but considerably simpler re-formulations of the moment basis and the associated LBE given in Eq. (8) will then form the foundation for the constructions of the 3D LB scheme on the cuboid lattice based on raw moments as well as its generalization in terms of central moments.

Efficient formulation of LBE on the cuboid lattice and its interpretation based on various matrices
First, we define the following moment basis Q involving just the linearly independent bare moments for the D3Q27 cuboid lattice (i.e., without any linear combinations of them as in T -see Eqs. (3) and (4)): Q = |1 , |e x , |e y , |e z , |e x e y , |e x e z , |e y e z , |e 2 x , |e 2 y , |e 2 z , |e x e 2 y , |e x e 2 z , |e 2 x e y , |e y e 2 z , |e 2 x e z , |e 2 y e z , |e x e y e z , |e 2 x e 2 y , |e 2 x e 2 z , |e 2 y e 2 z , |e 2 x e y e z , |e x e 2 y e z , |e x e y e 2 z , |e x e 2 y e 2 z , |e 2 x e y e 2 z , |e 2 x e 2 y e z , |e 2 x e 2 y e 2 z .
where the particle velocities |e x , |e x and |e z appearing in Eq. (50) are given in Eqs. (1a)-(1c) and thus Q depends on the grid aspect ratios r and s. We can then express a correspondence between this re-defined moment basis Q and the original one, i.e, T given Eqs. (3) and (4) using where B represents those operations that combine the various moments according to Eq. (4) and has the form of a block diagonal matrix. While an explicit specification of B is not necessary for the following discussion, Eq. (51) would be helpful in recasting the LBE in an equivalent form that leads to an efficient implementation.
Moreover, let us see how this moment basis Q for the cuboid lattice can be related to that for the cubic lattice 22 that involves the following set of particle velocities so as to construct a modular LB scheme: Such a moment basis for the cubic lattice, referred to as P in the following obviously does not depend on the grid aspect ratios, and can be written as P = |1 , |ē x , |ē y , |ē z , |ē xēy , |ē xēz , |ē yēz , |ē 2 x , |ē 2 y , |ē 2 z , |ē xē It follows from the above definitions that Q and P matrices are related via where S is a simple diagonal scaling matrix, and for the D3Q27 lattice it can be expressed as S = diag 1, 1, r, s, r, s, rs, 1, r 2 , s 2 , r 2 , s 2 , r, rs 2 , s, r 2 s, rs, r 2 , s 2 , r 2 s 2 , rs, r 2 s, rs 2 , r 2 s 2 , rs 2 , r 2 s, r 2 s 2 . (55) Moreover, the inverse mapping for the cuboid lattice Q −1 can be related to that for the cubic lattice P −1 using where S −1 is another diagonal matrix whose elements are the reciprocals of the corresponding elements of S. Thus, As we will see below, these considerations would facilitate the construction of the 3D cuboid LB schemes involving both forward and inverse transformations around collision with effort similar to that for the common cubic lattice that is augmented with some simple scalings of moments based on r and s rather than using lengthy formulas. Next, using the simpler moment basis Q we can relate the vectors of the bare moments (i.e., without involving any linear combinations) and the distribution functions via based on which we now list the countable raw moments for the D3Q27 lattice as follows: m = k 000 , k 100 , k 010 , k 001 , k 110 , k 101 , k 011 , k 200 , k 020 , k 002 , k 120 , k 102 , k 210 , k 012 , k 201 , k 021 , k 111 , k 220 , k 202 , k 022 , k 211 , k 121 , k 112 , k 122 , k 212 , k 221 , k 222 † .
Here, the components of the raw moments, their equilibria and the source moments k mnp and k eq mnp and σ mnp , respectively, are defined in Eq. (6).
With these preliminaries, we now rewrite the LBE given in Eq. (8) and using f = T −1 n in terms of the following collide-and-stream steps with the goal of constructing an efficient LB algorithm on the cuboid f (x + e∆t, t + ∆t) =f (x, t).
Note that here and in what follows, we use 'tilde' (·) to refer to the post-collision state of any variable. Then, we combine Eqs. (7) and (58) As a further simplification, in order to bring this above LB formulation for the cuboid lattice (Eq. (65)) closer to that for the common cubic lattice, we replace the transformations based on Q in terms P via Eq. (54), i.e., of moments only and not for any mappings to or from the velocity space, which are handled by operations similar to those for the cubic lattice augmented with simple scalings of moments based on the grid aspect ratios. These aspects greatly simplify the implementation of the cuboid LB algorithm. Noting that the terms within the brackets [· · · ] in Eq. (66) correspond to the post-collision raw moments, which we write as where the pre-collision raw moments m can be obtained from the distribution functions f via m = SPf , we can finally obtain the following equivalent but significantly more efficient cuboid LB formulation of that given in Eqs. (62)-(63) based on raw moments: Note that when we discuss the implementation details of each of the steps involved here in what follows, we will write down the explicit details of performing P and P −1 related to mappings between the raw moments (which can be readily constructed based on the usual cubic lattice), S and S −1 related to scalings of raw moments by the grid aspect ratios, and B and B −1 related to combining and segregating the moments around the collision steps involving the relaxation process and augmented by the source terms.
As discussed in the introduction earlier, schemes constructed using the central moments have been shown to be more robust than those based on raw moments. Hence, a natural generalization of the above is to consider next in using the central moments in lieu of the raw moments in performing the collision step. In m c,eq = k eq 000 , k eq 100 , k eq 010 , k eq 001 , k eq 110 , k eq 101 , k eq 011 , k eq 200 , k eq 020 , k eq 002 , k eq 120 , k eq 102 , k eq 210 , k eq 012 , k eq 201 , k eq 021 , k eq 111 , k eq 220 , k eq 202 , k eq 022 , k eq 211 , k eq 121 , k eq 112 , k eq 122 , k eq 212 , k eq 221 , k eq Φ c = σ 000 , σ 100 , σ 010 , σ 001 , σ 110 , σ 101 , σ 011 , σ 200 , σ 020 , σ 002 , σ 120 , σ 102 , σ 210 , σ 012 , σ 201 , σ 021 , σ 111 , σ 220 , σ 202 , σ 022 , σ 211 , σ 121 , σ 112 , σ 122 , σ 212 , σ 221 , σ 222 † , where the components of the central moment components k mnp and k eq mnp and σ mnp , respectively, are defined in Eq. (5). We note that the mappings between the central moments and raw moments can be accomplished via the frame transformation matrix F and its inverse F −1 as Here, F and F −1 are both lower triangular matrices dependent on the fluid velocity components (u x , u y , u z ) and are the same for the cuboid lattice as well as the cubic lattice. The elements of F can be readily obtained by listing the binomial transforms of all the linearly independent central moments for the D3Q27 lattice in terms of the raw moments and collecting them in a matrix-vector representation. As noted in our recent work on the rectangular central moment LB scheme [26],once F = F (u x , u y , u z ) is known, its inverse follows from it by exploiting an interesting property of the generating function representation of the binomial expansion, which involves simply reversing the signs of the velocity components appearing in F , Then, as an extension of Eq. (66), the changes under collision as well as the source term in the LBE can be prescribed in terms of central moments, and the resulting post-collision central moments are then mapped back to raw moments via F −1 , which leads to the following equation for the post-collision vector f : In view of Eq. (72), we can write the 3D cuboid central moment LB method (3DCCM-LBM) as a generalization of the previous raw moment formulation given in Eq. (67) as follows: The expressions associated with F and F −1 in addition to the other matrices will be written explicitly when we discuss the implementation details of the central moment LB scheme for the cuboid lattice next.

Algorithmic details of the 3DCCM-LBM
We will now discuss the implementation details of the 3DCCM-LBM based on Eq. (73). While the main steps will be identified systematically in the following, the formulas associated with each of steps will be presented in respective appendices.
• Compute pre-collision raw moments where P is given in Eq. (53), and the components of f , i.e., f α are at time level t, i.e., f α = f α (x, t). This provides intermediate values for all the pre-collision raw moments k mnp for the D3Q27 cubic lattice.
The implementation of this step is given in Appendix A.
• Perform scaling of pre-collision raw moments for cuboid lattice where S is given in Eq. (55) and this step yields the pre-collision raw moments for the cuboid lattice k mnp from those for the cubic lattice computed in the previous step. The equations representing this operation is presented in Appendix B.
• Compute pre-collision central moments where F transforms the pre-collision raw moments k mnp into central moments k mnp , and the associated details are listed in Appendix C.
• Compute post-collision central moments: Relaxation under collision using extended equilibria for corrections and sources for body forcẽ Here • Perform streaming of distribution functions f α (x, t + ∆t) = f α (x − e α ∆t, t).
• Update hydrodynamic fields Using f α (x, t + ∆t) at the new time level t + ∆t from the previous step, fluid density and velocities can be obtained as The following comments regarding the algorithm given above is now in order. (i) When compared to the central moment LB scheme for the common D3Q27 cubic lattice, we emphasize that the only extra computations incurred by the above discussed method for the D3Q27 cuboid lattice are those involving the corrections in the diagonal components of the second order moments in the collision step and the applications of scalings on the raw moments before and after the collision step (i.e., S and S −1 ); these result in a minor additional overhead, but as will be demonstrated in the simulations of an inhomogeneous and anisotropic shear flow case study later, it will be offset by far due to the overall significant computational advantages of using the cuboid lattice in lieu of the cubic lattice that accommodates the nature of the flow much more efficiently.
(ii) The algorithm given above based on central moments is modular in nature. It can be readily extended to other collision models. For example, a raw moment-based cuboid scheme as a modification of the above would simply involve the elimination of the mappings between the raw moments and central moments (i.e., those involving F and F −1 ) around collision and performing the collision step in terms of raw moments rather than central moments, but using the same corrections via the extended moment equilibria. However, as will be numerically shown later that the central moment formulation is more robust than that based on raw moments for the cuboid lattice. Schemes based on other collision models such as those based on cumulants [28] can be readily constructed from the above algorithm by introducing additional mappings between central moments and cumulants and performing the collision step in terms of relaxations of cumulants by using the same required corrections for the cuboid lattice as in the above algorithm. Such an extension can be considered in a future study. (iii) Finally, the present approach devised for the D3Q27 cuboid lattice readily extends for other lattice sets, such as D3Q15 and D3Q27 cuboid lattices, by using only a subset of the moment basis as necessary while incorporating the corrections to the equilibria derived in this work.

Boundary conditions on moving walls: Momentum-augmented bounce back scheme on cuboid lattice grids
Before discussing the numerical results based on the 3DCCM-LBM constructed above, we note that the use of cuboid lattice requires some modifications to the implementation of the standard link-based half-way bounce back boundary condition for moving boundaries. Since we will be investigating shear driven flows due to the motion of the walls in this work, we will now derive the necessary changes to such a boundary scheme.
If x f and x w represent the fluid node nearest to the boundary and the wall node, respectively, and α andᾱ denote, respectively, the outgoing and incoming particle directions, where eᾱ = −e α , then the general form of the momentum-augmented half-way bounce back scheme can be written as Here, the equilibrium distribution functions f eq α and f eq α associated with f eq at x w depend on the wall density and velocities. It can be obtained via an inverse mapping from the equilibrium moments m eq using f eq = Q −1 m eq , where the components of the equilibrium moments are evaluated using Eq. (11) from those based on the imposed wall conditions. Note that since Q −1 is dependent on the grid aspect ratios, the resulting formulas will be parameterized by r and s. Specifically, consider a moving boundary located in the x − z plane, whose outward normal is along the y direction and pointing into the plane of the paper (see Fig. 1). Let this boundary plane be moving along the x direction with an imposed velocity U ; thus based on the no-slip boundary conditions for the fluid, we can then write u x = U , u y = 0, u z = 0, and, as usual for the bounce back scheme, we approximate the wall node density to that based on the known fluid, i.e., ρ = ρ f .
For the case under consideration, referring to Fig. 1, the unknown or incoming distribution functions are fᾱ, whereᾱ = 4, 9,10,16,18,21,22,25,26, and can be written in terms of the known directions using Eq. (74) as Then, we determine the expressions of the known equilibrium distribution functions using f eq = Q −1 m eq , where we evaluate the components of the moment equilibria on the wall via Eq. (11) based on the imposed wall conditions mentioned above. Using these, we then finally obtain the momentum-augmented bounce-back scheme for the cuboid lattice as Clearly, the formulas given in Eq. (75), which will be utilized in the next section, depend on the grid aspect ratios r and s in addition to the specified wall velocity U .

Results and discussion: Numerical validation
In this section, we will perform numerical validations of our new 3D cuboid central moment LBM via simulations of an assortment of canonical fluid flow problems including flows in square duct, pulsatile flow in a square duct driven by a periodic body force, and lid-driven flow within a cubic cavity at various characteristic parameters and grid aspect ratios.

Flow through a square duct driven by a constant body force
First, we perform simulations of a fully developed flow through a square duct driven by a constant body force using the 3DCCM-LBM. The analytical solution of this flow problem for the velocity field u(y, z) can be written as [48] u(y, z) = 4L 2 F x π 3 ρν ∞ n=1,3,5,...
where L represents the side of the square duct, ν is the kinematic viscosity and F x is the body force imposed in the direction of the fluid flow, the x-direction. We resolve the computational domain using N x × N y × N z grid nodes and the grid aspect ratios in the y and z directions, i.e., r and s, respectively, are then calculated using r = L/N y and s = L/N z . Periodic boundary condition is applied along the flow direction (x-direction) and no slip boundary conditions at the four side walls are imposed using the half-way bounce back boundary conditions (see e.g., Eqs. (75)). We define the Reynolds number Re for this flow based on the maximum flow velocity occurring at the midpoint within the duct cross section and the duct side length as the characteristic velocity and length scales, respectively. Table 1  Excellent agreement can be seen between the computed results using the cuboid lattice and the analytical solution. Moreover, though not shown in this figure, the 3DCCM-LBM results were also found to be in similar agreements with the exact solution for all the other cases of the grid aspect ratios tested according to Table 1. We note here that due to the similarity of the governing equations for this Stokes flow problem with the anisotropic advection diffusion equation (ADE), it can be alternatively solved by the LBMs developed in this regard by Ginzburg. [49,50] via a coordinate grid transform, as in e.g., [51,52], for which the role of the associated anisotropy on the numerical stability was studied in [53]. More generally, the idea of the coordinate grid transform used for the ADE LBM may also be adopted for flow LBM.

Pulsatile flow in a square duct driven by a periodic body force
Next, we will assess the ability of the 3DCCM-LBM in accurately representing a time-dependent flow problem for which an analytical solution is available. In this regard, we simulate time-periodic flow through a square duct of side length 2a driven by a sinusoidally varying body force F x = F m cos ωt, where F m is its peak amplitude and ω = 2π/T is the angular frequency of the applied pulsations and with T being its time period. This case study admits an analytical solution for the spatially and temporally varying velocity field based on Fourier series and given by [54] u(y, (−1) n p n cosh(γ n y/a) cos(p n z/a) cosh(γ n ) + cosh(γ n z/a) cos(p n y/a) cosh(γ n ) e iωt ,(77) γ n = p 2 n + iWo 2 , p n = (2n + 1)π/2.
where −a ≤ y, z ≤ a and R[· · · ] represents the real part of the terms within the brackets. Here, Wo is  involves circulation patterns with 3D vortical structures whose details depend on the Reynolds number Re defined as Re = U H/ν [55]. In this paper, we will perform simulation of this flow at Re = 100, 400 and 1000 for which the reference numerical solutions for benchmarking can be found from various prior studies (see e.g. [56,57,58]). For applying the 3DCCM-LBM in this regard, we used the following three different grid aspect ratios r = {1.0, 0.5, 0.33} and s = 1, so that the resolutions are selectively varied only along the y direction which is normal to the shearing motion of the top lid. The choices made for the various parameters for each set of grid aspect ratios, including the number of grid nodes, lid velocity, speed of sound, and the relaxation time τ = 1.0/ω ν are presented in Table 2.  Table 2 earlier. Here, we compared them with the benchmark numerical solution based on the method of differential quadrature to solve the incompressible NS equation in Shu et al [58]. It is evident that our 3D central moment LBM based on the cuboid lattice is in very good agreement with the reference numerical data. Moreover, In addition, for visualizing the overall flow patterns, the streamlines along the midx − y plane, x − z plane and y − z plane of the cubic cavity at Re = 100, 400 and 1000 computed using r = 0.5 and s = 1 are presented in Fig. 8. The motion of the lid normal to the y direction is seen to set up a major vortex in x − y plane whose size and its center location appear to significantly change as the Reynolds number is varied. In particular, the eye of this vortex moves towards the center in the x − y plane as Re increases. Meanwhile, some secondary vortices start to form around the corners at higher Re and grow in size with increase in the inertial effects as the Reynolds number is increased. The three-dimensionality of the resulting flow field is evident from observing the other y − z and x − z planes in these figures. These features are in agreement with the prior numerical solutions [56,57,58]. (d,e,f) and Re=1000 (g,h,i) at z = 0.5H, y = 0.5H, and x = 0.5H, respectively.
6. Demonstration of computational advantages of using cuboid lattice over cubic lattice: 3D anisotropic shear flows in a lid-driven shallow cuboid cavity The results from the above case studies show the ability of our 3DCCM-LBM in computing the flow fields accurately for a variety of standard benchmark flow problems. Now, we will demonstrate the advantages of using the cuboid lattice over the cubic lattice in more efficiently simulating a flow case study characterized by an anisotropic and inhomogeneous shear flow where their spatial gradients in one or more direction are larger than those in the other directions. In this regard, we consider a 3D cuboid cavity of span length L in the x direction, height H in the y direction, and width W in the z direction enclosing a fluid of viscosity ν.
A flow is set up by the shearing motion of the lid located in the x − z plane at y = H at a velocity U (see  For performing flow simulations, we choose c 2 s = 1/3 for the cubic lattice and c 2 s = 0.04 for the cuboid lattice. Figure 10 presents comparisons of the velocity profiles of the u component along the y direction at x = 0.5L and z = 0.5W , and v component along the x direction at y = 0.5H and z = 0.5W obtained using the 3DCCM-LBM with the above three sets of grid aspect ratios for the cuboid lattice and the cubic lattice with (r, s) = (1.0, 1.0). Moreover, the streamlines along the three mid-planes of the 3D shallow cuboid cavity using both the cubic lattice and the cuboid lattice using the grid aspect ratios for the case (i) (corresponding to that with the lowest total number of grid points among the three cuboid lattice cases) are shown in Fig. 11. Clearly, the cuboid lattice results, while using significantly fewer total number of grid nodes, are in excellent agreement with those using the cubic lattice. In particular, to achieve practically similar accuracy for the numerical results, the flexibility accorded by the cuboid lattice yielded significant savings in the computer memory storage by a factor of 7.88, 5.33 and 7.11 for cases (i), (ii), and (iii), respectively, and about similar reductions in the computational cost for the simulation turnaround time when compared to the cubic lattice. While for the specific flow configuration considered here a reduction in the computational resource requirements by a factor between 5 to 7.5 was achieved for the choices made for the grid aspect ratios, if the flow geometry happens to be further skewed in the different coordinate directions, e.g., when H/L 1 and/or H/W 1, additional improvements with using the cuboid lattice can be expected for simulating such anisotropic shear flow problems. respectively. In both cases, we use c 2 s = 0.08, which is the minimum of the optimal values for the two grid aspect ratios, in order to have the same reference speed of sound in calculating the Mach number. All the relaxation parameters, other than that for the shear viscosity, are set to 1.0 for simplicity. The results are shown in Fig. 12. Clearly, the 3DCCM-LBM is able to support larger magnitudes of shear consistently than the 3DCRM-LBM for all the choices of the grid resolution used. Moreover, it was pointed out earlier that in our 3D cuboid LB algorithms, the trace of the diagonal components of the second order moments is evolved separately from the other moments with its own relaxation parameter ω ξ after applying B, and the individual components are subsequently recovered via applying B −1 , so that the bulk viscosity can be specified independently from shear viscosity. It is known that for the common LB schemes based on the cubic lattice sets that a higher value of the bulk viscosity (via decreasing ω ξ ) enhances stability by suppressing the spurious pressure waves. We will now verify if this also holds true for the cuboid LB formulation based on central moments. Repeating the stability test mentioned above by choosing ω ξ = 0.5, 1.0, and 1.4 at 30 × 60 × 30 and c 2 s = 0.08, Fig. 13 shows the resulting maximum threshold velocity of the lid U for numerical stability. This figure clarifies that the robustness of the 3DCCM-LBM is significantly improved by choosing lower ω ξ or higher bulk viscosity.

Summary and Conclusions
In this paper, we have presented a new 3D LB algorithm based on central moments for the D3Q27 lattice using a cuboid grid, which is parameterized by two grid aspect ratios that are related to the ratios of the particle speeds with respect to that along a reference coordinate direction. The additional degrees of freedom introduced by the flexible specification of the cuboid lattice grid enables the method to naturally and more efficiently compute flows having different characteristic length scales in different directions. It is constructed to simulate the Navier-Stokes equations consistently via introducing counteracting corrections to the second order moment equilibria that eliminate the associated grid anisotropy and the non-Galilean invariant third order velocity errors obtained via a Chapman-Enskog analysis. Such corrections are related to the diagonal components of the velocity gradient tensor, which are computed locally using non-equilibrium moments in our approach and depend on the grid aspect ratios. The implementation is shown to be compact and modular, with an interpretation based on special matrices, admitting ready extension of the standard algorithm for the cubic lattice to the cuboid lattice via appropriate scaling of moments based on grid aspect ratios before and after collision step and equilibria corrections. The resulting formulation is general in that the same grid corrections developed for the D3Q27 lattice for recovering the correct viscous stress tensor is applicable for other lattice subsets, and a variety of collision models, including those based on the relaxation of raw moments, central moments and cumulants, as well as their special case involving the distribution functions. The cuboid central moment LBM was validated against a variety of benchmark flows, and when used in lieu of the corresponding raw moment formulation for simulating shear flows, we show that it results in significant improvements in numerical stability. Finally, we demonstrated that our cuboid LB approach is efficient in simulating anisotropic shear flow problems with significant savings in computational cost and 43 memory storage when compared to that based on the cubic lattice. Future improvements are expected by extending the present formulation as a multiblock cuboid LB algorithm, which will be considered in a future study.