Local Vorticity Computation in Double Distribution Functions based Lattice Boltzmann Methods for Flow and Scalar Transport

TL;DR

This paper introduces a novel local vorticity computation method within double distribution function lattice Boltzmann frameworks, enabling accurate, fully local flow analysis suitable for parallel processing, demonstrated through benchmark flow simulations.

Contribution

It presents a new local vorticity computation approach using off-diagonal moments in DDF-LBMs, avoiding finite difference derivatives and supporting parallel computation.

Findings

Accurate vorticity fields with second order convergence.

Method is fully local and suitable for parallel implementation.

Demonstrated effectiveness on benchmark flow problems.

Abstract

Computation of vorticity in conjunction with the strain rate tensor, plays an important role in fluid mechanics in vortical structure identification and in the modeling of various complex fluids. For the simulation of flows accompanied by the advection-diffusion transport of a scalar field, double distribution functions (DDF) based lattice Boltzmann methods (LBMs) are commonly used. We present a new local vorticity computation approach by introducing an intensional anisotropy of the scalar flux in the third order, off-diagonal moment equilibria of the LB scheme for the scalar field, and then combining the second order non-equilibrium components of both the LBMs. As such, any pair of lattice sets in the DDF formulation that can independently support the third order off-diagonal moments would enable local determination of the complete flow kinematics, with the LBMs for the fluid motion and the transport of the passive scalar respectively providing the necessary moment relationships to determine the symmetric and skew-symmetric components of the velocity gradient tensor. Since the resulting formulation is completely local and do not rely on finite difference approximations for velocity derivatives, it is by design naturally suitable for parallel computation. As an illustration of our approach, we formulate a DDF-LB scheme for local vorticity computation using a pair of multiple relaxation times (MRT) based collision approaches on two-dimensional, nine velocity (D2Q9) lattices, where the necessary moment relationships to determine the velocity gradient tensor and the vorticity are established via a Chapman-Enskog analysis. Simulations of various benchmark flows demonstrate good accuracy of the predicted vorticity fields, with a second order convergence. Furthermore, extensions of our formulation for a variety of collision models to enable local vorticity computation are presented.

Equations

$${\partial }_j{u}_i=\frac{1}{2}({\partial }_j{u}_i+{\partial }_i{u}_j)+\frac{1}{2}({\partial }_j{u}_i-{\partial }_i{u}_j)=S_{ij}+{\Omega }_{ij},$$

$${\partial }_t\rho +\nabla \cdot (\rho \mathbf{u})$$

$${\partial }_t(\rho \mathbf{u})+\nabla \cdot (\rho \mathbf{u}\mathbf{u})$$

$${\partial }_t\varphi +\nabla \cdot (\varphi \mathbf{u})$$

$$|e_x⟩={\left(0,1,0,-1,0,1,-1,-1,1\right)}^{†},$$

$$|e_y⟩={\left(0,0,1,0,-1,1,1,-1,-1\right)}^{†},$$

$$|1⟩={\left(1,1,1,1,1,1,1,1,1\right)}^{†}.$$

$$T_0=|1⟩,T_1=|e_x⟩,T_2=|e_y⟩,T_3=|e_x^2+e_y^2⟩,T_4=|e_x^2-e_y^2⟩,$$

$$T_5=|e_x{e}_y⟩,T_6=|e_x^2{e}_y⟩,T_7=|e_x{e}_y^2⟩,T_8=|e_x^2{e}_y^2⟩.$$

$$\mathsf{T}=\left[T_0,T_1,T_2,T_3,T_4,T_5,T_6,T_7,T_8\right].$$

$$\displaystyle\left({\begin{array}[]{*{20}{l}}{{{\hat{\kappa}}_{{x^{m}}{y^{n}}}}}^{{}^{\prime}}\\ {\hat{\kappa}_{{x^{m}}{y^{n}}}^{eq^{\prime}}}\\ {\hat{\sigma}_{{x^{m}}{y^{n}}}^{eq^{\prime}}}\end{array}}\right)=\sum\limits_{\alpha=0}^{8}{\left({\begin{array}[]{*{20}{l}}{{f_{\alpha}}}\\ {f_{\alpha}^{eq}}\\ {S_{\alpha}}\end{array}}\right)}{e_{\alpha x}^{m}}{e_{\alpha y}^{n}}.$$

$$\hat{\mathbf{m}}=\mathsf{T}\mathbf{f},{\hat{\mathbf{m}}}^{eq}=\mathsf{T}{\mathbf{f}}^{eq},\hat{\mathbf{S}}=\mathsf{T}\mathbf{S},$$

$$\mathbf{f}$$

$$\mathbf{S}$$

$$\hat{\mathbf{m}}$$

$${\hat{\mathbf{m}}}^{eq}$$

$$\hat{\mathbf{S}}$$

$$\mathbf{f}\left(\mathbf{x}+{\mathbf{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)-\mathbf{f}\left(\mathbf{x},t\right)$$

$$\hat{\mathsf{\Lambda }}=\text{mbox}diag(0,0,0,{\omega }_3,{\omega }_4,{\omega }_5,{\omega }_6,{\omega }_7,{\omega }_8).$$

$$\bar{\mathbf{f}}\left(\mathbf{x}+{\mathbf{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)-\bar{\mathbf{f}}\left(\mathbf{x},t\right)={\mathsf{T}}^{-1}\left[-\hat{\mathsf{\Lambda }}\left(\widehat{\bar{\mathbf{m}}}-{\hat{\mathbf{m}}}^{eq}\right)\right]+{\mathsf{T}}^{-1}\left[\left(\mathsf{I}-\frac{1}{2}\hat{\mathsf{\Lambda }}\right)\hat{\mathbf{S}}\right]{\delta }_t,$$

$${\hat{\kappa }}_0^{e{q}^{\prime }}=\rho ,{\hat{\kappa }}_x^{e{q}^{\prime }}=\rho u_x,{\hat{\kappa }}_y^{e{q}^{\prime }}=\rho u_y,$$

$${\hat{\kappa }}_{xx}^{e{q}^{\prime }}=c_s^2\rho +\rho u_x^2,{\hat{\kappa }}_{yy}^{e{q}^{\prime }}=c_s^2\rho +\rho u_y^2,{\hat{\kappa }}_{xy}^{e{q}^{\prime }}=\rho u_x{u}_y,$$

$${\hat{\kappa }}_{xxy}^{e{q}^{\prime }}=c_s^2\rho u_y+\rho u_x^2{u}_y,{\hat{\kappa }}_{xyy}^{e{q}^{\prime }}=c_s^2\rho u_x+\rho u_x{u}_y^2,$$

$${\hat{\kappa }}_{xxyy}^{e{q}^{\prime }}=c_s^4\rho +c_s^2\rho (u_x^2+u_y^2)+\rho u_x^2{u}_y^2,$$

$${\hat{\sigma }}_0^{{}^{\prime }}=0,{\hat{\sigma }}_x^{{}^{\prime }}=F_x,{\hat{\sigma }}_y^{{}^{\prime }}=F_y,$$

$${\hat{\sigma }}_{xx}^{{}^{\prime }}=2F_x{u}_x,{\hat{\sigma }}_{yy}^{{}^{\prime }}=2F_y{u}_y,{\hat{\sigma }}_{xy}^{{}^{\prime }}=F_x{u}_y+F_y{u}_x,$$

$${\hat{\sigma }}_{xxy}^{{}^{\prime }}=F_y{u}_x^2+2F_x{u}_x{u}_y,{\hat{\sigma }}_{xyy}^{{}^{\prime }}=F_x{u}_y^2+2F_y{u}_y{u}_x,$$

$${\hat{\sigma }}_{xxyy}^{{}^{\prime }}=2(F_x{u}_x{u}_y^2+F_y{u}_y{u}_x^2),$$

$$\rho =\sum _{\alpha =0}^8{\bar{f}}_{\alpha },\rho \mathbf{u}=\sum _{\alpha =0}^8{\bar{f}}_{\alpha }{\mathbf{e}}_{\alpha }+\frac{1}{2}\mathbf{F}{\delta }_t,p=c_s^2\rho ,$$

$$\hat{\mathbf{m}}=\sum _{j=0}^{\infty }{ϵ}^j{\hat{\mathbf{m}}}^{(j)},{\partial }_t=\sum _{j=0}^{\infty }{ϵ}^j{\partial }_{t_j},$$

$$\bar{\mathbf{f}}(\mathbf{x}+{\mathbf{e}}_{\alpha }ϵ,t+ϵ)=\sum _{j=0}^{\infty }\frac{{ϵ}^j}{j!}({\partial }_t+{\mathbf{e}}_{\alpha }\cdot \nabla {)}^j\bar{\mathbf{f}}(\mathbf{x},t).$$

$$O({ϵ}^0):{\hat{\mathbf{m}}}^{(0)}={\hat{\mathbf{m}}}^{eq},$$

$$O({ϵ}^1):({\partial }_{t_0}+{\hat{\mathsf{E}}}_i{\partial }_i){\hat{\mathbf{m}}}^{(0)}=-\hat{\mathsf{\Lambda }}{\hat{\mathbf{m}}}^{(1)}+\hat{\mathbf{S}},$$

$$O({ϵ}^2):{\partial }_{t_1}{\hat{\mathbf{m}}}^{(0)}+({\partial }_{t_0}+{\hat{\mathsf{E}}}_i{\partial }_i)\left[\mathsf{I}-\frac{1}{2}\hat{\mathsf{\Lambda }}\right]{\hat{\mathbf{m}}}^{(1)}=-\hat{\mathsf{\Lambda }}{\hat{\mathbf{m}}}^{(2)},$$

$${\partial }_{t_0}\rho +{\partial }_x(\rho u_x)+{\partial }_y(\rho u_y)=0,$$

$${\partial }_{t_0}\left(\rho u_x\right)+{\partial }_x\left(c_s^2\rho +\rho u_x^2\right)+{\partial }_y\left(\rho u_x{u}_y\right)=F_x,$$

$${\partial }_{t_0}\left(\rho u_y\right)+{\partial }_x\left(\rho u_x{u}_y\right)+{\partial }_y\left(c_s^2\rho +\rho u_y^2\right)=F_y,$$

$${\partial }_{t_0}\left(2c_s^2\rho +\rho (u_x^2+u_y^2)\right)+{\partial }_x\left[(1+c_s^2)\rho u_x+\rho u_x{u}_y^2\right]+{\partial }_y\left[(1+c_s^2)\rho u_y+\rho u_x^2{u}_y\right]$$

$$=-{\omega }_3{\hat{m}}_3^{(1)}+2(F_x{u}_x+F_y{u}_y),$$

$${\partial }_{t_0}\left(\rho (u_x^2-u_y^2)\right)+{\partial }_x\left[(1-c_s^2)\rho u_x-\rho u_x{u}_y^2\right]+{\partial }_y\left[(-1+c_s^2)\rho u_y+\rho u_x^2{u}_y\right]$$

Full text

1 Introduction

Qualitative distribution and quantitative measures of vorticity is of fundamental interest in fluid mechanics. Indeed, fluid motions are often associated with vortical structures, which can be characterized by vorticity, and, more generally, by certain invariants of the velocity gradient tensor [1, 2]. The significance of the rigid-body like rotational component of the fluid element was first identified in a pioneering work by Helmholtz [3] and the subject has a long and rich history [4, 5]. This local rotational property of the flow, given by the curl of the velocity field, was termed vorticity by Lamb [6]. While there is no consensus on a rigorous definition of a vortex, various quantitative measures have been devised to identify regions associated with more rigid-body like rotations than stretching or shearing motions that aid in flow classification [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Such approaches are based on a complete knowledge of the velocity gradient tensor, and the local, Eulerian based methods for coherent structure identification are popular (see [20] for recent review).

In more detail, the velocity gradient tensor $A_{ij}\equiv\partial_{j}u_{i}$ of the velocity field $u_{i}$ can be decomposed into symmetric $S_{ij}$ and anti- or skew-symmetric parts $\Omega_{ij}$ as

[TABLE]

where $S_{ij}$ is the strain rate tensor and $\Omega_{ij}$ is the intrinsic rotation rate (spin) tensor, with $\Omega_{ij}=-\frac{1}{2}\epsilon_{ijk}\omega_{k}$. Here, $\omega_{k}$ is the Cartesian component of the vorticity and $\epsilon_{ijk}$ is the Levi-Civita (permutation) tensor, and the vorticity can be defined as $\omega_{i}=\epsilon_{ijk}\partial_{j}u_{k}$ or $\bm{\omega}=\bm{\nabla}\times\bm{u}$. Both $\omega_{i}$ and $S_{ij}$, or, in general, $\partial_{j}u_{i}$ play an important role in eduction techniques for vortex structure identification. In particular, many of these methods [20] are based on the second and third invariants of the velocity gradient tensor $\partial_{j}u_{i}$, i.e., $Q=-\frac{1}{2}S_{ij}S_{ij}+\frac{1}{4}\omega_{k}\omega_{k}$ and $R=\frac{1}{3}(S_{ij}S_{jk}S_{ki}+\frac{3}{4}\omega_{i}\omega_{j}S_{ij})$. Similarly, sometimes the Lamb vector $L_{i}=\epsilon_{ijk}\omega_{j}u_{k}$ plays a prominent role in the analysis of vortex dynamics [21]. Thus, a complete knowledge of the local velocity gradient tensor $\partial_{j}u_{i}$, or equivalently, $S_{ij}$ and $\Omega_{ij}$ or $\omega_{k}$ is of basic interest in structure identification and classification of flows. This also allows a local determination of the components of the convective acceleration of the fluid elements. In addition, the distribution of vorticity is related to the sound generation and propagation in flow generated acoustics [22]. Furthermore, many models for the representation of turbulence (e.g., [23]), rheological fluid flows such as those involving viscoelasticity, and complex fluid systems such as liquid crystals and polar fluids depend on the local measures of the complete velocity gradient tensor $\partial_{j}u_{i}$ [24, 25, 26, 27]. Moreover, molecular liquid flows under nanoscale confinement involves the relaxation of the intrinsic angular momentum to the vorticity of the fluid element, and its coupling to the linear momentum, which needs to be modeled [28, 29, 30, 31]. It is thus highly desirable for computational methods for fluid dynamics that allow especially local determination of all components of the velocity gradient tensor, including the skew symmetric part (i.e., the vorticity). Here, we emphasize that ‘local’ implies that such methods do not depend on finite difference approximations for velocity derivatives, but are entirely based on operations of suitable quantities available at the grid nodes, and hence are naturally suitable for parallel computing.

The lattice Boltzmann method (LBM) is a kinetic computational approach for a variety of fluid mechanics and transport problems  [32, 33, 34, 35, 36, 37, 38, 39]. Generally, the standard versions of the LB schemes can only represent the symmetric part of $\partial_{j}u_{i}$, i.e., the strain rate tensor $S_{ij}$ based on local algorithms via the second order non-equilibrium moments of the distribution function, which are, in turn, related to the spatial derivatives of the first and third order moment equilibria. The latter are constructed based on symmetry and isotropy considerations that respect the underlying isotropy of the viscous stress tensor of the fluid motion represented by the Navier-Stokes equations. It is known that such LB approaches can recover the strain rate tensor components locally with second order accuracy (see e.g.,  [40, 41, 42]). However, most of the existing LBMs are not constructed to recover the antisymmetric velocity gradient tensor $\Omega_{ij}$ locally and need to rely on finite difference computations. One notable exception is the recent and interesting work [43], which introduced an approach based on modifying the fifth order moment equilibria of the LB solver for fluid flow that enables vorticity computation. This approach is restricted to only lattices that can support fifth order independent moments and thus is applicable only to the three-dimensional, twenty seven velocity (D3Q27) lattice, and not for other standard lattice sets, including the common two-dimensional, nine velocity (D2Q9) lattice, and D3Q15 and D3Q19 lattices in 3D. Furthermore, since it is based on certain prescribed form of the higher order moment equilibria, it may be challenging to extend such LB scheme for thermal flows as well as those with significant compressibility effects that involve constraints on the higher moments of the single distribution function, and may also impact its Galilean invariance of solving the fluid motion. Also, since it involves combining second and fourth order non-equilibrium moments, which may be subjected to hyperviscosity effects [38], the attendant higher order moment equilibria for the solution of the fluid motion need to be constructed carefully.

Our approach is based on different considerations than the above mentioned work for vorticity computation. When the goal is to simulate the fluid motion along with an advection-diffusion transport of a scalar field, represented by the following Navier-Stokes equations (NSE) and the convection-diffusion equation (CDE), respectively:

[TABLE]

where $\rho$, $\bm{u}$ and $p$ are the fluid density, velocity, and pressure, respectively, $T_{ij}=2\rho\nu S_{ij}-\frac{2}{d}\rho\nu\partial_{k}u_{k}\delta_{ij}+\rho\zeta\partial_{k}u_{k}\delta_{ij}$ is the deviatoric stress tensor (with $\nu$ and $\zeta$ being the kinematic shear and bulk viscosities, respectively, and $d$ being the number of spatial dimensions), $\bm{F}$ is the local body force, and $\phi$ is the scalar field (with $D_{\phi}$ being its diffusivity), they can be solved by means of a double distribution functions (DDF) based approach using two LB schemes – one for the flow field and the other for the scalar field. Such situations related to solving the additional passive scalar field dynamics arise widely, including those related to the transport of energy or temperature field in thermal convection, and of the concentration field of a chemical species in reacting systems, as well as in the interface capturing using phase field models in multiphase flows. Indeed, the modeling of flow and scalar transport using DDF based LBEs is quite common and is a subject of a number of investigations (e.g., [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56]). In such cases, our essential philosophy is to use the additional degrees of freedom (DOF) available in the LBE for the solution of the CDE to construct a procedure for local vorticity computation [57]. This is possible because as the evolution of the scalar field $\phi$ is influenced by the local fluid velocity $\bm{u}$, its solution procedure can, in principle, contain the complete kinematics of the flow field, which can be obtained from the corresponding LBM with careful construction of its equilibria.

The basic idea behind our approach is as follows. Local vorticity computation in the DDF-LB schemes can be achieved by prescribing an intensional anisotropy of the scalar flux in the third order, off-diagonal moment equilibria of the LBM for the scalar field and then combining the second order, off-diagonal non-equilibrium moment components of both the LBMs. In essence, the LBM for the fluid flow provides local expressions for the strain rate tensor $S_{ij}$ and the LBM for the scalar field yields local relations for the skew-symmetric velocity gradient tensor $\Omega_{ij}$, and hence the vorticity $\omega_{k}$. This formulation leads to various advantages. The numerical characteristics of the LBM for the fluid motion are preserved as no additional modifications in terms of constraints on its equilibria are imposed (but only on those for the scalar field) and the resulting approach is non-invasive in representing the fluid flow. The freedom from the need to prescribing extra constraints for higher moments for the LB flow solver allows ready extension to construct LB schemes for complex flow physics. In addition, any pair of lattice sets, each supporting only lower (i.e., third) order independent moments, in this DDF-LBE approach can enable local vorticity computation. Thus, this method is applicable for all standard lattices (e.g., D2Q9, D3Q15, D3Q19 and D3Q27) and in different dimensions. Furthermore, the since method is based on two distribution functions which by themselves are generally solved with second order accuracy, the numerically predicted vorticity magnitudes are second order by construction, just like the computed strain rate tensor. Moreover, the local expressions for the vorticity field, which are not dependent on finite difference approximations of the velocity field, naturally lend themselves to parallel computation. Finally, it can be used to establish a LB framework for fully local modeling and computation of complex fluids (e.g., viscoelastic or polar fluids), which generally depend on both the symmetric and skew-symmetric velocity gradient tensor and are usually represented by the evolution of additional distribution functions to represent the attendant multiphysics effects.

For the purpose of illustration without losing generality, in this work, we will specialize our DDF approach by formulating in detail two LB schemes using natural (non-orthogonal) moment basis and multiple relaxation times (MRT) for the solution of flow and scalar transport using the standard D2Q9 lattice to locally compute the complete information about the flow kinematics, including the skew-symmetric velocity gradient tensor components. However, our method can be readily extended to LBM based on other collision models and various other lattice sets in different dimensions. For completeness, we will also present the extension of our approach for the single relaxation time (SRT) LBM and the cascaded central moment LB schemes for the D2Q9 lattice in the appendices. While the objective of this paper is to formulate, mathematically analyze and numerical validate our new LB approach in 2D, its extension to 3D lattices will be presented in a follow-up investigation.

This paper is organized as follows. The next section (Sec. 2) will present a MRT-LBM for computing the fluid motion, and its Chapman-Enskog (C-E) analysis to determine the symmetric components of the velocity gradient tensor. Section 3 will then discuss another MRT-LBM for representing the advection-diffusion transport of a scalar field with the required modifications as indicated earlier, and its C-E analysis to obtain the necessary relations for the skew-symmetric components of the velocity gradient tensor. The expression for the local computation of the vorticity field is derived in Sec. 4. Then, results and discussion of the comparisons of the computed vorticity fields against the analytical/numerical solutions for various representative fluid flow problems are given in Sec. 5. Finally, Sec. 6 presents a summary and conclusions of this work. In addition, A presents the system of non-equilibrium moments and spatial derivatives of various attendant components of moment equilibria needed in the derivation of our approach. B discusses a formulation to recover the skew-symmetric velocity gradient tensor for the SRT-LBM, while C and D present extensions of our idea for different versions of the LBM based on central moments.

2 MRT-LBM for Fluid Motion

In order to solve the fluid motion in two-dimensions (2D) represented by the mass and momentum conservation equations given in Eqs. (2) and (3), respectively, we will now present a MRT-LBM using a natural, non-orthogonal moment basis [58]. In this regard, a D2Q9 lattice is used, and whose particle velocities are given by the following:

[TABLE]

where ${\dagger}$ is the transpose operator and the standard Dirac’s bra-ket notation is used to represent the vectors. The Cartesian components for any particle direction $\alpha$ are represented by $e_{\alpha x}$ and $e_{\alpha y}$, where $\alpha=0,1,\ldots,8$. In addition, we need the following 9-dimensional vector whose inner product with the particle distribution function $f_{\alpha}$ yields its zeroth moment:

[TABLE]

The non-orthogonal basis vectors can then be written as

[TABLE]

In the above, symbols such as $\ket{e_{x}^{2}e_{y}}=\ket{e_{x}e_{x}e_{y}}$ denote a vector that arise from the elementwise vector multiplication of vectors $\ket{e_{x}}$, $\ket{e_{x}}$ and $\ket{e_{y}}$. In order to map changes of moments back to changes in the distribution function, we group the above set of vectors as a transformation matrix $\bm{\mathsf{T}}$, which reads as

[TABLE]

We then define the raw moments of order ($m+n$) of the distribution function $f_{\alpha}$, its equilibrium $f_{\alpha}^{eq}$, and the source terms $S_{\alpha}$ to represent the body force, respectively, as

[TABLE]

Here, and in what follows, the prime (′) symbols denote various raw moments. In terms of the nominal, nonorthogonal transformation matrix $\bm{\mathsf{T}}$ the relation between the various moments and their corresponding states in the velocity space can be written as

[TABLE]

where

[TABLE]

are the various quantities in the velocity space, and

[TABLE]

are the corresponding states in the moment space.

The MRT-LBM with trapezoidal rule to represent the source term with second order accuracy can be written as

[TABLE]

where the diagonal relaxation time matrix $\hat{\bm{\mathsf{\Lambda}}}$ can be represented as

[TABLE]

In order to obtain an effectively explicit scheme, we apply the transformation [59, 60] $\bar{f}_{\alpha}=f_{\alpha}-\frac{1}{2}S_{\alpha}\delta_{t}$, or equivalently $\widehat{\bar{\mathbf{m}}}=\widehat{\mathbf{m}}-\frac{1}{2}\widehat{\mathbf{S}}\delta_{t}$ and $\widehat{\bar{\kappa}}_{{x^{m}}{y^{n}}}^{{}^{\prime}}={\hat{\kappa}}_{{x^{m}}{y^{n}}}^{{}^{\prime}}-\frac{1}{2}{\hat{\sigma}}_{{x^{m}}{y^{n}}}^{{}^{\prime}}\delta_{t}$, and the MRT-LBE can be written as

[TABLE]

The moment equilibria $\hat{\kappa}_{{x^{m}}{y^{n}}}^{eq^{\prime}}$ at different orders can be written as [58]

[TABLE]

which are obtained from the discrete representation of the local Maxwellian by transforming back their central moments at a given order to their corresponding raw moments. Here, $c_{s}$ is the speed of sound, and in the present work, we typically set $c_{s}^{2}=1/3$. Also, moments of the source terms $\hat{\sigma}_{{x^{m}}{y^{n}}}^{eq^{\prime}}$ follows as [58]

[TABLE]

where $\bm{F}=(F_{x},F_{y})$. The hydrodynamic fields are given by

[TABLE]

where $\bm{u}=(u_{x},u_{y})$. The above represents the solution of the NSE (Eqs. (2) and (3)), with the kinematic bulk and shear viscosities related to the relaxation times via $\zeta=c_{s}^{2}\left(\frac{1}{\omega_{3}}-\frac{1}{2}\right)\delta_{t}$ and $\nu=c_{s}^{2}\left(\frac{1}{\omega_{j}}-\frac{1}{2}\right)\delta_{t}$, where $j=4,5$ respectively. The remaining relaxation times for the higher order moments, which influence the numerical stability, are set to unity in this work.

2.1 Moment relationships for the symmetric velocity gradient tensor: Chapman-Enskog Analysis

We will now perform a Chapman-Enskog analysis [61] to determine the expressions that relate the symmetric velocity gradient tensor to certain components of the local (non-equilibrium) moments. Expanding the moments about its equilibria as well as applying the standard multiscale expansion of the time derivatives in the MRT-LB scheme given in the previous section

[TABLE]

where $\epsilon$ is ${a}$ small bookkeeping perturbation parameter, and also performing a Taylor series expansion of the streaming operator in Eq. (20), i.e.,

[TABLE]

and converting all quantities in the velocity space to the moment space (via Eq. (16)) and using $\widehat{\bar{\mathbf{m}}}=\widehat{\mathbf{m}}-\frac{1}{2}\widehat{\mathbf{S}}\delta_{t}$, we obtain the following system of moment equations at consecutive order in $\epsilon$:

[TABLE]

where $\widehat{\bm{\mathsf{E}}}_{i}=\bm{\mathsf{T}}(\bm{e}_{i}\bm{\mathsf{I}})\bm{\mathsf{T}}^{-1},i\in\{x,y\}$. In order to obtain the hydrodynamic macroscopic equations, in the leading, i.e., $O(\epsilon)$ system (see Eq. (26b)), the equations representing the evolution of the moment components up to the second order are necessary, which read as (see A for details)

[TABLE]

Analogously, at the next, i.e., $O(\epsilon^{2})$ level (see Eq. (26c)), the relevant moment equations to recover the equations of the fluid motion written up to the first order as

[TABLE]

Here, the components of the second-order non-equilibrium moments $\widehat{m}_{3}^{(1)}$, $\widehat{m}_{4}^{(1)}$ and $\widehat{m}_{5}^{(1)}$ (which represent $\widehat{\kappa}_{xx}^{\prime(1)}+\widehat{\kappa}_{yy}^{\prime(1)}$, $\widehat{\kappa}_{xx}^{\prime(1)}-\widehat{\kappa}_{yy}^{\prime(1)}$ and $\widehat{\kappa}_{xy}^{\prime(1)}$, respectively) are unknowns. They can be obtained from Eqs. (27d), (27e) and (27f), respectively, where the time derivatives $\partial_{t_{0}}\left(2c_{s}^{2}\rho+\rho(u_{x}^{2}+u_{y}^{2})\right)$, $\partial_{t_{0}}\left(\rho(u_{x}^{2}-u_{y}^{2})\right)$ and $\partial_{t_{0}}\left(\rho u_{x}u_{y}\right)$ are eliminated in favor the spatial derivatives using the leading order mass and momentum equations (i.e., Eqs. (27a)–(27c), respectively). For details, see e.g., [58, 62]. Neglecting all terms of $O(u^{3})$ and higher, we can obtain the expressions for the various components of the non-equilibrium second order moments related to the symmetric part of the velocity gradient tensor $S_{ij}=\frac{1}{2}(\partial_{j}u_{i}+\partial_{i}u_{j})$ (i.e., $\partial_{x}u_{x}$, $\partial_{y}u_{y}$ and $\partial_{y}u_{x}+\partial_{y}u_{x}$), which read as [58, 62]

[TABLE]

When these expressions are substituted in Eqs. (28b) and (28c), and then combining the $O(\epsilon)$ and $O(\epsilon^{2})$ moment equations up to the first order, the NSE given Eqs. (2) and (3) follows. The non-equilibrium moment relations given in Eqs. (29a)–(29c) will be combined further with the developments given in the next section to develop a local computing approach for the vorticity field later in Sec. 4.

3 MRT-LBM for Transport of a Passive Scalar

The solution of the advection-diffusion of the passive scalar field $\phi$ given by the CDE in Eq. (4) will now be represented by using another MRT-LBM. Considering the D2Q9 lattice again, which, as required, supports the off-diagonal third order moment equilibria independently as noted in the Introduction, we use the same natural moment basis given in Eq. (7) as well as the resulting transformation matrix $\bm{\mathsf{T}}$ (see Eq. (8)). First, we define the relation between the various raw moments and the corresponding distribution function $g_{\alpha}$ and their equilibria $g_{\alpha}^{eq}$ for this MRT-LBE as

[TABLE]

where

[TABLE]

are given in the velocity space, and

[TABLE]

represent the equivalent states in the moment space. Here, the various sets of raw moments are defined as follows:

[TABLE]

Then the MRT-LBE using a non-orthogonal moment basis for the solution of the CDE can be written as

[TABLE]

where $\hat{\mathbf{\Lambda}}^{\phi}$ is the diagonal relaxation time matrix given by

[TABLE]

A key element in this work is the prescription of the moment equilibria $\hat{\mathbf{n}}^{eq}$ (Eq. (33)) used in Eq. (39) to enable a local computation of the antisymmetric velocity gradient tensor or the vorticity field. The passive scalar $\phi$ is advected by the local velocity field $\bm{u}$, and $\phi$ hence its solution procedure, in principle, has a complete information on the kinematics of the fluid elements undergoing a variety of motions when it is carefully designed. As such, most of the components of the moment equilibria $\hat{\mathbf{n}}^{eq}$ can be constructed in analogy with $\hat{\mathbf{m}}^{eq}$ given in Eq. (21), where the density $\rho$ is replaced by the scalar field $\phi$. On the other hand, in view of the above consideration, in order to extract the local intrinsic rotation rate of the fluid element related to the antisymmetric velocity gradient tensor, we prescribe anisotropy in the scalar flux $(\phi\bm{u})$ components used in the third order moment equilibria, which, as we shall see in the following, does not affect the recovery of the macroscopic CDE. Thus, we set

[TABLE]

where $c_{s\phi}$ is an independent parameter related to the diffusivity $D_{\phi}$ (see below), and we typically set $c_{s\phi}^{2}=1/3$ in this work. Here, $\beta_{1}$ and $\beta_{2}$ are free parameters that prescribe anisotropy on the scalar flux appearing in the third order moment equilibria. Typically, $\beta_{1}\approx 1$ and $\beta_{2}\approx 1$, but $\beta_{1}-\beta_{2}\neq 0$ , i.e., a small intentional anisotropy is introduced to locally recover the magnitude of the intrinsic rotation rate of the fluid motion (see the following section). The scalar field $\phi$ is then obtained as the zeroth moment of the distribution function $g_{\alpha}$, which evolves according to Eq. (39) in the form of the standard collide-and-steam steps:

[TABLE]

Then, the above represents the solution of the CDE (Eq. (4)), with the diffusivity related to the relaxation times via $D_{\phi}=c_{s\phi}^{2}\left(\frac{1}{\omega_{j}^{\phi}}-\frac{1}{2}\right)\delta_{t}$ where $j=1,2$. It may be noted that B–D present extensions of our approach to other collision models, including the SRT-LBM and central moment LB formulations.

3.1 Moment relationships for the scalar gradient vector and skew-symmetric velocity gradient tensor: Chapman-Enskog Analysis

We will now perform a C-E analysis of the MRT-LBE for the passive scalar field. Applying the moment expansion about its equilibria and a multiscale expansion of the time derivative to Eq. (39)

[TABLE]

where $\epsilon=\delta_{t}$ and also using a Taylor expansion of the streaming operator ${\mathbf{g}}_{\alpha}(\bm{x}+\bm{e}_{\alpha}\epsilon,t+\epsilon)=\sum_{j=0}^{\infty}\frac{\epsilon^{j}}{j!}(\partial_{t}+\bm{e}_{\alpha}\cdot\bm{\nabla})^{j}{\mathbf{g}}(\bm{x},t)$, the following moment equations at consecutive order in $\epsilon$ can be obtained:

[TABLE]

where $\widehat{\bm{\mathsf{E}}}_{i}$ is the same as that given earlier. Some of the relevant components at the leading order (i.e., $(O(\epsilon))$) of the moment system (see Eq. (44b)) are given as

[TABLE]

where the above can be obtained by replacing ${\widehat{\kappa}_{x^{m}y^{n}}^{eq^{\prime}}}$ in the corresponding C-E analysis for the fluid motion with $\hat{\eta}_{x^{m}y^{n}}^{eq^{\prime}}$ (see the previous section and A for details) and allowing for the relaxation of the first order moments, since only the scalar field $\phi$ is conserved in the present case. Similarly, the leading component (i.e., the zeroth order) of the moment system at the $O(\epsilon^{2})$ level to recover the CDE is obtained from Eqs. (44c) can be written as

[TABLE]

Now, in order to derive the CDE, we need to combine Eq. (45a) and $\epsilon$ times Eq. (46) by using $\partial_{t}=\partial_{t_{0}}+\epsilon\partial_{t_{1}}$, which requires $\hat{n}^{(1)}_{1}$ and $\hat{n}^{(1)}_{2}$. These first order non-equilibrium moments ($\hat{n}^{(1)}_{1}$ and $\hat{n}^{(1)}_{2}$) can be obtained from Eqs. (45b) and (45c), respectively, where the time derivatives are eliminated in favor of the spatial derivatives by using the leading order mass, momentum and scalar conservation equations (i.e., Eqs. (27a), (27b), (27c) and (45a)). Hence after some simplification, and neglecting terms of $O(u^{2})$ and higher, we get the components of the first order non-equilibrium moments in terms of the components of the scalar gradient vector $\partial_{i}\phi$ as

[TABLE]

It may be noted that in the derivation of these non-equilibrium moment components, only the spatial derivatives of the second order moment equilibria (i.e., $\hat{\eta}_{xx}^{eq^{\prime}}$, $\hat{\eta}_{yy}^{eq^{\prime}}$ and $\hat{\eta}_{xy}^{eq^{\prime}}$) are involved and do not involve the introduced anisotropy, which appears at a higher order level, i.e., for the third order moments of the equilibrium distribution via the factors $\beta_{1}$ and $\beta_{2}$ and hence the advection-diffusion of the passive scalar transport is correctly recovered.

As shown in the previous section, the symmetric components of the velocity gradient tensor $\partial_{x}u_{x},\partial_{y}u_{y}$ and $\partial_{x}u_{y}+\partial_{y}u_{x}$ can be obtained from the MRT-LBM for fluid flow. In order to obtain the skew-symmetric component, i.e., $\partial_{x}u_{y}-\partial_{y}u_{x}$, which would then provide a complete information about the velocity gradient tensor $\partial_{j}u_{i}$ and hence the vorticity field, we now exploit the additional degree of freedom available in the off-diagonal, second-order non-equilibrium moment equation resulting from the MRT-LBM for CDE, i.e., Eq. (45f). Simplifying this equation by eliminating the time derivative in favor of spatial derivatives and eliminating higher order terms (i.e., $(O(u^{2}))$ and above), we get

[TABLE]

which can be rewritten as

[TABLE]

Clearly, the anisotropy introduced into the scalar flux components in the third order moment equilibria results in an additional flexibility via an independent equation given above (Eq. (49)). In this equation, the gradients of the scalar field in the Cartesian coordinate directions $\partial_{x}\phi$ and $\partial_{y}\phi$ can be obtained locally from Eq. (47a) and Eq. (47b); and with the knowledge of the off-diagonal second-order non-equilibrium moment component $\hat{n}^{(1)}_{5}$, then Eq. (49) represents an additional independent equation to compute the antisymmetric velocity gradient tensor component, which will be exploited further in the next section.

4 Derivation of local expressions for the complete velocity gradient tensor and vorticity field

In order to independently determine the cross-derivative components of the velocity gradient tensor, i.e., $\partial_{y}u_{x}$ and $\partial_{x}u_{y}$, we combine the analysis presented in the two earlier sections. In particular, the Eq. (29c) resulting from the solution of the MRT-LBM for fluid flow and Eq. (49) from the MRT-LBM for CDE, can be rewritten as

[TABLE]

where, when $\phi\neq 0$,

[TABLE]

Solving Eqs. (50a) and (50b), we get following independent and local expressions for the off-diagonal components or the cross derivatives of the velocity field, which is one of the main results of this work:

[TABLE]

The diagonal components of the velocity gradient tensor, i.e., $\partial_{x}u_{x}$ and $\partial_{y}u_{y}$ follow from solving the Eqs. (29a) and (29b) resulting from the MRT-LBE for the fluid motion, which reads as

[TABLE]

and this completes the determination of all the components of the velocity gradient tensor. Finally, a local expression for the pseudo-vector, viz., the vorticity field $\bm{\omega}=\bm{\nabla}\times\bm{u}=(0,0,\omega_{z})$ can be obtained by combining Eqs. (52a) and (52b) as

[TABLE]

which is another key result arising from our analysis.

The terms $N_{xy}$ and $N_{xy}^{\phi}$ given in Eqs. (51a) and (51b), respectively, which are needed in Eqs. (52a), (52b) and (54) can be evaluated locally using

[TABLE]

and also since $\hat{n}^{(1)}_{1}={\widehat{\eta}_{x}^{{}^{\prime}}}-{\widehat{\eta}_{x}^{eq^{\prime}}}={\widehat{\eta}_{x}^{{}^{\prime}}}-\phi u_{x}$ and $\hat{n}^{(1)}_{2}={\widehat{\eta}_{y}^{{}^{\prime}}}-{\widehat{\eta}_{y}^{eq^{\prime}}}={\widehat{\eta}_{y}^{{}^{\prime}}}-\phi u_{y}$, and from Eqs. (47a) and (47b), we have the required local expressions for the derivatives of the scalar field, which read as

[TABLE]

Note that $\beta_{1}\approx 1$ and $\beta_{2}\approx 1$, but $\beta_{1}\neq\beta_{2}$ and are otherwise free parameters. We typically set $\beta_{1}=1,\beta_{2}=0.9$ in this work. In addition, the expressions for $\hat{m}_{3}^{(1)}$ and $\hat{m}_{4}^{(1)}$ needed in the diagonal components of the velocity gradient tensor, i.e., Eqs. (53a) and (53b) can be written as

[TABLE]

In the above, $\widehat{\kappa}_{xx}^{{}^{\prime}}$, $\widehat{\kappa}_{yy}^{{}^{\prime}}$, $\widehat{\kappa}_{xy}^{{}^{\prime}}$, $\widehat{\eta}_{x}^{{}^{\prime}}$, $\widehat{\eta}_{y}^{{}^{\prime}}$ and $\widehat{\eta}_{xy}^{{}^{\prime}}$ are the raw moment components of different orders of the respective distribution functions. The formulation presented above thus allows local computation of the complete velocity gradient tensor and hence the vorticity field without relying on any finite difference approximations of the velocity field.

5 Results and discussion

In this section, we will perform a numerical validation study of the new DDF MRT-LB scheme for vorticity computation. In this regard, we will consider a set of well-defined benchmark flow problems for which analytical solutions or numerical results for the vorticity field are available or can be derived. In the simulations results presented in the following, the relaxation times for the second order moments of the MRT-LBM for the flow field ($\omega_{4}=\omega_{5}=1/\tau$) are chosen to specify the desired fluid viscosity, while those for the first order moments of the MRT-LBM for the scalar field ($\omega_{1}^{\phi}=\omega_{2}^{\phi}=1/\tau^{\phi}$) are prescribed to select the diffusivity. The relaxation times of all the higher order moments for both the LB schemes are set to unity for simplicity. Unless otherwise specified, we consider the use of lattice units, i.e., $\delta_{x}=\delta_{t}=1.0$ typical for LB simulations and a reference density of unity is considered in this work. For all the benchmark problems reported in what follows, we set the coefficients for the scalar flux terms in the third order moment equilibria of the MRT-LBM for the scalar field to $\beta_{1}=1.0$ and $\beta_{2}=0.9$.

5.1 Poiseuille flow

As the first benchmark problem, a steady flow between two parallel plates with a width $2L$ driven by a constant body force $F_{x}$, i.e., the Poiseuille flow, is simulated. This flow problem has an analytical solution for the vorticity field as the linear profile $\omega_{z}(y)=2U_{max}y/L^{2}$, which can be obtained from the parabolic velocity profile $u_{x}(y)=U_{max}[1-\frac{y^{2}}{L^{2}}]$, where $U_{max}=\frac{F_{x}L^{2}}{2\rho\nu}$ is the maximum centerline velocity, $\nu$ and $\rho$ are fluid kinematic viscosity and density, respectively. Periodic boundary conditions are employed in the streamwise direction and no-slip condition for the velocity field are imposed using the half-way bounce back scheme. The computational domain is resolved using $3\times 151$ lattice nodes. For the scalar field, we consider fixed values at the bottom and top walls as $\phi_{L}=1.0$ and $\phi_{H}=2.0$, respectively, and its diffusivity is specified by choosing $\tau^{\phi}=0.57$. At a fixed body force $F_{x}=3\times 10^{-6}$, computations are carried by adjusting the fluid kinematic viscosity such that the following five sets of maximum centerline velocities are considered: $U_{max}=0.01,0.03,0.05$, and $0.08$. Figure 1 shows a comparison between numerical results for the vorticity profiles obtained using the DDF MRT-LB scheme and the analytical solutions for the above set of values for $U_{max}$. Excellent agreement is seen.

5.2 Four-rolls mill flow problem

In order to examine the validity of our approach for a flow problem with fully two-dimensional (2D) spatially varying distribution of the vorticity field, we consider next the four-rolls mill flow. It is a steady, rotational flow consisting of an array of counter-rotating vortices generated by the stirring action of a suitably specified local body force $F_{x}=F_{x}(x,y)$ and $F_{y}=F_{y}(x,y)$ in a periodic square domain of size $2\pi\times 2\pi$. It is a modified form of the Taylor-Green vortex flow. The spatially varying driving body force can be written as $F_{x}(x,y)=2\rho_{0}\nu u_{0}\sin x\sin y$ and $F_{y}(x,y)=2\rho_{0}\nu u_{0}\cos x\cos y$, where $\rho_{0}$ is the reference density, $\nu$ is kinematic viscosity and $u_{0}$ is the velocity scale and $0\leq x,y\leq 2\pi$. A solution of the simplified form of the Navier-Stokes equations with the above described body force yields the explicit form of the local velocity field, which reads as $u_{x}(x,y)=u_{0}\sin x\sin y$ and $u_{y}(x,y)=u_{0}\cos x\cos y$. Then, the analytical solution for the local vorticity field $\omega_{z}(x,y)$ can be derived by taking the curl of the above velocity field, which can be written as

[TABLE]

For the purpose of setting up simulations, the Reynolds number for this flow problem can be defined as $\mbox{Re}=u_{0}2\pi/\nu$ and the viscosity can be written as $\nu=\frac{1}{3}(\tau-\frac{1}{2})\Delta x$, where $\Delta x=\Delta t=2\pi/N$, where $N$ is the number of grid nodes in each direction. We consider a grid resolution of $84\times 84$ and a velocity scale $u_{0}=0.035$ to simulate four-rolls mill flow at $\mbox{Re}=40$. The scalar field is initialized to a uniform value of $2.0$ in this periodic domain with the relaxation time $\tau^{\phi}=0.57$. Figure 2 presents a comparison between the spatial distribution of the computed vorticity field obtained using the DDF MRT-LB scheme and the analytical solution. Due to the presence of a system of counter-rotating vortices, the vorticity field, represented by harmonic functions analytically, dramatically varies both in its magnitude and sign. Good agreement between the two results are evident.

Furthermore, in order to make a more head-on comparison, Fig. 3 shows the computed vorticity profiles $\omega_{z}(x,y)$ computed using our LB scheme along various horizontal sections at $y=0,\pi/4,\pi/2,\pi,5\pi/4$ along with results based on the analytical solution. It is evident that there is a very good agreement between our numerical results and the analytical solution.

5.2.1 Grid convergence study

We will now assess the order of accuracy of the convergence of the vorticity computation via our DDF MRT-LB scheme. In this regard, at a fixed viscosity of $\nu=0.00218$ with a velocity scale $u_{0}=0.045$, we consider the following sequence of four different resolutions: $24\times 24$, $48\times 48$, $96\times 96$ and $192\times 192$. For each case, we measure the following global relative error ($E_{g,\omega}$) between the vorticity field computed using the DDF MRT-LB scheme given by $\omega_{c}$ and the corresponding analytical solution denoted by $\omega_{a}$:

[TABLE]

where the summations in the above are for the whole computational domain. The rate of convergence of the global relative error is depicted using a log-log scale in Fig. 4. From this figure, it can be seen that the relative error exhibits a slope of -2.0, which demonstrates that the vorticity computation using our approach is second order accurate.

5.3 Womersley flow

In order to validate our approach for the calculation of the vorticity field in unsteady flows, a 2D pulsatile flow between two parallel plates separated by a width $2L$ driven by a sinusoidally time-dependent body force $F_{x}(t)$ is considered. This classical Womersley flow problem is subjected to a periodic body force given by $F_{x}=F_{m}\mbox{cos}(\Omega t)$, where $F_{m}$ is the maximum amplitude of the force and $\Omega=2\pi/T$ is the angular frequency and $T$ being the time period. Considering that this pulsatile flow is laminar and incompressible, the analytical solution for velocity field is given as [63]

[TABLE]

where $\gamma=\sqrt{i\mbox{Wo}^{2}}$ and $\mbox{Wo}=L\sqrt{(\Omega/\nu)}$ is the Womersley number. Here, and in the following $\mathcal{R}\{\cdot\}$ refers to the real part of the expression. Then, the analytical solution for the local time dependent vorticity field $\omega_{z}(y,t)$ can be readily obtained by taking the curl of the velocity field as

[TABLE]

We consider a grid resolution of $3\times 101$, maximum force amplitude $F_{m}=1.0\times 10^{-5}$ with a time period $T=10,000$ and two different values of the Womersley number, i.e., $\mbox{Wo}=4.0$ and $\mbox{Wo}=7.0$, which are specified by setting the relaxation times for the MRT-LBM for the flow field to be $\tau=0.781$ and $\tau=0.596$, respectively. Periodic boundary conditions and the no-slip boundary conditions are considered for the inlet/outlet in the streamwise direction and along the two parallel walls, respectively. The parameters and the boundary conditions for the scalar field are the same as those considered for the Poiseuille flow simulations discussed earlier. Figure 5 presents a comparison between the computed vorticity profiles obtained using the DDF MRT-LB scheme and the corresponding analytical solution at different time instants within a time period $T$. It is evident that the vorticity field is subjected to strong temporal and spatial variations, which are seen to increase with the Womersley number. These are very well reproduced quantitatively by our local computational approach.

5.4 Lid-driven cavity flow

As the final validation study, we consider simulation of a shear driven flow within a square cavity due to the motion of its top lid in order to compare the computed vorticity fields against those based on numerical results obtained by a finite-difference method. The lid-driven cavity flow is a classical benchmark problem characterized by complex flow features involving vortical patterns of different sizes which are strongly influenced by the nonlinear effects, i.e., the Reynolds number (see e.g., [64, 65, 66]). If $U_{0}$ is the velocity imposed on the top lid of a square cavity of side length $L$, its Reynolds number Re can be expressed as $\mbox{Re}=U_{0}L/\nu$. We perform numerical simulations of shear-driven flow within a cavity at $\mbox{Re}=400,1000$ and $3200$ by considering grid resolutions of $100\times 100$, $300\times 300$ and $450\times 450$. In this regard, the lid velocity $U_{0}$ is set to be $0.05$. The no-slip boundary conditions are prescribed on the walls via the standard half-way bounce-back condition, and including a momentum augmentation term for the moving top lid (see e.g., [42] for details). The scalar field $\phi$ is set be equal to $1.0$ on all the boundaries. Simulations are carried out until steady state is reached in each case. Figure 6 presents comparisons of the computed contours of the vorticity fields obtained using our DDF-LB scheme against numerical results based on the finite-difference method at $\mbox{Re}=400,1000$ and $3200$.

The vorticity patterns are found to agree very well with one another. It is worth emphasizing that our DDF-LB formulation used a completely local algorithm, while the finite-difference computations involved non-local operations to estimate the velocity derivatives. As Re increases, the flow becomes progressively more complex involving the clustering of finer vortical features near walls, which are well reproduced by our scheme.

6 Summary and Conclusions

A quantitative knowledge of the local skew-symmetric velocity gradient tensor, or equivalently the vorticity field, in conjunction with the symmetric velocity gradient tensor is crucial for various applications, including those related to techniques for the identification of flow structures and in the modeling of complex fluids. In many situations, it is required to compute the fluid motion coupled to the transport by advection and diffusion of a scalar field. In the mesoscopic LB methods, the hydrodynamics (i.e., the NSE) and the scalar transport (i.e., the CDE) are commonly computed via the evolution of a pair of distribution functions represented by means of two LBMs. In such double distribution functions (DDF) based LB approaches, we present a new strategy for computing the vorticity field locally via exploiting the additional degrees of freedom available in the construction of the higher order moment equilibria in the collision model for the representation of the scalar transport to obtain the necessary additional independent relations. In particular, we have shown that this can be achieved by introducing an intensional anisotropy in the scalar flux components in the third order, off-diagonal moment equilibria, and then combining the second-order, off-diagonal non-equilibrium moment components of both the LB schemes. This approach for local vorticity computation has several advantages, which include the following. Any pair of lattice sets in the DDF-LBMs that support the third order off diagonal moments independently, which includes the various standard lattice velocity models in different dimensions, can allow a local determination of the complete flow kinematics, including the skew-symmetric velocity gradient tensor. It imposes no additional constraints on the higher order equilibrium moments of the LBM for the flow field, which can be solved by using any standard formulation without modification thereby maintaining its numerical characteristics intact. Since the vorticity computation are based on distribution functions, which are generally solved to be second order accurate, the resulting mesoscopic and local computation of vorticity and the strain rate tensor are second order accurate as well. Moreover, the algorithm is completely local and do not depend on any finite difference approximations of the velocity derivatives, which is consistent with the general philosophy of the LBM as it well suited for implementation on parallel computers. The presented strategy is general and is applicable to a variety of collision models.

In the present work, for the purpose of demonstration, we formulate our approach by constructing in detail a DDF formulation using a MRT-LBM for the solution of the fluid motion and another MRT-LBM involving an anisotropy in the scalar flux components in the third order equilibria for the transport of a scalar field, each on a D2Q9 lattice. By means of a Chapman-Enskog analysis, we have shown that the former provides the necessary second order non-equilibrium moment equations to determine the symmetric velocity gradient tensor, while the latter yields additional corresponding moment relations to obtain the skew-symmetric velocity gradient tensor. For simplicity, the MRT-LBMs are constructed using natural, non-orthogonal moment bases. In order to validate our new approach, we have presented comparisons of the computed vorticity fields against the analytical and/or solutions for various benchmark problems such as the steady flow in a channel, four-rolls mill flow, time-dependent pulsatile (Womersley) flow in a channel, and lid-driven cavity flow at different Reynolds numbers, which demonstrate its good accuracy. In addition, an analysis of the method for various grid resolutions establishes its second order convergence for computing vorticity. While the focus of this work is on presenting and validating a new method for local computation of vorticity field using a DDF-LB scheme based on a MRT formulation, we have also discussed its deveopment for other collision models such as those based on SRT and cascaded central moments. Extensions our approach to 3D for various standard lattice sets (i.e., D3Q15, D3Q19 and D3Q27) and using a variety of collision models [67] will be reported in a future investigation. It may be noted that the method presented here can be extended to vectorial [68] and tensorial [69] forms of distribution functions to model and locally compute the skew-symmetric velocity gradient contributions in the general constitutive relations for complex fluids. Moreover, spin relaxation to the vorticity and the coupling of the intrinsic angular momentum to the linear momentum need to be accounted for in molecular liquid flows in nanoscale confined geometries [29, 30, 31], which can be modeled as generalization of the Cosserat theory for micropolar fluids [70, 71, 72, 28]. The approach presented here can also be used to construct LB schemes to locally represent such effects, which are subjects for future studies.

Acknowledgements

The authors would like to acknowledge the support of the US National Science Foundation (NSF) under Grant CBET-1705630.

Appendix A Relation between non-equilibrium moments and spatial derivatives of components of moment equilibria for D2Q9 lattice

For better clarity, the $O(\epsilon)$ moment system using a non-orthogonal moment basis given in Eq. (26b) in Sec. 2, i.e., $(\partial_{t_{0}}+\widehat{\bm{\mathsf{E}}}_{i}\partial_{i})\mathbf{\widehat{m}}^{(0)}=-\widehat{\bm{\mathsf{\Lambda}}}\mathbf{\widehat{m}}^{(1)}+\mathbf{\widehat{S}}$, which forms a main element in the derivation, can be expanded explicitly in terms of their various components as follows:

[TABLE]

In general, it can be seen that any non-equilibrium moment of order $n$ depends on the spatial derivatives of equilibrium moments of order $(n+1)$ and $(n-1)$. In particular, the diagonal components of the second order moment ($\hat{m}^{(1)}_{3}$ and $\hat{m}^{(1)}_{4}$) depend on the moment equilibria of first order ($\hat{\kappa}^{eq^{\prime}}_{x}$ and $\hat{\kappa}^{eq^{\prime}}_{y}$) and third order ($\hat{\kappa}^{eq^{\prime}}_{xxy}$ and $\hat{\kappa}^{eq^{\prime}}_{xyy}$), while the off-diagonal second order moment ($\hat{m}^{(1)}_{5}$) depends only on those of the third order equilibrium moments ($\hat{\kappa}^{eq^{\prime}}_{xxy}$ and $\hat{\kappa}^{eq^{\prime}}_{xyy}$). These considerations are important in establishing the relationship between the non-equilibrium second-order moments and the velocity gradient tensor components. In the case of the LBE for computing fluid flow, the symmetry of their moment equilibria to respect the isotropy of the viscous stress tensor limits the dependence of the corresponding non-equilibrium second order moments to only on the symmetric part of the velocity gradient tensor (i.e., the strain rate tensor). However, the construction of the LBE for computing the transport of a passive scalar represented by the CDE does not need to satisfy these restrictive constraints, and the additional degrees of freedom available for the higher order moments can be suitably exploited [57]. Indeed, since the diffusion term of the CDE need only to satisfy a lower degree of isotropy than that of the viscous term of the NSE, the third order moment equilibria for solving the former case can be specifically designed to locally represent the skew-symmetric part of the velocity gradient tensor via the respective off-diagonal non-equilibrium second-order moment (based on an equation analogous to the sixth equation in the above moment system with ${\hat{\kappa}_{x^{m}y^{n}}^{eq^{\prime}}}$ replaced by ${\hat{\eta}_{x^{m}y^{n}}^{eq^{\prime}}}$ and $\hat{m}^{(1)}_{j}$ by $\hat{n}^{(1)}_{j}$ – Sec. 3).

Appendix B SRT-LBM for solution of scalar transport to recover the skew-symmetric velocity gradient tensor

In this appendix, we will present a special case of the single-relaxation-time (SRT)-LBM for the solution of the convection-diffusion equation of a passive scalar field to recover the skew-symmetric velocity gradient tensor. This can be written as

[TABLE]

where the post-collision distribution functions $\tilde{g}_{\alpha}$ are prescribed by an update of $g_{\alpha}$ involving their relaxation to local equilibrium distribution functions $g_{\alpha}^{eq}$ at a single relaxation time $\tau_{\phi}$. A key aspect here is the construction of $g_{\alpha}^{eq}=g_{\alpha}^{eq}(\phi,\bm{u},\beta_{1},\beta_{2})$ that facilitates the recovery of $\Omega_{ij}=\frac{1}{2}(\partial_{j}u_{i}-\partial_{i}u_{j})$, which we achieve by mapping the various equilibrium moment components derived earlier for the D2Q9 lattice (see Sec. 3) to the velocity space. In this regard, defining

[TABLE]

we can relate it to ${\mathbf{g}}^{eq}=(g_{0}^{eq},g_{1}^{eq},g_{2}^{eq}\dots g_{8}^{eq})^{\dagger}$ via $\hat{\mathbf{q}}^{eq}=\bm{\mathsf{P}}\mathbf{g}^{eq}$ using the bare moment basis $\bm{\mathsf{P}}$ given by

[TABLE]

where

[TABLE]

Inverting, that is,

[TABLE]

we can then obtain the one-to-one mapping between the equilibrium moments and the equilibrium distribution functions. Thus, we have

[TABLE]

where $\hat{\eta}_{x^{m}y^{n}}^{eq^{\prime}}$ are given in Eq. (41). In particular, the third order moment equilibrium components $\hat{\eta}_{xxy}^{eq^{\prime}}$ and $\hat{\eta}_{xyy}^{eq^{\prime}}$ contain the intensional anisotropy needed for recovering the skew-symmetric velocity gradient tensor. Setting the relaxation time $\tau_{\phi}$ in terms of the relaxation parameter $\omega_{j}^{\phi}$ as $\omega_{j}^{\phi}=1/\tau_{\phi}$ and using the definitions of $\hat{\eta}_{x^{m}y^{n}}^{{}^{\prime}}$ and $\hat{\eta}_{x^{m}y^{n}}^{eq^{\prime}}$ given in Sec. 3, the local expressions derived earlier in Sec. 4 for $N_{xy}^{\phi}$, $\hat{n}_{5}^{(1)}$, $\partial_{x}\phi$, $\partial_{y}\phi$, $\partial_{x}u_{y}$ and $\partial_{y}u_{x}$ and $\omega_{z}$ in terms of the non-equilibrium moments are valid.

Appendix C Cascaded LBM based on central moments for solution of scalar transport to recover the skew-symmetric velocity gradient tensor

In this section, we will present further development of our formulation to a more general cascaded LBM based on central moments [73] extended for the solution of a scalar transport [54, 55, 56] capable of locally computing the skew-symmetric velocity gradient tensor. In this regard, we define the central moments of the distribution functions and their equilibrium as

[TABLE]

We prescribe the central moment equilibria based on those of the local Maxwellian, by replacing the density with the scalar field $\phi$ (see e.g., [54, 55, 56]). Usually, the third order central moment equilibria then become $\hat{\eta}^{eq}_{xxy}=\hat{\eta}^{eq}_{xyy}=0$ and the corresponding raw moment equilibria are $\hat{\eta}^{eq\prime}_{xxy}=c_{s\phi}^{2}\phi u_{y}+\phi u_{x}^{2}u_{y}$ and $\hat{\eta}^{eq\prime}_{xyy}=c_{s\phi}^{2}\phi u_{x}+\phi u_{x}u_{y}^{2}$ [54, 55, 56]. On the other hand, to enable local computation of the vorticity field, our derivation in Secs. 3 and 4 required the above raw moment components to be modified to $\hat{\eta}^{eq\prime}_{xxy}=\beta_{1}c_{s\phi}^{2}\phi u_{y}+\phi u_{x}^{2}u_{y}$ and $\hat{\eta}^{eq\prime}_{xyy}=\beta_{2}c_{s\phi}^{2}\phi u_{x}+\phi u_{x}u_{y}^{2}$. These are equivalent to modifying the central moment equilibria $\hat{\eta}^{eq}_{xxy}$ and $\hat{\eta}^{eq}_{xxy}$ as $\hat{\eta}^{eq}_{xxy}=(\beta_{1}-1)c_{s\phi}^{2}\phi u_{y}$ and $\hat{\eta}^{eq}_{xyy}=(\beta_{2}-1)c_{s\phi}^{2}\phi u_{x}$, where $(\beta_{1}-1)$ and $(\beta_{2}-1)$ represent the degree of anisotropy in the scalar flux components $\phi u_{y}$ and $\phi u_{x}$, respectively. Hence, we enumerate all the central moment equilibria for the D2Q9 lattice as

[TABLE]

The cascaded LBM then reads as

[TABLE]

where $\widehat{\mathbf{h}}=(\widehat{h}_{0},\widehat{h}_{1},\cdots,\widehat{h}_{8})^{\dagger}$ represents the changes in different moments due to collision via relaxation of central moments in a cascaded fashion. Here, $\bm{\mathsf{K}}=(K_{0},K_{1},\cdots,K_{8})^{\dagger}$ represents a matrix holding the orthogonal basis vectors given by

[TABLE]

To obtain the change in moments under collision $\widehat{\mathbf{h}}$, we need the following inner products:

[TABLE]

Then, we prescribe the relaxation of various central moments to their corresponding equilibria supported by the D2Q9 lattice as

[TABLE]

where $\mathbf{1}=\ket{1}$, $\mathbf{e}_{x}=\ket{e_{x}}$, $\mathbf{e}_{y}=\ket{e_{y}}$ and $\omega_{*}^{\phi}$ being the relaxation parameter of the central moment of order $(m+n)$. With the zeroth moment being conserved, i.e., a collision invariant, and evaluating Eq. (78) at various orders and then simplifying the resulting expressions, we obtain the changes in different moments due to cascaded collision as

[TABLE]

where $\omega_{1}^{\phi}=\omega_{2}^{\phi}$ controls the diffusivity $D_{\phi}=c_{s\phi}^{2}\left(\frac{1}{\omega_{j}^{\phi}}-\frac{1}{2}\right)\delta_{t}$, where $j=1,2$, while the relaxation parameters for the higher order moments $\omega_{3}^{\phi}=\omega_{4}^{\phi}$, $\omega_{5}^{\phi}$, $\omega_{6}^{\phi}$, $\omega_{7}^{\phi}$ and $\omega_{8}^{\phi}$ can be adjusted to improve numerical stability. Finally, expanding $(\bm{\mathsf{K}}\cdot\widehat{\mathbf{h}})_{\alpha}$ in Eq. (75a), the updates for the post-collision distribution functions read as

[TABLE]

Appendix D Non-cascaded central moment LBM for solution of scalar transport to recover the skew-symmetric velocity gradient tensor

For completeness, we will also present another version of a LBM based on central moments for solving the transport of the scalar field that allows local computation of the vorticity. Unlike C, the formulation given below is non-cascaded, i.e., the change of higher moments under collision do not depend on those of the lower moments. Rather, it is based on the relaxation of various central moments to their equilibria under collision, while involving systematic transformations between the distribution functions, raw moments and central moments before and after collision (similar to the algorithms presented in [38]). In this regard, we first enumerate the distribution functions, bare raw moments and central moments for the D2Q9 lattice, represented by vectors $\mathbf{g}$, $\hat{\mathbf{q}}$ and $\hat{\mathbf{q}}^{c}$, respectively, as

[TABLE]

[TABLE]

[TABLE]

Then, the mappings between the central moments, raw moments and distribution functions may be formally expressed in matrix-vector forms as

[TABLE]

where $\bm{\mathsf{P}}$ is a matrix representing the transformation from the distribution functions to the raw moments (see Eqs. (65) and (66)) and $\bm{\mathsf{\mathcal{F}}}$ is a frame transformation matrix that maps the raw moments to the central moments, i.e., containing the elements of $(\mathbf{e_{x}}-u_{x}\mathbf{1})^{m}(\mathbf{e_{y}}-u_{y}\mathbf{1})^{n}$. However, since $\bm{\mathsf{P}}$ and $\bm{\mathsf{P}}^{-1}$ are both sparse, while $\bm{\mathsf{\mathcal{F}}}$ as well as $\bm{\mathsf{\mathcal{F}}}^{-1}$ are of special lower triangular forms arising from the binomial expansions, it is neither necessary nor efficient to use them in matrix forms. Rather, we only list the resulting mapping expressions of the elements of each transformation before and after collision in the algorithm in what follows.

(a) Pre-collision raw moments

Expanding $\hat{\mathbf{q}}=\bm{\mathsf{P}}\mathbf{g}$, the raw moments before collision read as

[TABLE]

where

[TABLE]

(b) Pre-collision central moments

Based on $\hat{\mathbf{q}}^{c}=\bm{\mathsf{\mathcal{F}}}\hat{\mathbf{q}}$, the central moments from raw moments before collision follows. Hence, we obtain

[TABLE]

(c) Post-collision central moments: Relaxation of central moments under collision

We then prescribe the relaxation of various central moments to their equilibria at individual rates under collision, where the central moment equilibria that account for the anisotropy at the third order to recover the vorticity field are given in Eq. (74). Hence, the post-collision central moments can be written as

[TABLE]

The choices of the various relaxation times $\omega_{j}^{\phi}$, where $j=1,2,\ldots,8$ are the same as those given in C.

(d) Post-collision raw moments

The post-collision central moments can be mapped to those of raw moments via $\widetilde{\hat{\mathbf{q}}}=\bm{\mathsf{\mathcal{F}}}^{-1}\widetilde{\hat{\mathbf{q}}}^{c}$. It may be noted that the elements of $\bm{\mathsf{\mathcal{F}}}^{-1}$ (representing the inverse of binomial expansions) are the same of those of $\bm{\mathsf{\mathcal{F}}}$ (representing the binomial expansions) after making all the coefficients in the latter to be positive. Hence, we get

[TABLE]

(e) Post-collision distribution functions

Finally, the post-collision distribution functions can be obtained by simplifying $\mathbf{\widetilde{g}}=\bm{\mathsf{P}}^{-1}\widetilde{\hat{\mathbf{q}}}$, which yield

[TABLE]