The lattice Boltzmann (LB) method [1]}, [2]}, [3]} has been receiving a remarkable interest as a promising computational fluid dynamics (CFD) technique. It is a kinetic method that evolves the distribution functions due to the effects of collisions, which are often represented by relaxation model under certain symmetry, isotropy and conservation constraints (e.g., [4]}, [5]}, [6]}, [7]}) and due to their streaming along the particle characteristic directions. The features and applications of this mesoscopic computational technique have been discussed in various reviews [8]}, [9]}, [10]}, [11]}, [12]} and monographs [13]}, [14]}, [15]}. For simulating inhomogeneous fluid motion, such as those involving boundary layer flows with different scales for variations in different coordinate directions, or flows in domains where the spatial extent of one of the directions is considerably shorter than the others, i.e., characterized by geometric anisotropy such as in sheetlike porous media, it is highly desirable to use nonuniform grids for enabling efficient simulations. However, the symmetry constraints and the coupling of the particle velocity and coordinate space discretizations restrict the use of uniform grids, e.g., square lattice in two-dimensions (2D). To address this issue, broadly, the following two types of modifications to the LBM have been considered based on (a) the decoupling the discretizations of the velocity space and the spatial coordinate space and (b) the rectangular lattice grid to naturally accommodate the inhomogeneity in flows. In the first category, the resulting LBM does not maintain the lock-step advection during the streaming step, and either involves interpolations (e.g., [16]}) or the use of traditional discretizations such as finite volume or finite element schemes, which leads to a more complicated approach with attendant additional overhead (e.g., [17]}, [18]}, [19]}, [20]}). The second category maintains the perfect-shift advection during the streaming step that incurs relatively low numerical dissipation, an important numerical advantage, and is the focus of this work. However, to recover the inherent isotropy of the viscous stress tensor in LB simulations using such rectangular lattices require making certain modifications to the algorithm. Thus, the prior LB schemes on rectangular lattice grids, which was inspired from an early work [21]}, can be further classified according to the following modifications made: (i) designing the collision step with sufficient degrees of freedom and parametrization of the relaxation rates to maintain isotropy [22]}, [23]}, [24]}, (ii) extend the lattice with additional particle velocities [25]} (iii) the use of extended moment equilibria to correct for isotropy [26]}, [27]}, and (iv) use of coordinate and velocity transformations and counteracting source terms [28]}. Categorizing from a different consideration, such rectangular lattice-based LB algorithms use either the single relaxation time (SRT) collision model [21]}, [25]}, [26]}, [28]} or the multiple relaxation time (MRT) collision operator [22]}, [23]}, [24]}, [27]}, involving the relaxation of the distribution functions or the raw moments, respectively.
