1 Introduction↩︎

Isotropy is a central assumption in the derivation of classical lattice Boltzmann methods [1]. Not only do the position and velocity spaces have to be isotropic, the Maxwell-Boltzmann equilibrium distribution itself is isotropic.

Anisotropic fluid flow, however, is not uncommon in the rheology of complex fluids. One can think of polymer emulsions [2], often containing 1D or 2D components [3], or liquid crystals [4], where the anisotropy stems from the molecular shape itself. The hydrodynamics for molecular liquids with anisotropic molecules is also relevant for the interpretation of light scattering experiments [5], [6]. Another application is anisotropic advection-diffusion, which is a more common phenomenon since the diffusivity of many materials has a directional dependence. Therefore, anisotropic heat and mass transport can intrinsically not be described by standard theory. This limits the application of the method in the relevant fields of anisotropic heat conduction or diffusion [7].

A very important application of high technical relevance is flow in anisotropic porous media. The pore structure resulting for example from irregular particle shapes can introduce significant anisotropy in the resulting macroscopic flow field. In homogenized models this effective anisotropic flow behavior has to be treated accordingly to describe macroscopic flow patterns.

Other works related to anisotropic lattice Boltzmann models are often restricted by a very specific application, for example crystal growth [8]. An existing general approach [9] is limited to advection-diffusion problems. The presented work addresses this gap by deriving a universal lattice Boltzmann method for anisotropic equilibrium distributions. We fully derive the discretization in this new setting and compute the resulting macroscopic laws through Chapman-Enskog theory for a single-relaxation time collision operator. Consequently, we arrive at anisotropic Navier-Stokes equations. The methodology is described in detail allowing transfer of the scheme from flow to other applications such as advection-diffusion problems or heat transport. This is an important methodological development and provides new possibilities for the simulation of complex material behavior.

This article is structured as follows: In 2 we first give a short introduction to the Boltzmann equation, before introducing the anisotropic modification of the equilibrium distribution and performing the discretization of the velocity space. In 3 we perform the Chapman-Enskog analysis via perturbation in orders of Knudsen numbers yielding the governing macroscopic equations. We end with a short conclusion and an outlook on possible future extensions in 4. For ease of reading, most of the required mathematical background is given in the Appendix.

2 Methods↩︎

2.1 The Boltzmann Equation↩︎

Lattice-Boltzmann methods (LBM) aim to numerically solve the Boltzmann equation \[\frac{\mathop{}\!\mathrm{d}{f}}{\mathop{}\!\mathrm{d}{t}} = \Omega,\] which expresses the time evolution of a particle distribution \(f\) under the influence of interparticle collisions described by an operator \(\Omega\). In the most general sense a particle distribution function \(f\) maps the density distribution of a large number of non-interacting non-relativistic particles as a function of the microscopic particle velocities \(\boldsymbol{v}\), time \(t\), and spatial position \(\boldsymbol{x}\). The distribution is normalized such that the zeroth order moment yields the mass density \(\rho\) and the first order moment the momentum density \(\rho u_\alpha\) as time-dependent fields, i.e. \[\int_{-\infty}^{\infty}f(t,\boldsymbol{x},\boldsymbol{v}) \mathrm{d}^{d}v = \rho(t,\boldsymbol{x}), \label{eq:A1}\tag{1}\] \[\int_{-\infty}^{\infty}f(t,\boldsymbol{x},\boldsymbol{v}) v_\alpha\mathrm{d}^{d}v = \rho(t,\boldsymbol{x}) u_\alpha(t,\boldsymbol{x}). \label{eq:A2}\tag{2}\] Throughout this work, Greek indices run over the \(d\) spatial or velocity dimensions. Given enough time, inter-particle collisions will drive the system towards a local equilibrium state described by an equilibrium distribution \(f^\mathrm{eq}\) which depends only on the first two macroscopic moments and the microscopic particle velocities. The distribution will also yield the same moments to first order, i.e. \[\int_{-\infty}^{\infty}f^\mathrm{eq}(\rho,\boldsymbol{u},\boldsymbol{v}) \mathrm{d}^{d}v = \rho(t,\boldsymbol{x}), \label{eq:A3}\tag{3}\] \[\int_{-\infty}^{\infty}f^\mathrm{eq}(\rho,\boldsymbol{u},\boldsymbol{v}) v_\alpha\mathrm{d}^{d}v = \rho(t,\boldsymbol{x}) u_\alpha(t,\boldsymbol{x}). \label{eq:A4}\tag{4}\]

In local equilibrium, the space-time dependence enters only implicitly through the macroscopic fields. Here, and in the following, we only consider an isothermal system, where we can disregard the second order moment related to energy.

Usually, it is assumed that the interparticle collisions have no preferred direction and will hence even out any angular dependence of an initial state very rapidly, resulting in a isotropic and isothermal equilibrium described by a Maxwell-Boltzmann-type distribution \[\bar{f}^\mathrm{eq}(\rho,\boldsymbol{u},\boldsymbol{v}) = \rho(t,\boldsymbol{x}) \frac{1}{\sqrt{2\pi}^d} \exp\left(-\frac{[v_\alpha- u_\alpha(t,\boldsymbol{x})]\delta_{\alpha\beta}[v_\beta- u_\beta(t,\boldsymbol{x})]}{2} \right),\] where \(\delta_{\alpha\beta}\) is the Kronecker delta. Since we work in the isothermal model, we choose units in which the temperature constant is unity to simplify notation. This expression utilizes Einstein’s summation convention where we sum over repeated indices, which is applied throughout this work.

In the following, we derive an anisotropic LBM by skewing this distribution. We will develop the whole discretization framework based on this assumption and compute the governing macroscopic equations of a system described by this type of equilibrium.

2.2 Expansion of an Anisotropic Equilibrium Distribution↩︎

One of the simplest forms of a Maxwell-Boltzmann-like isothermal equilibrium distribution which introduces anisotropy can in suitable units and \(d\) dimensions be written as \[f^\mathrm{eq}(\rho,\boldsymbol{u},\boldsymbol{v}) = \rho(t,\boldsymbol{x}) \frac{1}{\sqrt{2\pi}^d} \exp \left(-\frac{[v_\alpha-u_\alpha(t,\boldsymbol{x}))] A_{\alpha\beta} [v_\beta-u_\beta(t,\boldsymbol{x})]}{2}\right). \label{eq:equil}\tag{5}\] Since we limit our investigation to the case where \(f^\mathrm{eq}\) remains Boltzmann-like, the expression in the exponential remains a quadratic form. Any quadratic form can be described by a unique symmetric matrix, which we call \(A\). To not pick up extra normalization factors in all expressions, we additionally require \(\det A=1\). Physically, \(A\) describes the strength of the directional dependence of the equilibrium velocity distribution.

The macroscopic moments of this new anisotropic equilibrium are still \[\int_{-\infty}^{\infty}f^\mathrm{eq}\mathrm{d}^{d}{v} = \rho(t,\boldsymbol{x}),\] \[\int_{-\infty}^{\infty}f^\mathrm{eq}v_\alpha\mathrm{d}^{d}{v} = \rho(t,\boldsymbol{x}) u_\alpha(t,\boldsymbol{x}),\] cf.5. The explicit treatment of a full exponential makes the theory laborious and any simulations computationally expensive. Since we only enforce to reproduce the first two moments, our physical observables, correctly, it is sufficient to approximate the equilibrium distribution. In this work, we limit the discussion to the single-relaxation time collision operator (BGK) [10], which recovers macroscopically the Navier-Stokes equations (cf.10) using an expansion up to second order in velocity for the isotropic case. An anisotropic equilibrium, however, benefits from third order expansion, as will we shown in 3. Using only second order, the anisotropic Navier-Stokes equations will include a discretization error depending on the anisotropy, see also 11.

To expand the exponential, LBM utilizes Hermite polynomial expansion because the resulting distribution has suitable properties for the following velocity space discretization. Details can be found in 6 and Ref.[1]. Expanded to third order, the equilibrium distribution can be approximated as \[f^\mathrm{eq}(\boldsymbol{v}) \simeq \omega(\boldsymbol{v}) \sum_{n=0}^3 \frac{1}{n!} \langle a^{(n)}(f^\mathrm{eq}) , H^{(n)} \rangle,\] where \(\omega\) is the generator of the Hermite polynomials \(H^{(n)}\), cf.6. The operation \(\langle\cdot,\cdot\rangle\) is the full contraction of the \(n\)-rank tensors \(H^{(n)}\) and \(a^{(n)}\). The expansion coefficients are projections onto the Hermite polynomials, i.e. \[a^{(n)}(f^\mathrm{eq}) = \int_{-\infty}^{\infty}f^\mathrm{eq}(\boldsymbol{v}) H^{(n)}(\boldsymbol{v}) \mathrm{d}^{d}v, \label{eq:expcoeff}\tag{6}\] which, to third order, results in \[\begin{align} \tilde{f}^\mathrm{eq}&(t,\boldsymbol{x},\boldsymbol{v}) = \omega(\boldsymbol{v}) \sum_{n=0}^3 \frac{1}{n!} \langle a^{(n)}(f^\mathrm{eq}) , H^{(n)} \rangle \\ &\qquad\quad\; = \omega(\boldsymbol{v}) \rho(t,\boldsymbol{x}) \left[ 1 + u_\alpha v_\alpha + \frac{1}{2} (u_\alpha u_\beta+ K_{\alpha\beta})(v_\alpha v_\beta- \delta_{\alpha\beta}) \right. \\ &\left. + \frac{1}{6}( u_\alpha u_\beta u_\gamma+ K_{\alpha\beta} u_\gamma+ K_{\beta\gamma} u_\alpha+ K_{\gamma\alpha} u_\beta) \left( v_\alpha v_\beta v_\gamma - v_\alpha\delta_{\beta\gamma} - v_\beta\delta_{\gamma\alpha} - v_\gamma\delta_{\alpha\beta} \right)\right], \end{align} \label{eq:3rd95order95equilibrium}\tag{7}\] where \(u_\alpha=u_\alpha(t,\boldsymbol{x})\), which we suppressed for the sake of readability. The matrix \(K\) is given by \[K_{\alpha\beta} = A^{-1}_{\alpha\beta} - \delta_{\alpha\beta}. \label{eq:K-term}\tag{8}\] Calculation of the first to moments confirms that truncation does not change relevant macroscopic properties. We find \[\int_{-\infty}^{\infty}\tilde{f}^\mathrm{eq}\mathrm{d}^{d}{v} = \int_{-\infty}^{\infty}\tilde{f}^\mathrm{eq}H^{(0)} \mathrm{d}^{d}{v} = a^{(0)}(f^\mathrm{eq}) = \rho, \label{eq:14}\tag{9}\] \[\int_{-\infty}^{\infty}\tilde{f}^\mathrm{eq}v_\alpha\mathrm{d}^{d}{v} = \int_{-\infty}^{\infty}\tilde{f}^\mathrm{eq}H^{(1)}_\alpha\mathrm{d}^{d}{v} = a^{(1)}_\alpha(f^\mathrm{eq}) = \rho u_\alpha, \label{eq:15}\tag{10}\] cf.@eq:eq:expcoeff and 39 in 6. It is not surprising that the relevant moments remain unchanged, since the Hermite polynomials are orthogonal. To a certain degree, they separate orders of velocity. Taking these integrals corresponds to projecting out the expansion coefficients \(a^{(n)}(f^\mathrm{eq})\).

2.3 Discretization of the Velocity Space↩︎

For calculation of the macroscopic quantities [eq:14,eq:15] we still have to compute the full velocity integrals to get the macroscopic quantities as moments of the distribution function. Just as in isotropic LBM, we can convert velocity integrals into finite sums through Gaussian quadrature. Here, we construct a weighted sum of orthogonal polynomials that are evaluated at specific points, referred to as nodes or abscissae. There are numerous variants of such quadratures and the choice mostly depends on the range of the interval to be evaluated. For the range \((-\infty,\infty)\) the required method is Gauss-Hermite quadrature. It shares its generator \(\omega\) with the Hermite polynomials and considerably simplifies the following treatment of the integrals. An overview of the method is available in 7 or Ref.[1].

In general, we encounter expressions including integrals of the form of 6 for the macroscopic moments. Note that in the following the function \(f^\mathrm{eq}\) can be considered the truncated one without affecting the accuracy of the scheme. Therefore, in the derivation below we refer to the truncated equilibrium distribution just as \(f^\mathrm{eq}\), if not otherwise specified.

We now utilize the relation of the Hermite expansion and the Gauss-Hermite quadrature and split \(f^\mathrm{eq}\) (7 ) into three parts: the normalization \(\rho\), the generator function \(\omega\), and the remaining polynomial \(\mathcal{P}\) which depends on \(\boldsymbol{v}\) and \(\boldsymbol{u}(t,\boldsymbol{x})\). This yields \[a^{(n)}(f^\mathrm{eq}) = \int_{-\infty}^{\infty}\omega(\boldsymbol{v})\rho(t,\boldsymbol{x}) \mathcal{P}(\boldsymbol{v},\boldsymbol{u}) H^{(n)}(\boldsymbol{v}) \mathrm{d}^{d}v.\] The quadrature procedure converts these integrals into sums by evaluating only a fixed set of points in the velocity space. Hence, the continuous velocity space essentially becomes discrete. The expansion coefficients are therefore given by \[a^{(n)}(f^\mathrm{eq}) = \sum_{i=1}^q w_i \rho(t,\boldsymbol{x})\mathcal{P}(\boldsymbol{v}_i)H^{(n)}(\boldsymbol{v}_i), \label{eq:quadrature}\tag{11}\] where \(q\) depends on the desired accuracy. Another convenient property of the Gauss-Hermite quadrature in contrast to other variants is that the expression remains exact for suitably large \(q\). The size of \(q\) depends on the order of the product \(\mathcal{P}H^{(n)}\).

To not stream off-lattice, spatial and velocity lattices need to overlap. But the nodes or abscissae, \(\boldsymbol{v}_i\), do not agree with a spatial standard lattice using \(\Delta x = 1\). It is hence standard practice in LBM to rescale the velocities to \[\boldsymbol{v}_i = \frac{\boldsymbol{c}_i}{c_\mathrm{s}},\quad \boldsymbol{u}= \frac{\boldsymbol{s}}{c_\mathrm{s}}.\] The constant factor \(c_\mathrm{s}\) could also be absorbed into the units. The updated equilibrium distribution reads \[\begin{align} f^\mathrm{eq}(\boldsymbol{c}_i) = \omega(\boldsymbol{c}_i) \rho \left[ 1 \right. &+ \left.\frac{c_{i\alpha} s_\alpha}{c_\mathrm{s}^2} + \frac{1}{2c_\mathrm{s}^4} (s_\alpha s_\beta+ c_\mathrm{s}^2 K_{\alpha\beta})(c_{i\alpha} c_{i\beta} - c_\mathrm{s}^2\delta_{\alpha\beta}) \right. \\ &+ \frac{1}{6c_\mathrm{s}^6} \left(s_\alpha s_\beta s_\gamma+ c_\mathrm{s}^2 K_{\alpha\beta} s_\gamma+ c_\mathrm{s}^2 K_{\beta\gamma} s_\alpha+ c_\mathrm{s}^2 K_{\gamma\alpha} s_\beta\right)\\ &\left.\times\left( c_{i\alpha} c_{i\beta} c_{i\gamma} - c_{i\alpha} c_\mathrm{s}^2 \delta_{\beta\gamma} - c_{i\beta} c_\mathrm{s}^2 \delta_{\gamma\alpha} - c_{i\gamma} c_\mathrm{s}^2 \delta_{\alpha\beta} \right) \right]. \end{align} \label{pnzkeidf}\tag{12}\] Aiming for a clean notation, we define the discrete equilibrium function \[f^\mathrm{eq}_i(t,\boldsymbol{x}) = \frac{w_i}{\omega(\boldsymbol{c}_i)} f^\mathrm{eq}(t,\boldsymbol{x},\boldsymbol{c}_i), \label{eq:feq95i}\tag{13}\] with the weights \(w_i\) included and perform the quadrature, 11 , yielding the physically relevant macroscopic moments as sums in the discretized velocity space as \[\rho = a^{(0)}(f^\mathrm{eq}) = \sum_{i=1}^q f^\mathrm{eq}_i, \label{eq:a0discrete}\tag{14}\] and \[\rho u_\alpha= a^{(1)}(f^\mathrm{eq}) = \sum_{i=1}^q f^\mathrm{eq}_i c_{i\alpha}. \label{eq:a1discrete}\tag{15}\]

Similarly, we can discretize the particle distribution function \(f\). Without repeating all the steps we just note that since we want to conserve mass and momentum, the distributions yield the same moments up to first order, cf.[eq:A1,eq:A2,eq:A3,eq:A4]. From this follows that the first two expansion coefficients must be identical, such that \[\rho = a^{(0)}(f^\mathrm{eq}) = a^{(0)}(f) = \sum_{i=1}^q f_i, \label{eq:a0discrete95f}\tag{16}\] and \[\rho u_\alpha= a^{(1)}(f^\mathrm{eq}) = a^{(1)}(f) = \sum_{i=1}^q f_i c_{i\alpha}, \label{eq:a1discrete95f}\tag{17}\] with \(f_i\) defined analogously to 13 .

Since the space-time discretization does not require modifications compared to the isotropic case, we will skip it here and refer the interested reader to Ref.[1] instead.

3 Chapman-Enskog Theory↩︎

The Chapman-Enskog [11] analysis will provide information about the macroscopic physical equations that are induced by the microscopic discrete formalism we developed above. A detailed but palatable derivation for the isotropic case can again be found in Ref.[1]. We will follow a similar derivation in the procedure below.

We want to derive the behavior of the anisotropic equilibrium function under single-relaxation time collision (BGK) [10], corresponding to the fully discrete Boltzmann equation \[f_i(t+\Delta t,x+\boldsymbol{c}_i \Delta t) - f_i(t,\boldsymbol{x}) = \frac{\Delta t}{\tau} [f^\mathrm{eq}_i(t,\boldsymbol{x}) - f_i(t,\boldsymbol{x})], \label{eq:LBE95BGK}\tag{18}\] where \(\tau\) is the relaxation time. Here, we performed the time discretization to first order only, since it is equivalent to second order assuming a suitable redefinition of \(f_i\). Some details about the expected form of the resulting Navier-Stokes equations for single-relaxation time collisions and their anisotropic modification are given in 10.

Performing the analysis on this collision operator, we first note that a general discrete distribution \(f_i\) deviates from the equilibrium by some amount \(\Delta f_i\), i.e. \[f_i = f^\mathrm{eq}_i + \Delta f_i \quad \Leftrightarrow \quad \Delta f_i = f_i - f^\mathrm{eq}_i. \label{eq:deltaf}\tag{19}\] This can be understood as the collection of particles consisting of an equilibrated part plus a (small) non-equilibrium perturbation. The macroscopic dynamic properties of the fluid fully depend on this perturbation. Similar to other perturbation theories, the perturbation can be expanded in orders of some ordering-parameter, in our case the Knudsen number \[\mathrm{Kn} = \frac{\lambda}{L},\] with \(\lambda\) and \(L\) being the mean free path of fluid particles and the characteristic physical length scale of the problem, respectively.

We expand the distribution as \[f_i = \sum_{n=0}^\infty \epsilon^n f_i^{(n)} = f_i^{(0)} + \sum_{n=1}^\infty \epsilon^n f_i^{(n)} = f^\mathrm{eq}_i + \sum_{n=1}^\infty \epsilon^n f_i^{(n)} \stackrel{(\ref{eq:deltaf})}{=} f^\mathrm{eq}_i + \Delta f_i, \label{eq:perturbation95expansion}\tag{20}\] where powers of expansion coefficients \(\epsilon\) determine the order in Kn, i.e. \[\frac{f^{(n)}}{f^\mathrm{eq}} = \mathcal{O}(\mathrm{Kn}^n).\]

By combining [eq:a0discrete,eq:a1discrete,eq:a0discrete_f,eq:a1discrete_f] it follows that \[\sum_i f^\mathrm{eq}_i = \sum_i f_i \quad \stackrel{(\ref{eq:deltaf})}{\Longrightarrow} \quad \sum_i \Delta f_i = 0 \label{eq:Dfzero0}\tag{21}\] and \[\sum_i f^\mathrm{eq}_i \boldsymbol{c}_i = \sum_i f_i \boldsymbol{c}_i \quad \stackrel{(\ref{eq:deltaf})}{\Longrightarrow} \quad \sum_i \Delta f_i \boldsymbol{c}_i = \boldsymbol{0}. \label{eq:Dfzero1}\tag{22}\] The expansion in 20 makes \(\Delta f_i\) a polynomial in \(\epsilon\). By comparison of coefficients [eq:Dfzero0,eq:Dfzero1] imply even stricter relations, namely \[\sum_i f_i^{(n)} = 0 \quad\forall n \geq 1, \label{eq:prop1}\tag{23}\] and \[\sum_i f_i^{(n)} \boldsymbol{c}_i = \boldsymbol{0} \quad\forall n \geq 1. \label{eq:prop2}\tag{24}\] Next, we expand the whole left hand side of 18 into a Taylor series up to second order, cf.@eq:eq:Taylor95expansion in 8, and apply 19 to the right hand side, which yields \[\Delta t(\partial_t + c_{i\alpha}\partial_\alpha) f_i + \frac{\Delta t^2}{2}(\partial_t + c_{i\alpha}\partial_\alpha)^2 f_i + \mathcal{O}(\Delta t^3) = -\frac{\Delta t}{\tau}\Delta f_i. \label{eq:step1}\tag{25}\] The treatment of a second derivative would be complicated. To eliminate it, we apply the operator \(\frac{\Delta t}{2}(\partial_t + c_{i\alpha}\partial_\alpha)\) to both sides of 25 , i.e. \[\frac{\Delta t^2}{2}(\partial_t + c_{i\alpha}\partial_\alpha)^2 f_i + \mathcal{O}(\Delta t^3) = -\frac{\Delta t^2}{2\tau}(\partial_t + c_{i\alpha}\partial_\alpha) \Delta f_i, \label{eq:step2}\tag{26}\] and subtract the resulting 26 from 25 , yielding \[\Delta t(\partial_t + c_{i\alpha}\partial_\alpha) f_i = -\frac{\Delta t}{\tau}\Delta f_i + \frac{\Delta t^2}{2\tau}(\partial_t + c_{i\alpha}\partial_\alpha) \Delta f_i, \label{eq:step3}\tag{27}\] where we neglected the \(\mathcal{O}(\Delta t^3)\) terms to simplify the expression.

The idea is to get a set of multiple equations, each corresponding to the next correction when considering higher and higher orders of Kn. The relevant order for the Navier-Stokes-regime is \(\mathcal{O}(\mathrm{Kn}^2)\).[1] Hence, we expand the time derivative in a similar fashion as \(\Delta f_i\) in 20 up to second order, i.e. \[\partial_t f_i = (\epsilon \partial_t^{(1)} + \epsilon^2 \partial_t^{(2)})f_i.\] The spatial derivative is sufficient at first order, but we want to introduce the corresponding notation \[c_{i\alpha} \partial_\alpha f_i = \epsilon c_{i\alpha} \partial_\alpha^{(1)}f_i.\] Applying these definitions to 27 yields \[\Delta t \left[ \left( \epsilon \partial_t^{(1)} + \epsilon^2 \partial_t^{(2)} \right) + \epsilon c_{i\alpha} \partial_\alpha^{(1)} \right] f_i = -\frac{\Delta t}{\tau}\Delta f_i + \frac{\Delta t^2}{2\tau} \left[ \left( \epsilon \partial_t^{(1)} + \epsilon^2 \partial_t^{(2)} \right) + \epsilon c_{i\alpha} \partial_\alpha^{(1)} \right] \Delta f_i.\] Using 20 we expand \(\Delta f_i\) and \(f_i\) to second order and split the result into two equations, one for each order of \(\epsilon\), yielding \[\begin{align} \epsilon \left( \partial_t^{(1)} + c_{i\alpha} \partial_\alpha^{(1)} \right) f^\mathrm{eq}_i &= - \frac{\epsilon}{\tau} f_i^{(1)}, \tag{28} \\ \epsilon^2 \partial_t^{(2)} f^\mathrm{eq}_i + \epsilon^2 \left(1 - \frac{\Delta t}{2\tau} \right) \left( \partial_t^{(1)} + c_{i\alpha}\partial_\alpha^{(1)} \right) f_i^{(1)} &= - \frac{\epsilon^2}{\tau} f_i^{(2)}. \tag{29} \end{align}\] Again, the macroscopic behavior is determined by the first two moments, zeroth and first order. And the second order moment is essentially a correction term. From the Gauss-Hermite quadrature we know that taking moments in the discrete velocity space corresponds to multiplying expressions by power of \(\boldsymbol{c}_i\) before summing over \(i\). We can take moments of the full equations by applying the operations to both sides. Considering [eq:prop1,eq:prop2], the zeroth moments of [eq:CE1,eq:CE2] are \[\epsilon \partial_t^{(1)} \rho + \epsilon \partial_\alpha^{(1)}(\rho u_\alpha) = 0, \label{eq:1st95order95mass}\tag{30}\] and \[\epsilon^2 \partial_t^{(2)} \rho = 0,\] respectively. Adding both equations and reversing the derivative expansions yields the governing equation for the zeroth order moment (mass), i.e.the continuity equation \[\partial_t \rho + \partial_\alpha(\rho u_\alpha) = 0.\]

With moment tensors \(M_{\alpha_1...\alpha_n}\) defined analogously to 9, the first order moments of [eq:CE1,eq:CE2] are \[\epsilon \partial_t^{(1)}(\rho u_\alpha) + \epsilon \partial_\beta^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) = 0 \label{eq:1st95order95momentum}\tag{31}\] and \[\epsilon^2 \partial_t^{(2)}(\rho u_\alpha) + \epsilon^2 \partial_\beta^{(1)} \left( 1 - \frac{\Delta t}{2\tau} \right) M_{\alpha\beta}(f_i^{(1)}) = 0,\] respectively. Again, adding these and recombining the expansions yields the governing equation for the first order moment (momentum), i.e. \[\partial_t(\rho u_\alpha) + \partial_\beta M_{\alpha\beta} (f^\mathrm{eq}_i) = -\epsilon \partial_\beta\left[ \left( 1 - \frac{\Delta t}{2\tau} \right) M_{\alpha\beta}(f^{(1)}_i) \right]. \label{eq:CE3}\tag{32}\] The higher order moments of the equilibrium encountered in this equation can be evaluated according to 9. But up to this point the second order moment of the perturbation \(f^{(1)}_i\) is unknown. However, taking the second order moment of 28 results in a defining equation for it, i.e. \[\epsilon \partial_t^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) + \epsilon\partial_\gamma^{(1)} M_{\alpha\beta\gamma}(f^\mathrm{eq}_i) = -\frac{\epsilon}{\tau} M_{\alpha\beta}(f_i^{(1)}). \label{eq:CE4}\tag{33}\] This equation provides us with the unknown moment in terms of second and third order moments of the equilibrium, which we can compute. Corresponding calculations can be found in 11. The needed equilibrium moments can again be found in 9. Using 33 , we find the perturbation moment \[M_{\alpha\beta}(f_i^{(1)}) = - \tau \rho c_\mathrm{s}^2 \left( A^{-1}_{\beta\gamma} \partial_\gamma^{(1)} s_\alpha + A^{-1}_{\alpha\gamma} \partial_\gamma^{(1)} s_\beta \right).\] Inserting the second order moment of the perturbation into 32 we arrive at \[\partial_t(\rho u_\alpha) + \partial_\beta M_{\alpha\beta} (f^\mathrm{eq}_i) = \epsilon \partial_\beta\left[ \rho c_\mathrm{s}^2 \left( \tau - \frac{\Delta t}{2} \right) \left( A^{-1}_{\beta\gamma} \partial_\gamma^{(1)} s_\alpha + A^{-1}_{\alpha\gamma} \partial_\gamma^{(1)} s_\beta \right) \right].\] After inserting the expression for the equilibrium moment and reversing the derivative expansion we finally arrive at \[\partial_t(\rho s_\alpha) + \partial_\beta\left( \rho c_\mathrm{s}^2 A^{-1}_{\alpha\beta} + \rho s_\alpha s_\beta \right) = \partial_\beta\left[ \rho c_\mathrm{s}^2 \left( \tau - \frac{\Delta t}{2} \right) \left( A^{-1}_{\beta\gamma} \partial_\gamma s_\alpha + A^{-1}_{\alpha\gamma} \partial_\gamma s_\beta \right) \right]. \label{eq:NSprefinal}\tag{34}\] By defining a shear viscosity \[\eta = \rho c_\mathrm{s}^2 \left( \tau - \frac{\Delta t}{2} \right)\] and applying the equation of state \[p = \rho c_\mathrm{s}^2,\] where \(p\) is the pressure of the fluid, we are left with a more familiar form of 34 , i.e. \[\partial_t(\rho s_\alpha) + \partial_\beta( \rho s_\alpha s_\beta) = - A^{-1}_{\alpha\beta} \partial_\beta p + \partial_\beta\left[ \eta \left( A^{-1}_{\beta\gamma} \partial_\gamma s_\alpha + A^{-1}_{\alpha\gamma} \partial_\gamma s_\beta \right) \right]. \label{eq:NSprefinal2}\tag{35}\] With our derivation, we imply that the bulk viscosity \[\zeta = \frac{2}{3} \eta,\] hence, two terms otherwise present in the Navier-Stokes equations cancel for us.

The reader might also have expected an anisotropic transformation acting on the macroscopic velocity \(\boldsymbol{s}\). This is in fact the case, although not explicitly visible. By defining the equilibrium as we have in 5 , \(A\) is already applied to \(\boldsymbol{s}\), hence \(\boldsymbol{s}\) is already anisotropically skewed by definition and requires no additional transformation.

4 Conclusion and Outlook↩︎

We were able to derive a consistent lattice Boltzmann framework for anisotropic flow problems using an anisotropic equilibrium distribution. We provided the proof through Chapman-Enskog theory that the governing macroscopic equations are correctly derived as being an anisotropic version of the Navier-Stokes equations. Important properties of LBM, like locality of the collision operator or isotropy of the space and velocity lattices are retained.

An interesting application of this model in future work is a volume averaged [12], [13] composite collision [14] where BGK and bounce back collision are combined in an anisotropic scheme to model fluid flow in unresolved anisotropic porous media. Additionally, the modification to anisotropic advection-diffusion would come with interesting applications and is relatively straight-forward to perform. A related topic is the anisotropic equilibrium distribution itself and how it might be measured or computed. It will be worth investigating how known LBM forcing schemes can be incorporated into the anisotropic model and if modifications are necessary.

Acknowledgements:

The presented research contributes to the Center for Electrochemical Energy Storage Ulm-Karlsruhe (CELEST).

Author Contributions:

Conceptualization: BK, JW. Methodology: BK. Validation: BK, JW. Formal analysis: BK. Investigation: BK. Resources: BK. Writing (Original Draft): BK. Writing (Review & Editing): BK, JW, AL, TD. Supervision: AL, TD. Project administration: AL, TD. Funding acquisition: AL, TD. All authors confirm that they have read the final manuscript and agree to its publication.

Conflicts of Interest:

The authors declare no conflicts of interest.

Abbreviations↩︎

We used the following abbreviations in the text:

Appendix↩︎

5 Useful Integral Identities↩︎

The integrals in 2 can be solved with these helpful identities: \[\int_{-\infty}^{\infty}\exp\left( -\frac{(v_\alpha- u_\alpha)A_{\alpha\beta}(v_\beta- u_\beta)}{2} \right) \mathrm{d}^{d}v = \sqrt{\frac{(2\pi)^d}{\det A}},\] \[\int_{-\infty}^{\infty}\exp\left( -\frac{(v_\alpha- u_\alpha)A_{\alpha\beta}(v_\beta- u_\beta)}{2} \right) v_\gamma\mathrm{d}^{d}v = \sqrt{\frac{(2\pi)^d}{\det A}} u_\gamma,\] \[\int_{-\infty}^{\infty}\exp\left( -\frac{(v_\alpha- u_\alpha)A_{\alpha\beta}(v_\beta- u_\beta)}{2} \right) v_\gamma v_\varepsilon\mathrm{d}^{d}v = \sqrt{\frac{(2\pi)^d}{\det A}} (A^{-1}_{\gamma\varepsilon} + u_\gamma u_\varepsilon),\] and finally \[\int_{-\infty}^{\infty}\exp\left( -\frac{(v_\alpha- u_\alpha)A_{\alpha\beta}(v_\beta- u_\beta)}{2} \right) v_\gamma v_\varepsilon v_\sigma\mathrm{d}^{d}v = \sqrt{\frac{(2\pi)^d}{\det A}} \left( A^{-1}_{\gamma\varepsilon} u_\sigma + A^{-1}_{\varepsilon\sigma} u_\gamma + A^{-1}_{\sigma\gamma} u_\varepsilon + u_\gamma u_\varepsilon u_\sigma \right).\]

6 Hermite polynomial expansion↩︎

Hermite polynomials of order \(n\) \[H^{(n)}(\boldsymbol{v}) = \frac{(-1)^n}{\omega(\boldsymbol{v})}\nabla^{(n)}\omega(\boldsymbol{v}). \label{eq:hermite95polynimals}\tag{36}\] are generated using the generating function¹ \[\omega(\boldsymbol{v}) = \frac{1}{\sqrt{2\pi}^d} \exp \left( -\frac{v_\alpha v_\alpha}{2} \right). \label{eq:generator}\tag{37}\]

The object \(\nabla^{(n)}\) is a tensor of rank \(n\) which applies \(n\) derivatives in the order of the indices, i. e. \[\nabla^{(n)} = \nabla_{\alpha_1...\alpha_n} = \frac{\partial}{\partial v_{\alpha_1}} ... \frac{\partial}{\partial v_{\alpha_n}},\] which makes \(H^{(n)}\) also a tensor of rank \(n\).

With this, any function \(g(\boldsymbol{v})\) can be expanded in the basis of Hermite polynomial as \[g(\boldsymbol{v}) = \omega(\boldsymbol{v}) \sum_{n=0}^\infty \frac{1}{n!} \langle a^{(n)}(g) , H^{(n)} \rangle.\] The expansion coefficients \(a^{(n)}(g)\) are found by projecting \(g(\boldsymbol{v})\) onto the \(n\)th basis polynomial, i.e. \[a^{(n)}(g) = \int_{-\infty}^{\infty}g(\boldsymbol{v}) H^{(n)}(\boldsymbol{v}) \mathrm{d}^{d}v. \label{eq:expansion95coefficients}\tag{38}\] The operation \(\langle\cdot,\cdot\rangle\) is the full contraction of the tensors, i. e. \[\langle a^{(n)}(g) , H^{(n)} \rangle = a^{(n)}_{\alpha_1...\alpha_n}(g) H^{(n)}_{\alpha_1...\alpha_n} .\]

From 36 we compute the first four Hermite polynomials \[\begin{align} &H^{(0)}(\boldsymbol{v}) = 1,\\ &H^{(1)}(\boldsymbol{v}) = v_\alpha,\\ &H^{(2)}(\boldsymbol{v}) = v_\alpha v_\beta- \delta_{\alpha\beta},\\ &H^{(3)}(\boldsymbol{v}) = v_\alpha v_\beta v_\gamma - \delta_{\alpha\beta} v_\gamma - \delta_{\beta\gamma} v_\alpha - \delta_{\gamma\alpha} v_\beta. \end{align} \label{eq:hermite95polynomials}\tag{39}\] Using the integral identities in 5, the expansion coefficients for \(f^\mathrm{eq}\) as defined in 5 are \[\begin{align} a^{(0)}(f^\mathrm{eq}) &= \rho, \\ a^{(1)}_\alpha(f^\mathrm{eq}) &= \rho u_\alpha, \\ a^{(2)}_{\alpha\beta}(f^\mathrm{eq}) &= \rho (u_\alpha u_\beta+ K_{\alpha\beta}), \\ a^{(3)}_{\alpha\beta\gamma}(f^\mathrm{eq}) &= \rho ( u_\alpha u_\beta u_\gamma + K_{\alpha\beta} u_\gamma + K_{\beta\gamma} u_\alpha + K_{\gamma\alpha} u_\beta), \end{align}\] where we define \[K_{\alpha\beta} = A^{-1}_{\alpha\beta} - \delta_{\alpha\beta}.\]

7 Gauss-Hermite Quadrature↩︎

Lets assume we need to evaluate an integral over the whole space in which the integrated function is a product of the Hermite generator \(\omega\) and an arbitrary 1D polynomial \(P^{(N)}(x)\) of order \(N\). Such integrals can be approximated using a weighted finite sum of the polynomial, such that \[\int_{-\infty}^{\infty}\omega(x) P^{(N)}(x) \mathrm{d}x = \sum_{i=1}^n w_i P^{(N)}(x_i).\] The specific value of \(n\) depends on the desired accuracy. In general \(N \leq 2n-1\). A nice property of the Gauss-Hermite quadrature – which sets it apart from most of the other quadrature variants – is that it becomes exact for sufficiently large \(n\), namely when \(N=2n-1\) or \(n=(N+1)/2\). The abscissae \(x_i\) are the roots of the \(n\)th order Hermite polynomial and the weights can be computed from one order lower via \[w_i = \frac{n!}{[n H^{(n-1)}(x_i)]^2}. \label{eq:weights}\tag{40}\] For example, a polynomial of order \(N=5\) requires \(n=3\) probing points. Hence, we need the third order Hermite polynomial \(H^{(3)}\) to compute the abscissae and the second order polynomial \(H^{(2)}\) for the weights. They can be derived from 36 by setting all the indices equal, such that \[H^{(3)}(x) = x^3 - 3x,\quad H^{(2)}(x) = x^2 -1.\] The roots of \(H^{(3)}\) are \[x_1 = 0, \quad x_{2,3} = \pm \sqrt{3},\] and using 40 yields the corresponding weights \[w_1 = \frac{2}{3}, \quad w_{2,3} = \frac{1}{6}.\]

The extension to \(d\) dimensions is straight-forward because of the two following properties of the generator and an arbitrary polynomial of order \(N\), \(P^{(N)}\): \[\begin{align} \omega(\boldsymbol{x}) = \prod_{\alpha=1}^d \omega(x_\alpha), \\ P^{(N)}(\boldsymbol{x}) = \sum_{ \{\boldsymbol{k}\} } a_{\boldsymbol{k}} \prod_{\alpha=1}^d x_\alpha^{N_\alpha}, \end{align}\] where the multi-index \(\boldsymbol{k}\) runs over the set \(\{\boldsymbol{k}\} = \{ (N_1,...,N_d) \in \mathbb{N}^d \mid \sum_\alpha N_\alpha\leq N \}\).

That way, the \(d\)-dimensional expression only includes 1D integrals, i. e. \[\int_{-\infty}^{\infty}\omega(\boldsymbol{x}) P^{(N)}(\boldsymbol{x}) \mathrm{d}^{d}\boldsymbol{x}= \sum_{\{\boldsymbol{k}\}} a_{\boldsymbol{k}} \prod_{\alpha=1}^d \int_{-\infty}^{\infty}\omega(x_\alpha) x_\alpha^{N_\alpha} \mathrm{d}x_\alpha, \label{eq:quadrature95general1}\tag{41}\] and the quadrature becomes \[\sum_{\{\boldsymbol{k}\}} a_{\boldsymbol{k}} \prod_{\alpha=1}^d \int_{-\infty}^{\infty}\omega(x_\alpha) x_\alpha^{N_\alpha} \mathrm{d}x_\alpha = \sum_{\{\boldsymbol{k}\}} a_{\boldsymbol{k}} \prod_{\alpha=1}^d \sum_{i_\alpha=1}^{n_\alpha} w_{i_\alpha} x_{i_\alpha}^{N_\alpha}. \label{eq:quadrature95general2}\tag{42}\] With the assumption that all dimensions are expanded to the same degree, i.e.\(n_\alpha=n\), one can immediately construct the D2Q9 and D3Q27 lattices. Through symmetry arguments, cf.Appendix in Ref.[15], one can reduce the lattice for \(d=3\) to D3Q19 without loss of precision in the Navier-Stokes regime. By mapping the products of weights to a single running index, \(w_{i1}...w_{id} \mapsto w_i\), one can rewrite 42 with a single sum \[\int_{-\infty}^{\infty}\omega(\boldsymbol{x}) P^{(N)}(\boldsymbol{x}) \mathrm{d}^{d}\boldsymbol{x}= \sum_{i=1}^q w_i P^{(N)}(\boldsymbol{x}_i)\] for the general lattice D\(d\)Q\(q\).

8 Series Expansion of the Discrete Boltzmann Equation↩︎

Let \(g:\mathbb{R}^d \to \mathbb{R}\) an \(N\)-times continuously differentiable function. Then \[g(x_\alpha+ h_\alpha) = g(x_\alpha) + \sum_{j=1}^N \frac{1}{j!} (h_\beta\partial_\beta)^j g(x_\alpha) + R,\] where \(R\) is a remainder of order \(N+1\). Note that the \(\beta\)-term is summed over because of the repeated index such that it corresponds to a inner product. With the definitions \(g(x_\alpha)=f_i(t,\boldsymbol{x})\), \(x_\alpha= (t,\boldsymbol{x})\), and \(h_\alpha=(\Delta t, \boldsymbol{c}_i\Delta t)\) we can write the discrete Boltzmann equation as \[f_i(t+\Delta t, \boldsymbol{x}+\boldsymbol{c}_i\Delta t) - f_i(t,\boldsymbol{x}) = \sum_{j=1}^N \frac{\Delta t^j}{j!}(\partial_t + c_{i\beta}\partial_\beta)^j f_i(t,\boldsymbol{x}) = \Omega_i \Delta t. \label{eq:Taylor95expansion}\tag{43}\]

9 Higher Order Moments of the Equilibrium↩︎

Although the physical observables are given by the zeroth and first order moments of the distribution function, the Chapman-Enskog analysis requires additionally the second and third order moments of the equilibrium function. For the continuous case, those can be computed using the integral identities in 5 and given by \[\mathcal{M}_{\alpha\beta}(f^\mathrm{eq}) = \int_{-\infty}^{\infty}f^\mathrm{eq}(\rho,\boldsymbol{u},\boldsymbol{v}) v_\alpha v_\beta\mathrm{d}^{d}{v} = \rho \left( A^{-1}_{\alpha\beta} + u_\alpha u_\beta \right),\] and \[\mathcal{M}_{\alpha\beta\gamma}(f^\mathrm{eq}) = \int_{-\infty}^{\infty}f^\mathrm{eq}(\rho,\boldsymbol{u},\boldsymbol{v}) v_\alpha v_\beta v_\gamma\mathrm{d}^{d}{v} = \rho \left( A^{-1}_{\alpha\beta} u_\gamma + A^{-1}_{\beta\gamma} u_\alpha + A^{-1}_{\gamma\alpha} u_\beta + u_\alpha u_\beta u_\gamma \right),\]

In the discrete velocity space, we have \[M_{\alpha\beta}(f^\mathrm{eq}_i) = \sum_i f^\mathrm{eq}_i v_{i\alpha} v_{i\beta} = \rho \left( A^{-1}_{\alpha\beta} + u_\alpha u_\beta \right),\] and \[M_{\alpha\beta\gamma}(f^\mathrm{eq}_i) = \sum_i f^\mathrm{eq}_i v_{i\alpha} v_{i\beta} v_{i\gamma} = \rho \left( A^{-1}_{\alpha\beta} u_\gamma + A^{-1}_{\beta\gamma} u_\alpha + A^{-1}_{\gamma\alpha} u_\beta + u_\alpha u_\beta u_\gamma \right).\] Using the compatible velocity lattice, where \(\boldsymbol{v}_i = \boldsymbol{c}_i / c_\mathrm{s}\) and \(\boldsymbol{u}= \boldsymbol{s}/ c_\mathrm{s}\), we find \[M_{\alpha\beta}(f^\mathrm{eq}_i) =\rho \left( c_\mathrm{s}^2 A^{-1}_{\alpha\beta} + s_\alpha s_\beta\right),\] and \[M_{\alpha\beta\gamma}(f^\mathrm{eq}_i) = \rho \left( c_\mathrm{s}^2 A^{-1}_{\alpha\beta} s_\gamma + c_\mathrm{s}^2 A^{-1}_{\beta\gamma} s_\alpha + c_\mathrm{s}^2 A^{-1}_{\gamma\alpha} s_\beta + s_\alpha s_\beta s_\gamma \right).\]

The reader should note that the moments of a full and truncated equilibrium distributions will be identical up to the order of expansion. Since we expanded \(f^\mathrm{eq}\) to third order, all the above moments remain identical. If we had expanded only to second order, the third moment of the second order equilibrium \(\hat{f^\mathrm{eq}}\) would lack the anisotropic correction and would therefore read \[\hat{M}_{\alpha\beta\gamma}(\hat{f^\mathrm{eq}}_i) = \rho \left( c_\mathrm{s}^2 \delta_{\alpha\beta} s_\gamma + c_\mathrm{s}^2 \delta_{\beta\gamma} s_\alpha + c_\mathrm{s}^2 \delta_{\gamma\alpha} s_\beta + s_\alpha s_\beta s_\gamma \right).\]

10 Navier-Stokes Equations for Single-relaxation Time Collisions↩︎

The force-free isotropic Navier-Stokes equations for a fluid of density \(\rho\) and a macroscopic velocity \(\boldsymbol{u}\) under pressure \(p\) can be written as \[\partial_t(\rho u_\alpha) + \partial_\beta( \rho u_\alpha u_\beta) = - \partial_\alpha p + \partial_\beta\left[ \eta \left( \partial_\beta u_\alpha + \partial_\alpha u_\beta \right) + \delta_{\alpha\beta} \left(\zeta - \frac{2}{3}\eta\right) \partial_\gamma u_\gamma \right]. \label{eq:appNS1}\tag{44}\] We introduced the shear viscosity \(\eta\) and the bulk or volume viscosity \(\zeta\). The Kronecker delta hints at the fact that the viscosities are actually not scalars but tensor quantities. Since we investigate the single-relaxation time (BGK) collision under isothermal assumption, the viscosities cannot be varied independently. The isentropic equation of state for pressure \(p\) and density \(\rho\) is \[p = p_0 \left(\frac{\rho}{\rho_0}\right)^\gamma,\] where we introduced the adiabatic index \(\gamma\). In general, the viscosities are related through \[\zeta = \eta \left( \frac{5}{3} - \gamma \right).\] The isothermal limit imposed by the equilibrium distribution 5 , however, is characterized by \(\gamma = 1\), such that \[\zeta = \frac{2}{3}\eta.\] This relation sets the last term in 44 to zero, and the isothermal Navier-Stokes equations read \[\partial_t(\rho u_\alpha) + \partial_\beta( \rho u_\alpha u_\beta) = - \partial_\alpha p + \partial_\beta\left[ \eta \left( \partial_\beta u_\alpha + \partial_\alpha u_\beta \right) \right]. \label{eq:appNS2}\tag{45}\]

The shear viscosity was already identified as a tensor. But the single-relaxation time collision leads to the relations \(\eta \propto \rho\) and \(p \propto \rho\), and hence, \(\eta \propto p\). We therefore assume that not only \(\eta\) but also \(p\) is implicitly a tensor that carries the anisotropy. We hence assume to arrive at Navier-Stokes equations where \(p\) and \(\eta\) both are anisotropic tensors, i.e. \[\partial_t(\rho u_\alpha) + \partial_\beta( \rho u_\alpha u_\beta) = - \partial_\beta\left(C_{\alpha\beta}p\right) + \partial_\beta\left[ \eta \left( C_{\beta\gamma}\partial_\gamma u_\alpha + C_{\alpha\gamma}\partial_\gamma u_\beta \right) \right], \label{akilvtzy}\tag{46}\] where \(C\) is a tensor related to \(A\). The Chapman-Enskog analysis in 3 shows that this is indeed the case, and that \(C \propto A^{-1}\).

More details on the thermodynamic relations can be found in [1] and the cited references therein.

11 Explicit Expression for the Perturbation Moment↩︎

During the Chapman-Enskog analysis we encounter the equation \[\partial_t(\rho u_\alpha) + \partial_\beta M_{\alpha\beta} (f^\mathrm{eq}_i) = -\epsilon \left( 1 - \frac{\Delta t}{2\tau} \right) \partial_\beta M_{\alpha\beta}(f^{(1)}_i). \label{eq:CE3b}\tag{47}\] and \[\epsilon \partial_t^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) + \epsilon \partial_\gamma^{(1)} M_{\alpha\beta\gamma}(f^\mathrm{eq}_i) = -\frac{\epsilon}{\tau} M_{\alpha\beta}(f_i^{(1)}), \label{eq:CE4b}\tag{48}\] which we want to compute here as an explicit expression of the observables \(\rho\) and \(\boldsymbol{u}=\boldsymbol{s}/c_\mathrm{s}\). The steps are provided in Ref.[1] and require only slight modifications.

Using the expression we found in 9, we write out the time derivative of the second order equilibrium moment, i.e. \[\partial_t^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) = \partial_t^{(1)} \left( c_\mathrm{s}^2 A^{-1}_{\alpha\beta} \rho + s_\alpha s_\beta\rho \right). \label{eq:F1}\tag{49}\] The second term on the right-hand side needs to be rewritten using a variation of the product rule. For some functions \(f\), \(g\), and \(h\) with a generic derivative of \(f\) denoted as \(f'\), we can write \[(fgh)' = f(gh)' + ghf',\] and since \[hf' = (fh)' - fh',\] we can write the first expression as \[(fgh)' = f(gh)' + g(fh)' - fgh'. \label{eq:product95rule}\tag{50}\] Applying the rule to 49 , we find \[\partial_t^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) = c_\mathrm{s}^2 A^{-1}_{\alpha\beta} \partial_t^{(1)} \rho + s_\alpha\partial_t^{(1)} (s_\beta\rho) + s_\beta\partial_t^{(1)} (s_\alpha\rho) - s_\alpha s_\beta\partial_t^{(1)} \rho \label{eq:F2}\tag{51}\] Next, we rewrite the first order balance equations for mass and momentum, [eq:1st_order_mass,eq:1st_order_momentum] as \[\partial_t^{(1)} \rho = - \partial_\alpha^{(1)}(s_\alpha\rho),\] and \[\partial_t^{(1)}(s_\alpha\rho) = - \partial_\beta^{(1)} \left( c_\mathrm{s}^2 A^{-1}_{\alpha\beta} \rho + s_\alpha s_\beta\rho \right),\] and plug them into 51 to eliminate some time derivatives, i.e. \[\begin{align} \begin{aligned} \partial_t^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) = &- c_\mathrm{s}^2 A^{-1}_{\alpha\beta} \partial_\gamma^{(1)}(s_\gamma\rho) - s_\alpha\partial_\gamma^{(1)} \left( c_\mathrm{s}^2 A^{-1}_{\beta\gamma} \rho + s_\beta s_\gamma\rho \right) \\ &- s_\beta\partial_\gamma^{(1)} \left( c_\mathrm{s}^2 A^{-1}_{\alpha\gamma} \rho + s_\alpha s_\gamma\rho \right) + s_\alpha s_\beta\partial_\gamma^{(1)}(s_\gamma\rho). \end{aligned} \label{eq:F3} \end{align}\tag{52}\] This can be rewritten as \[\begin{align} \begin{aligned} \partial_t^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) = &- c_\mathrm{s}^2 A^{-1}_{\alpha\beta} \partial_\gamma^{(1)}(s_\gamma\rho) - s_\alpha c_\mathrm{s}^2 A^{-1}_{\beta\gamma} \partial_\gamma^{(1)} \rho - s_\beta c_\mathrm{s}^2 A^{-1}_{\alpha\gamma} \partial_\gamma^{(1)} \rho \\ &+ s_\alpha s_\beta\partial_\gamma^{(1)}(s_\gamma\rho) - s_\alpha\partial_\gamma^{(1)} \left( s_\beta s_\gamma\rho \right) - s_\beta\partial_\gamma^{(1)} \left( s_\alpha s_\gamma\rho \right) , \end{aligned} \end{align}\] to apply the product rule in 50 in reverse to the second line, such that \[\partial_t^{(1)} M_{\alpha\beta}(f^\mathrm{eq}_i) = - c_\mathrm{s}^2 A^{-1}_{\alpha\beta} \partial_\gamma^{(1)}(s_\gamma\rho) - c_\mathrm{s}^2 A^{-1}_{\beta\gamma} s_\alpha\partial_\gamma^{(1)} \rho - c_\mathrm{s}^2 A^{-1}_{\alpha\gamma} s_\beta\partial_\gamma^{(1)} \rho -\partial_\gamma^{(1)} (s_\alpha s_\beta s_\gamma\rho).\]

Fortunately, the expression for the third order moment requires less effort and can be computed directly as \[\partial_\gamma^{(1)} M_{\alpha\beta\gamma}(f^\mathrm{eq}_i) = \partial_\gamma^{(1)} \left( c_\mathrm{s}^2 A^{-1}_{\alpha\beta} s_\gamma\rho + c_\mathrm{s}^2 A^{-1}_{\beta\gamma} s_\alpha\rho + c_\mathrm{s}^2 A^{-1}_{\alpha\gamma} s_\beta\rho + s_\alpha s_\beta s_\gamma\rho \right).\]

Going back to 48 we now found \[\begin{align} \begin{aligned} M_{\alpha\beta}(f_i^{(1)}) = \tau c_\mathrm{s}^2 A^{-1}_{\alpha\beta} \partial_\gamma^{(1)}(s_\gamma\rho) + &\tau c_\mathrm{s}^2 A^{-1}_{\beta\gamma} s_\alpha\partial_\gamma^{(1)} \rho + \tau c_\mathrm{s}^2 A^{-1}_{\alpha\gamma} s_\beta\partial_\gamma^{(1)} \rho + \tau \partial_\gamma^{(1)} (s_\alpha s_\beta s_\gamma\rho) \\ - & \tau \partial_\gamma^{(1)} \left( c_\mathrm{s}^2 A^{-1}_{\alpha\beta} s_\gamma\rho + c_\mathrm{s}^2 A^{-1}_{\beta\gamma} s_\alpha\rho + c_\mathrm{s}^2 A^{-1}_{\alpha\gamma} s_\beta\rho + s_\alpha s_\beta s_\gamma\rho \right), \end{aligned} \end{align}\] which, using the regular product rule, simplifies to \[\begin{align} \begin{aligned} M_{\alpha\beta}(f_i^{(1)}) &= - \tau \left[ c_\mathrm{s}^2 A^{-1}_{\beta\gamma} \left( \partial_\gamma^{(1)} (s_\alpha\rho) - s_\alpha\partial_\gamma^{(1)} \rho \right) + c_\mathrm{s}^2 A^{-1}_{\alpha\gamma} \left( \partial_\gamma^{(1)} (s_\beta\rho) - s_\beta\partial_\gamma^{(1)} \rho \right) \right] \\ &= - \tau \rho c_\mathrm{s}^2 \left[ A^{-1}_{\beta\gamma} \partial_\gamma^{(1)} s_\alpha + A^{-1}_{\alpha\gamma} \partial_\gamma^{(1)} s_\beta \right]. \end{aligned} \end{align}\] N. B.: This final expression is only that simple if we expand the equilibrium to the third order. Only then \(M_{\alpha\beta\gamma}\) also carries the \(A^{-1}\)-terms, cf.9. With a second order equilibrium the mixed order sum of equilibrium moments in 48 will include \(K\)-terms (8 ) as discretization errors. Additionally, we would see a \(\mathcal{O}(u^3)\) discretization error, which is also present in second order isotropic LBM and usually neglected.

References↩︎