# Curved boundary conditions of the lattice Boltzmann method for simulating
microgaseous flows in the slip regime

Liang Wang  Beijing Key Laboratory of Emission Surveillance and Control for
Thermal Power Generation, North China Electric Power University, Beijing
102206, China  School of Energy Power and Mechanical Engineering, North China
Electric Power University, Beijing 102206, China  Shi Tao  Guangdong
Provincial Key Laboratory of Distributed Energy Systems, Dongguan University
of Technology, Dongguan 523808, China  Junjie Hu  Faculty of Engineering,
China University of Geosciences, Wuhan 430074, China  Kai Zhang  Beijing Key
Laboratory of Emission Surveillance and Control for Thermal Power Generation,
North China Electric Power University, Beijing 102206, China  School of Energy
Power and Mechanical Engineering, North China Electric Power University,
Beijing 102206, China  Gui Lu<EMAIL_ADDRESS>Key Laboratory of Power
Station Energy Transfer Conversion and System of Ministry of Education, North
China Electric Power University, Beijing 102206, China

###### Abstract

The lattice Boltzmann method (LBM) has shown its promising capability in
simulating microscale gas flows. However, the suitable boundary condition is
still one of the critical issues for the LBM to model microgaseous flows
involving curved geometries. In this paper, a local boundary condition of the
LBM is proposed to treat curved solid walls of microgaseous flows. The
developed boundary treatment combines the Maxwellian diffuse reflection scheme
and a single-node boundary scheme which contains a free parameter as well as
the distance ratio. The curved boundary condition is analyzed within the
multiple-relaxation-time (MRT) model for a unidirectional microflow. It is
shown that the derived slip velocity depends on the free parameter as well as
the distance ratio and relaxation times. By virtue of the free parameter, the
combination parameter and the uniform relaxation time are theoretically
determined to realize the accurate slip boundary condition. In addition, it is
found that besides the halfway diffuse-bounce-back (DBB) scheme, previous
curved boundary schemes only containing the distance ratio cannot ensure
uniform relaxation times to realize the slip boundary condition. Some
numerical examples with planar and curved boundaries are carried out to
validate the present curved boundary scheme. The good and robust consistency
of numerical predictions with analytical solutions demonstrates our
theoretical analysis.

##  I Introduction

With the rapid progress of microscience and nanotechnology, microscale gaseous
flows have attracted increasing research attention in recent years Ho98 ;
Karniadakis02 ; Zhang12 . For such gas flows, the mean free path of gas
moleculars ($\lambda$) is usually comparable to the characteristic length
scale of the flow system ($H$). As such, the gas flows are far from the
thermodynamical equilibrium, and the classical Navier-Stokes equations based
on the continuum assumption are no longer valid. Generally, the Knudsen number
$\text{Kn}=\lambda/H$ is used to measure the departure degree of microscale
gaseous flows from the classic hydrodynamic theory. Based on the magnitude of
Kn, gas flows could be classified into four major regimes: the continuum flow
regime with $\text{Kn}\leq 0.001$, the slip flow regime with
$0.001<\text{Kn}\leq 0.1$, the transition flow regime with $0.1<\text{Kn}\leq
10$, and the free-molecular flow regime with $\text{Kn}>10$. It is widely
accepted that the Boltzmann equation can work for rarefied gas flows with
arbitrary Knudsen numbers Cercignani90 . Therefore, as a discrete scheme
derived from the Boltzmann equation HeX98 ; ShanX98 , the lattice Boltzmann
method (LBM) has been believed to have the potential for simulating microscale
gaseous flows. Since being applied to microscale gaseous flows first in 2002
Nie02 ; Lim02 , the LBM has received increasing interests over the past dozen
years Shu05 ; ZhangY05 ; ZhangJ05 ; Kima08 .

With the increase of Kn, the Navier-Stokes equations with the no-slip boundary
condition become invalid for gas flows, while the fluid slippage on a solid-
fluid boundary will arise due to the gaseous nonequilibrium effect Sone02 ;
Sharipov11 . Therefore, for the LBM to capture the gaseous slip phenomenon in
simulating microgaseous flows, the boundary condition is of critical
importance for effectively predicting gas-solid interactions. This important
issue has attracted substantial researches toward accurate slip boundary
treatments. Nie _et al_. Nie02  employed the standard bounce-back (BB)
boundary condition to predict the flows in microchannels, and found a nonzero
velocity proportional to the square of Kn on the channel wall. While, the slip
velocity was later revealed to be a numerical artifact Verhaeghe09  actually.
Lim _et al._ Lim02  investigated the specular reflection (SR) boundary
condition in microchannel flow simulations, and the slip velocity was not well
consistent with some existing analytical solutions. Ansumali and Karlin
Ansumali02  applied the Maxwellian diffusive (MD) boundary condition for the
Kramer’s problem, but the scheme tends to overpredict the slip velocity.
Therefore, as revealed in the subsequent researches, the above pure boundary
schemes cannot accurately capture the slip phenomenon.

The shortcomings in the above scenario have thus stimulated improving hybrid
schemes which combine the pure boundary conditions with an accommodation
coefficient Verhaeghe09 ; Succi02 ; Sbragaglia05 ; Tang05 . By mixing the BB
and the SR boundary conditions, Succi Succi02  proposed a hybrid scheme, i.e.,
the bounce-back and specular reflection (BBSR) scheme for simulating
microscale flows. A generalization of the BBSR scheme was subsequently
developed and analyzed by Sbragaglia and Succi Sbragaglia05 . Tang _et al._
Tang05  proposed the diffusive and the specular reflection (DSR) scheme, which
comes from the combination of the MD and SR boundary schemes. Thanks to the
tunable accommodation parameter, the degree of slip can be freely controlled
to recover different slip models in the two hybrid schemes. Another hybrid
boundary scheme is the diffusive and bounce-back (DBB) scheme Verhaeghe09 ;
Chai08 , which is a combination of the MD and the BB boundary conditions.
Owing to the advantage of local computation superior to the BBSR and DSR
schemes, the DBB scheme may possess more potential in simulating microscale
gaseous flows with complex geometries. For hybrid slip boundary conditions,
how to choose the combination parameter is not a convenient task. With the
Bhatnagar-Gross-Krook (BGK) model and the multiple-relaxation-time (MRT)
model, Guo _et al._ Guo07 ; Guo08  successively analyzed the discrete effects
of the BBSR and DSR boundary conditions, and found the interrelationship
between the two schemes. To realize the accurate slip boundary condition, they
proposed a strategy to determine the combination parameter. Verhaeghe _et al._
Verhaeghe09  and Chai _et al._ Chai08 ; Chai10  further mathematically
analyzed the discrete effect of the DBB scheme. Their results illustrated that
the combination parameter should be chosen carefully to impose the accurate
slip boundary condition. Noteworthily, these hybrid boundary schemes are
originally designed to implement the wall location with a definite distance
between lattice nodes. Szalm$\acute{a}$s Szal06  used an interpolated method
to combine the BB and the SR schemes, and it permits arbitrary locations of
the slip wall. However, this slip boundary condition is proposed for
microgaseous flows specially with straight walls. Therefore, for microscale
slip flows with curved boundaries, which bring variable wall locations in
relation to the underlying grid, the above-mentioned studies theoretically
cannot yield enough accurate results especially at a small grid resolution.

For a more accurate treatment of curved boundaries in microgaseous flows,
there have been developed a number of boundary conditions considering the
actual boundary shape in the literature. Suga Suga13  introduced an
interpolation method into the DBB boundary condition for microscale flow
simulations. The effect of curved surface is represented by its intersection
distance ratio between neighbouring lattice nodes. But, the combination
coefficient is still determined by that from the halfway bounce-back case.
Based on the non-equilibrium extrapolation method for curved boundaries and
the counter-extrapolation method for the velocity/temperature at curved
surfaces, Liu _et al._ LiuZ19  proposed a boundary condition involving the
distance ratio for thermal gaseous microflows with curved slip walls. However,
as revealed in the literature Tao15 ; Silva17 ; Silva18 , the discrete effects
also exist in curved boundary treatments for slip walls, which must be
minimized to capture correct microgaseous slip phenomenon. Within the MRT
model, Tao and Guo Tao15  incorporated the effect of distance ratio to analyze
the DBB scheme, and then developed a boundary scheme to realize the slip
boundary condition at curved boundaries. Silva and Semiao Silva17 ; Silva18
introduced the multireflection framework to gaseous slip flows, and put
forward curved slip boundary schemes within the two-relaxation-time (TRT)
model. The theoretical analysis on curved boundary schemes revealed that the
numerical accuracy is related with the combination parameter and the
relaxation times, which are functions of the wall cut-link distance ratio
Tao15 ; Silva17 ; Silva18 . Even for microscale binary gaseous flows, such
results are also exposed as the DBB scheme is analyzed to derive the slip
velocity at curved surfaces Ren19 . In the LBM for microflows with curved
walls, previous studies have indicated that the combination parameter can
change locally with the wall cut-link distance ratio. However, the relaxation
parameter should also vary with the distance ratio to realize the slip
boundary condition, which drives the anisotropic collision operator.
Unfortunately, such issue cannot be overcome with previous curved boundary
schemes which only contain the distance ratio.

To attain uniform relaxation parameters as noted above, one natural strategy
is to impose the halfway boundary scheme for curved slip walls, and the actual
locations between lattice nodes are approximated as halfway wall locations
Chai10 . In this way, the uniform relaxation parameter can be then determined
by the distance ratio fixed at $0.5$. Following this line of thought, there
have been some woks reported to treat the curved slip boundary by the halfway
DBB scheme Verhaeghe09 ; Chai10 ; Guo11 . However, the real curved geometry
with the halfway approximation will lost its fidelity under coarse grid
resolutions, and undesired errors may arise to contaminate the simulation
accuracy Silva17 ; Silva18 . Therefore, based on the above literature review,
it promotes us to resolve such a critical issue for microgaseous flows with
curved boundaries: how to retain the relaxation parameters unchanged to
realize the slip boundary condition. Accordingly, the present work has the
following twofold objectives: to develop a curved kinetic boundary condition
which involves additional parameters besides the distance ratio; and to
realize the accurate slip boundary condition while with invariant relaxation
parameters.

In this work, a kinetic boundary condition is first constructed by combining
an interpolation-based scheme and the diffuse reflection rule. Different from
previous curved boundary conditions, the present scheme is inspired by the
idea in Ref. Zhao19  to include a free parameter besides the distance ratio,
which can bring infinitely many curved boundary schemes for microgaseous
flows. As the subsequent step, the boundary scheme is mathematically analyzed
within the MRT model for a unidirectional flow. On the basis of theoretical
derivations, the combination parameter and a strategy to ensure uniform
relaxation parameters are given to realize the slip boundary condition at
curved slip walls. Numerical simulations are then carried out to validate the
developed method in the cases of aligned and inclined flat walls and curved
walls. As compared with the halfway DBB boundary scheme and other curved
boundary schemes with several specific free parameters, the present method can
predict the most accurate results consistent with the analytical solutions.

##  II MRT LBE for microscale gaseous flows

The LBM is derived from the discretization of continuous Boltzmann equation in
both time, space and velocity space. The discrete velocity distribution
functions evolve according to the following lattice Boltzmann equation (LBE),

$f_{i}(\bm{x}+\bm{c}_{i}\delta_{t},t+\delta_{t})-f_{i}(\bm{x},t)=\Omega_{i}(f)(\bm{x},t)+\delta_{t}F_{i}(\bm{x},t),\quad
i=0,1,\cdots,b-1,$ (1)

where $f_{i}(\bm{x},t)$ is the distribution function associated with the
discrete velocity $\bm{c}_{i}$ at position $\bm{x}$ and time $t$, $\delta_{t}$
is the time step, $\Omega_{i}(f)$ is the discrete collision operator, and
$F_{i}$ is the discrete forcing term, and $b$ is the number of discrete
velocities.

The BGK or single-relaxation-time model is the most widely used collision
operator in the LBM. However, the slip velocity derived within the BGK model
depends on the relaxation time, which is grid resolution dependent for a given
Kn Verhaeghe09 ; Guo08 . This means that some unphysical numerical artifacts
besides the physical part exist in the slip velocity. To avoid this problem,
in this work we turn to employ the MRT collision operator adhered to the LBE,
which is written as

$\Omega_{i}(f)=-\sum_{j}(\bm{M}^{-1}\bm{S}\bm{M})_{ij}\left[f_{j}-f_{j}^{(\text{eq})}\right],$
(2)

where $\bm{M}$ is a $b\times b$ transformation matrix, which maps $f_{i}$ onto
the moment space via $\bm{m}=\bm{M}\bm{f}$ with
$\bm{f}=(f_{0},f_{1},\cdots,f_{b-1})^{\text{T}}$,
$\bm{S}=\text{diag}(\tau_{0},\tau_{1},\cdots,\tau_{b-1})^{-1}$ is a diagonal
relaxation matrix with its non-negative element $\tau_{i}$ being the
relaxation time for the $i$-th moment. $f_{j}^{(\text{eq})}$ is the
equilibrium distribution function which is dependent on the gas density
$\rho$, velocity $\bm{u}$ and temperature $T$,

$f_{j}^{(\text{eq})}=\omega_{j}\rho\left[1+\frac{\bm{c}_{j}\cdot\bm{u}}{c_{s}^{2}}+\frac{(\bm{c}_{j}\cdot\bm{u})^{2}}{2c_{s}^{4}}-\frac{\bm{u}^{2}}{2c_{s}^{2}}\right],\quad
j=0,1,\cdots,b-1,$ (3)

where $\omega_{j}$ is the weight coefficient, $c_{s}=\sqrt{RT}$ ($R$ is the
gas constant) is the lattice sound speed. For isothermal flows, $c_{s}$ is
determined by the lattice speed $c=\delta_{x}/\delta_{t}$ with $\delta_{x}$
being the lattice spacing. For the discrete forcing term $F_{i}$ in Eq. (1),
it should be taken as Guob13

$\bm{F}=\bm{M}^{-1}\left(\bm{I}-\frac{\bm{S}}{2}\right)\bm{M}\bm{\overline{F}},$
(4)

where $\bm{I}$ is the identity matrix,
$\bm{F}=(F_{0},F_{1},\cdots,F_{b-1})^{\text{T}}$, and
$\bm{\overline{F}}=(\overline{F}_{0},\overline{F}_{1},\cdots,\overline{F}_{b-1})^{\text{T}}$
is expressed as

$\overline{F}_{i}=\omega_{i}\rho\left[\frac{\bm{c}_{j}\cdot\bm{a}}{c_{s}^{2}}+\frac{\bm{u}\bm{a}:(\bm{c}_{i}\bm{c}_{i}-c_{s}^{2}\bm{I})}{c_{s}^{4}}\right],$
(5)

where $\bm{G}=\rho\bm{a}$ is the external force.

In this work, we employed the two-dimensional nine-velocity (D2Q9) model,
where the discrete velocities $\bm{c}_{i}$ are defined by

$\bm{c}_{i}:=c\bm{e}_{i}=\begin{cases}c(0,0),&i=0,\\\
c\bigl{(}\cos\bigl{[}(i-1)\pi/2\bigr{]},\sin\bigl{[}(i-1)\pi/2\bigr{]}\bigr{)},&i=1,2,3,4,\\\
\sqrt{2}c\bigl{(}\cos\bigl{[}(i-1)\pi/2+\pi/4\bigr{]},\sin\bigl{[}(i-1)\pi/2+\pi/4\bigr{]}\bigr{)},&i=5,6,7,8,\end{cases}$
(6)

where $c=\delta_{x}/\delta_{t}$ with $\delta_{x}$ denoting the lattice
spacing. Correspondingly, the sound speed $c_{s}=\sqrt{RT}=c/\sqrt{3}$, and
the weight coefficients are given by $\omega_{0}=4/9$, $\omega_{1-4}=1/9$ and
$\omega_{5-8}=1/36$. Via the Gram-Schmidt orthogonalization procedure on the
discrete velocities $\bm{c}_{i}$, there are different versions of
transformation matrices $\bm{M}$, and one form of $\bm{M}$ as $c=1$ is given
by Lallemand00

$\displaystyle\bm{M}=\left(\begin{array}[]{ccccccccc}1&1&1&1&1&1&1&1&1\\\
-4&-1&-1&-1&-1&2&2&2&2\\\ 4&-2&-2&-2&-2&1&1&1&1\\\ 0&1&0&-1&0&1&-1&-1&1\\\
0&-2&0&2&0&1&-1&-1&1\\\ 0&0&1&0&-1&1&1&-1&-1\\\ 0&0&-2&0&2&1&1&-1&-1\\\
0&1&-1&1&-1&0&0&0&0\\\ 0&0&0&0&0&1&-1&1&-1\end{array}\right).$ (16)

As a result, the nine discrete velocity moments $m_{i}$ from the distribution
functions $f_{i}$ are expressed as

$\bm{m}=\bm{M}\bm{f}=(\rho,e,\varepsilon,j_{x},q_{x},j_{y},q_{y},p_{xx},p_{xy})^{\text{T}}.$
(17)

The corresponding relaxation matrix $\bm{S}$ for the nine moments is written
as

$\bm{S}=\text{diag}(\tau_{\rho},\tau_{e},\tau_{\varepsilon},\tau_{j},\tau_{q},\tau_{j},\tau_{q},\tau_{s},\tau_{s})^{-1}.$
(18)

The fluid density $\rho$ and velocity $\bm{u}=(u,~{}v)$ are respectively
defined as the zeroth and first-order moments of $f_{i}$

$\rho=\sum_{i}f_{i},\quad\rho\bm{u}=\sum_{i}\bm{c}_{i}f_{i}+\frac{\delta_{t}}{2}\rho\bm{a}.$
(19)

Through the Chapman-Enskog or linear analysis, the Navier-Stokes equations can
be derived from the MRT-LBE model (1). The fluid pressure $p$ is determined by
$p=c_{s}^{2}\rho=\rho RT$, and the shear and bulk viscosities are respectively
given by

$\nu=c_{s}^{2}\left(\tau_{s}-\frac{1}{2}\right)\delta_{t},\qquad\zeta=c_{s}^{2}\left(\tau_{e}-\frac{1}{2}\right)\delta_{t}.$
(20)

The numerical implementation of Eq. (1) is divided into two steps:

$\displaystyle\text{Collision}:\quad
f_{i}^{\prime}(\bm{x},t)=f_{i}(\bm{x},t)-\sum_{j}(\bm{M}^{-1}\bm{S}\bm{M})_{ij}\left[f_{j}(\bm{x},t)-f_{j}^{(\text{eq})}(\bm{x},t)\right]+\delta_{t}F_{i}(\bm{x},t),$
$\displaystyle\text{Streaming}:\quad
f_{i}(\bm{x}+\bm{c}_{i}\delta_{t},t+\delta_{t})=f_{i}^{\prime}(\bm{x},t),$
(21)

where $f_{i}^{\prime}(\bm{x},t)$ is the post-collision distribution function.
In the framework of MRT model, the collision step is usually implemented in
the moment space

$\begin{split}\bm{m}^{\prime}(\bm{x},t)&=\bm{m}(\bm{x},t)-\bm{S}\left[\bm{m}(\bm{x},t)-\bm{m}^{(\text{eq})}(\bm{x},t)\right]+\delta_{t}\bm{\hat{F}}(\bm{x},t),\\\
\bm{f}^{\prime}(\bm{x},t)&=\bm{M}^{-1}\bm{m}^{\prime}(\bm{x},t),\end{split}$
(22)

where $\bm{m}^{\prime}:=\bm{M}\bm{f}^{\prime}$ with
$\bm{f}^{\prime}=(f_{0}^{\prime},f_{1}^{\prime},\cdots,f_{8}^{\prime})^{\text{T}}$
is the post-collision moment,
$\bm{m}^{(\text{eq})}:=\bm{M}\bm{f}^{(\text{eq})}$ and
$\bm{\hat{F}}:=\bm{M}\bm{F}=\left(\bm{I}-\frac{\bm{S}}{2}\right)\bm{M}\bm{\overline{F}}$
denote the equilibria and forcing term in the moment space, respectively.
After the collision step is completed, the streaming step is then executed by
transforming the moments back to the velocity space.

In contrast to continuum flows, the LBM for microscale flows should consider
the rarefaction effects, which are commonly characterized by the dimensionless
Knudsen number. The most important characteristic parameter in continuum flows
is the Reynolds number Re, which determines the relaxation time $\tau_{s}$ in
the simulations. However, for microscale gas flows, the relationship between
$\tau_{s}$ and Kn should be carefully established to ensure the consistency
criterion Guo06 , which is read as

$\tau_{s}=\frac{1}{2}+\sqrt{\frac{6}{\pi}}\text{Kn}\frac{H}{\delta_{x}}.$ (23)

Another important issue for the LBM applied to microscale flows is the kinetic
boundary condition. As reviewed previously, there have been some boundary
conditions proposed for the LBE to simulate microgaseous flows, such as the
BBSR scheme, DSR scheme and DBB scheme. These kinetic boundary conditions are
originally designed for the case of flat walls with definite locations between
lattice nodes. Furthermore, it has been exposed in the literature Guo07 ;
Guo08  that some discrete effects exist in the boundary schemes and should be
corrected to realize the slip boundary condition. For a more accurate
treatment of curved walls, it is _common_ to include the distance ratio in the
boundary condition to preserve the actual shape of curved geometries Suga13 ;
LiuZ19 . However, when these curved boundary schemes are implemented to
realize the slip boundary condition at curved walls, the relaxation time in
the LBE will depend on the distance ratio and thus varies with link directions
at different boundary nodes Tao15 ; Silva17 ; Silva18 . Noteworthily, such a
problem cannot be overcome in previous curved boundary schemes only containing
the distance ratio.

##  III Kinetic boundary condition for curved slip walls

To remedy the above issue of nonuniform relaxation parameters, one natural and
promising way is to introduce additional parameters besides the distance ratio
into a kinetic boundary condition. As a matter of fact, by resorting to a
single-node boundary scheme Zhao19  for a free parameter besides the distance
ratio, we recently overcome the numerical slip at curved no-slip walls WangL20
with invariant relaxation parameters. Thus, it motivates and inspires us to
construct such a parameterized boundary condition for microflows with curved
walls, and then generalize the strategy to ensure uniform relaxation
parameters from continuum flows to microgaseous flows.

###  III.1 Curved boundary scheme for microscale gaseous flows

We now propose a kinetic boundary condition that involves an additional
parameter besides the distance ratio as noted above. Because a curved boundary
in simulations can be modeled separately for each lattice direction, we
consider a single direction $\bm{e}_{i}$ pointing from the boundary surface at
$\bm{x}_{b}$ to the boundary node $\bm{x}_{f}$, as shown in Fig. 1.

Figure 1: Schematic of a curved-wall boundary along one single lattice
direction. The thin solid line is the grid line, and the thick curved one
represents the boundary surface. White circles ($\circ$): the fluid nodes;
Black circle ($\bullet$): the intersection point of the boundary with the grid
line; Square box ($\square$): the solid node outside the computational domain.

The surface point $\bm{x}_{b}$ is intersected by $\bm{x}_{f}$ and the solid
node $\bm{x}_{r}$, and then the actual location of curved boundary can be
depicted by the distance ratio $\gamma$ as
$\gamma=|\bm{x}_{b}-\bm{x}_{f}|/|\bm{x}_{r}-\bm{x}_{f}|$. At the boundary node
$\bm{x}_{f}$, the unknown distribution function
$f_{i}(\bm{x}_{f},t+\delta_{t})$ in the LBM is specified according to the
kinetic boundary condition, which aims to realize the slip velocity condition
at physical walls. Noteworthily, for microscale flows with curved boundaries
or in complex geological porous media, a local boundary condition that
involves the current information at boundary nodes is desired especially. In
the present work, our construction aims to a local curved boundary condition
for microscale gas flows.

Recall that the DBB boundary condition, which has been developed for both
planar and curved slip walls Verhaeghe09 ; Chai08 ; Tao15 , pertains to a
local boundary scheme. This is due to its combination of the BB and the MD
schemes, which is formulated for the unknown distribution function as follows:

$f_{i}(\bm{x}_{f},t+\delta_{t})=r\left[f_{\bar{i}}^{\prime}(\bm{x}_{f},t)+2\omega_{i}\rho_{f}\frac{\bm{c}_{i}\cdot\bm{u}_{b}}{c_{s}^{2}}\right]+(1-r)f_{i}^{(\text{eq})}(\rho_{f},\bm{u}_{b}),$
(24)

where the subscript $\bar{i}$ indicates $\bm{e}_{\bar{i}}=-\bm{e}_{i}$,
$\rho_{f}=\rho(\bm{x}_{f},t)$ and $\bm{u}_{b}=\bm{u}(\bm{x}_{b},t)$ is the
wall velocity. The combination parameter $r$ represents the bounce-back
fraction and ranges in the region of $0\leq r\leq 1$. Only the lattice
directions of $\bm{e}_{i}$ and $\bm{e}_{\bar{i}}$ are involved at the current
node $\bm{x}_{f}$, and thus the DBB scheme shares the perfect feature of local
computation for curved walls. For this local scheme, the bounce-back part
accounts for the no-slip boundary condition, while the Maxwell diffuse part is
responsible for the gas slippage at solid walls. On the other hand, Zhao _et
al_. Zhao19  recently proposed a single-node boundary scheme, which contains a
free parameter besides $\gamma$, for the no-slip boundary condition. Later,
its ability to overcome the discrete effect with uniform parameters has been
revealed in our recent work WangL20 . Along with the underlying structure of
the DBB scheme, we thus replace the bounce-back part in Eq. (24) by the scheme
in Ref. Zhao19 , and derive the following boundary condition

$\displaystyle f_{i}(\bm{x}_{f},t+\delta_{t})$
$\displaystyle=r\left[\frac{1+l-2\gamma}{1+l}f_{\bar{i}}(\bm{x}_{f},t)+\frac{l}{1+l}f_{i}^{\prime}(\bm{x}_{f},t)+\frac{2\gamma-l}{1+l}f_{\bar{i}}^{\prime}(\bm{x}_{f},t)+\frac{2}{1+l}\omega_{i}\rho_{f}\frac{\bm{c}_{i}\cdot\bm{u}_{b}}{c_{s}^{2}}\right]$
$\displaystyle\qquad+(1-r)f_{i}^{(\text{eq})}(\rho_{f},\bm{u}_{b}),$ (25)

where $l$ is a free parameter besides $\gamma$ such that
$\bm{x}_{b}=(\bm{x}_{1}+\bm{x}_{2})/2$ (see Fig. 1). Clearly, the boundary
condition (III.1) is a local boundary scheme, and can preserve the geometry
fidelity of curved walls by the distance ratio $\gamma$. Since the no-slip
boundary scheme in Ref. Zhao19  is incorporated, the present boundary
condition is inherently associated with the diffusive scaling $\delta_{t}=\eta
h^{2},~{}h=\delta_{x}$ ($\eta$ is an adjustable parameter). Because $l$ takes
values along the same lattice direction as $\gamma$ does (cf Fig. 1), it can
be considered as a function of $\gamma$. In addition, the parameter $l$ ranges
in $\text{max}\\{0,2\gamma-1\\}\leq l\leq 2\gamma$ to ensure the convex
combination of distribution functions for the no-slip part.

In the LBM for microscale flows, previous curved boundary schemes only contain
the distance ratio $\gamma$, while the present boundary scheme (III.1)
introduces an adjustable parameter $l$ besides $\gamma$. With the variable
parameter $l$, the boundary scheme (III.1) can hence bring numerous boundary
conditions for microgaseous flows. More importantly, we will show later that
owing to the free parameter $l$, uniform relaxation parameters can be
fulfilled to realize an exact prescribed slip boundary condition. Some remarks
about the boundary curved scheme are given as follows:

######  Remark 1.

When $r=1$, the boundary scheme (III.1) will degenerate to the single-node
scheme proposed in Ref. Zhao19  for continuum flows with no-slip walls. From
this viewpoint, the present scheme is the generalized version of that in Ref.
Zhao19 .

######  Remark 2.

When $l=0$ and $\gamma=\frac{1}{2}$, the boundary scheme (III.1) then
degenerates to the DBB scheme (24). This also confirms the intrinsic halfway
consumption of the bounce-back part adhered to the boundary scheme (24).

######  Remark 3.

In the present work, the effect of curved wall geometry is directly embodied
in the boundary condition (III.1) by the distance ratio $\gamma$. While in
Ref. Tao15 , it is absorbed in the DBB scheme (24) through the combination
parameter to realize the slip boundary condition at curved walls.

As noted before, the combination parameter $r$ in hybrid boundary schemes
plays the vital role and directly affects the degree of slip at physical
boundaries Guo07 ; Guo08 . Therefore, it should be carefully chosen to derive
reasonable simulation results. Previous studies with the MRT model have shown
that $r$ is related with sever factors, including the relaxation time, the
gas-solid interaction parameter and the distance ratio Guo08 ; Tao15 . For the
present boundary scheme (III.1) with the free parameter $l$, the influence of
$r$ is expected to be more complicated in treating curved slip walls. Next, we
will conduct a theoretical analysis to investigate how to specify $r$ here,
and simultaneously ensure the above-mentioned uniform relaxation parameters.

###  III.2 Analysis of the curved boundary scheme

To simplify the analysis, we consider the steady incompressible Poiseuille
flow in a microchannel, which has been extensively employed for theoretical
analysis in previous studies Shu05 ; Verhaeghe09 ; Chai08 ; Guo08 ; Guo11 ;
Tao15 .

Figure 2: Schematic of the flow and lattice arrangement with a distance ratio
$\gamma$. The wall boundary corresponding to the halfway DBB boundary
condition is placed with $\gamma=1/2$.

As sketched in Fig. 2, the flow is driven by a constant force
$\rho\bm{a}=\rho(a,~{}0)$ along the $x$ direction, and is subjected to the
following assumptions:

$\rho=\text{const},\quad v=0,\quad\partial_{x}\phi=0,\quad\partial_{t}\phi=0,$
(26)

where $\phi$ is an arbitrary flow variable. The lattice node at the layer
$j=0$ is located with an arbitrary distance $\gamma\delta_{x}$ away from the
solid wall. Here $j$ is the index of the grid line at
$y_{j}=(j+\gamma)\delta_{x}$. After the streaming step, the unknown
distribution functions, $f_{2}$, $f_{5}$ and $f_{6}$ at $j=0$ are determined
according to the boundary scheme (III.1):

$f_{2}^{0}=r\left[\frac{1+l-2\gamma}{1+l}f_{4}^{0}+\frac{l}{1+l}f_{2}^{{}^{\prime}0}+\frac{2\gamma-l}{1+l}f_{4}^{{}^{\prime}0}+\frac{2}{1+l}\omega_{2}\rho\frac{\bm{c}_{2}\cdot\bm{u}_{b}}{c_{s}^{2}}\right]+(1-r)f_{2}^{(\text{eq})}(\rho,u_{w}),$
(27a)
$f_{5}^{0}=r\left[\frac{1+l-2\gamma}{1+l}f_{7}^{0}+\frac{l}{1+l}f_{5}^{{}^{\prime}0}+\frac{2\gamma-l}{1+l}f_{7}^{{}^{\prime}0}+\frac{2}{1+l}\omega_{5}\rho\frac{\bm{c}_{5}\cdot\bm{u}_{b}}{c_{s}^{2}}\right]+(1-r)f_{5}^{(\text{eq})}(\rho,u_{w}),$
(27b)
$f_{6}^{0}=r\left[\frac{1+l-2\gamma}{1+l}f_{8}^{0}+\frac{l}{1+l}f_{6}^{{}^{\prime}0}+\frac{2\gamma-l}{1+l}f_{8}^{{}^{\prime}0}+\frac{2}{1+l}\omega_{6}\rho\frac{\bm{c}_{6}\cdot\bm{u}_{b}}{c_{s}^{2}}\right]+(1-r)f_{6}^{(\text{eq})}(\rho,u_{w}),$
(27c)

where $u_{w}$ is the wall velocity, and $f_{i}^{0}=f_{i}(y_{0})$ and
$f_{i}^{{}^{\prime}0}=f_{i}^{\prime}(y_{0})$ with $y_{0}=\gamma\delta_{x}$.

Following the procedures exhibited in Refs. Guo08 ; WangL20 ; Guo08P , we can
obtain the relationship between the velocities $u_{0}$ and $u_{1}$
respectively at $j=0$ and $j=1$:

$u_{1}=\mathcal{A}u_{0}+\mathcal{B}a\delta_{t}+\mathcal{C}u_{w},$ (28)

where

$\mathcal{A}=\frac{(1+l)(1-r)\tau_{s}+2r(1+\gamma)}{(1+l)(1-r)\tau_{s}+2r\gamma},\qquad\mathcal{C}=-\frac{2r}{(1+l)(1-r)\tau_{s}+2r\gamma},$
(29a)
$\mathcal{B}=\frac{(1+r+l-lr)(1+4\tau_{q}-8\tau_{s}-8\tau_{q}\tau_{s})+12r[(2\gamma-l)(2\tau_{s}-1)+(1-2\gamma)\tau_{s}]}{2\left[(1+l)(1-r)\tau_{s}+2r\gamma\right](2\tau_{s}-1)}.$
(29b)

For the Poiseuille flow between two plates located at $y=0$ and $y=H$, the
analytical solution can be expressed as

$u_{j}=4u_{c}\frac{y_{j}}{H}\left(1-\frac{y_{j}}{H}\right)+u_{w}+u_{s},$ (30)

where $u_{c}=aH^{2}/8\nu$, and $u_{s}$ is the slip velocity at the physical
wall. Substituting Eq. (30) into Eq. (28) to replace $u_{0}$ and $u_{1}$, we
can derive the dimensionless slip velocity

$U_{s}:=\frac{u_{s}}{u_{c}}=\frac{2(1+l)(1-r)}{r}\tau_{s}\frac{\delta_{x}}{H}-\frac{1}{3r}\left[(1+l)(1+4\tau_{q}-2\tau_{s}+12\gamma\tau_{s}-8\tau_{s}\tau_{q})+r\mathcal{D}\right]\frac{\delta_{x}^{2}}{H^{2}}.$
(31)

where
$\mathcal{D}=12\gamma^{2}+12\gamma(\tau_{s}-1)+(2\tau_{s}-1)(4l\tau_{q}-1-4\tau_{q})-l[11(2\tau_{s}-1)+12\gamma\tau_{s}]$.
Further invoking the $\tau_{s}$-Kn relation given by Eq. (23), Equation (31)
can be rewritten in terms of Kn as

$U_{s}=4\frac{(1-r)(1+l)\tau_{s}}{r(2\tau_{s}-1)}\sqrt{\frac{6}{\pi}}\text{Kn}+\frac{8}{\pi}\frac{\mathcal{K}}{r(2\tau_{s}-1)^{2}}\text{Kn}^{2}$
(32)

where $\mathcal{K}=(2\tau_{s}-1)(1+r+l+11lr)+12\gamma[r+\tau_{s}(lr-
l-r-1)]-12r\gamma^{2}-4\tau_{q}(2\tau_{s}-1)(lr-l-r-1)$. The expression of Eq.
(32) is the slip velocity derived from the boundary scheme (III.1). Compared
with previous similar studies Szal06 ; Tao15 , the most striking feature here
is that the slip velocity $U_{s}$ includes an additional free parameter $l$
besides $\gamma$. Owing to the degree of freedom from $l$, we now could obtain
invariant relaxation parameters to realize the adopted slip boundary condition
by the boundary scheme (III.1).

For non-continuum gaseous flows, a second-order slip boundary condition is
widely used and read as Suga13 ; Hadji06

$u_{s}=\mathcal{L}_{1}\lambda\frac{\partial u}{\partial
n}-\mathcal{L}_{2}\lambda^{2}\frac{\partial^{2}u}{\partial n^{2}},$ (33)

where $\bm{n}$ is the unit outer vector normal to the wall, and
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ refer to the slippage coefficients
related with the gas-wall interactions. Under the slip boundary condition
(33), the Poiseuille flow exhibits the following slip velocity at the wall

$U_{s}=4\mathcal{L}_{1}\text{Kn}+8\mathcal{L}_{2}\text{Kn}^{2}.$ (34)

Therefore, to realize the slip boundary condition (34) with the boundary
scheme (III.1), $r$ and $\tau_{q}$ must be determined by comparing Eq. (32)
with Eq. (34), which are given by

$\displaystyle r$
$\displaystyle=\frac{\sqrt{\frac{6}{\pi}}\tau_{s}(1+l)}{\mathcal{L}_{1}(2\tau_{s}-1)+\sqrt{\frac{6}{\pi}}\tau_{s}(1+l)},$
(35) $\displaystyle\tau_{q}$
$\displaystyle=\frac{\mathcal{L}_{1}(2\tau_{s}-1)\left[2\tau_{s}(6\gamma-1)+1\right]+\sqrt{\frac{6}{\pi}}\tau_{s}\bigl{[}\mathcal{L}_{2}\pi(2\tau_{s}-1)^{2}+12\gamma^{2}+2(2\tau_{s}-1)(6\gamma-6l-1)\bigr{]}}{4(2\tau_{s}-1)\left[2\sqrt{\frac{6}{\pi}}\tau_{s}+\mathcal{L}_{1}(2\tau_{s}-1)\right]}.$
(36)

As noted previously, the parameter $l$ is related to the distance ratio
$\gamma$. Then, it can be seen from Eqs. (35) and (36) that $r$ and $\tau_{q}$
are dependent on $\gamma$ as well as the gas-wall interaction parameters and
the relaxation time $\tau_{s}$. For microgaseous flows with curved walls, it
is noted that different boundary nodes bring variable values of $\gamma$ along
different lattice directions. It is reasonable for the parameter $r$, which
represents the fraction of the no-slip velocity condition in the combination,
to change locally with boundary nodes. However, for previous curved boundary
schemes which only contains $\gamma$, the relaxation time $\tau_{q}$ therein
is inevitably changeable with $\gamma$ to realize the slip boundary condition.
In the following, we will fix the choice of $l$ to address the issue of
nonuniform relaxation parameters.

From Eq. (23), the relaxation time $\tau_{s}$ is determined by Kn, and should
be constant for a concrete flow problem. Noting that $l$ is a function of
$\gamma$, we denote
$P(\gamma)=\mathcal{L}_{1}(2\tau_{s}-1)\left[2\tau_{s}(6\gamma-1)+1\right]+\sqrt{\frac{6}{\pi}}\tau_{s}\bigl{[}\mathcal{L}_{2}\pi(2\tau_{s}-1)^{2}+12\gamma^{2}+2(2\tau_{s}-1)(6\gamma-6l-1)\bigr{]}$.
In order to achieve the uniform $\tau_{q}$, the numerator of Eq. (36) must be
irrelevant to $\gamma$, which leads to

$\frac{\text{d}P(\gamma)}{\text{d}\gamma}=0\Rightarrow(1-2\tau_{s})\sqrt{\frac{6}{\pi}}\frac{\text{d}l}{\text{d}\gamma}+2\sqrt{\frac{6}{\pi}}\gamma+\bigl{(}\mathcal{L}_{1}+\sqrt{\frac{6}{\pi}}\bigr{)}(2\tau_{s}-1)=0.$
(37)

For the above ordinary differential equation, the solution of $l$ is

$l=\frac{\sqrt{\frac{6}{\pi}}\gamma^{2}+\gamma(2\tau_{s}-1)\bigl{(}\mathcal{L}_{1}+\sqrt{\frac{6}{\pi}}\bigr{)}}{\sqrt{\frac{6}{\pi}}(2\tau_{s}-1)}+\mathcal{E},$
(38)

where $\mathcal{E}$ is a constant independent of $\gamma$. With the solved $l$
substituted into Eqs. (35) and (36), we can obtain the new expressions of $r$
and $\tau_{s}$ as

$\displaystyle r$
$\displaystyle=\frac{\mathcal{L}_{1}\tau_{s}\gamma(2\tau_{s}-1)+\sqrt{\frac{6}{\pi}}\tau_{s}[\gamma^{2}+\gamma(2\tau_{s}-1)+(1+\mathcal{E})(2\tau_{s}-1)]}{\mathcal{L}_{1}(2\tau_{s}-1)[(2+\gamma)\tau_{s}-1]+\sqrt{\frac{6}{\pi}}\tau_{s}[\gamma^{2}+\gamma(2\tau_{s}-1)+(1+\mathcal{E})(2\tau_{s}-1)]},$
(39) $\displaystyle\tau_{q}$
$\displaystyle=\frac{(2\tau_{s}-1)\left(\tau_{s}\sqrt{6\pi}\mathcal{L}_{2}-\mathcal{L}_{1}\right)-2\sqrt{\frac{6}{\pi}}\tau_{s}(1+6\mathcal{E})}{8\sqrt{\frac{6}{\pi}}\tau_{s}+4\mathcal{L}_{1}(2\tau_{s}-1)}.$
(40)

It is clear that the relaxation time $\tau_{q}$ is irrelevant to $\gamma$ now.
In summary, our approach to resolve the uniform relaxation time $\tau_{q}$ is
prescribed as follow: When $l$ is given by Eq. (38), the relaxation time
$\tau_{q}$ can take uniform value as Eq. (40), and simultaneously the slip
boundary condition (33) can be correctly realized by the boundary scheme
(III.1) with $r$ determined by Eq. (39).

We note that the above treatment cannot be accomplished in previous studies
for curved slip walls. For the sake of clarification, we refer to the boundary
scheme (III.1) at $l=0$ temporarily as one example of previous curved boundary
conditions. Correspondingly, the values of $r$ and $\tau_{q}$ expressed in
Eqs. (35) and (36) turn to

$\displaystyle r$
$\displaystyle=\frac{\sqrt{\frac{6}{\pi}}\tau_{s}}{\mathcal{L}_{1}(2\tau_{s}-1)+\sqrt{\frac{6}{\pi}}\tau_{s}},$
(41) $\displaystyle\tau_{q}$
$\displaystyle=\frac{\mathcal{L}_{1}(2\tau_{s}-1)\left[2\tau_{s}(6\gamma-1)+1\right]+\sqrt{\frac{6}{\pi}}\tau_{s}\bigl{[}\mathcal{L}_{2}\pi(2\tau_{s}-1)^{2}+12\gamma^{2}+2(2\tau_{s}-1)(6\gamma-1)\bigr{]}}{4(2\tau_{s}-1)\left[2\sqrt{\frac{6}{\pi}}\tau_{s}+\mathcal{L}_{1}(2\tau_{s}-1)\right]}.$
(42)

This clearly indicates that without the parameter $l$ in previous studies, the
relaxation time $\tau_{q}$ should be changeable with $\gamma$ to realize the
prescribed slip boundary condition, as noted before. For the case of flat
walls, $\tau_{q}$ can be invariant because the distance ratio is fixed for all
boundary nodes. While for the case of curved wall geometries, approximated
measurements for $\gamma$ would not be avoided to obtain the uniform
$\tau_{q}$. However, as will shown later in the numerical examples, this
hinders us to correctly realize the slip boundary condition (33) Silva17 .
Specifically, corresponding to the halfway DBB scheme, Eqs. (41) and (42) at
$\gamma=1/2$ give the parameters of $r$ and $\tau_{q}$ that are required to
realize the prescribed slip boundary condition

$\displaystyle r$
$\displaystyle=\frac{\sqrt{\frac{6}{\pi}}\tau_{s}}{\mathcal{L}_{1}(2\tau_{s}-1)+\sqrt{\frac{6}{\pi}}\tau_{s}},$
(43) $\displaystyle\tau_{q}$
$\displaystyle=\frac{\mathcal{L}_{1}(2\tau_{s}-1)\left(4\tau_{s}+1\right)+\sqrt{\frac{6}{\pi}}\tau_{s}\bigl{[}\mathcal{L}_{2}\pi(2\tau_{s}-1)^{2}+4(2\tau_{s}-1)+3\bigr{]}}{4(2\tau_{s}-1)\left[2\sqrt{\frac{6}{\pi}}\tau_{s}+\mathcal{L}_{1}(2\tau_{s}-1)\right]}.$
(44)

Now some comments on the above derivation results are given in order. First,
if we set $\mathcal{L}_{1}=\mathcal{L}_{2}=0$ in Eq. (33), the slip velocity
$u_{s}$ degenerates to $u_{s}=0$, which means the no-slip boundary condition
at solid walls. As a consequence, the combination parameter $r$ in Eq. (39)
becomes to $r=1$, and the boundary scheme (III.1) reduces to that for the no-
slip boundary condition Zhao19 , as pointed out in Remark 1. Then, the values
of $l$ and $\tau_{q}$ in Eqs. (38) and (40) will follow the same method as
those given in our recent work WangL20  for curved no-slip walls. Second, to
treat curved slip walls with the uniform $\tau_{q}$, previous studies usually
approximate the distance ratio $\gamma$ as $\gamma=\frac{1}{2}$ (i.e., the
halfway boundary scheme) or by an artificial measurement. These inaccurate
values of $\gamma$ degrade the fidelity of curved geometry. Further based on
Eq. (42), the actual $\tau_{q}$ cannot be obtained to realize the slip
boundary condition under coarse grid resolutions. Third, although the above
derivations are based on the planar Poiseuille flow, they may be applicable to
general cases with curved walls if the flows in the near-wall region can be
assumed to have a second-order polynomial profile locally. This assumption has
been adopted in many previous studies, and its reasonability has been
demonstrated for continuum flows and microscale gaseous flows Verhaeghe09 ;
Chai08 ; Chai08 ; Guo07 ; Guo08 ; Tao15 ; Guo11 ; Silva17 ; WangL20 ; He97 ;
Ginzburg03 ; Pan06 . Finally, there are many choices of $l$ to derive
infinitely specific curved boundary conditions from Eq. (III.1). The present
theoretical analysis clearly indicates that in order to correctly realize a
certain slip boundary condition at curved walls, the free parameter $l$
besides $r$ and $\tau_{q}$ must be also chosen carefully to ensure invariable
relaxation parameters. Eqs. (38)-(40) give the theoretical formulae to
determine $l$ and $r$ in the boundary scheme (III.1) together with the unform
$\tau_{q}$.

Within the framework of BGK model, the corresponding $l$, $r$ and $\tau_{q}$
in the above equations can be determined by taking $\tau_{s}=\tau_{q}=\tau$.
After a direct comparison, one can find that the relaxation time $\tau$ from
Eq. (40) remains invariant to $\gamma$, however, $\tau$ from Eq. (42) without
$l$ still changes with $\gamma$ as its MRT counterpart. Based on this fact, we
can conclude that only by seeking more relaxation parameters in the LBE (e.g.,
extending the BGK model to the MRT model), the uniform relaxation time cannot
be accomplished to realize the slip boundary condition at curved walls. To
resolve such problem, one more feasible method based on our analysis is by
adding free parameters to the curved boundary scheme.

Finally, we discuss the parameter range in Eqs. (38)-(40) adapted to the
curved boundary scheme (III.1). To ensure the non-negativity of $l$ in Eq.
(38), the constant $\mathcal{E}$ should satisfy

$\mathcal{E}\geq-\frac{\sqrt{\frac{6}{\pi}}\gamma^{2}+\gamma(2\tau_{s}-1)\bigl{(}\mathcal{L}_{1}+\sqrt{\frac{6}{\pi}}\bigr{)}}{\sqrt{\frac{6}{\pi}}(2\tau_{s}-1)}.$
(45)

With the above requirement together with $\tau_{s}>0.5$, one can derive that
the numerator and denominator in Eq. (39) are both positive, and further the
combination parameter $r$ locates in $0<r<1$. On the other hand, because of
the stability condition that $\tau_{q}>0.5$, another requirement of
$\mathcal{E}$ from Eq. (40) can be obtained as

$\mathcal{E}<\frac{(\pi\tau_{s}\mathcal{L}_{2}-3\sqrt{\frac{\pi}{6}}\mathcal{L}_{1})(2\tau_{s}-1)}{12\tau_{s}}-\frac{1}{2}.$
(46)

For the case that $\gamma$ is close to zero, the parameter $l$ from Eq. (38)
approximates to $\mathcal{E}$ and hence should obey Eq. (46) as well. When the
right-hand-side term is smaller than zero, the parameter $l$ in the
simulations would be negative. In this case, as we have done for continuum
flows WangL20 , to balance the numerical stability and the accurate
implementation of the boundary scheme (III.1), the parameter $l$ is
compulsively set as zero if the actual $l<-0.4$ obtained from Eq. (38) at very
small $\gamma$. As will be shown in the subsequent simulations, this
compromised treatment can bring good predictions with the analytical
solutions.

##  IV Numerical results and discussions

To validate the proposed curved boundary condition [Eq. (III.1)] and the
theoretical derivations, some well-established microscale flows are simulated
in this section. The problems under consideration include the microscale
aligned and inclined plane Poiseuille flow and the microcylindrical Couette
flow. The slip boundary condition at curved walls is modeled by the boundary
scheme (III.1), where $l$ and $r$ are determined from Eqs. (38) and (39). For
comparison with previous studies, some other choices of $l$ as well as the
halfway boundary scheme (corresponds to $\gamma=1/2$) are also tested to
predict the microslip velocity at curved walls. To clearly expose the
difference between them, coarse grid resolutions are used in the numerical
simulations of all considered problems.

In the simulations, the slip coefficients $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$ in Eq. (33) are taken as
$\mathcal{L}_{1}=(2-\sigma)(1-0.1817\sigma)/\sigma$, and
$\mathcal{L}_{2}=\pi^{-1}+\mathcal{L}_{1}^{2}/2$ Guo08P , where $\sigma$ is
the wall accommodation coefficient. The value of $\sigma$ in this work is set
as $\sigma=1$, which means that the wall is fully diffusive. The relaxation
times $\tau_{\rho}$ and $\tau_{j}$ are specified as
$\tau_{\rho}=\tau_{j}=1.0$, and $\tau_{s}$ is determined via the Knudsen
number as Eq. (23). To realize the slip boundary condition by the proposed
boundary scheme (III.1), $\tau_{q}$ is chosen according to Eq. (40) as derived
above. The remaining relaxation times $\tau_{e}$ and $\tau_{\varepsilon}$ are
given by $\tau_{e}=1.1$ and $\tau_{\varepsilon}=1.2$. Actually, the effects of
$\tau_{\rho}$, $\tau_{j}$, $\tau_{e}$ and $\tau_{\varepsilon}$ are negligible
on the numerical results. For the diffusive scaling used in the boundary
scheme (III.1), $\delta_{x}$ and $\delta_{t}$ obey the relation of
$\delta_{t}=\eta\delta_{x}^{2}$, and in the simulations, they are determined
as $\delta_{t}=\eta\delta_{x}^{2},~{}\eta=\frac{\tau_{s}-\frac{1}{2}}{3\nu}$.

###  IV.1 Force-driven microchannel flow

We first consider the Poiseuille flow driven by a constant force in a
microchannel. This microflow, which has an analytical solution, has been
recognized as a benchmark problem in the LBM. In the following, the aligned
and inclined plates with respect to the computational grid are successively
considered in the simulations.

####  IV.1.1 Aligned channel case

The force-driven Poiseuille flow between two parallel plates is first
simulated. As schematically shown in Fig. 3, the flows in the channel with
width $H$ are driven by a constant acceleration $\bm{a}=(a,0)$ along the
$x$-direction.

Figure 3: Schematic of the Poiseuille flow in a microscale channel with
aligned straight walls.

Periodic boundary conditions are executed at the entrance and outlet of the
channel. With the wall slip velocity given by Eq. (34), the dimensionless
velocity has the following analytical solution

$U:=\frac{u(y)}{u_{c}}=4\frac{y}{H}\left(1-\frac{y}{H}\right)+4\mathcal{L}_{1}\text{Kn}+8\mathcal{L}_{2}\text{Kn}^{2},\quad
V:=\frac{v(y)}{u_{c}}=0,$ (47)

where $0\leq y\leq H$, and $u_{c}=aH^{2}/8\nu$ is the maximum streamwise
velocity.

In the simulations, the lower and upper plates are placed with distance
$\gamma\delta_{x}$ away from boundary lattice nodes (as illustrated in Fig.
2). The grid number spanning in the vertical direction is $M$, and this gives
the lattice spacing as $\delta_{x}=H/(M+2\gamma)$. To ensure the low Mach
number for a finite Knudsen number, the driven acceleration $a$ is set to be
$10^{-4}$. We performed some simulations with the boundary scheme (III.1) to
measure the dimensionless slip velocities under different choices of $l$.

---
(a) $\text{Kn}=0.02$ (b) $\text{Kn}=0.1$
Figure 4: Slip velocity against $r$ for the boundary scheme (III.1) with
different $l$ at $M=16$ and [(a) $\text{Kn}=0.02$; (b) $\text{Kn}=0.1$]. As a
specific but representative case, the distance ratio is $\gamma=0.25$, and
$\tau_{q}=\tau_{s}$.

The numerical results against $r$ at two Knudsen numbers
($\text{Kn}=0.02,~{}0.1$) are shown in Fig. 4 where $M=16$, $\nu=0.01$, and
$\tau_{q}=\tau_{s}$. The value of $\gamma=0.25$ is taken as a representative
case for $\gamma$ ranging in $0\leq\gamma\leq 1$ . As clearly seen from the
figure, the numerical predictions agree well with the theoretical derivations
given by Eq. (32). In addition, similar excellent consistency results can be
also obtained for some other values of $\gamma$ and $\tau_{q}$, which confirms
the derivation result for the slip velocity $U_{s}$.

For the case of aligned channel flows, the distance ratios possess the same
value at boundary nodes along different lattice directions. Thus, as noted
previously, the issue of nonuniform relaxation times will not occur for the
halfway DBB and curved boundary schemes to realize the slip boundary
condition. We next investigate the discrete effects only for the proposed
boundary condition (III.1). Simulations with different lattice sizes are
carried out for $\text{Kn}=0.02$ and $\text{Kn}=0.2$. The predicted velocity
profiles are shown in Fig. 5, where $\gamma=0.5$ and $l=\gamma$ are used as a
representative value.

---
Figure 5: Velocity profiles predicted by the present boundary scheme with
different lattice sizes [(a) $\text{Kn}=0.02$; (b) $\text{Kn}=0.2$]. Filled
shapes denote the results obtained by Eqs. (35) and (36). Empty shapes denote
the results obtained under (a) $r=0.5$ and $\tau_{q}$ given by Eq. (36); (b)
$r$ given by Eq. (35) while $\tau_{q}=\tau_{s}$. The analytical velocity
solutions (solid lines) are also included.

It is clearly seen that if $r$ and $\tau_{q}$ is unitedly given by Eqs. (35)
and (36), the numerical predictions agree well with the analytical solution
even with only four grid points. Otherwise, apparent grid-independent results
are observed to deviate from the analytical velocity profiles. Furthermore, as
Kn increases or $\delta_{x}$ decreases, we found that such difference from the
analytical solutions become more pronounced in the simulations. These results
demonstrate that the parameters $r$ and $\tau_{q}$ must be carefully chosen to
realize the desired slip boundary condition, as have been revealed in many
published studies.

####  IV.1.2 Inclined channel case

In contrast to the above aligned microchannel case, the Poiseuille flow in an
inclined channel is more complex and further considered. As shown in Fig. 6,
the flat walls are inclined at an inclination angle $\theta$ with respect to
the underlying grid.

Figure 6: Schematic illustration of the Poiseuille flow in an inclined
microchannel.

From the entrance to the exit, the flow in the inclined channel is periodic
along the direction of driven force. The microgaseous flows in the channel
have the analytical velocity profiles as

$U:=\frac{u(y^{\prime})}{u_{c}}=4\frac{y^{\prime}}{H}\left(1-\frac{y^{\prime}}{H}\right)+4\mathcal{L}_{1}\text{Kn}+8\mathcal{L}_{2}\text{Kn}^{2},$
(48)

where $x^{\prime}$ and $y^{\prime}$ are the coordinates respectively parallel
and perpendicular to the inclined channel wall, $0\leq y^{\prime}\leq H$, and
$u_{c}=aH^{2}/8\nu$.

Denote $N_{x}$ and $N_{AB}$ as the grid number in the $x$-direction and for
the segment of $AB$. It follows that the total grid number in the
$y$-direction is $N_{y}=N_{AB}+N_{x}\text{tan}\theta$. It should be noted that
different from the aligned case, the present skew boundary geometry cannot
bring uniform values of $\gamma$ at boundary nodes any more. The curved
boundary scheme (III.1) with the derived Eqs. (38)-(40) under
$\mathcal{E}=-0.65$ is then employed for the slip boundary condition at the
inclined plates. In the simulations, if the computation from Eq. (38) gives
$l<-0.4$ at a boundary node with very small $\gamma$, the parameter $l$ will
be set to zero as note before. Figure 7 presents the velocity profiles in the
inclined microchannel at $\text{Kn}=0.0194$. Three inclination angles,
$\text{tan}\theta=0.2,~{}1.2$, and $2.0$, are considered under $N_{x}=120$ but
different values of $N_{AB}$.

---
Figure 7: Normalized velocity profiles of the inclined microchannel flow with
$\text{Kn}=0.0194$ for different inclination angles [(a)
$\text{tan}\theta=0.2$; (b) $\text{tan}\theta=1.2$; (c)
$\text{tan}\theta=2.0$]. Filled shapes denote the predicted results by the
boundary scheme (III.1) with Eqs. (38)-(40). Empty shapes represent the
results obtained by the halfway DBB scheme with Eqs. (43) and (44).

Clearly, good agreement with the analytical solution is achieved for the
boundary scheme (III.1) with Eqs. (38)-(40), even with the artificial
treatment of $l$ at very small $\gamma$. For comparisons, the results
predicted by the halfway DBB scheme with Eqs. (43) and (44) are also shown.
Grid-dependent derivations from the analytical velocity profile are clearly
observed. This confirms that notwithstanding the uniform $\tau_{q}$ determined
by Eq. (44), the halfway DBB scheme cannot realize the accurate slip boundary
condition at curved walls theoretically. As noted previously, the numerical
error is induced by the insufficient accuracy of discrete zigzag ghost
boundary to match the real curved wall. In contrast, owing to the local
$\gamma$ handling the actual curved geometry and the free parameter $l$ [Eq.
(38)], the present boundary scheme (III.1) ($r$ given by Eq. (39)) can
excellently capture the analytical solution with a uniform $\tau_{q}$ [Eq.
(40)] even under a low grid resolution.

In Fig. 8, the velocity profiles predicted by the boundary scheme (III.1) with
$l=0$ and $r$ given by Eq. (41) are shown against different $\gamma$ at
$\text{Kn}=0.0194$ and $N_{AB}=16$. For each $\gamma$, the relaxation time
$\tau_{q}$ is computed according to Eq. (42). In this case, this invariable
$\tau_{q}$ is artificially obtained by approximating the same value of
$\gamma$ at all boundary nodes. As can be clearly seen, the simulated velocity
profiles deviate from the analytical solution in all cases. Furthermore, as
also exposed in Fig. 7, such deviations aggravate as the inclination angle
$\theta$ increases. Similar deviation results are also found in the
simulations for other choices of $l$ that dissatisfy Eq. (38).

---
Figure 8: Velocity profiles of the inclined microchannel flow at
$\text{Kn}=0.0194$ and $N_{AB}=16$ for different inclination angles [(a)
$\text{tan}\theta=0.2$; (b) $\text{tan}\theta=1.2$; (c)
$\text{tan}\theta=2.0$]. Filled square shapes denote the predicted results by
the boundary scheme (III.1) with Eqs. (38)-(40). Empty shapes represent the
results obtained by the boundary scheme (III.1) with $l=0$ and $r$ given by
Eq. (43), while $\gamma$ is fixed in Eq. (44) to obtain the uniform
$\tau_{q}$.

These results demonstrate that for the halfway DBB scheme and previous curved
boundary schemes only including $\gamma$, the slip boundary condition at
curved walls cannot be successfully realized with uniform relaxation times. In
contrast, when $l$ is determined by Eq. (38), the numerical outcome can
produce good predictions consistent with the analytical solution.

From the theoretical analysis presented in Sec. III, it is found that the
parameter $l$ in the boundary scheme (III.1) affects the numerical slip
velocity to match the physical one. To investigate the effect from the choice
of $l$, some simulations are carried out with four cases of $l$
($l=\gamma,~{}2\gamma,~{}\gamma^{2},~{}\gamma^{2}+\gamma$) besides the case of
$l$ given by Eq. (38), while $r$ and $\tau_{q}$ are fixedly determined by Eqs.
(39) and (40) with $\mathcal{E}=-0.65$.

---
Figure 9: Velocity profiles of the inclined microchannel for different choices
of $l$ at $\text{tan}\theta=1.2$. (a) $\text{Kn}=0.0194,~{}N_{AB}=4$; (b)
$\text{Kn}=0.1,~{}N_{AB}=6$. The analytical solution is included as the
reference result for comparison.

In Fig. 9, the results predicted by the five cases of $l$ are shown and
compared with the analytical solutions. Clearly, the results with the
parameter $l$ conforming to Eq. (38) exhibit the best agreement with the
analytical solutions. However, clear deviations from the analytical solutions
are observed for the other four cases of $l$. This is in expectation because
the four choices of $l$ cannot generate uniform relaxation times $\tau_{q}$ to
realize the slip boundary condition. These observations demonstrate the
superiority of the present boundary scheme, as mentioned previously, over
previous curved ones in capturing micro flows with curved walls.

###  IV.2 Couette flow between two concentric cylinders

The proposed boundary scheme is further applied to a microgaseous flow with
more complex geometries, i.e., the microcylindrical Couette flow between two
cylinders, to validate the theoretical analysis. This classical problem has
been studied by many researchers as a benchmark case in fluid dynamics Guo11 ;
LiuZ19 ; Tao15 ; Lockerb04 ; Yuhong05 .

Figure 10: Schematic illustration of the microcylindrical Couette flow.

As shown in Fig. 10, the two cylinders of the problem are concentric with
respective radii of $R_{1}$ and $R_{2}$ ($R_{1}<R_{2}$), and rotate at angular
velocities $\omega_{1}$ and $\omega_{2}$, respectively. As the flow between
the cylinders reaches at steady state, they can be described by the following
reduced Navier-Stokes equations in cylinder polar coordinates $(r,\theta)$

$\frac{\text{d}^{2}u_{\theta}}{\text{d}r^{2}}+\frac{\text{d}}{\text{d}r}\left(\frac{u_{\theta}}{r}\right)=0,$
(49)

where $u_{\theta}$ is the tangential velocity and $r$ is the radial distance.
Based on the Maxwell’s diffusive boundary condition for gas-wall interactions,
the slip boundary condition at the inner and outer cylinder walls can be
expressed as

$u_{\theta}\big{|}_{r=R_{1}}=\omega_{1}R_{1}+\frac{2-\sigma_{1}}{\sigma_{1}}\lambda\left(\frac{\text{d}u_{\theta}}{\text{d}r}-\frac{u_{\theta}}{r}\right)\bigg{|}_{r=R_{1}},\quad
u_{\theta}\big{|}_{r=R_{2}}=\omega_{2}R_{2}-\frac{2-\sigma_{2}}{\sigma_{2}}\lambda\left(\frac{\text{d}u_{\theta}}{\text{d}r}-\frac{u_{\theta}}{r}\right)\bigg{|}_{r=R_{2}},$
(50)

where $\sigma_{1}$ and $\sigma_{2}$ are the accommodation coefficients of the
inner and outer cylinder surface, respectively. Under the boundary condition
(50), the solution of Eq. (49) for the velocity profile can be analytical
solved as Yuhong05

$u_{\theta}=\frac{\mathcal{M}\omega_{1}-\mathcal{N}\omega_{2}}{\mathcal{M}-\mathcal{N}}r+\frac{\omega_{1}-\omega_{2}}{\mathcal{N}-\mathcal{M}}\frac{1}{r},$
(51)

where

$\mathcal{M}=\frac{1}{R_{2}^{2}}\left(1-\frac{2-\sigma_{2}}{\sigma_{2}}\frac{2\lambda}{R_{2}}\right),\quad\mathcal{N}=\frac{1}{R_{1}^{2}}\left(1+\frac{2-\sigma_{1}}{\sigma_{1}}\frac{2\lambda}{R_{1}}\right).$
(52)

In the simulations, the two cylinder surfaces are both assumed to be fully
diffusive, i.e., $\sigma_{1}=\sigma_{2}=\sigma=1$. To mimic the slip boundary
condition (50) at the two cylinders’ surfaces, we set
$\mathcal{L}_{1}=(2-\sigma)/2,~{}\mathcal{L}_{2}=0$ in the boundary condition
(III.1) implemented with Eqs. (38)-(40). The Knudsen number for the flow is
defined as $\text{Kn}=\lambda/(R_{2}-R_{1})$. The radius ratio of the two
cylinders, $\beta=R_{1}/R_{2}$, is obtained by changing $R_{1}$ under the
fixed $R_{2}=1.0$. The two cylinder’s center is placed at that of a square
domain covered by $M$ grid cells. Considering the small grid sizes used, we
set $\text{Kn}=0.01$ and $\omega_{1}=\omega_{2}=\omega=0.001$ in the
simulations. As done before, the value of $l$ is taken as zero in case
$l<-0.4$ computed from Eq. (38), in which the constant $\mathcal{E}$ is
assigned in the range of Eq. (46). In Fig. 11, the tangential velocities with
different lattice sizes are shown for two cases of rotating cylinders.

---
Figure 11: Velocity profiles of the microcylindrical Couette flow for (a)
$\beta=0.2,~{}\omega_{1}=0.001,~{}\omega_{2}=0$; (b)
$\beta=0.6,~{}\omega_{1}=0,~{}\omega_{2}=0.001$. Filled shapes denote the
results predicted by the boundary scheme (III.1) with Eqs. (38)-(40)
($\mathcal{E}=-0.68$). Empty shapes are the results obtained by the halfway
DBB scheme with Eqs. (43) and (44).

Clearly, the velocity profiles predicted by the boundary scheme (III.1) with
Eqs. (38)-(40) are in good agreement with the analytical solutions. Whereas
for the halfway DBB scheme, Eqs. (43) and (44) cannot yield grid-independent
results consistent with the analytical solutions. This is ascribed to the fact
that the midway approximation for the cylinder surfaces distort the actual
curved geometries. Our theoretical analyses on curved boundary conditions are
demonstrated again here.

Figure 12 displays the velocity profiles between the cylinders predicted under
different cases of $l$. Four additional choices of $l$, i.e.,
$l=\gamma,~{}2\gamma,~{}\gamma^{2},~{}\gamma^{2}+\gamma$ besides $l$ given by
Eq. (38) are considered in the simulations. In all these cases, the
combination parameter $r$ is determined by Eq. (39), and $\tau_{q}$ by Eq.
(40) with $\mathcal{E}=-0.65$.

---
Figure 12: Velocity profiles of the microcylindrical Couette flow for (a)
$\beta=0.2,~{}\omega_{1}=0.001,~{}\omega_{2}=0$; (b)
$\beta=0.6,~{}\omega_{1}=0,~{}\omega_{2}=0.001$ at $M=16$.

From the figure, it is observed that compared with the other choices of $l$,
the optimal consistency with analytical solutions is obtained for the case
that $l$ conforms to Eq. (38), even with the artificial value $l=0$ at several
boundary points. This also indicates that the uniform relaxation time
$\tau_{q}$ cannot be achieved for the other cases of $l$ to realize the slip
boundary condition. In addition, the tangential velocity predicted with
$l=\gamma$ in Eqs. (35) and (36) are investigated as $\gamma$ is artificially
fixed to derive the uniform $\tau_{q}$. The simulated results are delineated
in Fig. 13, where $\text{Kn}=0.045$ is set to ensure $\tau_{q}>0.5$.

---
Figure 13: Velocity profiles of the microcylindrical Couette flow at
$\text{Kn}=0.0045$ and $M=16$ for (a)
$\beta=0.2,~{}\omega_{1}=0.001,~{}\omega_{2}=0$; (b)
$\beta=0.6,~{}\omega_{1}=0,~{}\omega_{2}=0.001$. Filled square shapes denote
the predicted results by the boundary scheme (III.1) with Eqs. (38)-(40).
Empty shapes represent the results obtained by the boundary scheme (III.1)
with $l=\gamma$ in Eq. (35) for $r$, and in Eq. (36) to obtain the uniform
$\tau_{q}$ by fixed values of $\gamma$.

Again, apparent deviations from the analytical solutions are observed for the
artificial approximation of $\gamma$. While with Eqs. (38)-(40), the boundary
scheme (III.1) brings robust agreement results with the analytical solutions.
These observations strengthen and confirm the capability of the present
boundary scheme to realize the slip boundary condition at curved walls with
uniform relaxation parameters.

##  V Conclusions

In this work, a kinetic boundary condition has been developed for the LBM
simulating microgaseous flows with curved geometries. This curved boundary
scheme is a combination of the Maxwellian diffuse reflection scheme and a
single-node boundary scheme for curved no-slip walls. In addition to the
distance ratio $\gamma$, which is only involved in previous curved boundary
treatments, an additional free parameter $l$ is incorporated in the present
local boundary scheme. Based on the theoretical analysis within the framework
of MRT model, it is shown that the free parameter $l$ as well as the distance
ratio $\gamma$ and the relaxation times ($\tau_{s}$ and $\tau_{q}$) unitedly
affect the slip velocity derived from the boundary scheme. Thanks to the free
parameter $l$, an available strategy to determine the uniform $\tau_{q}$
together with the combination parameter $r$ is proposed to realize the slip
boundary condition at solid walls. Furthermore, it is found that without the
free parameter $l$, the accurate slip boundary condition cannot be realized
with an invariable $\tau_{q}$ for the halfway DBB scheme and previous curved
boundary schemes.

The proposed curved boundary condition is applied to some benchmark problems
with planar and curved walls, including the aligned and inclined microchannel
flows and the microcylindrical Couette flow. To avoid the instability at very
small $\gamma$ in the simulations, the free parameter $l$ is artificially
assigned with zero value at several boundary nodes. Good and robust
predictions from the present derivations are obtained to match the analytical
solutions satisfactorily even under a small lattice size. The numerical
results also show that for microgaseous flows with nonplanar and curved walls,
the halfway DBB scheme and previous curved boundary conditions which only
contains $\gamma$ bring clear grid-dependent deviations from the analytical
solutions.

On the basis of the present study, a noteworthy point is that to ensure
uniform relaxation parameters in realizing a prescribed slip boundary
condition at curved walls, adding free parameters to the boundary scheme would
be more efficient than seeking more relaxation parameters in the LBE. It
should be noted that this work on curved boundary conditions focused on
microscale gas flows in the slip regime. However, the present analysis is also
instructive for our extension work to the transition regime, where the effect
of the Knudsen layer must be incorporated. In addition, we would like to point
out that by replacing the Maxwellian diffusive part by the specular reflection
scheme in Eq. (III.1), it can generate another hybrid but nonlocal boundary
scheme for microscale flows. Following the method presented in this work, a
similar strategy for uniform relaxation parameters can be obtained to realize
the prescribed slip boundary condition at curved walls. Furthermore, such idea
of introducing a free parameter can be also applied to combining the
Maxwellian diffusive boundary scheme and the specular reflection scheme.
Finally, the present work can be extended to the three-dimensional case
without much difficulty. These interesting topics will be left for our future
work.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(No. 51776068 and No. 51906044) and the Fundamental Research Funds for the
Central Universities (No. 2018MS060). L. Wang would like to thank Profs. Wen-
An Yong and Zhaoli Guo and Dr. Weifeng Zhao for their helpful discussions.

## References

  * (1) C. M. Ho and Y.C. Tai, Annu. Rev. Fluid Mech. 30, 579 (1998).
  * (2) G. Karniadakis, A. Beskok, N. Aluru, _Microflows: fundamentals and simulation_ (Springer, NewYork, 2002).
  * (3) W.-M. Zhang, G. Meng, and X. Y. Wei, Microfluid. Nanofluid. 13, 845 (2012).
  * (4) C. Cercignani, _Mathematical Methods in Kinetic Theory_ (Plenum, New York, 1990).
  * (5) X. He, L. S. Luo, Phys. Rev. E 55, R6333 (1997).
  * (6) X. Shan, X. He, Phys. Rev. Lett. 80, 65 (1998).
  * (7) X. Nie, G. D. Doolen, and S. Chen, J. Stat. Phys. 107, 279 (2002).
  * (8) C. Y. Lim, C. Shu, X. D. Niu, and Y. T. Chew, Phys. Fluids 14, 2299 (2002).
  * (9) C. Shu, X. D. Niu, and Y. T. Chew, J. Stat. Phys. 121, 239 (2005).
  * (10) Y. H. Zhang, R. Qin, and D. R. Emerson, Phys. Rev. E 71, 047702 (2005).
  * (11) J. F. Zhang, Microfluid. Nanofluid. 10, 1 (2005).
  * (12) S. H. Kima, H. Pitsch, I. D. Boyd, J. Comput. Phys. 227, 8655 (2008).
  * (13) Y. Sone, _Kinetic Theory and Fluid Dynamics_ (Birkhäuser, Boston, 2002).
  * (14) F. Sharipov, J. Phys. Chem. Ref. Data 40, 023101 (2011).
  * (15) F. Verhaeghe, L. S. Luo, and B. Blanpain, J. Comput. Phys. 228, 147 (2009).
  * (16) C. Y. Lim, C. Shu, X. D. Niu, and Y. T. Chew, Phys. Fluids 14, 2299 (2002).
  * (17) S. Ansumali and I. V. Karlin, Phys. Rev. E 66, 026311 (2002).
  * (18) S. Succi, Phys. Rev. Lett. 89, 064502 (2002).
  * (19) M. Sbragaglia and S. Succi, Phys. Fluids 17, 093602 (2005).
  * (20) G. H. Tang, W. Q. Tao, and Y. L. He, Phys. Fluids 17, 058101 (2005).
  * (21) Z. H. Chai, Z. L. Guo, L. Zheng, and B. C. Shi, J. Appl. Phys. 104, 014902 (2008).
  * (22) Z. Guo, B. Shi, T.S. Zhao, and C. Zheng, Phys. Rev. E 76, 056704 (2007).
  * (23) Z. Guo and C. Zheng, Int. J. Comput. Fluid Dyn. 22, 465 (2008).
  * (24) Z. H. Chai, B. C. Shi, Z. L. Guo, and J. H. Lu, Commun. Comput. Phys. 8, 1052 (2010).
  * (25) L. Szalm$\acute{a}$s, Phys. Rev. E 73, 066710 (2006).
  * (26) K. Suga, Fluid Dyn. Res. 45, 034501 (2013).
  * (27) Z. Liu, Z. Mu, H. Wu, Microfluid. Nanofluid. 23 27 (2019).
  * (28) S. Tao and Z. Guo, Phys. Rev. E 91 043305 (2015).
  * (29) G. Silva and V. Semiao, Phys. Rev. E 96 013311 (2017).
  * (30) G. Silva, Phys. Rev. E 98 023302 (2018).
  * (31) J. J. Ren, Q. Zheng, and Y. L. Li, Physica A 550 124181 (2020).
  * (32) Z. L. Guo, B. C. Shi, and C. G. Zheng, Comput. Math. Appl. 61, 3519 (2011).
  * (33) W. F. Zhao, J. T. Huang, and W.-A. Yong, Multiscale. Model. Simul. 17, 854 (2019).
  * (34) Z. L. Guo and C. Shu, _Lattice Boltzmann Method and its Applications in Engineering_(World Scientific Press, Singapore, 2013).
  * (35) P. Lallemand and L. S. Luo, Phys. Rev. E 61, 6546 (2000).
  * (36) Z. L. Guo, T. S. Zhao, and Y. Shi, J. Appl. Phys. 99, 074903 (2006).
  * (37) L. Wang, S. Tao, X. H. Meng, K. Zhang, and G. Lu, Phys. Rev. E 101, 0633070 (2020).
  * (38) Z. L. Guo, C. G. Zheng, and B. C. Shi, Phys. Rev. E 77, 036707 (2008).
  * (39) N. G. Hadjiconstantinou, Phys. Fluids 18, 111301 (2006).
  * (40) X. He, Q. Zou, L.-S. Luo, and M. Dembo, J. Stat. Phys. 87, 115 (1997).
  * (41) I. Ginzburg and D. d’Humi$\grave{e}$res, Phys. Rev. E 68, 066614 (2003).
  * (42) C. Pan, L.-S. Luo, and C. T. Miller, Comput. Fluids 35 898 (2006).
  * (43) D. A. Lockerby, J. M. Reese, D. R. Emerson, and R. W. Barber, Phys. Rev. E 70, 017303 (2004).
  * (44) S. Yuhong, R. W. Barber, and D. R. Emerson, Phys. Fluids 17, 047102 (2005).