Curved boundary conditions of the lattice Boltzmann method for simulating microgaseous flows in the slip flow regime
aa r X i v : . [ phy s i c s . c o m p - ph ] D ec Curved boundary conditions of the lattice Boltzmann method forsimulating microgaseous flows in the slip regime
Liang Wang,
1, 2
Shi Tao, Junjie Hu, Kai Zhang,
1, 2 and Gui Lu ∗ Beijing Key Laboratory of Emission Surveillanceand 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 Guangdong Provincial Key Laboratory of Distributed Energy Systems,Dongguan University of Technology, Dongguan 523808, China Faculty of Engineering, China University of Geosciences, Wuhan 430074, China Key Laboratory of Power Station Energy TransferConversion and System of Ministry of Education,North China Electric Power University, Beijing 102206, China (Dated: January 1, 2021) bstract The lattice Boltzmann method (LBM) has shown its promising capability in simulating mi-croscale gas flows. However, the suitable boundary condition is still one of the critical issues forthe LBM to model microgaseous flows involving curved geometries. In this paper, a local boundarycondition of the LBM is proposed to treat curved solid walls of microgaseous flows. The developedboundary treatment combines the Maxwellian diffuse reflection scheme and a single-node boundaryscheme which contains a free parameter as well as the distance ratio. The curved boundary condi-tion 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 ratioand relaxation times. By virtue of the free parameter, the combination parameter and the uniformrelaxation time are theoretically determined to realize the accurate slip boundary condition. Inaddition, it is found that besides the halfway diffuse-bounce-back (DBB) scheme, previous curvedboundary schemes only containing the distance ratio cannot ensure uniform relaxation times torealize the slip boundary condition. Some numerical examples with planar and curved boundariesare carried out to validate the present curved boundary scheme. The good and robust consistencyof numerical predictions with analytical solutions demonstrates our theoretical analysis. ∗ Corresponding author: [email protected] . INTRODUCTION With the rapid progress of microscience and nanotechnology, microscale gaseous flowshave attracted increasing research attention in recent years [1–3]. For such gas flows, themean free path of gas moleculars ( λ ) is usually comparable to the characteristic length scaleof 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 longervalid. Generally, the Knudsen number Kn = λ/H is used to measure the departure degreeof microscale gaseous flows from the classic hydrodynamic theory. Based on the magnitudeof Kn, gas flows could be classified into four major regimes: the continuum flow regime withKn ≤ . . < Kn ≤ .
1, the transition flow regime with0 . < Kn ≤
10, and the free-molecular flow regime with Kn >
10. It is widely accepted thatthe Boltzmann equation can work for rarefied gas flows with arbitrary Knudsen numbers[4]. Therefore, as a discrete scheme derived from the Boltzmann equation [5, 6], the latticeBoltzmann method (LBM) has been believed to have the potential for simulating microscalegaseous flows. Since being applied to microscale gaseous flows first in 2002 [7, 16], the LBMhas received increasing interests over the past dozen years [9–12].With the increase of Kn, the Navier-Stokes equations with the no-slip boundary conditionbecome invalid for gas flows, while the fluid slippage on a solid-fluid boundary will arise dueto the gaseous nonequilibrium effect [13, 14]. Therefore, for the LBM to capture the gaseousslip phenomenon in simulating microgaseous flows, the boundary condition is of critical im-portance for effectively predicting gas-solid interactions. This important issue has attractedsubstantial researches toward accurate slip boundary treatments. Nie et al . [7] employedthe 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 [15] actually. Lim et al. [16]investigated the specular reflection (SR) boundary condition in microchannel flow simula-tions, and the slip velocity was not well consistent with some existing analytical solutions.Ansumali and Karlin [17] applied the Maxwellian diffusive (MD) boundary condition forthe Kramer’s problem, but the scheme tends to overpredict the slip velocity. Therefore, asrevealed in the subsequent researches, the above pure boundary schemes cannot accuratelycapture the slip phenomenon. 3he shortcomings in the above scenario have thus stimulated improving hybrid schemeswhich combine the pure boundary conditions with an accommodation coefficient [15, 18–20].By mixing the BB and the SR boundary conditions, Succi [18] proposed a hybrid scheme, i.e.,the bounce-back and specular reflection (BBSR) scheme for simulating microscale flows. Ageneralization of the BBSR scheme was subsequently developed and analyzed by Sbragagliaand Succi [19]. Tang et al. [20] proposed the diffusive and the specular reflection (DSR)scheme, which comes from the combination of the MD and SR boundary schemes. Thanks tothe tunable accommodation parameter, the degree of slip can be freely controlled to recoverdifferent slip models in the two hybrid schemes. Another hybrid boundary scheme is thediffusive and bounce-back (DBB) scheme [15, 21], which is a combination of the MD and theBB boundary conditions. Owing to the advantage of local computation superior to the BBSRand DSR schemes, the DBB scheme may possess more potential in simulating microscalegaseous flows with complex geometries. For hybrid slip boundary conditions, how to choosethe combination parameter is not a convenient task. With the Bhatnagar-Gross-Krook(BGK) model and the multiple-relaxation-time (MRT) model, Guo et al. [22, 23] successivelyanalyzed the discrete effects of the BBSR and DSR boundary conditions, and found theinterrelationship between the two schemes. To realize the accurate slip boundary condition,they proposed a strategy to determine the combination parameter. Verhaeghe et al. [15] andChai et al. [21, 24] further mathematically analyzed the discrete effect of the DBB scheme.Their results illustrated that the combination parameter should be chosen carefully to imposethe accurate slip boundary condition. Noteworthily, these hybrid boundary schemes areoriginally designed to implement the wall location with a definite distance between latticenodes. Szalm´ a s [25] 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 conditionis proposed for microgaseous flows specially with straight walls. Therefore, for microscaleslip flows with curved boundaries, which bring variable wall locations in relation to theunderlying grid, the above-mentioned studies theoretically cannot yield enough accurateresults especially at a small grid resolution.For a more accurate treatment of curved boundaries in microgaseous flows, there havebeen developed a number of boundary conditions considering the actual boundary shapein the literature. Suga [26] introduced an interpolation method into the DBB boundarycondition for microscale flow simulations. The effect of curved surface is represented by4ts intersection distance ratio between neighbouring lattice nodes. But, the combinationcoefficient is still determined by that from the halfway bounce-back case. Based on thenon-equilibrium extrapolation method for curved boundaries and the counter-extrapolationmethod for the velocity/temperature at curved surfaces, Liu et al. [27] proposed a boundarycondition involving the distance ratio for thermal gaseous microflows with curved slip walls.However, as revealed in the literature [28–30], the discrete effects also exist in curved bound-ary treatments for slip walls, which must be minimized to capture correct microgaseous slipphenomenon. Within the MRT model, Tao and Guo [28] incorporated the effect of distanceratio to analyze the DBB scheme, and then developed a boundary scheme to realize the slipboundary condition at curved boundaries. Silva and Semiao [29, 30] introduced the mul-tireflection framework to gaseous slip flows, and put forward curved slip boundary schemeswithin the two-relaxation-time (TRT) model. The theoretical analysis on curved boundaryschemes revealed that the numerical accuracy is related with the combination parameter andthe relaxation times, which are functions of the wall cut-link distance ratio [28–30]. Evenfor microscale binary gaseous flows, such results are also exposed as the DBB scheme isanalyzed to derive the slip velocity at curved surfaces [31]. In the LBM for microflows withcurved walls, previous studies have indicated that the combination parameter can changelocally with the wall cut-link distance ratio. However, the relaxation parameter shouldalso vary with the distance ratio to realize the slip boundary condition, which drives theanisotropic collision operator. Unfortunately, such issue cannot be overcome with previouscurved boundary schemes which only contain the distance ratio.To attain uniform relaxation parameters as noted above, one natural strategy is to imposethe halfway boundary scheme for curved slip walls, and the actual locations between latticenodes are approximated as halfway wall locations [24]. In this way, the uniform relaxationparameter can be then determined by the distance ratio fixed at 0 .
5. Following this lineof thought, there have been some woks reported to treat the curved slip boundary by thehalfway DBB scheme [15, 24, 32]. However, the real curved geometry with the halfwayapproximation will lost its fidelity under coarse grid resolutions, and undesired errors mayarise to contaminate the simulation accuracy [29, 30]. Therefore, based on the above lit-erature review, it promotes us to resolve such a critical issue for microgaseous flows withcurved boundaries: how to retain the relaxation parameters unchanged to realize the slipboundary condition. Accordingly, the present work has the following twofold objectives: to5evelop a curved kinetic boundary condition which involves additional parameters besidesthe distance ratio; and to realize the accurate slip boundary condition while with invariantrelaxation 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 con-ditions, the present scheme is inspired by the idea in Ref. [33] to include a free parameterbesides the distance ratio, which can bring infinitely many curved boundary schemes for mi-crogaseous flows. As the subsequent step, the boundary scheme is mathematically analyzedwithin the MRT model for a unidirectional flow. On the basis of theoretical derivations, thecombination parameter and a strategy to ensure uniform relaxation parameters are givento realize the slip boundary condition at curved slip walls. Numerical simulations are thencarried out to validate the developed method in the cases of aligned and inclined flat wallsand curved walls. As compared with the halfway DBB boundary scheme and other curvedboundary schemes with several specific free parameters, the present method can predict themost 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 bothtime, space and velocity space. The discrete velocity distribution functions evolve accordingto the following lattice Boltzmann equation (LBE), f i ( x + c i δ t , t + δ t ) − f i ( x , t ) = Ω i ( f )( x , t ) + δ t F i ( x , t ) , i = 0 , , · · · , b − , (1)where f i ( x , t ) is the distribution function associated with the discrete velocity c i at position x and time t , δ t is the time step, Ω i ( f ) is the discrete collision operator, and F i is the discreteforcing term, and b is the number of discrete velocities.The BGK or single-relaxation-time model is the most widely used collision operator in theLBM. However, the slip velocity derived within the BGK model depends on the relaxationtime, which is grid resolution dependent for a given Kn [15, 23]. This means that someunphysical numerical artifacts besides the physical part exist in the slip velocity. To avoidthis problem, in this work we turn to employ the MRT collision operator adhered to the6BE, which is written as Ω i ( f ) = − X j ( M − SM ) ij h f j − f (eq) j i , (2)where M is a b × b transformation matrix, which maps f i onto the moment space via m = M f with f = ( f , f , · · · , f b − ) T , S = diag( τ , τ , · · · , τ b − ) − is a diagonal relaxationmatrix with its non-negative element τ i being the relaxation time for the i -th moment. f (eq) j is the equilibrium distribution function which is dependent on the gas density ρ , velocity u and temperature T , f (eq) j = ω j ρ (cid:20) c j · u c s + ( c j · u ) c s − u c s (cid:21) , j = 0 , , · · · , b − , (3)where ω j is the weight coefficient, c s = √ RT ( R is the gas constant) is the lattice soundspeed. For isothermal flows, c s is determined by the lattice speed c = δ x /δ t with δ x beingthe lattice spacing. For the discrete forcing term F i in Eq. (1), it should be taken as [34] F = M − (cid:18) I − S (cid:19) M F , (4)where I is the identity matrix, F = ( F , F , · · · , F b − ) T , and F = ( F , F , · · · , F b − ) T isexpressed as F i = ω i ρ (cid:20) c j · a c s + ua : ( c i c i − c s I ) c s (cid:21) , (5)where G = ρ a is the external force.In this work, we employed the two-dimensional nine-velocity (D2Q9) model, where thediscrete velocities c i are defined by c i := c e i = c (0 , , i = 0 ,c (cid:0) cos (cid:2) ( i − π/ (cid:3) , sin (cid:2) ( i − π/ (cid:3)(cid:1) , i = 1 , , , , √ c (cid:0) cos (cid:2) ( i − π/ π/ (cid:3) , sin (cid:2) ( i − π/ π/ (cid:3)(cid:1) , i = 5 , , , , (6)where c = δ x /δ t with δ x denoting the lattice spacing. Correspondingly, the sound speed c s = √ RT = c/ √
3, and the weight coefficients are given by ω = 4 / ω − = 1 / ω − = 1 /
36. Via the Gram-Schmidt orthogonalization procedure on the discrete velocities c i , there are different versions of transformation matrices M , and one form of M as c = 17s given by [35] M = − − − − − − − − − − − − − − − − − −
10 0 − − −
10 1 − − − − . (7)As a result, the nine discrete velocity moments m i from the distribution functions f i areexpressed as m = M f = ( ρ, e, ε, j x , q x , j y , q y , p xx , p xy ) T . (8)The corresponding relaxation matrix S for the nine moments is written as S = diag( τ ρ , τ e , τ ε , τ j , τ q , τ j , τ q , τ s , τ s ) − . (9)The fluid density ρ and velocity u = ( u, v ) are respectively defined as the zeroth andfirst-order moments of f i ρ = X i f i , ρ u = X i c i f i + δ t ρ a . (10)Through the Chapman-Enskog or linear analysis, the Navier-Stokes equations can be derivedfrom the MRT-LBE model (1). The fluid pressure p is determined by p = c s ρ = ρRT , andthe shear and bulk viscosities are respectively given by ν = c s (cid:18) τ s − (cid:19) δ t , ζ = c s (cid:18) τ e − (cid:19) δ t . (11)The numerical implementation of Eq. (1) is divided into two steps:Collision : f ′ i ( x , t ) = f i ( x , t ) − X j ( M − SM ) ij h f j ( x , t ) − f (eq) j ( x , t ) i + δ t F i ( x , t ) , Streaming : f i ( x + c i δ t , t + δ t ) = f ′ i ( x , t ) , (12)8here f ′ i ( x , t ) is the post-collision distribution function. In the framework of MRT model,the collision step is usually implemented in the moment space m ′ ( x , t ) = m ( x , t ) − S (cid:2) m ( x , t ) − m (eq) ( x , t ) (cid:3) + δ t ˆ F ( x , t ) , f ′ ( x , t ) = M − m ′ ( x , t ) , (13)where m ′ := M f ′ with f ′ = ( f ′ , f ′ , · · · , f ′ ) T is the post-collision moment, m (eq) := M f (eq) and ˆ F := M F = (cid:0) I − S (cid:1) M F denote the equilibria and forcing term in the moment space,respectively. After the collision step is completed, the streaming step is then executed bytransforming the moments back to the velocity space.In contrast to continuum flows, the LBM for microscale flows should consider the rar-efaction effects, which are commonly characterized by the dimensionless Knudsen number.The most important characteristic parameter in continuum flows is the Reynolds numberRe, which determines the relaxation time τ s in the simulations. However, for microscalegas flows, the relationship between τ s and Kn should be carefully established to ensure theconsistency criterion [36], which is read as τ s = 12 + r π Kn Hδ x . (14)Another important issue for the LBM applied to microscale flows is the kinetic boundarycondition. As reviewed previously, there have been some boundary conditions proposedfor the LBE to simulate microgaseous flows, such as the BBSR scheme, DSR scheme andDBB scheme. These kinetic boundary conditions are originally designed for the case of flatwalls with definite locations between lattice nodes. Furthermore, it has been exposed in theliterature [22, 23] that some discrete effects exist in the boundary schemes and should becorrected to realize the slip boundary condition. For a more accurate treatment of curvedwalls, it is common to include the distance ratio in the boundary condition to preserve theactual shape of curved geometries [26, 27]. However, when these curved boundary schemesare implemented to realize the slip boundary condition at curved walls, the relaxation timein the LBE will depend on the distance ratio and thus varies with link directions at differentboundary nodes [28–30]. Noteworthily, such a problem cannot be overcome in previouscurved boundary schemes only containing the distance ratio.9 l x f x b x r γ h | e i | e ¯ i e i lh | e i | x x FIG. 1. Schematic of a curved-wall boundary along one single lattice direction. The thin solid lineis the grid line, and the thick curved one represents the boundary surface. White circles ( ◦ ): thefluid nodes; Black circle ( • ): the intersection point of the boundary with the grid line; Square box( (cid:3) ): the solid node outside the computational domain. III. KINETIC BOUNDARY CONDITION FOR CURVED SLIP WALLS
To remedy the above issue of nonuniform relaxation parameters, one natural and promis-ing way is to introduce additional parameters besides the distance ratio into a kinetic bound-ary condition. As a matter of fact, by resorting to a single-node boundary scheme [33] for afree parameter besides the distance ratio, we recently overcome the numerical slip at curvedno-slip walls [37] with invariant relaxation parameters. Thus, it motivates and inspires usto construct such a parameterized boundary condition for microflows with curved walls, andthen generalize the strategy to ensure uniform relaxation parameters from continuum flowsto microgaseous flows.
A. Curved boundary scheme for microscale gaseous flows
We now propose a kinetic boundary condition that involves an additional parameterbesides the distance ratio as noted above. Because a curved boundary in simulations canbe modeled separately for each lattice direction, we consider a single direction e i pointingfrom the boundary surface at x b to the boundary node x f , as shown in Fig. 1. The surfacepoint x b is intersected by x f and the solid node x r , and then the actual location of curved10oundary can be depicted by the distance ratio γ as γ = | x b − x f | / | x r − x f | . At theboundary node x f , the unknown distribution function f i ( x f , t + δ t ) in the LBM is specifiedaccording to the kinetic boundary condition, which aims to realize the slip velocity conditionat physical walls. Noteworthily, for microscale flows with curved boundaries or in complexgeological porous media, a local boundary condition that involves the current informationat boundary nodes is desired especially. In the present work, our construction aims to alocal curved boundary condition for microscale gas flows.Recall that the DBB boundary condition, which has been developed for both planarand curved slip walls [15, 21, 28], pertains to a local boundary scheme. This is due toits combination of the BB and the MD schemes, which is formulated for the unknowndistribution function as follows: f i ( x f , t + δ t ) = r (cid:20) f ′ ¯ i ( x f , t ) + 2 ω i ρ f c i · u b c s (cid:21) + (1 − r ) f (eq) i ( ρ f , u b ) , (15)where the subscript ¯ i indicates e ¯ i = − e i , ρ f = ρ ( x f , t ) and u b = u ( x b , t ) is the wall velocity.The combination parameter r represents the bounce-back fraction and ranges in the regionof 0 ≤ r ≤
1. Only the lattice directions of e i and e ¯ i are involved at the current node 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 otherhand, Zhao et al . [33] recently proposed a single-node boundary scheme, which contains afree parameter besides γ , for the no-slip boundary condition. Later, its ability to overcomethe discrete effect with uniform parameters has been revealed in our recent work [37]. Alongwith the underlying structure of the DBB scheme, we thus replace the bounce-back part inEq. (15) by the scheme in Ref. [33], and derive the following boundary condition f i ( x f , t + δ t ) = r (cid:20) l − γ l f ¯ i ( x f , t ) + l l f ′ i ( x f , t ) + 2 γ − l l f ′ ¯ i ( x f , t ) + 21 + l ω i ρ f c i · u b c s (cid:21) + (1 − r ) f (eq) i ( ρ f , u b ) , (16)where l is a free parameter besides γ such that x b = ( x + x ) / γ . Since the no-slip boundary scheme in Ref. [33]is incorporated, the present boundary condition is inherently associated with the diffusivescaling δ t = ηh , h = δ x ( η is an adjustable parameter). Because l takes values along11he same lattice direction as γ does (cf Fig. 1), it can be considered as a function of γ .In addition, the parameter l ranges in max { , γ − } ≤ l ≤ γ to ensure the convexcombination of distribution functions for the no-slip part.In the LBM for microscale flows, previous curved boundary schemes only contain thedistance ratio γ , while the present boundary scheme (16) introduces an adjustable parameter l besides γ . With the variable parameter l , the boundary scheme (16) can hence bringnumerous boundary conditions for microgaseous flows. More importantly, we will showlater that owing to the free parameter l , uniform relaxation parameters can be fulfilled torealize an exact prescribed slip boundary condition. Some remarks about the boundarycurved scheme are given as follows: Remark 1.
When r = 1 , the boundary scheme (16) will degenerate to the single-nodescheme proposed in Ref. [33] for continuum flows with no-slip walls. From this viewpoint,the present scheme is the generalized version of that in Ref. [33]. Remark 2.
When l = 0 and γ = , the boundary scheme (16) then degenerates to the DBBscheme (15) . This also confirms the intrinsic halfway consumption of the bounce-back partadhered to the boundary scheme (15) . Remark 3.
In the present work, the effect of curved wall geometry is directly embodied inthe boundary condition (16) by the distance ratio γ . While in Ref. [28], it is absorbed in theDBB scheme (15) through the combination parameter to realize the slip boundary conditionat curved walls. As noted before, the combination parameter r in hybrid boundary schemes plays thevital role and directly affects the degree of slip at physical boundaries [22, 23]. Therefore,it should be carefully chosen to derive reasonable simulation results. Previous studies withthe MRT model have shown that r is related with sever factors, including the relaxationtime, the gas-solid interaction parameter and the distance ratio [23, 28]. For the presentboundary scheme (16) with the free parameter l , the influence of r is expected to be morecomplicated in treating curved slip walls. Next, we will conduct a theoretical analysis toinvestigate how to specify r here, and simultaneously ensure the above-mentioned uniformrelaxation parameters. 12 all FIG. 2. Schematic of the flow and lattice arrangement with a distance ratio γ . The wall boundarycorresponding to the halfway DBB boundary condition is placed with γ = 1 / B. Analysis of the curved boundary scheme
To simplify the analysis, we consider the steady incompressible Poiseuille flow in a mi-crochannel, which has been extensively employed for theoretical analysis in previous stud-ies [9, 15, 21, 23, 28, 32]. As sketched in Fig. 2, the flow is driven by a constant force ρ a = ρ ( a,
0) along the x direction, and is subjected to the following assumptions: ρ = const , v = 0 , ∂ x φ = 0 , ∂ t φ = 0 , (17)where φ is an arbitrary flow variable. The lattice node at the layer j = 0 is located withan arbitrary distance γδ x away from the solid wall. Here j is the index of the grid line at y j = ( j + γ ) δ x . After the streaming step, the unknown distribution functions, f , f and f at j = 0 are determined according to the boundary scheme (16): f = r (cid:20) l − γ l f + l l f ′ + 2 γ − l l f ′ + 21 + l ω ρ c · u b c s (cid:21) + (1 − r ) f (eq)2 ( ρ, u w ) , (18a) f = r (cid:20) l − γ l f + l l f ′ + 2 γ − l l f ′ + 21 + l ω ρ c · u b c s (cid:21) + (1 − r ) f (eq)5 ( ρ, u w ) , (18b) f = r (cid:20) l − γ l f + l l f ′ + 2 γ − l l f ′ + 21 + l ω ρ c · u b c s (cid:21) + (1 − r ) f (eq)6 ( ρ, u w ) , (18c)where u w is the wall velocity, and f i = f i ( y ) and f ′ i = f ′ i ( y ) with y = γδ x .13ollowing the procedures exhibited in Refs. [23, 37, 38], we can obtain the relationshipbetween the velocities u and u respectively at j = 0 and j = 1: u = A u + B aδ t + C u w , (19)where A = (1 + l )(1 − r ) τ s + 2 r (1 + γ )(1 + l )(1 − r ) τ s + 2 rγ , C = − r (1 + l )(1 − r ) τ s + 2 rγ , (20a) B = (1 + r + l − lr )(1 + 4 τ q − τ s − τ q τ s ) + 12 r [(2 γ − l )(2 τ s −
1) + (1 − γ ) τ s ]2 [(1 + l )(1 − r ) τ s + 2 rγ ] (2 τ s − . (20b)For the Poiseuille flow between two plates located at y = 0 and y = H , the analyticalsolution can be expressed as u j = 4 u c y j H (cid:16) − y j H (cid:17) + u w + u s , (21)where u c = aH / ν , and u s is the slip velocity at the physical wall. Substituting Eq. (21)into Eq. (19) to replace u and u , we can derive the dimensionless slip velocity U s := u s u c = 2(1 + l )(1 − r ) r τ s δ x H − r [(1 + l )(1 + 4 τ q − τ s + 12 γτ s − τ s τ q ) + r D ] δ x H . (22)where D = 12 γ + 12 γ ( τ s −
1) + (2 τ s − lτ q − − τ q ) − l [11(2 τ s −
1) + 12 γτ s ]. Furtherinvoking the τ s -Kn relation given by Eq. (14), Equation (22) can be rewritten in terms ofKn as U s = 4 (1 − r )(1 + l ) τ s r (2 τ s − r π Kn + 8 π K r (2 τ s − Kn (23)where K = (2 τ s − r + l +11 lr )+12 γ [ r + τ s ( lr − l − r − − rγ − τ q (2 τ s − lr − l − r − U s includes an additional free parameter l besides γ . Owing to the degree offreedom from l , we now could obtain invariant relaxation parameters to realize the adoptedslip boundary condition by the boundary scheme (16).For non-continuum gaseous flows, a second-order slip boundary condition is widely usedand read as [26, 39] u s = L λ ∂u∂n − L λ ∂ u∂n , (24)where n is the unit outer vector normal to the wall, and L and L refer to the slippagecoefficients related with the gas-wall interactions. Under the slip boundary condition (24),14he Poiseuille flow exhibits the following slip velocity at the wall U s = 4 L Kn + 8 L Kn . (25)Therefore, to realize the slip boundary condition (25) with the boundary scheme (16), r and τ q must be determined by comparing Eq. (23) with Eq. (25), which are given by r = q π τ s (1 + l ) L (2 τ s −
1) + q π τ s (1 + l ) , (26) τ q = L (2 τ s −
1) [2 τ s (6 γ −
1) + 1] + q π τ s (cid:2) L π (2 τ s − + 12 γ + 2(2 τ s − γ − l − (cid:3) τ s − h q π τ s + L (2 τ s − i . (27)As noted previously, the parameter l is related to the distance ratio γ . Then, it canbe seen from Eqs. (26) and (27) that r and τ q are dependent on γ as well as the gas-wallinteraction parameters and the relaxation time τ s . For microgaseous flows with curved walls,it is noted that different boundary nodes bring variable values of γ along different latticedirections. It is reasonable for the parameter r , which represents the fraction of the no-slipvelocity condition in the combination, to change locally with boundary nodes. However, forprevious curved boundary schemes which only contains γ , the relaxation time τ q therein isinevitably changeable with γ to realize the slip boundary condition. In the following, wewill fix the choice of l to address the issue of nonuniform relaxation parameters.From Eq. (14), the relaxation time τ s is determined by Kn, and should be constant fora concrete flow problem. Noting that l is a function of γ , we denote P ( γ ) = L (2 τ s −
1) [2 τ s (6 γ −
1) + 1] + q π τ s (cid:2) L π (2 τ s − + 12 γ + 2(2 τ s − γ − l − (cid:3) . In order toachieve the uniform τ q , the numerator of Eq. (27) must be irrelevant to γ , which leads tod P ( γ )d γ = 0 ⇒ (1 − τ s ) r π d l d γ + 2 r π γ + (cid:0) L + r π (cid:1) (2 τ s −
1) = 0 . (28)For the above ordinary differential equation, the solution of l is l = q π γ + γ (2 τ s − (cid:0) L + q π (cid:1)q π (2 τ s −
1) + E , (29)where E is a constant independent of γ . With the solved l substituted into Eqs. (26) and1527), we can obtain the new expressions of r and τ s as r = L τ s γ (2 τ s −
1) + q π τ s [ γ + γ (2 τ s −
1) + (1 + E )(2 τ s − L (2 τ s − γ ) τ s −
1] + q π τ s [ γ + γ (2 τ s −
1) + (1 + E )(2 τ s − , (30) τ q = (2 τ s − (cid:0) τ s √ π L − L (cid:1) − q π τ s (1 + 6 E )8 q π τ s + 4 L (2 τ s − . (31)It is clear that the relaxation time τ q is irrelevant to γ now. In summary, our approachto resolve the uniform relaxation time τ q is prescribed as follow: When l is given by Eq.(29), the relaxation time τ q can take uniform value as Eq. (31), and simultaneously theslip boundary condition (24) can be correctly realized by the boundary scheme (16) with r determined by Eq. (30).We note that the above treatment cannot be accomplished in previous studies for curvedslip walls. For the sake of clarification, we refer to the boundary scheme (16) at l = 0temporarily as one example of previous curved boundary conditions. Correspondingly, thevalues of r and τ q expressed in Eqs. (26) and (27) turn to r = q π τ s L (2 τ s −
1) + q π τ s , (32) τ q = L (2 τ s −
1) [2 τ s (6 γ −
1) + 1] + q π τ s (cid:2) L π (2 τ s − + 12 γ + 2(2 τ s − γ − (cid:3) τ s − h q π τ s + L (2 τ s − i . (33)This clearly indicates that without the parameter l in previous studies, the relaxation time τ q should be changeable with γ to realize the prescribed slip boundary condition, as notedbefore. For the case of flat walls, τ q can be invariant because the distance ratio is fixed for allboundary nodes. While for the case of curved wall geometries, approximated measurementsfor γ would not be avoided to obtain the uniform τ q . However, as will shown later in thenumerical examples, this hinders us to correctly realize the slip boundary condition (24)[29]. Specifically, corresponding to the halfway DBB scheme, Eqs. (32) and (33) at γ = 1 / r and τ q that are required to realize the prescribed slip boundary16ondition r = q π τ s L (2 τ s −
1) + q π τ s , (34) τ q = L (2 τ s −
1) (4 τ s + 1) + q π τ s (cid:2) L π (2 τ s − + 4(2 τ s −
1) + 3 (cid:3) τ s − h q π τ s + L (2 τ s − i . (35)Now some comments on the above derivation results are given in order. First, if we set L = L = 0 in Eq. (24), the slip velocity u s degenerates to u s = 0, which means theno-slip boundary condition at solid walls. As a consequence, the combination parameter r in Eq. (30) becomes to r = 1, and the boundary scheme (16) reduces to that for the no-slipboundary condition [33], as pointed out in Remark 1. Then, the values of l and τ q in Eqs.(29) and (31) will follow the same method as those given in our recent work [37] for curvedno-slip walls. Second, to treat curved slip walls with the uniform τ q , previous studies usuallyapproximate the distance ratio γ as γ = (i.e., the halfway boundary scheme) or by anartificial measurement. These inaccurate values of γ degrade the fidelity of curved geometry.Further based on Eq. (33), the actual τ q cannot be obtained to realize the slip boundarycondition under coarse grid resolutions. Third, although the above derivations are based onthe planar Poiseuille flow, they may be applicable to general cases with curved walls if theflows 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 beendemonstrated for continuum flows and microscale gaseous flows [15, 21–23, 28, 29, 32, 37,40–42]. Finally, there are many choices of l to derive infinitely specific curved boundaryconditions from Eq. (16). The present theoretical analysis clearly indicates that in orderto correctly realize a certain slip boundary condition at curved walls, the free parameter l besides r and τ q must be also chosen carefully to ensure invariable relaxation parameters.Eqs. (29)-(31) give the theoretical formulae to determine l and r in the boundary scheme(16) together with the unform τ q .Within the framework of BGK model, the corresponding l , r and τ q in the above equationscan be determined by taking τ s = τ q = τ . After a direct comparison, one can find that therelaxation time τ from Eq. (31) remains invariant to γ , however, τ from Eq. (33) without l still changes with γ as its MRT counterpart. Based on this fact, we can conclude thatonly by seeking more relaxation parameters in the LBE (e.g., extending the BGK model to17he MRT model), the uniform relaxation time cannot be accomplished to realize the slipboundary condition at curved walls. To resolve such problem, one more feasible methodbased on our analysis is by adding free parameters to the curved boundary scheme.Finally, we discuss the parameter range in Eqs. (29)-(31) adapted to the curved boundaryscheme (16). To ensure the non-negativity of l in Eq. (29), the constant E should satisfy E ≥ − q π γ + γ (2 τ s − (cid:0) L + q π (cid:1)q π (2 τ s − . (36)With the above requirement together with τ s > .
5, one can derive that the numerator anddenominator in Eq. (30) are both positive, and further the combination parameter r locatesin 0 < r <
1. On the other hand, because of the stability condition that τ q > .
5, anotherrequirement of E from Eq. (31) can be obtained as E < ( πτ s L − p π L )(2 τ s − τ s − . (37)For the case that γ is close to zero, the parameter l from Eq. (29) approximates to E andhence should obey Eq. (37) 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 forcontinuum flows [37], to balance the numerical stability and the accurate implementation ofthe boundary scheme (16), the parameter l is compulsively set as zero if the actual l < − . γ . 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. (16)] and the theoretical deriva-tions, some well-established microscale flows are simulated in this section. The problemsunder consideration include the microscale aligned and inclined plane Poiseuille flow andthe microcylindrical Couette flow. The slip boundary condition at curved walls is modeledby the boundary scheme (16), where l and r are determined from Eqs. (29) and (30). Forcomparison with previous studies, some other choices of l as well as the halfway boundaryscheme (corresponds to γ = 1 /
2) are also tested to predict the microslip velocity at curvedwalls. To clearly expose the difference between them, coarse grid resolutions are used in thenumerical simulations of all considered problems.18
IG. 3. Schematic of the Poiseuille flow in a microscale channel with aligned straight walls.
In the simulations, the slip coefficients L and L in Eq. (24) are taken as L = (2 − σ )(1 − . σ ) /σ , and L = π − + L / σ is the wall accommodation coefficient.The value of σ in this work is set as σ = 1, which means that the wall is fully diffusive.The relaxation times τ ρ and τ j are specified as τ ρ = τ j = 1 .
0, and τ s is determined viathe Knudsen number as Eq. (14). To realize the slip boundary condition by the proposedboundary scheme (16), τ q is chosen according to Eq. (31) as derived above. The remainingrelaxation times τ e and τ ε are given by τ e = 1 . τ ε = 1 .
2. Actually, the effects of τ ρ , τ j , τ e and τ ε are negligible on the numerical results. For the diffusive scaling used in theboundary scheme (16), δ x and δ t obey the relation of δ t = ηδ x , and in the simulations, theyare determined as δ t = ηδ x , η = τ s − ν . A. Force-driven microchannel flow
We first consider the Poiseuille flow driven by a constant force in a microchannel. Thismicroflow, which has an analytical solution, has been recognized as a benchmark problem inthe LBM. In the following, the aligned and inclined plates with respect to the computationalgrid are successively considered in the simulations.
1. Aligned channel case
The force-driven Poiseuille flow between two parallel plates is first simulated. As schemat-ically shown in Fig. 3, the flows in the channel with width H are driven by a constantacceleration a = ( a,
0) along the x -direction. Periodic boundary conditions are executed at19 r U s Analytical, Eq.(23)l=0l=l=2l= r U s Analytical, Eq. (23)l=0l=l=2l= (a) Kn = 0 .
02 (b) Kn = 0 . r for the boundary scheme (16) with different l at M = 16 and [(a)Kn = 0 .
02; (b) Kn = 0 . γ = 0 .
25, and τ q = τ s . the entrance and outlet of the channel. With the wall slip velocity given by Eq. (25), thedimensionless velocity has the following analytical solution U := u ( y ) u c = 4 yH (cid:16) − yH (cid:17) + 4 L Kn + 8 L Kn , V := v ( y ) u c = 0 , (38)where 0 ≤ y ≤ H , and u c = aH / ν is the maximum streamwise velocity.In the simulations, the lower and upper plates are placed with distance γδ x away fromboundary lattice nodes (as illustrated in Fig. 2). The grid number spanning in the verticaldirection is M , and this gives the lattice spacing as δ x = H/ ( M + 2 γ ). To ensure the lowMach number for a finite Knudsen number, the driven acceleration a is set to be 10 − . Weperformed some simulations with the boundary scheme (16) to measure the dimensionlessslip velocities under different choices of l . The numerical results against r at two Knudsennumbers (Kn = 0 . , .
1) are shown in Fig. 4 where M = 16, ν = 0 .
01, and τ q = τ s . Thevalue of γ = 0 .
25 is taken as a representative case for γ ranging in 0 ≤ γ ≤ γ and τ 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 bound-ary nodes along different lattice directions. Thus, as noted previously, the issue of nonuni-20 AnalyticalEqs.(26)-(27) , M=4Eqs.(26)-(27) , M=8Eqs.(26)-(27) , M=16r=0.5, Eq.(27), M=4r=0.5, Eq.(27), M=8r=0.5, Eq.(27), M=16 q = s , M=4Eq. (26), q = s , M=8Eq. (26), q = s , M=16 FIG. 5. Velocity profiles predicted by the present boundary scheme with different lattice sizes [(a)Kn = 0 .
02; (b) Kn = 0 . r = 0 . τ q given by Eq. (27); (b) r given by Eq.(26) while τ q = τ s . The analytical velocity solutions (solid lines) are also included. form relaxation times will not occur for the halfway DBB and curved boundary schemesto realize the slip boundary condition. We next investigate the discrete effects only for theproposed boundary condition (16). Simulations with different lattice sizes are carried outfor Kn = 0 .
02 and Kn = 0 .
2. The predicted velocity profiles are shown in Fig. 5, where γ = 0 . l = γ are used as a representative value. It is clearly seen that if r and τ q isunitedly given by Eqs. (26) and (27), the numerical predictions agree well with the analyt-ical solution even with only four grid points. Otherwise, apparent grid-independent resultsare observed to deviate from the analytical velocity profiles. Furthermore, as Kn increasesor δ x decreases, we found that such difference from the analytical solutions become morepronounced in the simulations. These results demonstrate that the parameters r and τ q mustbe carefully chosen to realize the desired slip boundary condition, as have been revealed inmany published studies.
2. Inclined channel case
In contrast to the above aligned microchannel case, the Poiseuille flow in an inclinedchannel is more complex and further considered. As shown in Fig. 6, the flat walls areinclined at an inclination angle θ with respect to the underlying grid. From the entrance to21 B FIG. 6. Schematic illustration of the Poiseuille flow in an inclined microchannel. the exit, the flow in the inclined channel is periodic along the direction of driven force. Themicrogaseous flows in the channel have the analytical velocity profiles as U := u ( y ′ ) u c = 4 y ′ H (cid:18) − y ′ H (cid:19) + 4 L Kn + 8 L Kn , (39)where x ′ and y ′ are the coordinates respectively parallel and perpendicular to the inclinedchannel wall, 0 ≤ y ′ ≤ H , and u c = aH / ν .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 tan θ . It shouldbe noted that different from the aligned case, the present skew boundary geometry cannotbring uniform values of γ at boundary nodes any more. The curved boundary scheme (16)with the derived Eqs. (29)-(31) under E = − .
65 is then employed for the slip boundarycondition at the inclined plates. In the simulations, if the computation from Eq. (29) gives l < − . γ , the parameter l will be set to zero as notebefore. Figure 7 presents the velocity profiles in the inclined microchannel at Kn = 0 . θ = 0 . , .
2, and 2 .
0, are considered under N x = 120 butdifferent values of N AB . Clearly, good agreement with the analytical solution is achievedfor the boundary scheme (16) with Eqs. (29)-(31), even with the artificial treatment of l at very small γ . For comparisons, the results predicted by the halfway DBB scheme withEqs. (34) and (35) are also shown. Grid-dependent derivations from the analytical velocityprofile are clearly observed. This confirms that notwithstanding the uniform τ q determined22 AnalyticalEqs. (29)-(31), N AB =4Eqs. (29)-(31), N AB =8Eqs. (29)-(31), N AB =16Half-way DBB, N AB =4Half-way DBB, N AB =8Half-way DBB, N AB =16 AnalyticalEqs. (29)-(31), N AB =4Eqs. (29)-(31), N AB =8Eqs. (29)-(31), N AB =16Half-way DBB, N AB =4Half-way DBB, N AB =8Half-way DBB, N AB =16 AnalyticalEqs. (29)-(31), N AB =4Eqs. (29)-(31), N AB =8Eqs. (29)-(31), N AB =16Half-way DBB, N AB =4Half-way DBB, N AB =8Half-way DBB, N AB =16 FIG. 7. Normalized velocity profiles of the inclined microchannel flow with Kn = 0 . θ = 0 .
2; (b) tan θ = 1 .
2; (c) tan θ = 2 . by Eq. (35), the halfway DBB scheme cannot realize the accurate slip boundary conditionat curved walls theoretically. As noted previously, the numerical error is induced by theinsufficient accuracy of discrete zigzag ghost boundary to match the real curved wall. Incontrast, owing to the local γ handling the actual curved geometry and the free parameter l [Eq. (29)], the present boundary scheme (16) ( r given by Eq. (30)) can excellently capturethe analytical solution with a uniform τ q [Eq. (31)] even under a low grid resolution.In Fig. 8, the velocity profiles predicted by the boundary scheme (16) with l = 0 and r given by Eq. (32) are shown against different γ at Kn = 0 . N AB = 16. For each γ ,the relaxation time τ q is computed according to Eq. (33). In this case, this invariable τ q isartificially obtained by approximating the same value of γ at all boundary nodes. As can beclearly 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 θ increases. Similar deviation results are also found in the simulations for other choices of l that dissatisfy Eq. (29). These results demonstrate that for the halfway DBB scheme andprevious curved boundary schemes only including γ , the slip boundary condition at curvedwalls cannot be successfully realized with uniform relaxation times. In contrast, when l isdetermined by Eq. (29), the numerical outcome can produce good predictions consistentwith the analytical solution.From the theoretical analysis presented in Sec. III, it is found that the parameter l inthe boundary scheme (16) affects the numerical slip velocity to match the physical one. To23 AnalyticalEqs. (29)-(31)l=0, Eq. (32), Eq.(33), =0.3l=0, Eq. (32), Eq.(33), =0.5l=0, Eq. (32), Eq.(33), =0.7l=0, Eq. (32), Eq.(33), =0.9
AnalyticalEqs. (29)-(31)l=0, Eq. (32), Eq.(33), =0.3l=0, Eq. (32), Eq.(33), =0.5l=0, Eq. (32), Eq.(33), =0.7l=0, Eq. (32), Eq.(33), =0.9
AnalyticalEqs. (29)-(31)l=0, Eq.(32), Eq. (33), =0.3l=0, Eq.(32), Eq. (33), =0.5l=0, Eq.(32), Eq. (33), =0.7l=0, Eq.(32), Eq. (33), =0.9
FIG. 8. Velocity profiles of the inclined microchannel flow at Kn = 0 . N AB = 16 fordifferent inclination angles [(a) tan θ = 0 .
2; (b) tan θ = 1 .
2; (c) tan θ = 2 . l = 0 and r given by Eq. (34),while γ is fixed in Eq. (35) to obtain the uniform τ q . AnalyticalEq. (29)l=l=2l= l= + AnalyticalEq. (29)l=l=2l= l= + FIG. 9. Velocity profiles of the inclined microchannel for different choices of l at tan θ = 1 .
2. (a)Kn = 0 . , N AB = 4; (b) Kn = 0 . , N AB = 6. The analytical solution is included as thereference result for comparison. investigate the effect from the choice of l , some simulations are carried out with four casesof l ( l = γ, γ, γ , γ + γ ) besides the case of l given by Eq. (29), while r and τ q are fixedlydetermined by Eqs. (30) and (31) with E = − .
65. In Fig. 9, the results predicted by thefive cases of l are shown and compared with the analytical solutions. Clearly, the resultswith the parameter l conforming to Eq. (29) exhibit the best agreement with the analyticalsolutions. However, clear deviations from the analytical solutions are observed for the otherfour cases of l . This is in expectation because the four choices of l cannot generate uniform24 ω R R FIG. 10. Schematic illustration of the microcylindrical Couette flow. relaxation times τ q to realize the slip boundary condition. These observations demonstratethe superiority of the present boundary scheme, as mentioned previously, over previouscurved ones in capturing micro flows with curved walls. B. Couette flow between two concentric cylinders
The proposed boundary scheme is further applied to a microgaseous flow with morecomplex geometries, i.e., the microcylindrical Couette flow between two cylinders, to validatethe theoretical analysis. This classical problem has been studied by many researchers as abenchmark case in fluid dynamics [27, 28, 32, 43, 44]. As shown in Fig. 10, the two cylindersof the problem are concentric with respective radii of R and R ( R < R ), and rotate atangular velocities ω and ω , respectively. As the flow between the cylinders reaches atsteady state, they can be described by the following reduced Navier-Stokes equations incylinder polar coordinates ( r, θ ) d u θ d r + dd r (cid:16) u θ r (cid:17) = 0 , (40)where u θ is the tangential velocity and r is the radial distance. Based on the Maxwell’sdiffusive boundary condition for gas-wall interactions, the slip boundary condition at the25nner and outer cylinder walls can be expressed as u θ (cid:12)(cid:12) r = R = ω R + 2 − σ σ λ (cid:18) d u θ d r − u θ r (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) r = R , u θ (cid:12)(cid:12) r = R = ω R − − σ σ λ (cid:18) d u θ d r − u θ r (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) r = R , (41)where σ and σ are the accommodation coefficients of the inner and outer cylinder surface,respectively. Under the boundary condition (41), the solution of Eq. (40) for the velocityprofile can be analytical solved as [44] u θ = M ω − N ω M − N r + ω − ω N − M r , (42)where M = 1 R (cid:18) − − σ σ λR (cid:19) , N = 1 R (cid:18) − σ σ λR (cid:19) . (43)In the simulations, the two cylinder surfaces are both assumed to be fully diffusive, i.e., σ = σ = σ = 1. To mimic the slip boundary condition (41) at the two cylinders’ surfaces,we set L = (2 − σ ) / , L = 0 in the boundary condition (16) implemented with Eqs.(29)-(31). The Knudsen number for the flow is defined as Kn = λ/ ( R − R ). The radiusratio of the two cylinders, β = R /R , is obtained by changing R under the fixed R = 1 . M grid cells.Considering the small grid sizes used, we set Kn = 0 .
01 and ω = ω = ω = 0 .
001 in thesimulations. As done before, the value of l is taken as zero in case l < − . E is assigned in the range of Eq. (37). In Fig. 11, thetangential velocities with different lattice sizes are shown for two cases of rotating cylinders.Clearly, the velocity profiles predicted by the boundary scheme (16) with Eqs. (29)-(31) arein good agreement with the analytical solutions. Whereas for the halfway DBB scheme, Eqs.(34) and (35) cannot yield grid-independent results consistent with the analytical solutions.This is ascribed to the fact that the midway approximation for the cylinder surfaces distortthe actual curved geometries. Our theoretical analyses on curved boundary conditions aredemonstrated again here.Figure 12 displays the velocity profiles between the cylinders predicted under differentcases of l . Four additional choices of l , i.e., l = γ, γ, γ , γ + γ besides l given by Eq.(29) are considered in the simulations. In all these cases, the combination parameter r isdetermined by Eq. (30), and τ q by Eq. (31) with E = − .
65. From the figure, it is observedthat compared with the other choices of l , the optimal consistency with analytical solutions26 AnalyticalEqs. (29)-(31), M=16Eqs. (29)-(31), M=24Eqs. (29)-(31), M=30Half-way DBB, M=16Half-way DBB, M=24Half-way DBB, M=30
FIG. 11. Velocity profiles of the microcylindrical Couette flow for (a) β = 0 . , ω = 0 . , ω = 0;(b) β = 0 . , ω = 0 , ω = 0 . E = − . AnalyticalEq. (29)l=l=2l= l= + AnalyticalEq. (29)l=l=2l= l= + FIG. 12. Velocity profiles of the microcylindrical Couette flow for (a) β = 0 . , ω = 0 . , ω = 0;(b) β = 0 . , ω = 0 , ω = 0 .
001 at M = 16. is obtained for the case that l conforms to Eq. (29), even with the artificial value l = 0at several boundary points. This also indicates that the uniform relaxation time τ q cannotbe achieved for the other cases of l to realize the slip boundary condition. In addition,the tangential velocity predicted with l = γ in Eqs. (26) and (27) are investigated as γ isartificially fixed to derive the uniform τ q . The simulated results are delineated in Fig. 13,where Kn = 0 .
045 is set to ensure τ q > .
5. Again, apparent deviations from the analyticalsolutions are observed for the artificial approximation of γ . While with Eqs. (29)-(31), the27 AnalyticalEqs. (29)-(31)l= , Eq. (26), Eq. (27), =0.3l= , Eq. (26), Eq. (27), =0.5l= , Eq. (26), Eq. (27), =0.7l= , Eq. (26), Eq. (27), =0.9
FIG. 13. Velocity profiles of the microcylindrical Couette flow at Kn = 0 . M = 16 for(a) β = 0 . , ω = 0 . , ω = 0; (b) β = 0 . , ω = 0 , ω = 0 . l = γ in Eq. (26) for r , and in Eq. (27) toobtain the uniform τ q by fixed values of γ . boundary scheme (16) brings robust agreement results with the analytical solutions. Theseobservations strengthen and confirm the capability of the present boundary scheme to realizethe 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 simulatingmicrogaseous flows with curved geometries. This curved boundary scheme is a combinationof the Maxwellian diffuse reflection scheme and a single-node boundary scheme for curvedno-slip walls. In addition to the distance ratio γ , which is only involved in previous curvedboundary treatments, an additional free parameter l is incorporated in the present localboundary 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 γ and the relaxation times( τ s and τ q ) unitedly affect the slip velocity derived from the boundary scheme. Thanks tothe free parameter l , an available strategy to determine the uniform τ q together with thecombination 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 boundarycondition cannot be realized with an invariable τ q for the halfway DBB scheme and previous28urved boundary schemes.The proposed curved boundary condition is applied to some benchmark problems withplanar and curved walls, including the aligned and inclined microchannel flows and themicrocylindrical Couette flow. To avoid the instability at very small γ in the simulations,the free parameter l is artificially assigned with zero value at several boundary nodes. Goodand robust predictions from the present derivations are obtained to match the analyticalsolutions satisfactorily even under a small lattice size. The numerical results also showthat for microgaseous flows with nonplanar and curved walls, the halfway DBB schemeand previous curved boundary conditions which only contains γ bring clear grid-dependentdeviations from the analytical solutions.On the basis of the present study, a noteworthy point is that to ensure uniform relaxationparameters in realizing a prescribed slip boundary condition at curved walls, adding freeparameters to the boundary scheme would be more efficient than seeking more relaxationparameters in the LBE. It should be noted that this work on curved boundary conditionsfocused on microscale gas flows in the slip regime. However, the present analysis is alsoinstructive for our extension work to the transition regime, where the effect of the Knudsenlayer must be incorporated. In addition, we would like to point out that by replacing theMaxwellian diffusive part by the specular reflection scheme in Eq. (16), it can generateanother hybrid but nonlocal boundary scheme for microscale flows. Following the methodpresented in this work, a similar strategy for uniform relaxation parameters can be obtainedto realize the prescribed slip boundary condition at curved walls. Furthermore, such ideaof introducing a free parameter can be also applied to combining the Maxwellian diffusiveboundary scheme and the specular reflection scheme. Finally, the present work can beextended to the three-dimensional case without much difficulty. These interesting topicswill be left for our future work. ACKNOWLEDGMENTS
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 Uni-versities (No. 2018MS060). L. Wang would like to thank Profs. Wen-An Yong and Zhaoli29uo and Dr. Weifeng Zhao for their helpful discussions. [1] C. M. Ho and Y.C. Tai, Annu. Rev. Fluid Mech. , 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. , 845 (2012).[4] C. Cercignani, Mathematical Methods in Kinetic Theory (Plenum, New York, 1990).[5] X. He, L. S. Luo, Phys. Rev. E , R6333 (1997).[6] X. Shan, X. He, Phys. Rev. Lett. , 65 (1998).[7] X. Nie, G. D. Doolen, and S. Chen, J. Stat. Phys. , 279 (2002).[8] C. Y. Lim, C. Shu, X. D. Niu, and Y. T. Chew, Phys. Fluids , 2299 (2002).[9] C. Shu, X. D. Niu, and Y. T. Chew, J. Stat. Phys. , 239 (2005).[10] Y. H. Zhang, R. Qin, and D. R. Emerson, Phys. Rev. E , 047702 (2005).[11] J. F. Zhang, Microfluid. Nanofluid. , 1 (2005).[12] S. H. Kima, H. Pitsch, I. D. Boyd, J. Comput. Phys. , 8655 (2008).[13] Y. Sone, Kinetic Theory and Fluid Dynamics (Birkh¨auser, Boston, 2002).[14] F. Sharipov, J. Phys. Chem. Ref. Data , 023101 (2011).[15] F. Verhaeghe, L. S. Luo, and B. Blanpain, J. Comput. Phys. , 147 (2009).[16] C. Y. Lim, C. Shu, X. D. Niu, and Y. T. Chew, Phys. Fluids , 2299 (2002).[17] S. Ansumali and I. V. Karlin, Phys. Rev. E , 026311 (2002).[18] S. Succi, Phys. Rev. Lett. , 064502 (2002).[19] M. Sbragaglia and S. Succi, Phys. Fluids , 093602 (2005).[20] G. H. Tang, W. Q. Tao, and Y. L. He, Phys. Fluids , 058101 (2005).[21] Z. H. Chai, Z. L. Guo, L. Zheng, and B. C. Shi, J. Appl. Phys. , 014902 (2008).[22] Z. Guo, B. Shi, T.S. Zhao, and C. Zheng, Phys. Rev. E , 056704 (2007).[23] Z. Guo and C. Zheng, Int. J. Comput. Fluid Dyn. , 465 (2008).[24] Z. H. Chai, B. C. Shi, Z. L. Guo, and J. H. Lu, Commun. Comput. Phys. , 1052 (2010).[25] L. Szalm´ a s, Phys. Rev. E , 066710 (2006).[26] K. Suga, Fluid Dyn. Res. , 034501 (2013).[27] Z. Liu, Z. Mu, H. Wu, Microfluid. Nanofluid.
27 (2019).
28] S. Tao and Z. Guo, Phys. Rev. E , 3519 (2011).[33] W. F. Zhao, J. T. Huang, and W.-A. Yong, Multiscale. Model. Simul. , 854 (2019).[34] Z. L. Guo and C. Shu, Lattice Boltzmann Method and its Applications in Engineering (WorldScientific Press, Singapore, 2013).[35] P. Lallemand and L. S. Luo, Phys. Rev. E , 6546 (2000).[36] Z. L. Guo, T. S. Zhao, and Y. Shi, J. Appl. Phys. , 074903 (2006).[37] L. Wang, S. Tao, X. H. Meng, K. Zhang, and G. Lu, Phys. Rev. E , 0633070 (2020).[38] Z. L. Guo, C. G. Zheng, and B. C. Shi, Phys. Rev. E , 036707 (2008).[39] N. G. Hadjiconstantinou, Phys. Fluids , 111301 (2006).[40] X. He, Q. Zou, L.-S. Luo, and M. Dembo, J. Stat. Phys. , 115 (1997).[41] I. Ginzburg and D. d’Humi` e res, Phys. Rev. E , 066614 (2003).[42] C. Pan, L.-S. Luo, and C. T. Miller, Comput. Fluids
898 (2006).[43] D. A. Lockerby, J. M. Reese, D. R. Emerson, and R. W. Barber, Phys. Rev. E , 017303(2004).[44] S. Yuhong, R. W. Barber, and D. R. Emerson, Phys. Fluids , 047102 (2005)., 047102 (2005).