# Cascaded Lattice Boltzmann Modeling and Simulations of Three-Dimensional
Non-Newtonian Fluid Flows

Saad Adam<EMAIL_ADDRESS>Farzaneh Hajabdollahi
<EMAIL_ADDRESS>Kannan N. Premnath
<EMAIL_ADDRESS>Department of Mechanical Engineering, University
of Colorado Denver, 1200 Larimer street, Denver, Colorado 80217 , U.S.A

###### Abstract

Non-Newtonian fluid flows, especially in three dimensions (3D), arise in
numerous settings of interest to physics. Prior studies using the lattice
Boltzmann method (LBM) of such flows have so far been limited to mainly to two
dimensions and used less robust collision models. In this paper, we develop a
new 3D cascaded LBM based on central moments and multiple relaxation times on
a three-dimensional, nineteen velocity (D3Q19) lattice for simulation of
generalized Newtonian (power law) fluid flows. The relaxation times of the
second order moments are varied locally based on the local shear rate and
parameterized by the consistency coefficient and the power law index of the
nonlinear constitutive relation of the power law fluid. Numerical validation
study of the 3D cascaded LBM for various benchmark problems, including the
complex 3D non-Newtonian flow in a cubic cavity at different Reynolds numbers
and power law index magnitudes encompassing shear thinning and shear
thickening fluids, are presented. Furthermore, numerical stability comparisons
of the proposed advanced LBM scheme against the LBM based on a single
relaxation time model are made. Numerical results demonstrate the accuracy,
second order grid convergence and improvements in stability of the 3D cascaded
LBM for simulation of 3D non-Newtonian flows of power law fluids.

###### keywords:

Lattice Boltzmann method, Central moments, Multiple relaxation times, Three
dimensional (3D) non-Newtonian flows, Power law fluids,

††journal: Journal Name

##  1 Introduction

In the last few decades, the lattice Boltzmann Method (LBM) has become a
preferred method for simulating complicated physical, chemical, and fluid
mechanics problems [1, 2]. It is a kinetic-based approach for fluid flow
computations. This is why it is especially useful for computing fluid flows
with multiple components involving interfacial dynamics, nonlinear
constitutive models and complex boundaries. Because the LBM lies on the scale
of mesoscopic level, challenges that have been encountered with using
conventional CFD methods are not found with the LBM [3, 4]. The lattice
Boltzmann Equation (LBE) can be obtained through several ways. One can derive
the LBE through discrete velocity models or the Boltzmann kinetic equation. In
addition, there are several approaches to derive Navier-Stokes equations from
the LBE. With Chapman-Enskog expansion being more popular among the rest,
other approaches include the asymptotic expansion, extended Taylor series
expansion and order of magnitude analysis.

The lattice Boltzmann methods are comprised of two fundamental steps, which
are the streaming step and collision step. The streaming step is the same in
various models of the LBM. However, because the collision step is more
complicated, researchers have been devoting considerable efforts into finding
the most suitable collision model for the lattice Boltzmann method. Among
those different collision models, the simplest and commonly used is the so-
called single-relaxation-time (SRT) model [5], which is based on the
Bhatnagar-Gross-Krook (BGK) approximation [6]. On the other hand, a multiple-
relaxation-time (MRT) model has been developed in order to improve numerical
stability [7]. MRT has been confirmed as a more stable collision model in
various problems. In the MRT collision model, various moments are relaxed to
their equilibrium states at different rates during the collision step. More
recently, a new class of collision model, the so-called cascaded LBM, has been
introduced by Geier et al. [8]. Later on, Asinari [9] reinterpreted this
approach based on relaxation to a generalized equilibrium. Premnath and
Banerjee [10] incorporated forcing terms in the cascaded LBM by the method of
central moments and systematically derived it and demonstrated its consistency
to the Navier-Stokes equations via a Chapman-Enskog expansion. In such an
approach, Galilean invariance is naturally enforced to the lattice Boltzmann
equation (LBE) based on the relaxation central moments. This involves
computing moments which are shifted by the macroscopic fluid velocity. In
another words, the moments are prescribed in a moving frame of reference.
Comparatively, the moments in the prior approaches are computed in a rest
frame of reference, which are termed as the raw moments.

In various recent studies [11, 12], the cascaded LBM based on central moments
and multiple relaxation times has been shown to be significantly more stable
when compared to the SRT collision model based LBM for the simulation of
Newtonian fluids flows. Three-dimensional central moment LBMs including the
forcing terms for D3Q15 and D3Q27 lattices were systematically derived in
[13]. A variant of such a formulation involving discrete moment equilibria
rather than the continuous Maxwellian equilibria is presented in [14]. A
preconditioned cascaded LBM to accelerate steady flow simulations with
improved accuracy was recently developed by [15]. Furthermore, the cascaded
LBM has recently been extended to simulate thermal convective flows in two-
dimensions (2D) [16, 17] and three-dimensions (3D) [18]. In addition, a
modified forcing formulation and implementation of the central moment LBM has
been presented in [19]. A systematic survey of various forcing approaches,
which were categorized according to their being either split or unsplit
formulations were made recently in [20], which also developed an efficient and
second order accurate method to include sources in the cascaded LBM based on
symmetric operator or Strang splitting.

Complex non-Newtonian fluid flows with nonlinear constitutive relations
represent an important class of flows with numerous applications in various
scenarios including in engineering, materials and food processing, and
geophysical processes [21, 22]. The extension of LBM to non-Newtonian fluid
flows has received significant but limited attention so far, which has been
reviewed in [23], and some more recent examples of such studies include the
works of [24, 25, 26]. While most of the prior studies have demonstrated the
applications of different LBM to 2D non-Newtonian fluid flows, very little
focus has been given to the development and validation of LBM, particularly
with using advanced collision models, for 3D non-Newtonian fluid flows
involving well defined benchmark problems. Indeed, due to complexity of
handling nonlinear rheological behavior via more general constitutive
relations with attendant numerical stability issues and the need for high grid
resolutions with relatively higher computational demands, there are only very
few studies in 3D even in the context of computational fluid dynamics (CFD)
using conventional numerical schemes. For the well-defined problem of 3D non-
Newtonian cubic cavity flow, Refs. [27, 28] presented some results at coarse
grid resolutions and, recently, Ref. [29] reported benchmark quality results
using a fractional step based finite volume method.

In the present work, we will present a new cascaded LBM in 3D for the
simulation of non-Newtonian flows represented by power law fluid based
rheology. In this regard, we construct a collision model based on central
moments and multiple relaxation times, on a three-dimensional nineteen
velocity (D3Q19) lattice. Prior versions of the ‘cascaded’ formulations of the
central moment LBM, where the changes in the higher moments depend
successively on those preceding lower moments, were developed for Newtonian
flows for the D3Q27 or D3Q15 lattice [8, 13]. In this work, we will present a
derivation of the cascaded LBM for the D3Q19 lattice, which is a compromise
between stability and efficiency, to handle nonlinear constitutive relations.
Here, the relaxation times for the second order moments that emulate the
nonlinear viscous behavior of the fluid are adjusted by the local shear rate
and parametrized by the consistency coefficient and the power law index of the
generalized Newtonian (i.e., power law) fluid. Expressions for computing the
shear rates based on estimating the components of the strain rate tensor
locally via the non-equilibrium moments will be provided for this formulation.
The implementation strategy for the sources corresponding to the applied body
forces for our 3D cascaded LBM will be according to the operator splitting
strategy presented in [20], and will not be discussed further here as the
present focus is on developing a modified collision term based on central
moments for non-Newtonian flows. Compared to the conventional schemes for CFD,
our 3D LB algorithm is local and is thus naturally suitable for efficient
implementation on large scale parallel computers. Moreover, our approach is
based on advanced formulation of the collision terms, which is expected to be
more stable when compared to other LB model for 3D non-Newtonian flows. We
will present a validation study of our 3D cascaded LB scheme for the non-
Newtonian flow in a channel, duct and a cubic cavity. In particular, for the
latter case involving complex nonlinear fluid motion with power-law rheology
in a cubic cavity, we will make direct comparisons of the velocity profiles
against the recent benchmark solution [29] for various Reynolds numbers and
power law index magnitudes encompassing both shear thinning and shear
thickening fluid behavior. We will also assess the order of accuracy of grid
convergence of our new 3D scheme for non-Newtonian flow simulations. Finally,
as an additional objective, we will present a comparative study of our new 3D
cascaded LBM formulation against the SRT formulation for simulation of non-
Newtonian flows and demonstrate an improvement in numerical stability achieved
by the former when compared to the latter approach.

This paper is organized as follows. In the next section (Sec. 2), we present
an overview of the macroscopic governing equations for generalized non-
Newtonian fluid flows and list the attendant non-linear constitutive relation.
A detailed exposition of the construction of the 3D cascaded LBM for the D3Q19
lattice is given in Sec. 3. This includes the changes of different moments
under cascaded collision and the specification of the local relaxation times
for the second order moments parameterized for the power law fluid, where the
shear rate is related to the non-equilibrium moments. Validation studies
involving various canonical non-Newtonian fluid flow problems, including those
in a cubic cavity are presented and discussed in Sec. 4. Section 5 makes a
numerical stability comparison of our 3D central moment LBM against a SRT-LBM
for the simulation of power law fluid flows in a shear driven problem. The
conclusions are finally summarized in Sec. 6.

##  2 Governing equations for three-dimensional generalized non-Newtonian
fluid flows

In the present work, we consider the simulation of three-dimensional (3D)
generalized Newtonian flows (GNF). The mathematical model, which is given in
terms of the continuity and momentum equations, can be written as

$\displaystyle\frac{\partial{\rho}}{\partial{t}}+\frac{\partial}{\partial{x_{j}}}{(\rho
u_{j})}$ $\displaystyle=0,$ (1a)
$\displaystyle\frac{\partial}{\partial{t}}{(\rho
u_{i})}+\frac{\partial}{\partial{x_{j}}}{(\rho u_{i}u_{j})}$
$\displaystyle=-\frac{\partial{P}}{\partial{x_{i}}}+\frac{\partial{\sigma_{ij}}}{\partial{x_{j}}}+F_{i},$
(1b)

where $\rho$ and $u_{i}$ are the local fluid density and the velocity field,
respectively, with $i\in{(x,y,z)}$. Here, $P$, $F_{i}$ and $\sigma_{ij}$
represent the pressure, Cartesian component of the imposed body force, and
deviatoric viscous stress tensor induced within the GNF, respectively. The
viscous stress tensor can be represented as

$\sigma_{ij}=\mu{(|\dot{\gamma}_{ij}|)}\dot{\gamma}_{ij},$ (2)

where $\dot{\gamma}_{ij}$ is the shear rate tensor. It is related to the
strain rate tensor $S_{ij}$ by

$\dot{\gamma}_{ij}=2S_{ij},\qquad
S_{ij}=\frac{1}{2}\left(\partial_{j}u_{i}+\partial_{i}u_{j}\right),$ (3)

and $\mu{(|\dot{\gamma}_{ij}|)}$ is assumed to be an effective or apparent
dynamic viscosity of the GNF, where the magnitude of the shear rate
$|\dot{\gamma}_{ij}|$ is related to the second invariant of the symmetric
strain rate tensor $S_{ij}$ as

$|\dot{\gamma}_{ij}|=\sqrt{2S_{ij}S_{ij}}.$ (4)

Here, the usual convention of the summation of the repeated indices is
assumed. Thus,

$|\dot{\gamma}_{ij}|=\sqrt{2\left[S^{2}_{xx}+S^{2}_{yy}+S^{2}_{zz}+2(S^{2}_{xy}+S^{2}_{yz}+S^{2}_{xz})\right]}.$
(5)

In this study, we consider the constitutive relations for the power law fluid
as a representative of a class of GNF for their numerical solution using the
3D cascaded LB method presented in the next section. The constitutive relation
for the power law fluid can be described mathematically as follows:

$\sigma_{ij}=\mu_{p}{|\dot{\gamma}_{ij}|}^{n-1}\dot{\gamma}_{ij},$ (6)

where the model parameter $\mu_{p}$ and the exponent $n$ are known as the
consistency coefficient and power law index of the fluid, respectively. Based
on the values of the power-law index $n$, the following three different types
of non-Newtonian fluids can be represented: $n<1$ corresponds to shear
thinning or pseudo-plastic fluid, whereas $n>1$ models shear thickening
fluids, and $n=1$ reduces to the Newtonian constitutive relation. Comparing
Eqs. (2) and (6), the effective viscosity for the power-law model is given by
the following expression

$\mu_{power}{|\dot{\gamma}_{ij}|}=\mu_{p}{|\dot{\gamma}_{ij}|}^{n-1}.$ (7)

##  3 Three-dimensional cascaded LBM for non-Newtonian fluid flows

We will formulate our non-Newtonian fluid flow solver based on a cascaded B
approach using the three-dimensional nineteen velocity (D3Q19) lattice, which
is shown in Fig. 1. The corresponding particle velocity
$\overrightarrow{e_{\alpha}}$ is represented as

$\bm{e}_{\alpha}=\left\\{\begin{array}[]{ll}{(0,0,0),}&{\alpha=0}\\\ {(\pm
1,0,0),(0,\pm 1,0),(0,0,\pm 1),}&{\alpha=1,\cdots,6}\\\ {(\pm 1,\pm 1,0),(\pm
1,0,\pm 1),(0,\pm 1,\pm 1),}&{\alpha=7,\cdots,18}\end{array}\right.$ (8)

Figure 1:  Three-dimensional, nineteen particle velocity (D3Q19) lattice.

In this work, the Greek and Latin subscripts are used for the particle
velocity directions and Cartesian coordinate directions, respectively. The
bare raw moments in the LBM are defined in terms of the distribution function
$f_{\alpha}$ as $\sum_{\alpha=0}^{18}e_{\alpha x}^{m}e_{\alpha y}^{n}e_{\alpha
z}^{p}f_{\alpha}$, where $m$, $n$ and $p$ are integers, and $(m+n+p)$
represents the order of the corresponding moment component. We use the Dirac’s
bra-ket notation in this paper to denote the basis vectors. For example, if
$\bra{\mathbf{a}}$ and $\ket{\mathbf{b}}$ represents a row vector and a column
vector, respectively, then their inner product is given by
$\braket{\mathbf{a}}{\mathbf{b}}$. Now, the nineteen non-orthogonal basis
vectors based on the combination of the monomials $e_{\alpha x}^{m}e_{\alpha
y}^{n}e_{\alpha z}^{p}$ in an ascending order supported by the D3Q19 lattice
can be listed as follows:

$\displaystyle\ket{T_{0}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{1}},\quad\ket{T_{1}}=\ket{\mathbf{e_{x}}},\quad\ket{T_{2}}=\ket{\mathbf{e_{y}}},\quad\ket{T_{3}}=\ket{\mathbf{e_{z}}},$
$\displaystyle\ket{T_{4}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}\mathbf{e_{y}}},\quad\ket{T_{5}}=\ket{\mathbf{e_{x}}\mathbf{e_{z}}},\quad\ket{T_{6}}=\ket{\mathbf{e_{y}}\mathbf{e_{z}}},$
$\displaystyle\ket{T_{7}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}^{2}-\mathbf{e_{y}}^{2}},\quad\ket{T_{8}}=\ket{\mathbf{e_{x}}^{2}-\mathbf{e_{z}}^{2}},\quad\ket{T_{9}}=\ket{\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}},$
$\displaystyle\ket{T_{10}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}\left(\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}\right)},\quad\ket{T_{11}}=\ket{\mathbf{e_{y}}\left(\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}\right)},$
$\displaystyle\ket{T_{12}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{z}}\left(\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}\right)},$
$\displaystyle\ket{T_{13}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}\left(\mathbf{e_{y}}^{2}-\mathbf{e_{z}}^{2}\right)},\quad\ket{T_{14}}=\ket{\mathbf{e_{y}}\left(\mathbf{e_{z}}^{2}-\mathbf{e_{x}}^{2}\right)},\quad\ket{T_{15}}=\ket{\mathbf{e_{z}}\left(\mathbf{e_{x}}^{2}-\mathbf{e_{y}}^{2}\right)},$
$\displaystyle\ket{T_{16}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}^{2}\mathbf{e_{y}}^{2}+\mathbf{e_{x}}^{2}\mathbf{e_{z}}^{2}+\mathbf{e_{y}}^{2}\mathbf{e_{z}}^{2}},\quad\ket{T_{17}}=\ket{\mathbf{e_{x}}^{2}\mathbf{e_{y}}^{2}+\mathbf{e_{x}}^{2}\mathbf{e_{z}}^{2}-\mathbf{e_{y}}^{2}\mathbf{e_{z}}^{2}},$
$\displaystyle\ket{T_{18}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}^{2}\mathbf{e_{y}}^{2}-\mathbf{e_{x}}^{2}\mathbf{e_{z}}^{2}}.$

Here, the components of the basis vectors
$\ket{\mathbf{1}},\ket{\mathbf{e_{x}}},\ket{\mathbf{e_{y}}}$ and
$\ket{\mathbf{e_{z}}}$ corresponding to the conserved moments based on which
the moment basis can be constructed are given by

$\displaystyle\ket{\mathbf{1}}$ $\displaystyle=$
$\displaystyle\left(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1\right)^{\dagger},$
$\displaystyle\ket{\mathbf{e_{x}}}$ $\displaystyle=$
$\displaystyle\left(0,1,-1,0,0,0,0,1,-1,1,-1,1,-1,1,-1,0,0,0,0\right)^{\dagger},$
$\displaystyle\ket{\mathbf{e_{y}}}$ $\displaystyle=$
$\displaystyle\left(0,0,0,1,-1,0,0,1,1,-1,-1,0,0,0,0,1,-1,1,-1\right)^{\dagger},$
$\displaystyle\ket{\mathbf{e_{z}}}$ $\displaystyle=$
$\displaystyle\left(0,0,0,0,0,1,-1,0,0,0,0,1,1,-1,-1,1,1,-1,-1\right)^{\dagger}.$

As in our previous work [13], in order to formulate a ‘cascaded’ LB
formulation, where the changes in higher order moments successively depend on
those of the preceding lower order ones, for the D3Q19 lattice, the above
$\ket{T_{j}}$ vectors, where $j=0,1,2,..,18,$ can be transformed into the
following equivalent set of orthogonal basis vectors via the standard Gram-
Schmidt procedure, which then read as

$\displaystyle\ket{K_{0}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{1}},\quad\ket{K_{1}}=\ket{\mathbf{e_{x}}},\quad\ket{K_{2}}=\ket{\mathbf{e_{y}}},\quad\ket{K_{3}}=\ket{\mathbf{e_{z}}},$
$\displaystyle\ket{K_{4}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}\mathbf{e_{y}}},\quad\ket{K_{5}}=\ket{\mathbf{e_{x}}\mathbf{e_{z}}},\quad\ket{K_{6}}=\ket{\mathbf{e_{y}}\mathbf{e_{z}}},$
$\displaystyle\ket{K_{7}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}^{2}-\mathbf{e_{y}}^{2}},\quad\ket{K_{8}}=\ket{\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}}-3\ket{\mathbf{e_{z}}^{2}},$
$\displaystyle\ket{K_{9}}$ $\displaystyle=$ $\displaystyle
19\ket{\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}}-30\ket{\mathbf{1}},$
$\displaystyle\ket{K_{10}}$ $\displaystyle=$ $\displaystyle
5\ket{\mathbf{e_{x}}\left(\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}\right)}-9\ket{\mathbf{e_{x}}},\quad\ket{K_{11}}=5\ket{\mathbf{e_{y}}\left(\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}\right)}-9\ket{\mathbf{e_{y}}},$
$\displaystyle\ket{K_{12}}$ $\displaystyle=$ $\displaystyle
5\ket{\mathbf{e_{z}}\left(\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}\right)}-9\ket{\mathbf{e_{z}}},$
$\displaystyle\ket{K_{13}}$ $\displaystyle=$
$\displaystyle\ket{\mathbf{e_{x}}\mathbf{e_{y}}^{2}-\mathbf{e_{x}}\mathbf{e_{z}}^{2}},\quad\ket{K_{14}}=\ket{\mathbf{e_{z}}^{2}\mathbf{e_{y}}-\mathbf{e_{y}}\mathbf{e_{x}^{2}}},\quad\ket{K_{15}}=\ket{\mathbf{e_{x}}^{2}\mathbf{e_{z}}-\mathbf{e_{y}}^{2}\mathbf{e_{z}}},$
$\displaystyle\ket{K_{16}}$ $\displaystyle=$ $\displaystyle
21\ket{\mathbf{e_{x}}^{2}\mathbf{e_{y}}^{2}+\mathbf{e_{x}}^{2}\mathbf{e_{z}}^{2}+\mathbf{e_{y}}^{2}\mathbf{e_{z}}^{2}}-16\ket{\mathbf{e_{x}}^{2}+\mathbf{e_{y}}^{2}+\mathbf{e_{z}}^{2}}+12\ket{\mathbf{1}},$
$\displaystyle\ket{K_{17}}$ $\displaystyle=$ $\displaystyle
3\ket{\mathbf{e_{x}}^{2}\mathbf{e_{y}}^{2}+\mathbf{e_{x}}^{2}\mathbf{e_{z}}^{2}-2\mathbf{e_{y}}^{2}\mathbf{e_{z}}^{2}}-2\ket{2\mathbf{e_{x}}^{2}-\mathbf{e_{y}}^{2}-\mathbf{e_{z}}^{2}},$
$\displaystyle\ket{K_{18}}$ $\displaystyle=$ $\displaystyle
3\ket{\mathbf{e_{x}}^{2}\mathbf{e_{y}}^{2}-\mathbf{e_{x}}^{2}\mathbf{e_{z}}^{2}}-2\ket{\mathbf{e_{y}}^{2}-\mathbf{e_{z}}^{2}}.$

These vectors can be used to form an orthogonal matrix $\bm{\mathsf{K}}$ that
maps changes in moments under cascaded collision as

$\bm{\mathsf{K}}=\left[\ket{K_{0}},\ket{K_{1}},\ket{K_{2}},\ldots,\ket{K_{18}}\right]$
(9)

whose elements are given by
$\bm{\mathsf{K}}=$
$\left(\begin{array}[]{rrrrrrrrrrrrrrrrrrrrrrrrrrr}1&0&0&0&0&0&0&0&0&-30&0&0&0&0&0&0&12&0&0\\\
1&1&0&0&0&0&0&1&1&-11&-4&0&0&0&0&0&-4&-4&0\\\
1&-1&0&0&0&0&0&1&1&-11&4&0&0&0&0&0&-4&-4&0\\\
1&0&1&0&0&0&0&-1&1&-11&0&-4&0&0&0&0&-4&2&-2\\\
1&0&-1&0&0&0&0&-1&1&-11&0&4&0&0&0&0&-4&2&-2\\\
1&0&0&1&0&0&0&0&-2&-11&0&0&-4&0&0&0&-4&2&2\\\
1&0&0&-1&0&0&0&0&-2&-11&0&0&4&0&0&0&-4&2&2\\\
1&1&1&0&1&0&0&0&2&8&1&1&0&1&-1&0&1&1&1\\\
1&-1&1&0&-1&0&0&0&2&8&-1&1&0&-1&-1&0&1&1&1\\\
1&1&-1&0&-1&0&0&0&2&8&1&-1&0&1&1&0&1&1&1\\\
1&-1&-1&0&1&0&0&0&2&8&-1&-1&0&-1&1&0&1&1&1\\\
1&1&0&1&0&1&0&1&-1&8&1&0&1&-1&0&1&1&1&-1\\\
1&-1&0&1&0&-1&0&1&-1&8&-1&0&1&1&0&1&1&1&-1\\\
1&1&0&-1&0&-1&0&1&-1&8&1&0&-1&-1&0&-1&1&1&-1\\\
1&-1&0&-1&0&1&0&1&-1&8&-1&0&-1&1&0&-1&1&1&-1\\\
1&0&1&1&0&0&1&-1&-1&8&0&1&1&0&-1&-1&1&-2&0\\\
1&0&-1&1&0&0&-1&-1&-1&8&0&-1&1&0&-1&-1&1&-2&0\\\
1&0&1&-1&0&0&-1&-1&-1&8&0&1&-1&0&1&1&1&-2&0\\\
1&0&-1&-1&0&0&1&-1&-1&8&0&-1&-1&0&-1&1&1&-2&0\\\ \end{array}\right)$

Next, we define the continuous central moment of the particle distribution
function $f$ and its attractor or the local equilibriuam $f^{at}$ shifted by
the local fluid velocity of order $(m+n+p)$ as

$\displaystyle\left[\begin{array}[]{l}{{\hat{\Pi}}_{{x^{m}}{y^{n}}{z^{p}}}}\\\
{{\hat{\Pi}}^{at}_{{x^{m}}{y^{n}}{z^{p}}}}\end{array}\right]=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{\left[\begin{array}[]{l}{f}\\\
f^{at}\\\
\end{array}\right]}(\xi_{x}-u_{x})^{m}(\xi_{y}-u_{y})^{n}(\xi_{z}-u_{z})^{p}d\xi_{x}d\xi_{y}d\xi_{z}.$
(14)

Here, and in the rest of this paper, we employ "hat" over a symbol to
represent the values in the space of moments. One possibility is to consider
the local Maxwell-Boltzmann distribution function in the continuous particle
velocity space ($\bm{\xi}=(\xi_{x}$,$\xi_{y}$,$\xi_{z})$) as the attractor

$f^{\mathcal{M}}\equiv f^{\mathcal{M}}(\rho,\bm{u};\bm{\xi})=\frac{\rho}{2\pi
c_{s}^{3/2}}\exp\left[-\left(\bm{\xi}-\bm{u}\right)^{2}/(2c_{s}^{2})\right],$
(15)

where we typically choose $c_{s}^{2}=1/3$, and its central moments follow via
the definition given in Eq. (14), which can be written as follows:

$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{0}$ $\displaystyle=$
$\displaystyle\rho,$ $\displaystyle\widehat{\Pi}^{\mathcal{M}}_{i}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{ii}$ $\displaystyle=$
$\displaystyle c_{s}^{2}\rho,$ $\displaystyle\widehat{\Pi}^{\mathcal{M}}_{ij}$
$\displaystyle=$ $\displaystyle 0,\quad i\neq j,$
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{ijj}$ $\displaystyle=$
$\displaystyle 0,\quad i\neq j,$
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{ijk}$ $\displaystyle=$
$\displaystyle 0,\quad i\neq j\neq k,$
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{iijj}$ $\displaystyle=$
$\displaystyle c_{s}^{4}\rho,\quad i\neq j.$ (16)

More generally, as suggested in Ref. [30], we consider the factorized form of
attractors that guide the collision process, i.e.,

$\displaystyle\widehat{\Pi}^{at}_{i}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{i}=0,$
$\displaystyle\widehat{\Pi}^{at}_{ij}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{i}\widetilde{\widehat{\Pi}}_{j}=0,\quad
i\neq j$ $\displaystyle\widehat{\Pi}^{at}_{iij}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{ii}\widetilde{\widehat{\Pi}}_{j}=0,$
$\displaystyle\widehat{\Pi}^{at}_{ijk}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{i}\widetilde{\widehat{\Pi}}_{j}\widetilde{\widehat{\Pi}}_{k}=0,$
$\displaystyle\widehat{\Pi}^{at}_{iijj}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{ii}\widetilde{\widehat{\Pi}}_{jj},$
(17)

where ‘tilde’ denote the post-collision values and
$\widetilde{\widehat{\Pi}}_{x^{m}y^{n}z^{p}}$ represents the central moments
of order $(m+n+p)$ of the post-collision distribution function
$\widetilde{f}$:

$\widetilde{\widehat{\Pi}}_{x^{m}y^{n}z^{p}}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\widetilde{f}(\xi_{x}-u_{x})^{m}(\xi_{y}-u_{y})^{n}(\xi_{z}-u_{z})^{p}d\xi_{x}d\xi_{y}d\xi_{z}.$
(18)

Nevertheless, the diagonal components of the lower second-order central
moments of the attractor should be obtained from the Maxwellian central
moments, i.e., $\widehat{\Pi}^{\mathcal{M}}_{ii}=c_{s}^{2}\rho$ in order to
correctly recover the pressure and momentum flux in the macroscopic flow
equations [30]. Then, the $19$ independent components of the central moment
attractors based on the factorized formulation can be written as

$\displaystyle\widehat{\Pi}^{at}_{0}$ $\displaystyle=$ $\displaystyle
0,\widehat{\Pi}^{at}_{x}=\widehat{\Pi}^{at}_{y}=\widehat{\Pi}^{at}_{z}=0,$
$\displaystyle\widehat{\Pi}^{at}_{xx}$ $\displaystyle=$
$\displaystyle\widehat{\Pi}^{at}_{yy}=\widehat{\Pi}^{at}_{zz}=c_{s}^{2}\rho,$
$\displaystyle\widehat{\Pi}^{at}_{xy}$ $\displaystyle=$
$\displaystyle\widehat{\Pi}^{at}_{xz}=\widehat{\Pi}^{at}_{yz}=0,$
$\displaystyle\widehat{\Pi}^{at}_{xyy}$ $\displaystyle=$
$\displaystyle\widehat{\Pi}^{at}_{xzz}=\widehat{\Pi}^{at}_{xxy}=\widehat{\Pi}^{at}_{yzz}=\widehat{\Pi}^{at}_{xxz}=\widehat{\Pi}^{at}_{yyz}=\widehat{\Pi}^{at}_{xyz}=0,$
$\displaystyle\widehat{\Pi}^{at}_{xxyy}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{xx}\widetilde{\widehat{\Pi}}_{yy},$
$\displaystyle\widehat{\Pi}^{at}_{xxzz}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{xx}\widetilde{\widehat{\Pi}}_{zz},$
$\displaystyle\widehat{\Pi}^{at}_{yyzz}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\Pi}}_{yy}\widetilde{\widehat{\Pi}}_{zz}.$
(19)

Based on the above, we will now construct a 3D cascaded LBM by defining a
discrete distribution function supported by the discrete particle velocity set
$\bm{e}_{\alpha}$ as $\mathbf{f}=(f_{0},f_{1},f_{2},\ldots,f_{18})^{\dagger}$
and its collision term as
$\mathbf{\Omega}^{c}=\ket{\Omega_{\alpha}^{c}}=(\Omega_{0}^{c},\Omega_{1}^{c},\Omega_{2}^{c},\ldots,\Omega_{18}^{c})^{\dagger}$,
which will be based on the relaxation of various central moments to their
attractors. The cascaded LBM results from a special discretization of the
continuous Boltzmann equation by integrating along the particle
characteristics as

$f_{\alpha}(\bm{x}+\bm{e}_{\alpha},t+1)=f_{\alpha}(\bm{x},t)+\Omega^{c}_{\alpha}(\bm{x},t)$
(20)

Here, the collision term $\Omega_{\alpha}^{c}$ is a function of the unknown
changes in various moments under collision
$\mathbf{\widehat{g}}=(\widehat{g}_{0},\widehat{g}_{1},\widehat{g}_{2},\ldots,\widehat{g}_{18})^{\dagger}$
resulting from the state $\mathbf{f}$, which are then mapped to the changes in
the distribution functions via the orthogonal mapping matrix $\bm{\mathsf{K}}$
as [8]

$\Omega_{\alpha}^{c}\equiv\Omega_{\alpha}^{c}(\mathbf{f},\mathbf{\widehat{g}})=(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})_{\alpha},$
(21)

For convenience, the above cascaded LBM (Eq. (20)) can reformulated as an
algorithm in the form of the usual collision and streaming steps,
respectively, as

$\displaystyle\widetilde{f}_{\alpha}(\bm{x},t)$ $\displaystyle=$
$\displaystyle f_{\alpha}(\bm{x},t)+\Omega^{c}_{\alpha}{(\bm{x},t)},$ (22)
$\displaystyle f_{\alpha}(\bm{x}+\bm{e}_{\alpha},t+1)$ $\displaystyle=$
$\displaystyle\widetilde{f}_{\alpha}(\bm{x},t).$ (23)

The fluid dynamical variables, i.e., the local density and momentum, and the
pressure field $P$ are determined from the updated distribution functions as

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{18}f_{\alpha}=\braket{\mathbf{f}}{\mathbf{1}},\quad
P=\rho c_{s}^{2}$ (24a) $\displaystyle\rho u_{i}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{18}f_{\alpha}e_{\alpha
i}=\braket{\mathbf{f}}{\mathbf{e}_{i}},i\in{x,y,z}.$ (24b)

In order to derive the expressions for the change of moments under collision
$\widehat{\mathbf{g}}$ for the cascaded LBM, we need to also define the
following _discrete_ central moments of the distribution functions
$\mathbf{f}=(f_{0},f_{1},f_{2},\ldots,f_{18})^{\dagger}$, their attractors
$\mathbf{f^{at}}=(f_{0}^{at},f_{1}^{at},f_{2}^{at},\ldots,f_{18}^{at})^{\dagger}$,
and the post-collision distribution functions
$\mathbf{\widetilde{f}}=(\widetilde{f}_{0},\widetilde{f}_{1},\widetilde{f}_{2},\ldots,\widetilde{f}_{18})^{\dagger}$:

$\displaystyle\widehat{\kappa}_{x^{m}y^{n}z^{p}}$ $\displaystyle=$
$\displaystyle\braket{(\mathbf{e_{x}}-u_{x}\mathbf{1})^{m}(\mathbf{e_{y}}-u_{y}\mathbf{1})^{n}(\mathbf{e_{z}}-u_{z}\mathbf{1})^{p}}{\mathbf{f}},$
$\displaystyle\widehat{\kappa}_{x^{m}y^{n}z^{p}}^{at}$ $\displaystyle=$
$\displaystyle\braket{(\mathbf{e_{x}}-u_{x}\mathbf{1})^{m}(\mathbf{e_{y}}-u_{y}\mathbf{1})^{n}(\mathbf{e_{z}}-u_{z}\mathbf{1})^{p}}{\mathbf{f^{at}}},$
$\displaystyle\widetilde{\widehat{\kappa}}_{x^{m}y^{n}z^{p}}$ $\displaystyle=$
$\displaystyle\braket{(\mathbf{e_{x}}-u_{x}\mathbf{1})^{m}(\mathbf{e_{y}}-u_{y}\mathbf{1})^{n}(\mathbf{e_{z}}-u_{z}\mathbf{1})^{p}}{\mathbf{\widetilde{f}}},.$
(25)

For achieving highest possible accuracy, we match the discrete central moments
of the attractors supported by the D3Q19 lattice with the corresponding
continuous central moments, i.e.,
$\widehat{\kappa}_{x^{m}y^{n}z^{p}}^{at}=\widehat{\Pi}^{at}_{x^{m}y^{n}z^{p}}$.
As a result, the discrete central moments of the attractors can be listed as
follows:

$\displaystyle\widehat{\kappa}^{at}_{0}$ $\displaystyle=$ $\displaystyle
0,\widehat{\kappa}^{at}_{x}=\widehat{\kappa}^{at}_{y}=\widehat{\kappa}^{at}_{z}=0,$
$\displaystyle\widehat{\kappa}^{at}_{xx}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}^{at}_{yy}=\widehat{\kappa}^{at}_{zz}=c_{s}^{2}\rho,$
$\displaystyle\widehat{\kappa}^{at}_{xy}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}^{at}_{xz}=\widehat{\kappa}^{at}_{yz}=0,$
$\displaystyle\widehat{\kappa}^{at}_{xyy}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}^{at}_{xzz}=\widehat{\kappa}^{at}_{xxy}=\widehat{\kappa}^{at}_{yzz}=\widehat{\kappa}^{at}_{xxz}=\widehat{\kappa}^{at}_{yyz}=\widehat{\kappa}^{at}_{xyz}=0,$
$\displaystyle\widehat{\kappa}^{at}_{xxyy}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{yy},$
$\displaystyle\widehat{\kappa}^{at}_{xxzz}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{zz},$
$\displaystyle\widehat{\kappa}^{at}_{yyzz}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{yy}\widetilde{\widehat{\kappa}}_{zz}.$
(26)

We also require the raw moments in order to relate to the results obtained in
terms of the central moments in the construction of the cascaded collision
term in what follows, which we define as

$\widehat{\kappa}_{x^{m}y^{n}z^{p}}^{{}^{\prime}}=\braket{\mathbf{e_{x}}^{m}\mathbf{e_{y}}^{n}\mathbf{e_{z}}^{p}}{\mathbf{f}},$
(27)

where the superscript ‘prime’ (′) is used for the various raw moments in order
to distinguish the raw moments from the central moments that are designated
without the primes.

###  3.1 Change of moments under cascaded collision for D3Q19 lattice

We now need the various moments of the cascaded collision term, i.e.,

$\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}^{m}\mathbf{e_{y}}^{n}\mathbf{e_{z}}^{p}}$
(28)

to proceed further. Since the zeroth and first order moments are collision
invariants, it follows that the corresponding moment changes under collision
are zero, i.e.,
$\widehat{g}_{0}=\widehat{g}_{1}=\widehat{g}_{2}=\widehat{g}_{3}=0$. In view
of the fact that $\bm{\mathsf{K}}$ is orthogonal, the expressions for the
various non-conserved components can be readily evaluated, which read as

$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}\mathbf{e_{y}}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{4},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}\mathbf{e_{z}}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{5},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{y}}\mathbf{e_{z}}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{6},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}^{2}}$
$\displaystyle=$ $\displaystyle
6\widehat{g}_{7}+6\widehat{g}_{8}+42\widehat{g}_{9},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{y}}^{2}}$
$\displaystyle=$
$\displaystyle-6\widehat{g}_{7}+6\widehat{g}_{8}+42\widehat{g}_{9},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{z}}^{2}}$
$\displaystyle=$ $\displaystyle-12\widehat{g}_{8}+42\widehat{g}_{9},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}\mathbf{e_{y}}^{2}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{10}+4\widehat{g}_{13},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}\mathbf{e_{z}}^{2}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{10}-4\widehat{g}_{13},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}^{2}\mathbf{e_{y}}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{11}-4\widehat{g}_{14},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{y}}\mathbf{e_{z}}^{2}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{11}+4\widehat{g}_{14},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}^{2}\mathbf{e_{z}}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{12}+4\widehat{g}_{15},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{y}}^{2}\mathbf{e_{z}}}$
$\displaystyle=$ $\displaystyle 4\widehat{g}_{12}-4\widehat{g}_{15},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}^{2}\mathbf{e_{y}}^{2}}$
$\displaystyle=$ $\displaystyle
8\widehat{g}_{8}+32\widehat{g}_{9}+4\widehat{g}_{16}+4\widehat{g}_{17}+4\widehat{g}_{18},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{x}}^{2}\mathbf{e_{z}}^{2}}$
$\displaystyle=$ $\displaystyle
4\widehat{g}_{7}-4\widehat{g}_{8}+32\widehat{g}_{9}+4\widehat{g}_{16}+4\widehat{g}_{17}-4\widehat{g}_{18},$
$\displaystyle\braket{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}{\mathbf{e_{y}}^{2}\mathbf{e_{z}}^{2}}$
$\displaystyle=$
$\displaystyle-4\widehat{g}_{7}-4\widehat{g}_{8}+32\widehat{g}_{9}+4\widehat{g}_{16}-8\widehat{g}_{17}.$

Then, the expressions for the change of different moments under collision
$\mathbf{\widehat{g}}$ of the 3D cascaded LBM appearing in the term
$\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}}$ can be derived as follows. In
essence, the procedure starts from the lowest order, non-conserved, central
moments (i.e., $\widehat{\kappa}_{xy}$, $\widehat{\kappa}_{xz}$,
$\widehat{\kappa}_{yz}$ and higher), and their post-collision central moments
(i.e., $\widetilde{\widehat{\kappa}}_{xy},\widetilde{\widehat{\kappa}}_{xz}$,
$\widetilde{\widehat{\kappa}}_{yz}$ and higher) are successively set equal to
the corresponding attractors given in Eq. (26) (i.e.,
$\widehat{\kappa}_{xy}^{at},\widehat{\kappa}_{xz}^{at}$ and
$\widehat{\kappa}_{yz}^{at}$, and higher). This intermediate step can provide
tentative expressions for $\mathbf{\widehat{g}}$ based on an equilibrium
assumption, which are then modified to allow for collision as a relaxation
process. They are multiplied with a corresponding relaxation parameter that
results in the final expressions for the change of moments under collision
$\widehat{g}_{\alpha}$ for a given order [8]. In this step, the relaxation
parameter needs to multiply only with those terms that are not yet in post-
collision states. In other words, the relaxation of different central moments
to their corresponding attractors may be formally represented as

$\braket{(\mathbf{e_{x}}-u_{x}\mathbf{1})^{m}(\mathbf{e_{y}}-u_{y}\mathbf{1})^{n}(\mathbf{e_{z}}-u_{z}\mathbf{1})^{p}}{(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})}=\omega_{*}\left(\widehat{\kappa}_{x^{m}y^{n}z^{p}}^{at}-\widehat{\kappa}_{x^{m}y^{n}z^{p}}\right),$
(29)

where $\omega_{*}$ represents the relaxation parameter for the central moment
of order $(m+n+p)$. In simplifying the above equation (Eq. (29)) for different
combinations of $m$, $n$ and $p$ representing the independent moments
supported by the D3Q19 lattice, we use Eq. (3.1), which then lead to a
cascaded structure, i.e., the change of moments under cascaded collision
depend successively on those of the preceding lower order moment changes. The
final results, after transforming the central moments to raw moments (see Eq.
(27) for definition) for convenience, can then be summarized as follows:

$\widehat{g}_{4}=\frac{\omega_{4}}{4}\left[-\widehat{\kappa}_{xy}^{{}^{\prime}}+\rho
u_{x}u_{y}\right],$ (30)
$\widehat{g}_{5}=\frac{\omega_{5}}{4}\left[-\widehat{\kappa}_{xz}^{{}^{\prime}}+\rho
u_{x}u_{z}\right],$ (31)
$\widehat{g}_{6}=\frac{\omega_{6}}{4}\left[-\widehat{\kappa}_{yz}^{{}^{\prime}}+\rho
u_{y}u_{z}\right].$ (32)
$\widehat{g}_{7}=\frac{\omega_{7}}{12}\left[-(\widehat{\kappa}_{xx}^{{}^{\prime}}-\widehat{\kappa}_{yy}^{{}^{\prime}})+\rho(u_{x}^{2}-u_{y}^{2})\right],$
(33) $\displaystyle\widehat{g}_{8}$ $\displaystyle=$
$\displaystyle\frac{\omega_{8}}{36}\left[-(\widehat{\kappa}_{xx}^{{}^{\prime}}+\widehat{\kappa}_{yy}^{{}^{\prime}}-2\widehat{\kappa}_{zz}^{{}^{\prime}})+\rho(u_{x}^{2}+u_{y}^{2}-2u_{z}^{2})\right],$
(34) $\displaystyle\widehat{g}_{9}$ $\displaystyle=$
$\displaystyle\frac{\omega_{9}}{126}\left[-(\widehat{\kappa}_{xx}^{{}^{\prime}}+\widehat{\kappa}_{yy}^{{}^{\prime}}+\widehat{\kappa}_{zz}^{{}^{\prime}})+\rho(u_{x}^{2}+u_{y}^{2}+u_{z}^{2})+\rho\right].$
(35) $\displaystyle\widehat{g}_{10}$ $\displaystyle=$
$\displaystyle\frac{\omega_{10}}{8}\left[-(\widehat{\kappa}_{xyy}^{{}^{\prime}}+\widehat{\kappa}_{xzz}^{{}^{\prime}})+2(u_{y}\widehat{\kappa}_{xy}^{{}^{\prime}}+u_{z}\widehat{\kappa}_{xz}^{{}^{\prime}})+u_{x}(\widehat{\kappa}_{yy}^{{}^{\prime}}+\widehat{\kappa}_{zz}^{{}^{\prime}})\right.$
(36) $\displaystyle\left.-2\rho
u_{x}(u_{y}^{2}+u_{z}^{2})\right]+(u_{y}\widehat{g}_{4}+u_{z}\widehat{g}_{5})+\frac{3}{4}u_{x}(-(\widehat{g}_{7}+\widehat{g}_{8})+14\widehat{g}_{9}),$
$\displaystyle\widehat{g}_{11}$ $\displaystyle=$
$\displaystyle\frac{\omega_{11}}{8}\left[-(\widehat{\kappa}_{yzz}^{{}^{\prime}}+\widehat{\kappa}_{xxy}^{{}^{\prime}})+2(u_{x}\widehat{\kappa}_{xy}^{{}^{\prime}}+u_{z}\widehat{\kappa}_{yz}^{{}^{\prime}})+u_{y}(\widehat{\kappa}_{xx}^{{}^{\prime}}+\widehat{\kappa}_{zz}^{{}^{\prime}})\right.$
(37) $\displaystyle\left.-2\rho
u_{y}(u_{x}^{2}+u_{z}^{2})\right]+(u_{x}\widehat{g}_{4}+u_{z}\widehat{g}_{6})+\frac{3}{4}u_{y}(\widehat{g}_{7}-\widehat{g}_{8}+14\widehat{g}_{9}),$
$\displaystyle\widehat{g}_{12}$ $\displaystyle=$
$\displaystyle\frac{\omega_{12}}{8}\left[-(\widehat{\kappa}_{xxz}^{{}^{\prime}}+\widehat{\kappa}_{yyz}^{{}^{\prime}})+2(u_{x}\widehat{\kappa}_{xz}^{{}^{\prime}}+u_{y}\widehat{\kappa}_{yz}^{{}^{\prime}})+u_{z}(\widehat{\kappa}_{xx}^{{}^{\prime}}+\widehat{\kappa}_{yy}^{{}^{\prime}})\right.$
(38) $\displaystyle\left.-2\rho
u_{z}(u_{x}^{2}+u_{y}^{2})\right]+(u_{x}\widehat{g}_{5}+u_{y}\widehat{g}_{6})+\frac{3}{2}u_{z}(\widehat{g}_{8}+7\widehat{g}_{9}),$
$\displaystyle\widehat{g}_{13}$ $\displaystyle=$
$\displaystyle\frac{\omega_{13}}{8}\left[-(\widehat{\kappa}_{xyy}^{{}^{\prime}}-\widehat{\kappa}_{xzz}^{{}^{\prime}})+2(u_{y}\widehat{\kappa}_{xy}^{{}^{\prime}}-u_{z}\widehat{\kappa}_{xz}^{{}^{\prime}})+u_{x}(\widehat{\kappa}_{yy}^{{}^{\prime}}-\widehat{\kappa}_{zz}^{{}^{\prime}})\right.$
(39) $\displaystyle\left.-2\rho
u_{x}(u_{y}^{2}-u_{z}^{2})\right]+(u_{y}\widehat{g}_{4}-u_{z}\widehat{g}_{5})+\frac{3}{4}u_{x}(-\widehat{g}_{7}+3\widehat{g}_{8}),$
$\displaystyle\widehat{g}_{14}$ $\displaystyle=$
$\displaystyle\frac{\omega_{14}}{8}\left[-(\widehat{\kappa}_{yzz}^{{}^{\prime}}-\widehat{\kappa}_{xxy}^{{}^{\prime}})+2(u_{z}\widehat{\kappa}_{yz}^{{}^{\prime}}-u_{x}\widehat{\kappa}_{xy}^{{}^{\prime}})+u_{y}(\widehat{\kappa}_{zz}^{{}^{\prime}}-\widehat{\kappa}_{xx}^{{}^{\prime}})\right.$
(40) $\displaystyle\left.-2\rho
u_{y}(u_{z}^{2}-u_{x}^{2})\right]+(-u_{x}\widehat{g}_{4}+u_{z}\widehat{g}_{6})+\frac{3}{4}u_{y}(-\widehat{g}_{7}-3\widehat{g}_{8}),$
$\displaystyle\widehat{g}_{15}$ $\displaystyle=$
$\displaystyle\frac{\omega_{15}}{8}\left[-(\widehat{\kappa}_{xxz}^{{}^{\prime}}-\widehat{\kappa}_{yyz}^{{}^{\prime}})+2(u_{x}\widehat{\kappa}_{xz}^{{}^{\prime}}-u_{y}\widehat{\kappa}_{yz}^{{}^{\prime}})+u_{z}(\widehat{\kappa}_{xx}^{{}^{\prime}}-\widehat{\kappa}_{yy}^{{}^{\prime}})\right.$
(41) $\displaystyle\left.-2\rho
u_{z}(u_{x}^{2}-u_{y}^{2})\right.+(u_{x}\widehat{g}_{5}-u_{y}\widehat{g}_{6})+\frac{3}{2}u_{z}\widehat{g}_{7},$
$\displaystyle\widehat{g}_{16}$ $\displaystyle=$
$\displaystyle\frac{\omega_{16}}{12}\left[-(\widehat{\kappa}_{xxyy}^{{}^{\prime}}+\widehat{\kappa}_{xxzz}^{{}^{\prime}}+\widehat{\kappa}_{yyzz}^{{}^{\prime}})+2u_{x}(\widehat{\kappa}_{xyy}^{{}^{\prime}}+\widehat{\kappa}_{xzz}^{{}^{\prime}})+2u_{y}(\widehat{\kappa}_{xxy}^{{}^{\prime}}\right.$
(42)
$\displaystyle\left.+\widehat{\kappa}_{yzz}^{{}^{\prime}})+2u_{z}(\widehat{\kappa}_{xxz}^{{}^{\prime}}+\widehat{\kappa}_{yyz}^{{}^{\prime}})-u_{x}^{2}(\widehat{\kappa}_{yy}^{{}^{\prime}}+\widehat{\kappa}_{zz}^{{}^{\prime}}))-u_{y}^{2}(\widehat{\kappa}_{zz}^{{}^{\prime}}+\widehat{\kappa}_{xx}^{{}^{\prime}})\right.$
$\displaystyle\left.-u_{z}^{2}(\widehat{\kappa}_{xx}^{{}^{\prime}}+\widehat{\kappa}_{yy}^{{}^{\prime}})-4(u_{x}u_{y}\widehat{\kappa}_{xy}^{{}^{\prime}}+u_{x}u_{z}\widehat{\kappa}_{xz}^{{}^{\prime}}+u_{y}u_{z}\widehat{\kappa}_{yz}^{{}^{\prime}})\right.$
$\displaystyle\left.+3\rho(u_{x}^{2}u_{y}^{2}+u_{x}^{2}u_{z}^{2}+u_{y}^{2}u_{z}^{2})+(\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{yy}+\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{zz}+\widetilde{\widehat{\kappa}}_{yy}\widetilde{\widehat{\kappa}}_{zz})\right]$
$\displaystyle-\frac{4}{3}u_{x}u_{y}\widehat{g}_{4}-\frac{4}{3}u_{x}u_{z}\widehat{g}_{5}-\frac{4}{3}u_{y}u_{z}\widehat{g}_{6}+\frac{1}{2}(u_{x}^{2}-u_{y}^{2})\widehat{g}_{7}+\frac{1}{2}(u_{x}^{2}+u_{y}^{2}-2u_{z}^{2})\widehat{g}_{8}$
$\displaystyle+(-7(u_{x}^{2}+u_{y}^{2}+u_{z}^{2})-8)\widehat{g}_{9}+\frac{4}{3}u_{x}\widehat{g}_{10}+\frac{4}{3}u_{y}\widehat{g}_{11}+\frac{4}{3}u_{z}\widehat{g}_{12},$
$\displaystyle\widehat{g}_{17}$ $\displaystyle=$
$\displaystyle\frac{\omega_{17}}{24}\left[-(\widehat{\kappa}_{xxyy}^{{}^{\prime}}+\widehat{\kappa}_{xxzz}^{{}^{\prime}}-2\widehat{\kappa}_{yyzz}^{{}^{\prime}})+2\left(u_{x}\widehat{\kappa}_{xyy}^{{}^{\prime}}+u_{x}\widehat{\kappa}_{xzz}^{{}^{\prime}}+u_{y}\widehat{\kappa}_{xxy}^{{}^{\prime}}+u_{z}\widehat{\kappa}_{xxz}^{{}^{\prime}})-\right.\right.$
(43)
$\displaystyle\left.\left.2(u_{y}\widehat{\kappa}_{yzz}^{{}^{\prime}}+u_{z}\widehat{\kappa}_{yyz}^{{}^{\prime}})\right)-u_{x}^{2}(\widehat{\kappa}_{yy}^{{}^{\prime}}+\widehat{\kappa}_{zz}^{{}^{\prime}})-u_{y}^{2}(\widehat{\kappa}_{xx}^{{}^{\prime}}-2\widehat{\kappa}_{zz}^{{}^{\prime}})-u_{z}^{2}(\widehat{\kappa}_{xx}^{{}^{\prime}}-2\widehat{\kappa}_{yy}^{{}^{\prime}})\right.$
$\displaystyle\left.-4(u_{x}u_{y}\widehat{\kappa}_{xy}^{{}^{\prime}}+u_{x}u_{z}\widehat{\kappa}_{xz}^{{}^{\prime}}-2u_{y}u_{z}\widehat{\kappa}_{yz}^{{}^{\prime}})+(\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{yy}+\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{zz}-2\widetilde{\widehat{\kappa}}_{yy}\widetilde{\widehat{\kappa}}_{zz})\right.$
$\displaystyle\left.+3\rho(u_{x}^{2}u_{y}^{2}+u_{x}^{2}u_{z}^{2}-2u_{y}^{2}u_{z}^{2})\right]-\frac{2}{3}u_{x}u_{y}\widehat{g}_{4}-\frac{2}{3}u_{x}u_{z}\widehat{g}_{5}+\frac{4}{3}u_{y}u_{z}\widehat{g}_{6}+\frac{1}{4}(u_{x}^{2}$
$\displaystyle-
u_{y}^{2}-3u_{z}^{2}-2)\widehat{g}_{7}+\frac{1}{4}(u_{x}^{2}-5u_{y}^{2}+u_{z}^{2}-2)\widehat{g}_{8}+\frac{7}{4}(-2u_{x}^{2}+u_{y}^{2}+u_{z}^{2})\widehat{g}_{9}$
$\displaystyle+\frac{2}{3}u_{x}\widehat{g}_{10}-\frac{1}{3}u_{y}\widehat{g}_{11}-\frac{1}{3}u_{z}\widehat{g}_{12}-u_{y}\widehat{g}_{14}+u_{z}\widehat{g}_{15},$
$\displaystyle\widehat{g}_{18}$ $\displaystyle=$
$\displaystyle\frac{\omega_{18}}{8}\left[-(\widehat{\kappa}_{xxyy}^{{}^{\prime}}-\widehat{\kappa}_{xxzz}^{{}^{\prime}})+2\left(u_{x}\widehat{\kappa}_{xyy}^{{}^{\prime}}-u_{x}\widehat{\kappa}_{xzz}^{{}^{\prime}}+u_{y}\widehat{\kappa}_{xxy}^{{}^{\prime}}-u_{z}\widehat{\kappa}_{xxz}^{{}^{\prime}})\right.\right.$
(44)
$\displaystyle\left.-u_{x}^{2}(\widehat{\kappa}_{yy}^{{}^{\prime}}-\widehat{\kappa}_{zz}^{{}^{\prime}})-u_{y}^{2}\widehat{\kappa}_{xx}^{{}^{\prime}}+u_{z}^{2}\widehat{\kappa}_{xx}^{{}^{\prime}}-4(u_{x}u_{y}\widehat{\kappa}_{xy}^{{}^{\prime}}-u_{x}u_{z}\widehat{\kappa}_{xz}^{{}^{\prime}})\right.$
$\displaystyle+(\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{yy}-\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{zz})\left.+3\rho(u_{x}^{2}u_{y}^{2}-u_{x}^{2}u_{z}^{2})\right]-2u_{x}u_{y}\widehat{g}_{4}+2u_{x}u_{z}\widehat{g}_{5}$
$\displaystyle+\frac{3}{4}(u_{x}^{2}-u_{y}^{2}+u_{z}^{2}+\frac{2}{3})\widehat{g}_{7}+\frac{3}{4}(-3u_{x}^{2}-u_{y}^{2}+u_{z}^{2}-2)\widehat{g}_{8}+\frac{21}{4}(-u_{y}^{2}+u_{z}^{2})\widehat{g}_{9}$
$\displaystyle+u_{y}\widehat{g}_{11}-u_{z}\widehat{g}_{12}+2u_{x}\widehat{g}_{13}-u_{y}\widehat{g}_{14}-u_{z}\widehat{g}_{15}.$

In the above, $\omega_{4}$, $\omega_{5}\dots$ and $\omega_{18}$ are relaxation
parameters of the central moments of different orders. Similar to the 2D
central moment LBM [10], we can apply the Chapman-Enskog expansion [31] to the
above 3D cascaded LBM to show that it represents the Navier-Stokes equations.
Some of the relaxation parameters in the collision model can be related to the
transport coefficients of the momentum transfer. For example, those
corresponding to the second-order moments control the shear viscosity $\nu$
and those related to the trace (isotropic part) of such moments determine the
bulk viscosity $\zeta$ of the fluid. That is,
$\nu=c_{s}^{2}\left(\frac{1}{\omega^{\nu}}-\frac{1}{2}\right)$ where
$\omega^{\nu}=\omega_{j}$ where $j=4,5,6,7,8$ and
$\zeta=\frac{2}{3}c_{s}^{2}\left(\frac{1}{\omega^{\zeta}}-\frac{1}{2}\right)$
where $\omega^{\zeta}=\omega_{9}$. The rest of the parameters can be set
either to $1$ (i.e., equilibration) or adjusted independently to carefully
control and improve numerical stability by means of a linear stability
analysis, while all satisfying the usual bounds $0<\omega_{\beta}<2$. In this
work, for simplicity, we set the relaxation parameters for higher order
moments to unity, i.e., $\omega_{j}=1$, where $j=9,10,11,...,18$.

###  3.2 Nonlinear constitutive relations via local variations of relaxation
times and local computation of strain rate tensor

Then, for the simulation of non-Newtonian flows of power law fluids using the
above 3D cascaded LBM, where the nonlinear constitutive relation prescribe the
apparent viscosity dependent on the shear rate $|\dot{\gamma}_{ij}|$ , i.e.,
$\mu(|\dot{\gamma}_{ij}|)=\rho\nu(|\dot{\gamma}_{ij}|)$ (see Sec. 2), the
relaxation times of the second order moments
$\omega_{4},\omega_{5},\omega_{6},\omega_{7}$ and $\omega_{8}$ are locally
adjusted as follows:

$\omega^{\nu}(\bm{x},t)=\omega_{4}=\omega_{5}=\omega_{6}=\omega_{7}=\omega_{8}=\left[\frac{\mu(|\dot{\gamma}_{ij}|)}{\rho
c_{s}^{2}}+\frac{1}{2}\right]^{-1},$ (45)

where $\mu(|\dot{\gamma}_{ij}|)$ can be obtained from Eq. (7), which is
parametrized by the consistency coefficient $\mu_{p}$ and the power law index
$n$. The magnitude of the shear rate $|\dot{\gamma}_{ij}|$ is related to the
second invariant of the strain rate tensor $S_{ij}$ via Eq. (5). The strain
rate tensor components $S_{xx},S_{yy},S_{zz},S_{xy},S_{yz}$ and $S_{xz}$ can
be obtained locally using non-equilibrium moments as follows.

A Chapman-Enskog analysis [31] of the 3D cascaded LBM, whose details are
omitted here for brevity but follows as an extension of that presented in [10,
15], relates the components of the second order non-equilibrium moments to the
components of the symmetric part of the velocity gradient tensor. The results
can be summarized as follows:

$\displaystyle\widehat{m}^{(1)}_{4}$ $\displaystyle\approx$
$\displaystyle-\frac{1}{3\omega_{4}}\rho(\partial_{y}u_{x}+\partial_{x}u_{y}),\quad\widehat{m}^{(1)}_{5}\approx-\frac{1}{3\omega_{5}}\rho(\partial_{z}u_{x}+\partial_{x}u_{z}),$
$\displaystyle\widehat{m}^{(1)}_{6}$ $\displaystyle\approx$
$\displaystyle-\frac{1}{3\omega_{6}}\rho(\partial_{z}u_{y}+\partial_{y}u_{z}),\quad\widehat{m}^{(1)}_{7}\approx-\frac{2}{3\omega_{7}}\rho(\partial_{x}u_{x}-\partial_{y}u_{y}),$
$\displaystyle\widehat{m}^{(1)}_{8}$ $\displaystyle\approx$
$\displaystyle-\frac{2}{3\omega_{8}}\rho(\partial_{x}u_{x}-\partial_{z}u_{z}),\quad\widehat{m}^{(1)}_{9}\approx-\frac{2}{3\omega_{9}}\rho(\partial_{x}u_{x}+\partial_{y}u_{y}++\partial_{z}u_{z}),$

where the non-equilibrium moments
$\widehat{m}^{(1)}_{4},\widehat{m}^{(1)}_{5},\ldots,\widehat{m}^{(1)}_{9}$ are
given by

$\displaystyle\widehat{m}^{(1)}_{4}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}^{{}^{\prime}}_{xy}-\rho u_{x}u_{y},$ (46a)
$\displaystyle\widehat{m}^{(1)}_{5}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}^{{}^{\prime}}_{xz}-\rho u_{x}u_{z},$ (46b)
$\displaystyle\widehat{m}^{(1)}_{6}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}^{{}^{\prime}}_{yz}-\rho u_{y}u_{z},$ (46c)
$\displaystyle\widehat{m}^{(1)}_{7}$ $\displaystyle=$
$\displaystyle(\widehat{\kappa}^{{}^{\prime}}_{xx}-\widehat{\kappa}^{{}^{\prime}}_{yy})-\rho(u_{x}^{2}-u_{y}^{2}),$
(46d) $\displaystyle\widehat{m}^{(1)}_{8}$ $\displaystyle=$
$\displaystyle(\widehat{\kappa}^{{}^{\prime}}_{xx}-\widehat{\kappa}^{{}^{\prime}}_{zz})-\rho(u_{x}^{2}-u_{z}^{2}),$
(46e) $\displaystyle\widehat{m}^{(1)}_{9}$ $\displaystyle=$
$\displaystyle(\widehat{\kappa}^{{}^{\prime}}_{xx}+\widehat{\kappa}^{{}^{\prime}}_{yy}+\widehat{\kappa}^{{}^{\prime}}_{zz})-(\rho(u_{x}^{2}+u_{y}^{2}+u_{z}^{2})+\rho).$
(46f)

From the above, and considering Eq. (45) and invoking the definition of the
strain rate tensor $S_{ij}=\frac{1}{2}(\partial_{j}u_{i}+\partial_{i}u_{j})$,
its components can be solved in terms of the non-equilibrium moments
$\widehat{m}^{(1)}_{4},\widehat{m}^{(1)}_{5},\ldots,\widehat{m}^{(1)}_{9}$ as

$\displaystyle S_{xy}$ $\displaystyle=$
$\displaystyle-\frac{3}{2\rho}\omega^{\nu}(\bm{x},t)\widehat{m}^{(1)}_{4},$
(47a) $\displaystyle S_{xz}$ $\displaystyle=$
$\displaystyle-\frac{3}{2\rho}\omega^{\nu}(\bm{x},t)\widehat{m}^{(1)}_{5},$
(47b) $\displaystyle S_{yz}$ $\displaystyle=$
$\displaystyle-\frac{3}{2\rho}\omega^{\nu}(\bm{x},t)\widehat{m}^{(1)}_{6},$
(47c) $\displaystyle S_{xx}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\rho}\left\\{\omega^{\nu}(\bm{x},t)\left[\widehat{m}^{(1)}_{7}+\widehat{m}^{(1)}_{8}\right]+\omega^{\zeta}\widehat{m}^{(1)}_{9}\right\\},$
(47d) $\displaystyle S_{yy}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\rho}\left\\{\omega^{\nu}(\bm{x},t)\left[-2\widehat{m}^{(1)}_{7}+\widehat{m}^{(1)}_{8}\right]+\omega^{\zeta}\widehat{m}^{(1)}_{9}\right\\},$
(47e) $\displaystyle S_{zz}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\rho}\left\\{\omega^{\nu}(\bm{x},t)\left[\widehat{m}^{(1)}_{7}-2\widehat{m}^{(1)}_{8}\right]+\omega^{\zeta}\widehat{m}^{(1)}_{9}\right\\}.$
(47f)

The above two sets of equations enable local computation of the strain rate
tensor in simulating non-Newtonian flows and they are naturally suitable for
parallel computing, as opposed to using finite difference approximations.

###  3.3 Post-collision distribution functions

Finally, evaluating the elements of
$(\bm{\mathsf{K}}\cdot\mathbf{\widehat{g}})_{\alpha}$ in the collision step of
the 3D cascaded LBM (see Eqs. (21) and (22)), the post collision distribution
functions can be obtained. These are summarized as follows:

$\displaystyle\widetilde{f}_{0}$ $\displaystyle=$ $\displaystyle
f_{0}+\left[\widehat{g}_{0}-30\widehat{g}_{9}+12\widehat{g}_{16}\right],$
$\displaystyle\widetilde{f}_{1}$ $\displaystyle=$ $\displaystyle
f_{1}+\left[\widehat{g}_{0}+\widehat{g}_{1}+\widehat{g}_{7}+\widehat{g}_{8}-11\widehat{g}_{9}-4\widehat{g}_{10}-4\widehat{g}_{16}-4\widehat{g}_{17}\right],$
$\displaystyle\widetilde{f}_{2}$ $\displaystyle=$ $\displaystyle
f_{2}+\left[\widehat{g}_{0}-\widehat{g}_{1}+\widehat{g}_{7}+\widehat{g}_{8}-11\widehat{g}_{9}+4\widehat{g}_{10}-4\widehat{g}_{16}-4\widehat{g}_{17}\right],$
$\displaystyle\widetilde{f}_{3}$ $\displaystyle=$ $\displaystyle
f_{3}+\left[\widehat{g}_{0}+\widehat{g}_{2}-\widehat{g}_{7}+\widehat{g}_{8}-11\widehat{g}_{9}-4\widehat{g}_{11}-4\widehat{g}_{16}+2\widehat{g}_{17}-2\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{4}$ $\displaystyle=$ $\displaystyle
f_{4}+\left[\widehat{g}_{0}-\widehat{g}_{2}-\widehat{g}_{7}+\widehat{g}_{8}-11\widehat{g}_{9}+4\widehat{g}_{11}-4\widehat{g}_{16}+2\widehat{g}_{17}-2\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{5}$ $\displaystyle=$ $\displaystyle
f_{5}+\left[\widehat{g}_{0}+\widehat{g}_{3}-2\widehat{g}_{8}-11\widehat{g}_{9}-4\widehat{g}_{12}-4\widehat{g}_{16}+2\widehat{g}_{17}+2\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{6}$ $\displaystyle=$ $\displaystyle
f_{6}+\left[\widehat{g}_{0}-\widehat{g}_{3}-2\widehat{g}_{8}-11\widehat{g}_{9}+4\widehat{g}_{12}-4\widehat{g}_{16}+2\widehat{g}_{17}+2\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{7}$ $\displaystyle=$ $\displaystyle
f_{7}+\left[\widehat{g}_{0}+\widehat{g}_{1}+\widehat{g}_{2}+\widehat{g}_{4}+2\widehat{g}_{8}+8\widehat{g}_{9}+\widehat{g}_{10}+\widehat{g}_{11}+\widehat{g}_{13}-\widehat{g}_{14}+\widehat{g}_{16}+\widehat{g}_{17}+\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{8}$ $\displaystyle=$ $\displaystyle
f_{8}+\left[\widehat{g}_{0}-\widehat{g}_{1}+\widehat{g}_{2}-\widehat{g}_{4}+2\widehat{g}_{8}+8\widehat{g}_{9}-\widehat{g}_{10}+\widehat{g}_{11}-\widehat{g}_{13}-\widehat{g}_{14}+\widehat{g}_{16}+\widehat{g}_{17}+\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{9}$ $\displaystyle=$ $\displaystyle
f_{9}+\left[\widehat{g}_{0}+\widehat{g}_{1}-\widehat{g}_{2}-\widehat{g}_{4}+2\widehat{g}_{8}+8\widehat{g}_{9}+\widehat{g}_{10}-\widehat{g}_{11}+\widehat{g}_{13}+\widehat{g}_{14}+\widehat{g}_{16}+\widehat{g}_{17}+\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{10}$ $\displaystyle=$ $\displaystyle
f_{10}+\left[\widehat{g}_{0}-\widehat{g}_{1}-\widehat{g}_{2}+\widehat{g}_{4}+2\widehat{g}_{8}+8\widehat{g}_{9}-\widehat{g}_{10}-\widehat{g}_{11}-\widehat{g}_{13}+\widehat{g}_{14}+\widehat{g}_{16}+\widehat{g}_{17}+\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{11}$ $\displaystyle=$ $\displaystyle
f_{11}+\left[\widehat{g}_{0}+\widehat{g}_{1}+\widehat{g}_{3}+\widehat{g}_{5}+\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}+\widehat{g}_{10}+\widehat{g}_{12}-\widehat{g}_{13}+\widehat{g}_{15}+\widehat{g}_{16}+\widehat{g}_{17}-\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{12}$ $\displaystyle=$ $\displaystyle
f_{12}+\left[\widehat{g}_{0}-\widehat{g}_{1}+\widehat{g}_{3}-\widehat{g}_{5}+\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}-\widehat{g}_{10}+\widehat{g}_{12}+\widehat{g}_{13}+\widehat{g}_{15}+\widehat{g}_{16}+\widehat{g}_{17}-\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{13}$ $\displaystyle=$ $\displaystyle
f_{13}+\left[\widehat{g}_{0}+\widehat{g}_{1}-\widehat{g}_{3}-\widehat{g}_{5}+\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}+\widehat{g}_{10}-\widehat{g}_{12}-\widehat{g}_{13}-\widehat{g}_{15}+\widehat{g}_{16}+\widehat{g}_{17}-\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{14}$ $\displaystyle=$ $\displaystyle
f_{14}+\left[\widehat{g}_{0}-\widehat{g}_{1}-\widehat{g}_{3}+\widehat{g}_{5}+\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}-\widehat{g}_{10}-\widehat{g}_{12}+\widehat{g}_{13}-\widehat{g}_{15}+\widehat{g}_{16}+\widehat{g}_{17}-\widehat{g}_{18}\right],$
$\displaystyle\widetilde{f}_{15}$ $\displaystyle=$ $\displaystyle
f_{15}+\left[\widehat{g}_{0}+\widehat{g}_{2}+\widehat{g}_{3}-\widehat{g}_{6}-\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}+\widehat{g}_{11}+\widehat{g}_{12}+\widehat{g}_{14}-\widehat{g}_{15}+\widehat{g}_{16}-2\widehat{g}_{17}\right],$
$\displaystyle\widetilde{f}_{16}$ $\displaystyle=$ $\displaystyle
f_{16}+\left[\widehat{g}_{0}-\widehat{g}_{2}+\widehat{g}_{3}-\widehat{g}_{6}-\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}-\widehat{g}_{11}+\widehat{g}_{12}-\widehat{g}_{14}-\widehat{g}_{15}+\widehat{g}_{16}-2\widehat{g}_{17}\right],$
$\displaystyle\widetilde{f}_{17}$ $\displaystyle=$ $\displaystyle
f_{17}+\left[\widehat{g}_{0}+\widehat{g}_{2}-\widehat{g}_{3}-\widehat{g}_{6}-\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}+\widehat{g}_{11}-\widehat{g}_{12}+\widehat{g}_{14}+\widehat{g}_{15}+\widehat{g}_{16}-2\widehat{g}_{17}\right],$
$\displaystyle\widetilde{f}_{18}$ $\displaystyle=$ $\displaystyle
f_{18}+\left[\widehat{g}_{0}-\widehat{g}_{2}-\widehat{g}_{3}+\widehat{g}_{6}-\widehat{g}_{7}-\widehat{g}_{8}+8\widehat{g}_{9}-\widehat{g}_{11}-\widehat{g}_{12}-\widehat{g}_{14}+\widehat{g}_{15}+\widehat{g}_{16}-2\widehat{g}_{17}\right].$

This is then followed by the streaming step via Eq. (23), after which the
hydrodynamic fields of the 3D non-Newtonian fluid flow can be obtained from
Eqs. (24a) and (24b).

##  4 Results and discussion

To validate our 3D cascaded LBM for non-Newtonian flows, we consider the
following benchmark problems for which either the analytical solutions and/or
numerical benchmark results available: (i) non-Newtonian Poiseuille flow of a
power law fluid, (ii) three-dimensional fully developed flow in a square duct,
and (iii) power-law fluid flow in a cubic lid-driven cavity (see Fig. 2). For
the first two problems, there exist analytical solutions, while for the third
problem, we use the recent benchmark numerical data of [29] for comparison.

###  4.1 Poiseuille flow of power-law fluid

First, we consider the flow of power-law fluids between two parallel plates
separated by a height $H$ and driven by a constant body force. Periodic
boundary conditions along the streamwise direction and spanwise directions and
a half-way bounce-back boundary condition is applied on the fixed walls. The
flow is driven by a body force of $F_{x}=1.0\times 10^{-6}$, which represents
the pressure gradient $-\frac{\partial P}{\partial x}$. The characteristic
dimensionless Reynolds number Re for this problem may be written as
$\mbox{Re}=\rho(\frac{H}{2})^{n}u_{max}^{2-n}/\mu_{p}$, where $u_{max}$ is the
maximum velocity of the flow occurring midway between the two plates (see
below). We consider a grid resolution of $3\times 3\times 101$ for resolving
the domain of the flow of power law fluids with $\mbox{Re}=100$ and
$n=0.8,1.0$ and $1.5$, thereby encompassing both shear thinning and shear
thickening fluids. The analytical solution for this non-Newtonian flow problem
is given by

$u(z)=\frac{n}{n+1}\left(-\frac{1}{\mu_{p}}\frac{\partial P}{\partial
x}\right)^{1/n}\left[\left(\frac{H}{2}\right)^{(n+1)/n}-\left|z\right|^{(n+1)/n}\right],$

where the maximum fluid velocity is obtained using

$u_{max}=\frac{n}{n+1}\left(-\frac{1}{\mu_{p}}\frac{\partial P}{\partial
x}\right)^{1/n}\left(\frac{H}{2}\right)^{(n+1)/n}.$

Figure 3 presents a comparison between the computed velocity profiles for
different values of the power law index obtained using the 3D cascaded LBM
implemented on a D3Q19 lattice with the analytical solution. It is clear that
the numerical simulation results for different power law fluids are in very
good agreement with the analytical solution.

In addition, Fig. 4 demonstrates that our 3D cascaded LBM is second order
accurate for the simulation of power law fluids. In this grid convergence
test, we have used diffusive scaling to establish asymptotic convergence of
our scheme to the incompressible macroscopic non-Newtonian fluid flow
equations of shear thinning, Newtonian and shear thickening fluids.

###  4.2 Three-dimensional fully developed flow in a square duct

The fluid flow through a 3D square duct (see Fig. 2b) is considered as the
next example. While an analytical solution for such a problem does not exist
for power law fluid at all values of the power law index, an exact solution
for the Newtonian fluid flow case ($n=1$) exists (see below), which could be
used to test our 3D cascaded LB formulation on a D3Q19 lattice. The cross
section of the square duct is defined by $-a\leq y\leq a$ and $-a\leq z\leq
a$, where $a$ is the half width of the duct, and $x$ is the streamwise flow
direction. Periodic boundary conditions are applied at the inlet and outlet,
whereas a no-slip boundary condition is adopted on the wall boundaries by
using the standard half-way bounce back approach. The Reynolds number Re for
this problem is defined based on the maximum fluid velocity and the duct half
width. The flow is driven by applying a body force of $F_{x}=1.0\times
10^{-6}$, and a grid resolution of $3\times 45\times 45$ is employed to
represent the flow domain. The body force is implementation is discussed in
Ref. [20], who is adopted for the D3Q19 lattice. The analytical solution for
the velocity field of the fully developed square duct flow can be expressed in
terms of the following infinite series based on harmonic functions [32]:

$u(y,z)=\frac{16a^{2}F_{x}}{\rho\nu\pi^{3}}\sum_{n=1}^{\infty}(-1)^{(n-1)}\left[1-\frac{\cosh\left(\frac{(2n-1)\pi
z}{2a}\right)}{\cosh\left(\frac{(2n-1)\pi}{2}\right)}\right]\frac{\cos\left(\frac{(2n-1)\pi
y}{2a}\right)}{(2n-1)^{3}}.$

Figures 5a and 5b show the comparison of the surface contours of the velocity
field computed using the 3D cascaded LBM using a D3Q19 lattice for a Reynolds
number of 80 and the analytical solution (see the equation above). It is clear
that the computed results are in good agreement with the exact solution, which
shows a paraboloid distribution for the velocity profile. In order to get a
better perspective for comparison, Figs. 6 and 8 show the comparisons of the
computed velocity profiles at different locations in the duct for
$\mbox{Re}=80$ and $20$, respectively, which show excellent agreement of our
3D cascaded LBM with the analytical solution. In addition, the iso-speed
contours of the velocity field at these two Re are presented in Figs. 7 and 9,
respectively, which are consistent with expectations.

(a)

(b)

Figure 2: Schematics of the non-Newtonian benchmark flow problems. (a) Power
law fluid flow in a lid-driven cubic cavity, (b) 3D fully developed flow in a
square duct.  (a)
Figure 3: Comparison of the normalized velocity profiles of power law fluids
in a channel computed using 3D Cascaded LBM against the analytical solution
for power law index $n=0.8,1.0,$ and $1.5$. at $\mbox{Re}=100$. (a)
Figure 4: Relative global error against the variation in the grid spacing,
plotted in log-log scales, for the simulation of power law fluid flow in a
channel at $n=0.8$, $n=1.0$ and $1.5$ obtained using the 3D cascaded LBM. The
solid line (black) represents an ideal slope of -2.0

(a)

(b)

Figure 5:  Flow through a square duct with side length 2a subjected to a
constant body force: Comparison of surface contours of the velocity field for
Reynolds number Re = 80 (a) computed by the D3Q19 formulation of the cascaded
LBM with forcing term with (b) analytical solution. (a)  Figure 6: Flow
through a square duct with side length 2a subjected to a constant body force:
Comparison of velocity profiles computed by the D3Q19 formulation of the
cascaded LBM with forcing term (symbols) with analytical solution (lines) at
different locations in the duct cross-section for Reynolds number Re = 80. (a)
Figure 7:  Iso-speed contours of the duct velocity field for Reynolds number
Re = 80. (a)  Figure 8: Flow through a square duct with side length 2a
subjected to a constant body force: Comparison of velocity profiles computed
by the D3Q19 formulation of the cascaded LBM with forcing term (symbols) with
analytical solution (lines) at different locations in the duct cross-section
for Reynolds number Re = 20. (a)
Figure 9:  Iso-speed contours of the duct velocity field for Reynolds number
Re = 20.

(a) Re=100

(b) Re=400

(c) Re=1000

Figure 10:  Computed velocity component $u/U_{p}$ along the vertical
centerline and component $v/U_{p}$ along the horizontal centerline obtained
using the cascaded LB method (lines) and compared with the benchmark solution
of Jin et al (2017) (symbols) for lid-driven cubic cavity flow of power law
fluids. Reynolds number $Re=100,400,1000$ and the power law index $n=1.5$.

(a) Re=100

(b) Re=400

(c) Re=1000

Figure 11:  Computed velocity component $u/U_{p}$ along the vertical
centerline and component $v/U_{p}$ along the horizontal centerline obtained
using the cascaded LB method (lines) and compared with the benchmark solution
of Jin et al (2017) (symbols) for lid-driven cubic cavity flow of power law
fluids. Reynolds number $Re=100,400,1000$ and the power law index $n=1.0$.

(a) Re=100

(b) Re=400

(c) Re=1000

Figure 12:  Computed velocity component $u/U_{p}$ along the vertical
centerline and component $v/U_{p}$ along the horizontal centerline obtained
using the cascaded LB method (lines) and compared with the benchmark solution
of Jin et al (2017) (symbols) for lid-driven cubic cavity flow of power law
fluids. Reynolds number $Re=100,400,1000$ and the power law index $n=0.8$.

###  4.3 3D Power-law fluid flow in a lid-driven cubic cavity

The 3D lid-driven flow of non-Newtonian fluids in a cubic cavity is a classic
benchmark problem with complex flow features and is a stringent test to
validate numerical methods. The geometric configuration of this problem is
shown in Fig. 2a, where a power-law fluid is enclosed within a cubic cavity of
side length $H$. The upper surface of the cavity or the lid is subjected to a
uniform velocity $U_{p}$, while all the remaining walls of the cavity are
maintained stationary. As noted in the introduction, only very few studies
exist on the 3D non-Newtonian flows of power law fluids in cubic cavities
driven by its top lid. A recent work reported by [29] presented high quality
results based on fine enough grid resolutions using a finite volume method,
which will be used as a benchmark solution for comparison.

We employ a grid resolution of $128\times 128\times 128$ in our simulations,
which was found to lead to grid independent results for the Reynolds number
cases considered in this regard. The no-slip boundary conditions on the walls
are imposed using the standard half-way bounce back scheme, with the motion of
the lid represented by applying a momentum correction to the respective bounce
back condition. The characteristic Reynolds number for this problem is defined
by $\mbox{Re}=\rho H^{n}U_{p}^{2-n}/\mu_{p}$. We consider the simulation of
cubic cavity flow at three different Reynolds numbers of $\mbox{Re}=100,400$
and $1000$. In each case, we consider three types of non-Newtonian fluids with
power law index $n=1.5,1.0$ and $0.8$, thereby encompassing both shear
thickening and shear thinning cases. In all these cases the flow was found to
be stationary and hence all the simulations were run until steady state was
reached.

The numerical results for the computed profiles of the velocity components
along the horizontal and vertical centerlines in the plane $z=0.5H$ at the
above three Re for the power law $n=1.5,1.0$ and $0.8$ are shown in Figs. 10,
11 and 12, respectively. While for the Newtonian fluid, the viscosity is a
constant, the local viscosity varies appreciably depending on the strain rate
tensor components that changes due to variations in the vortical flow
structures inside the cubic cavity. As a result, in the case of shear
thickening fluid, as the flow encounters high shear zones near the moving lid,
its effective viscosity increases in these regions, diffusing its momentum
(and conversely for the shear thinning case) and altering the velocity
profiles in the respective cases. Hence, the boundary layer becomes thicker
for fluid with $n>1$ when compared to the case where $n<1$. In addition, it is
found that as the power law index $n$ increases, the peak magnitude of the
velocity component $v$ along the horizontal centerline increases. These
observations are consistent with the recent findings of [29]. Overall, the
computed results obtained using the present 3D central moment LBM are in good
quantitative agreement with the benchmark results of [29].

##  5 Comparison of numerical stability of different collision models for 3D
non-Newtonian flow simulations

We now make a direct comparison of the numerical stability between the popular
single relaxation time (SRT) model and the cascaded MRT model presented in
this work, both using the D3Q19 lattice, for the simulation of 3D lid-driven
cubic cavity flow simulations. Due to it being a shear driven flow and
accompanied by geometric singularity around the corners, this flow
configuration serves as a good benchmark problem to test the stability of the
numerical schemes. For both these approaches, for a given resolution and a
fixed lid velocity, the relaxation parameter $\tau$ for the second order
moments that controls viscosity was decreased gradually until the computation
becomes unstable. Hence, the numerical instability is deemed to occur when the
global error of the velocity field becomes exponentially large. Figure 13
shows the maximum Reynolds number that could be attained before the
computations become unstable for both the collision models. Results are
presented for three different values of the power law index ($n=0.8,1.0$ and
$1.5$) encompassing both shear thinning and shear thickening fluids at various
set of grid resolutions ($48\times 48\times 48$, $64\times 64\times 64$,
$96\times 96\times 96$ and $128\times 128\times 128$). It is clear that the
cascaded LBM computations achieve significantly higher Reynolds numbers than
the SRT model for the same grid resolution and for a given type of fluid or a
power law index. Thus, we have demonstrated that the 3D cascaded LB
formulation is more stable than the SRT-LBM to simulate 3D non-Newtonian
flows.

(a) n=0.8

(b) n=1.0

(c) n=1.5

Figure 13: Comparison of the maximum Reynolds number for numerical stability
of the 3D SRT-LBM and the 3D cascaded LBM based on central moments for
simulation of non-Newtonian lid driven cubic cavity flow with $n=0.8,1.0,1.5$.

##  6 Summary and Conclusions

Complex fluid flows satisfying nonlinear constitutive relationships arise in
numerous engineering, chemical and materials processing applications and
geophysical situations. Simulations of such flows pose numerical challenges
especially in 3D, whose adequate resolution are associated with high
computational demands. Lattice Boltzmann methods (LBM) are inherently parallel
approaches for flow simulations. Prior efforts in using such kinetic schemes
have mainly focused on computing non-Newtonian flows in 2D, and generally
utilizing less robust underlying collision models. In this work, we present a
new 3D cascaded LBM for the simulation of non-Newtonian power law fluids on a
D3Q19 lattice. The effect of collisions are prescribed in terms of changes in
different central moments using multiple relaxation times, whose derivation is
given in detail. In particular, the relaxation times of the second order
moments are related to the local shear rate and parameterized by the model
coefficient and exponent representing the power law fluids. We performed a
validation of our 3D central moment LBM for various benchmark problems
including the complex non-Newtonian flow inside a cubic cavity driven by its
top lid at Reynolds numbers of 100, 400 and 1000, and with the power law index
specified to be 0.8, 1.0 and 1.5. Comparisons against recent numerical
benchmark solutions demonstrate very good accuracy. In addition, our 3D
cascaded LBM is shown to have superior numerical stability when compared to
the SRT-LBM and exhibits second order accuracy under grid refinement for
simulation of 3D shear driven flows of both shear thinning and shear
thickening fluids. Thus, the advanced cascaded LB formulation based on central
moments appears to be promising for simulations of complex 3D non-Newtonian
fluid flows.

## Acknowledgements

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

## References

## References

  * [1] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual review of fluid mechanics 30 (1) (1998) 329–364.
  * [2] S. Succi, S. Succi, The lattice Boltzmann equation: for fluid dynamics and beyond, Oxford university press, 2001.
  * [3] C. Aidun, J. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472.
  * [4] T. Krüger, et al., The lattice Boltzmann method: Principles and practice. sl: Springer international publishing (2017).
  * [5] Y. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, EPL (Europhysics Letters) 17 (6) (1992) 479.
  * [6] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems, Physical review 94 (3) (1954) 511.
  * [7] D. d’Humières, Multiple–relaxation–time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 360 (1792) (2002) 437–451.
  * [8] M. Geier, J. Greiner, F. Korvink, Cascaded digital lattice Boltzmann automata for high Reynolds number flow, Phys. Rev. E 73 (2006) 066705.
  * [9] P. Asinari, Generalized local equilibrium in the cascaded lattice Boltzmann method, Phys. Rev. E 78 (2008) 016701.
  * [10] K. N. Premnath, S. Banerjee, Incorporating forcing terms in cascaded lattice Boltzmann approach by method of central moments, Physical Review E 80 (3) (2009) 036702.
  * [11] M. Geier, M. Schonherr, A. Pasquali, M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory and validation, Comp. Math. Appl. 704 (2015) 507–547.
  * [12] Y. Ning, K. Premnath, D. Patil, Numerical study of the properties of the central moment lattice Boltzmann method, Int. J. Numer Meth Fluids 82 (2015) 59–90.
  * [13] K. N. Premnath, S. Banerjee, On the three-dimensional central moment lattice Boltzmann method, Journal of Statistical Physics 143 (4) (2011) 747–794.
  * [14] A. D. Rosis, Non-orthogonal central moments relaxing to a discrete equilibrium, Europhys. Lett. 116 (2016) 44003.
  * [15] F. Hajabdollahi, K. N. Premnath, Galilean-invariant preconditioned central-moment lattice Boltzmann method without cubic velocity errors for efficient steady flow simulations, Physical Review E 97 (5) (2018) 053303.
  * [16] L. Fei, K. H. Luo, C. Lin, Q. Li, Modeling incompressible thermal flows using a central-moments-based lattice Boltzmann method, International Journal of Heat and Mass Transfer 120 (2018) 624–634.
  * [17] F. M. Elseid, S. W. Welch, K. N. Premnath, A cascaded lattice Boltzmann model for thermal convective flows with local heat sources, International Journal of Heat and Fluid Flow 70 (2018) 279–298.
  * [18] F. Hajabdollahi, K. N. Premnath, Central moments-based cascaded lattice Boltzmann method for thermal convective flows in three-dimensions, Int. J. Heat Mass Transfer 120 (2018) 838 – 850.
  * [19] L. Fei, K. H. Luo, Q. Li, Three-dimensional cascaded lattice Boltzmann method: Improved implementation and consistent forcing scheme, Physical Review E 97 (5) (2018) 053309.
  * [20] F. Hajabdollahi, K. N. Premnath, Symmetrized operator split schemes for force and source modeling in cascaded lattice Boltzmann methods for flow and scalar transport, Physical Review E 97 (6) (2018) 063303.
  * [21] R. I. Tanner, Engineering rheology, Vol. 52, Oxford University Press, Oxford, 2000\.
  * [22] M. Deville, T. B. Gatski, Mathematical modeling for complex fluids and flows, Springer Science & Business Media, 2012.
  * [23] T. Phillips, G. W Roberts, Lattice Boltzmann models for non-Newtonian flows 76 (2011) 790–816.
  * [24] S. Gabbanelli, G. Drazer, J. Koplik, Lattice Boltzmann method for non-Newtonian (power-law) fluids, Physical review E 72 (4) (2005) 046312.
  * [25] Z. Chai, B. Shi, Z. Guo, F. Rong, Multiple-relaxation-time lattice Boltzmann model for generalized Newtonian fluid flows, Journal of Non-Newtonian Fluid Mechanics 166 (5-6) (2011) 332–342.
  * [26] D. Conrad, A. Schneider, M. Böhle, Accuracy of non-Newtonian lattice Boltzmann simulations, Journal of Computational Physics 301 (2015) 218–229.
  * [27] M. Reddy, J. Reddy, Finite-element analysis of flows of non-Newtonian fluids in three-dimensional enclosures, International journal of non-linear mechanics 27 (1) (1992) 9–26.
  * [28] R. N. Elias, M. A. Martins, A. L. Coutinho, Parallel edge-based solution of viscoplastic flows with the SUPG/PSPG formulation, Computational Mechanics 38 (4-5) (2006) 365–381.
  * [29] K. Jin, S. Vanka, R. Agarwal, B. Thomas, GPU accelerated simulations of three-dimensional flow of power-law fluids in a driven cube, International Journal of Computational Fluid Dynamics 31 (1) (2017) 36–56.
  * [30] M. Geier, A. Greiner, J. Korvink, A factorized central moment lattice Boltzmann method, The European Physical Journal Special Topics 171 (1) (2009) 55–61.
  * [31] S. Chapman, T. G. Cowling, The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases, Cambridge university press, 1990.
  * [32] F. M. White, I. Corfield, Viscous fluid flow, Vol. 3, McGraw-Hill New York, 2006\.