## Abstract

In the present paper, we provide an analytical expression for the first- and second-order velocity slip coefficients by means of a variational technique that applies to the integrodifferential form of the Boltzmann equation based on the true linearized collision operator and the Cercignani–Lampis scattering kernel of the gas–surface interaction. The polynomial form of the Knudsen number obtained for the Poiseuille mass flow rate and the values of the velocity slip coefficients are analysed in the frame of potential applications of the lattice Boltzmann methods in simulations of microscale flows.

## 1. Introduction

The use of the Lattice Boltzmann (LB) method to describe flows in nano- and microfluidic devices has been an area of intense research in recent years. Since, historically, the LB method was derived as a Navier–Stokes solver, criticism has been raised about the validity of applying the LB method to finite-Knudsen-number flows. However, it is now generally accepted that the classical hydrodynamical equations can supply realistic results well beyond the slip region, provided that higher order boundary conditions are employed. Therefore, one of the critical issues with the LB method is to specify suitable slip boundary conditions to ensure microscale simulation accuracy. Very recently, Sbragaglia & Succi [1] demonstrated how, in the presence of a non-zero slip coefficient, the LB develops a physical slip flow component at the wall that can be matched exactly to experimental data up to second order in the Knudsen number, provided the free parameters appearing in the LB model are tuned in a correct way. Since the LB equation has its roots in kinetic theory, in the current investigation, a variational technique has been used to obtain asymptotic near-continuum solutions of the linearized Boltzmann equation for hard-sphere molecules and general models of boundary conditions, in order to provide analytical expressions for the first- and second-order velocity slip coefficients that enable the mesoscopic approach presented in Sbragaglia & Succi [1] to become predictive.

## 2. Mathematical formulation

Let us consider two plates separated by a distance *d* and a gas flowing parallel to them, in the *z*-direction, owing to a pressure gradient. Both boundaries, fixed at *x*=±*d*/2, are held at a constant temperature *T*_{0}. For small values of the pressure gradient, it can be assumed that the velocity distribution of the flow is nearly the same as that occurring in an equilibrium state. Under these conditions, the Boltzmann equation can be linearized [2],
2.1
where *h*(*x*,**c**) is the small perturbation on the basic equilibrium state, **c** is the molecular velocity expressed in units of (2*RT*)^{1/2} (with *R* being the gas constant and *T* being the temperature), *k*=(1/*p*)(∂*p*/∂*z*)=(1/*ρ*)(∂*ρ*/∂*z*) (with *p* and *ρ* being the gas pressure and density, respectively) and *Lh* is the linearized collision operator. In the current investigation, we analyse the Poiseuille problem on the basis of the linearized collision operator for hard-sphere molecules in order to obtain a better approximation of real-gas behaviour for isothermal flows. Appropriate boundary conditions on the two plates must be supplied for the Boltzmann equation (2.1) to be solved. Once the deviation *h* from the equilibrium distribution is known, the bulk velocity of the gas can be written as
2.2
and the flow rate *F* (per unit time through unit thickness) as
2.3
To apply the variational formulation, we shall rewrite equation (2.1) in the following symbolic form:
2.4
where *Dh*=*c*_{x}(∂*h*/∂*x*), *S*=−*kc*_{z}. The boundary conditions to be matched with equation (2.4) have the general expression [3]: *h*^{+}=*Kh*^{−}, where the explicit form of the operator *K* depends on the scattering kernel used. Here, *h*^{±} are the restrictions of the function *h*, defined on the boundary, to positive (negative) values of *c*_{x}. Most of the works in rarefied gas dynamics are based on the use of the classical Maxwell gas–surface interaction law, characterized by a single accommodation coefficient *α*, while in practice, every physical property has its own accommodation coefficient. In order to include better physics, in the current investigation, we will focus on the Cercignani–Lampis (CL) scattering kernel [3], which is based on two different adjustable parameters. In this case, the boundary conditions can be written as
2.5
where
2.6
with **c**_{t}=(*c*_{y},*c*_{z}) being the two-dimensional vector of the tangential molecular velocity and *I*_{o} the modified Bessel function of first kind and zeroth order. Physically, the parameters *α*_{t} and *α*_{n} appearing in equation (2.6) represent the accommodation coefficient of the tangential momentum and the accommodation of the kinetic energy owing to the velocity normal to the bounding walls, respectively. The CL model recovers, as limiting cases, the specular reflection (for *α*_{t}=*α*_{n}=0) and the diffuse re-emission (for *α*_{t}=*α*_{n}=1). Using the variational principle described in Cercignani & Lorenzani [2], we introduce the following functional *J* of the test function :
2.7
where *P* is the parity operator in velocity space and ((,)), (,)_{B} denote two scalar products defined in Cercignani & Lorenzani [2]. The functional attains its minimum value when solves equation (2.4) with appropriate boundary conditions. If we let , equation (2.7) gives: *J*(*h*)=−((*PS*,*h*)). Looking at the definitions (2.2) and (2.3), it can be easily shown that the stationary value of *J* is related to a quantity of physical interest, the flow rate of the gas, through the relation: *F*=−*ρ*/*kJ*(*h*). Since the purpose of the present investigation is to provide an analytical expression for the first- and second-order velocity slip coefficients, it is sufficient to consider asymptotic results (near-continuum) for mass flow rates. Therefore, the following simplified test function, rescaled by the relative pressure gradient *k* and the length parameter *θ*, has been used to evaluate equation (2.7):
2.8
which is the same trial function that, for the Bhatnagar–Gross–Krook (BGK) kinetic model, has been proved to produce very accurate results in the near-continuum flow limit [2]. In equation (2.8), *A* and *B* are adjustable constants to be varied in order to obtain the best value of . Substituting , given by equation (2.8), into equation (2.7), we obtain the following polynomial of the second order with respect to the constants *A* and *B*, which are to be determined:
2.9
where the same symbol has been used to denote the dimensionless quantity . The coefficients in non-dimensional form are given by
2.10
2.11
2.12
where and . In equations (2.10)–(2.12), *δ*=*d*/*θ* is the rarefaction parameter (inverse Knudsen number) and the symbol *J*_{i} stands for integral expressions defined by using the brackets [*ϕ*,*ψ*],
2.13
with *Lψ* being the linearized Boltzmann collision operator. For hard spheres of diameter *σ*, the length parameter *θ* is given by and the mean free path *λ* reads as (where *n* is the number density). Therefore,
2.14
where *ψ* is a function of **c**, while *ψ*_{1} refers to **c**_{1}. *V* is the relative velocity |**c**−**c**_{1}|. *ψ*′≡*ψ*(**c**′) and *ψ*′_{1}≡*ψ*(**c**′_{1}), where **c**′ and **c**′_{1} are the velocities after collision of two molecules with velocities **c** and **c**_{1}. The collision geometry, in conjunction with the conservation laws, relates the velocities after collision to the velocities before collision [2]. In equation (2.14), *Θ* is the angle through which the relative velocity has turned and *ϵ* is the azimuthal angle the plane containing the relative velocities before and after collision makes with a fixed reference plane. The integrals *J*_{1} and *J*_{2} appearing in equation (2.10) are eight-fold integrals: *J*_{1}=[*c*_{x}*c*_{z},*c*_{x}*c*_{z}] and *J*_{2}=−[*c*_{x}^{2}*c*_{z},*c*_{x}^{2}*c*_{z}], where . The accuracy of our variational analysis solution is mostly related to the precision in the computation of such integrals, which have been numerically evaluated using a Monte Carlo integration on an eight-dimensional space (see equations (2.13)–(2.14)) [2]. The derivatives of with respect to *A* and *B* vanish in correspondence with the optimal values of these constants. The resulting expression for the minimum of is
2.15
Thus, the computation of the optimal value of the functional *J*(*h*) (equation (2.15)) will lead to an accurate estimate of the flow rate of the gas, which in non-dimensional form is
2.16

## 3. Near-continuum solution and slip coefficients

A notable advantage of the variational approach is that it permits one to write down simple approximate equations, computed via kinetic theory, to be used in classical hydrodynamic numerical tools in order to extend the continuum description (that is significantly more efficient compared with molecular-based approaches), even beyond the slip regime. In recent years, there has been considerable success in the implementation of second-order slip boundary conditions to extend the Navier–Stokes equations into the transition regime. Assuming a second-order boundary condition at a flat wall, in the isothermal case, the slip velocity is [2]
3.1
where the gas-velocity gradients are evaluated at the wall. Here, *λ* is the mean free path of the molecules, , where *μ* is the gas viscosity, *P* is the pressure, *T* is the absolute temperature and *R* is the universal gas constant. In equation (3.1), *A*_{1} and are the first- and second-order slip coefficients, respectively. Recent experimental studies [2] have revealed large discrepancies between the experimentally determined values of the second-order slip coefficient and the theoretical values listed in the literature, most of which include empirical parameters suitably *a posteriori* calibrated and, therefore, not useful in a predictive sense. The lack of a well-founded value of makes it difficult to extend slip-flow predictions into the transition regime. The asymptotic near-continuum solution for the Poiseuille mass flux obtained by means of our variational technique can be used to predict slip coefficients. When the linearized Boltzmann equation for hard-sphere molecules is considered and the CL scattering kernel is used to describe the gas–wall interaction, equations (2.15) and (2.16) give, in the limit *δ*≫1,
3.2
where
3.3
3.4
3.5
3.6
3.7
and
3.8
with and . Rewriting equation (3.2) in terms of deviations from the continuum solution, one obtains
3.9
where the Knudsen number (*Kn*) is given by . A comparison with the solution of the Navier–Stokes equations obtained by using the boundary condition (3.1),
3.10
gives
3.11
where the term *O*(1) indicates that a kinetic boundary layer, of the order of *Kn*^{2}, needs to be superposed to the Navier–Stokes component of the flow field in order to capture the true solution of the Boltzmann equation [2]. *A*_{2}, evaluated starting from mass flow rate computations, can be considered a sort of ‘effective’ second-order slip coefficient since it includes the contribution of the Knudsen layer. Once the effective second-order slip coefficients are evaluated by means of our variational computation of the Poiseuille flow rate, one can filter out the Knudsen layer effects, following the theoretical considerations presented below. In the BGK approximation, the Poiseuille flow rate can be split as follows:
3.12
where represents the Navier–Stokes component of the true flow field, while the second term in the integral (3.12) accounts for the contribution of the Knudsen layer. According to Cercignani & Lorenzani [2], the Knudsen layer contribution to the flow rate is largely independent of the molecular model considered and, therefore, also for a hard-sphere gas, *ξ* retains approximately the same value it keeps for the BGK model equation, *ξ*≃0.296. Thus, the second-order slip coefficient that should be used in the boundary conditions (3.1) is given by . The latter highlights that a successful second-order slip model is one that does not agree with the Boltzmann equation solutions in the Knudsen layers, in order to be able to predict within the same framework all the physical quantities of interest: not only the flow rate but also the stress field, which is not subject to a Knudsen layer correction.

## 4. Results and discussion

The first-order (*A*_{1}) and second-order (*A*_{2}) slip coefficients, obtained on the basis of our variational solution of the linearized Boltzmann equation for hard-sphere molecules, are given in tables 1 and 2, for several values of the accommodation coefficients *α*_{t} and *α*_{n}. To test the reliability of our variational approach, we list in table 1, the values of the first-order slip coefficient *A*_{1} computed by using equation (3.11), along with the highly accurate numerical results reported by Siewert [4]. We have rescaled the values of the viscous slip coefficient given by Siewert [4] in order to use the same definitions of mean free path and units as in the present investigation. As shown in the table, the agreement is fairly good with the relative error lying within 1.2 per cent. This comparison reveals that our variational datas are more accurate than those reported in McCormick [5], where analytical equations for the first-order viscous slip have been obtained by means of a variational principle applied to the linear Boltzmann equation for a rigid-sphere gas and the CL scattering kernel. We have also included in tables 1 and 2 results based on a variational solution of the linearized BGK model equation in order to provide evidence that the BGK first-order slip coefficients are similar to those determined by the solution of the Boltzmann equation for hard-sphere molecules, while the second-order slip coefficients are found to be significantly dependent on the interaction model. As stated in §2, the Maxwell gas–surface interaction law is characterized by a single accommodation coefficient *α* for every molecular property. Thus, one may expect that assuming the same accommodation coefficients, the Maxwell kernel and the CL model will provide results close to each other.

In table 3, a comparison of the slip coefficients, obtained for both kernels under the condition *α*=*α*_{t}=*α*_{n}, is performed. One can see that the dependence of the first-order slip coefficient on the scattering kernel used is weak, while the second-order slip coefficient depends significantly on the gas–surface interaction model. In spite of that, whatever model of boundary conditions is chosen, the estimates we obtain for *A*_{2}, in the case of hard-sphere molecules, seem inconsistent with most of the available theoretical models derived on heuristic or empirical grounds, while they are very close to the values obtained in recent experiments (see [2] for a detailed discussion). This means that our analytical expressions for the first- and second-order velocity slip coefficients can be used to adjust the free parameters appearing in the mesoscopic approach of the LB boundary conditions proposed in Sbragaglia & Succi [1], in order to predict experimental data. In recent years, the LB method has been recognized as one of the most promising approaches for the simulation of microscale gas flows provided that boundary conditions derived from continuous kinetic theory are employed. Two basis types of kinetic boundary condition have been proposed: one is the discrete Maxwellian boundary condition, which is a straightforward discretization of Maxwell’s diffuse-specular-reflection boundary condition in kinetic theory [6–8]; the other is the combination boundary condition (CBC), which combines the no-slip bounce-back and the free-slip specular-reflection schemes [1, 8]. Recently, it was shown that both schemes are virtually equivalent in principle [8] and that they correspond to a second-order slip boundary condition [1, 8]. In particular, the CBC approach provides free tunable parameters to model slip flow phenomena. Therefore, comparing the velocity slip coefficients obtained through the continuous Boltzmann equation with those derived by using the LB method with the combination boundary conditions, written in terms of adjustable parameters, one can fix the values of these free parameters, in order to recover quantitative agreement up to second order in the Knudsen number, and then use the LB method as a valuable tool for predictive simulations.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society