Discrete Boltzmann method with Maxwell-type boundary condition for slip flow

TL;DR

This paper develops a discrete Boltzmann model with Maxwell-type boundary conditions to accurately simulate slip flow in microchannels, capturing velocity slip and Knudsen layer effects across various conditions.

Contribution

It introduces a Maxwell-type boundary condition with tangential momentum accommodation coefficient into the discrete Boltzmann model for slip flow simulation.

Findings

Successfully describes velocity slip and Knudsen layer effects.

Validates the model with Couette and Poiseuille flow simulations.

Clarifies the relation between different Knudsen number definitions.

Abstract

The rarefied effect of gas flow in microchannel is significant and cannot be well described by traditional hydrodynamic models. It has been know that discrete Boltzmann model (DBM) has the potential to investigate flows in a relatively wider range of Knudsen number because of its intrinsic kinetic nature inherited from Boltzmann equation. It is crucial to have a proper kinetic boundary condition for DBM to capture the velocity slip and the flow characteristics in the Knudsen layer. In this paper, we present a DBM combined with Maxwell-type boundary condition model for slip flow. The tangential momentum accommodation coefficient is introduced to implement a gas-surface interaction model. Both the velocity slip and the Knudsen layer under various Knudsen numbers and accommodation coefficients can be well described. Two kinds of slip flows, including Couette flow and Poiseuille flow, are simulated to verify the model. To dynamically compare results from different models, the relation between the definition of Knudsen number in hard sphere model and that in BGK model is clarified.

Tables

Table 1. Table 1: Discrete velocity model, where λ = 1 3 𝜆 1 3 \lambda=\frac{1}{\sqrt{3}} , φ = 1 + 5 2 𝜑 1 5 2 \varphi=\frac{1+\sqrt{5}}{2} , and ϕ = 2 5 + 5 italic-ϕ 2 5 5 \phi=\frac{\sqrt{2}}{\sqrt{5+\sqrt{5}}}

$i$ direaction	Unit vector($v_{ix},v_{iy},v_{iz}$)
$i=1-8$	$\lambda (\pm 1,\pm 1,\pm 1)$
$i=9-12$	$\lambda (0,\pm {\phi }^{-1},\pm \phi )$
$i=13-16$	$\lambda (\pm \phi ,0,\pm {\phi }^{-1})$
$i=17-20$	$\lambda (\pm {\phi }^{-1},\pm \phi ,0)$
$i=21-24$	$\varphi (0,\pm \phi ,\pm 1)$
$i=25-28$	$\varphi (\pm 1,0,\pm \phi )$
$i=29-32$	$\varphi (\pm \phi ,\pm 1,0)$

Table 2. Table 2: Coefficients for several kinds of values of α 𝛼 \alpha .

$\alpha$	$A_0$	$A_1$	$A_2$	$A_3$	$A_4$	$A_5$	$A_6$	$A_7$	$k_s$
0.2	-0.8601	2.8585	-10.0275	18.7261	-19.2088	10.4579	-2.8826	0.3036	$8.2248$
0.5	-0.6698	2.0540	-7.0263	12.9588	-13.2973	7.2675	-2.0206	0.2141	$2.8612$
0.8	-0.5008	1.4022	-4.6499	8.4399	-8.6668	4.7619	-1.3393	0.1431	$1.4877$
1.0	-0.3989	1.0437	-3.3750	6.0435	-6.2111	3.4289	-0.9740	0.1047	$1.0162$

Equations

$$Kn=\frac{\lambda }{L}.$$

$${\lambda }_{HS}=\frac{1}{\sqrt{2}n\pi d^2},$$

$$\mu =\frac{5}{16}\frac{\sqrt{mkT/\pi }}{d^2},$$

$${\lambda }_{HS}=\frac{4}{5}\frac{\mu }{p}\sqrt{\frac{8RT}{\pi }}.$$

$${\lambda }_{BGK}=\tau \overline{c},$$

$$\overline{c}=\sqrt{\frac{8RT}{\pi }}$$

$$\mu =\tau p.$$

$${\lambda }_{BGK}=\frac{\mu }{p}\sqrt{\frac{8RT}{\pi }}.$$

$${\lambda }_{HS}=\frac{4}{5}{\lambda }_{BGK}.$$

$$\frac{\partial f_{ki}}{\partial t}+v_{ki\alpha }\frac{\partial f_{ki}}{\partial r_{\alpha }}-\frac{a_{\alpha }(v_{ki\alpha }-u_{\alpha })}{T}{f}_{ki}^{eq}=-\frac{1}{\tau }(f_{ki}-f_{ki}^{eq}),$$

$$f_{ki}^{eq}=\rho F_k[(1-\frac{u^2}{2T}+\frac{u^4}{8T^2})+\frac{1}{T}(1-\frac{u^2}{2T})v_{ki\xi }{u}_{\xi }$$

$$+\frac{1}{2T^2}(1-\frac{u^2}{2T})v_{ki\xi }{v}_{ki\eta }{u}_{\xi }{u}_{\eta }$$

$$+\frac{1}{6T^3}{v}_{ki\xi }{v}_{ki\eta }{v}_{ki\tau }{u}_{\xi }{u}_{\eta }{u}_{\tau }$$

$$+\frac{1}{24T^4}{v}_{ki\xi }{v}_{ki\eta }{v}_{ki\tau }{v}_{ki\chi }{u}_{\xi }{u}_{\eta }{u}_{\tau }{u}_{\chi }],$$

$$F_0=1-32(F_1+F_2+F_3+F_4),$$

$$F_1=\frac{1}{c_1^2(c_1^2-c_2^2)(c_1^2-c_3^2)(c_1^2-c_4^2)}$$

$$\times [\frac{945}{32}{T}^4-\frac{105}{32}(c_2^2+c_3^2+c_4^2)T^3$$

$$+\frac{15}{32}(c_2^2{c}_3^2+c_3^2{c}_4^2+c_4^2{c}_2^2)T^2-\frac{3}{32}{c}_2^2{c}_3^2{c}_4^{2}T],$$

$$F_2=\frac{1}{c_2^2(c_2^2-c_3^2)(c_2^2-c_4^2)(c_2^2-c_1^2)}$$

$$\times [\frac{945}{32}{T}^4-\frac{105}{32}(c_3^2+c_4^2+c_1^2)T^3$$

$$+\frac{15}{32}(c_3^2{c}_4^2+c_4^2{c}_1^2+c_1^2{c}_3^2)T^2-\frac{3}{32}{c}_3^2{c}_4^2{c}_1^{2}T],$$

$$F_3=\frac{1}{c_3^2(c_3^2-c_4^2)(c_3^2-c_1^2)(c_3^2-c_2^2)}$$

$$\times [\frac{945}{32}{T}^4-\frac{105}{32}(c_4^2+c_1^2+c_2^2)T^3$$

$$+\frac{15}{32}(c_4^2{c}_1^2+c_1^2{c}_2^2+c_2^2{c}_4^2)T^2-\frac{3}{32}{c}_4^2{c}_1^2{c}_2^{2}T],$$

$$F_4=\frac{1}{c_4^2(c_4^2-c_1^2)(c_4^2-c_2^2)(c_4^2-c_3^2)}$$

$$\times [\frac{945}{32}{T}^4-\frac{105}{32}(c_1^2+c_2^2+c_3^2)T^3$$

$$+\frac{15}{32}(c_1^2{c}_2^2+c_2^2{c}_3^2+c_3^2{c}_1^2)T^2-\frac{3}{32}{c}_1^2{c}_2^2{c}_3^{2}T].$$

$$f_{ki}^{t+\Delta t}=f_{ki}^t-v_{ki\alpha }\frac{\partial f_{ki}}{\partial r_{\alpha }}\Delta t-\frac{1}{\tau }(f_{ki}-f_{ki}^{eq})\Delta t$$

$$+\frac{a_{\alpha }(v_{ki\alpha }-u_{\alpha })}{T}{f}_{ki}^{eq}\Delta t.$$

$$\frac{{\partial{f_{ki}}}}{{\partial{r_{\alpha}}}}=\left\{\begin{array}[]{l}\frac{{3{f_{ki,I}}-4{f_{ki,I-1}}+{f_{ki,I-2}}}}{{2\Delta{r_{\alpha}}}}{\kern 14.0pt}{\rm{if}}{\kern 5.0pt}{v_{ki\alpha}}\geq 0,\\ \frac{{3{f_{ki,I}}-4{f_{ki,I+1}}+{f_{ki,I+2}}}}{{-2\Delta{r_{\alpha}}}}{\kern 14.0pt}{\rm{if}}{\kern 5.0pt}{v_{ki\alpha}}<0.\\ \end{array}\right.$$

$$\frac{\partial f_{ki}}{\partial r_{\alpha }}=\frac{f_{ki,2}-f_{ki,1}}{\Delta r_{\alpha }},$$

$$\frac{\partial f_{ki}}{\partial r_{\alpha }}=\frac{f_{ki,N-1}-f_{ki,N}}{-\Delta r_{\alpha }},$$

$$f_{ki,N}={f_{ki}}^{eq}({\rho }_w^R,u_w^R,T_w^R),v_{ki\alpha }<0,$$

$$f_{ki,1}={f_{ki}}^{eq}({\rho }_w^L,u_w^L,T_w^L),v_{ki\alpha }>0.$$

$$\sum _{c_{ki\alpha }>0}{f}_{ki,N}{c}_{ki\alpha }+{\rho }_w^R\sum _{c_{ki\alpha }<0}{f}_{ki}^{eq}(1.0,v_w^R,e_w^R)c_{ki\alpha }=0,$$

$$\sum _{c_{ki\alpha }<0}{f}_{ki,1}{c}_{ki\alpha }+{\rho }_w^L\sum _{c_{ki\alpha }>0}{f}_{ki}^{eq}(1.0,{v_w}^L,e_w^L)c_{ki\alpha }=0.$$

Full text

PACS numbers

51.10.+y, 47.11.-j, 47.45.-n.

Key words

Rarefied gas, Discrete Boltzmann method, Maxwell-type boundary, slip flow.

Suggested keywords

††preprint: APS/123-CTP

I Introduction

In recent years, the development of natural science and engineering technology has moved towards miniaturization. One of the most typical examples is Micro-Electro- Mechanical System (MEMS)Gad2001Flow ; Zhang2005Gas ; Bao2005Linear ; ChingShen2005Rarefied . It is extremely important to investigate the underlying physics of unconventional phenomena at the micro-scale. Those unconventional phenomena cannot be explained by the traditional macro-model and has become a key bottleneck limiting the further development of MEMS. Among these unconventional physical problems, the gas flow and heat transfer characteristics at the mico-scale are especially critical.

Due to the reduction of the geometric scale, the mean free path of gas molecules may be comparable to the length scale of the device. The Knudsen number ($Kn$), a dimensionless parameter used to measure the degree of rarefaction of the flow and defined as the ratio of the mean free path of molecules to characteristic length of the device, may be larger than 0.001 and reaches slip-flow regime(0.001 $<$ $Kn$ $<$ 0.1 ) or even transition-flow regime (0.1 $<$ $Kn$ $<$ 10). As we know, Navier-Stokes (NS) equations are applicable to continuum flow ($Kn$ $<$ 0.001) where the continuum hypothesis is acceptable. The slip boundary models based on kinetic theory should be adopted to describe the slip-flow regime. While in the transition-flow region, the NS equations totally fail to calculate viscous stress and heat flow accurately and higher order equations such as Burnett equationBurnett1935The ; Chenweifang2016 and Grad’s 13-moment equationsgrad1949kinetic are needed.

In practical, gas flow in a microchannel may encounter continuum, slip and transition regimes simultaneously. Traditionally macro-scale models cannot apply to such a broad range of Knudsen numbers by only one set of equations. Besides, the numerical solutions of higher order macro-equations are difficult to obtain because of the numerical stability problem Bao2005Linear . It is commonly accepted that Direct Simulation Monte Carlo (DSMC)Bird2003Molecular ; ChingShen2005Rarefied is a accurate method for rarefied gas flow which has been verified by experimental data. However, the computation cost in numerical simulations is too expensive for low speed gas flow. To reduce the huge ratio of the noise to the useful information, extremely large sample size is needed. Although, the information preservation (IP) method was presented to treat this problemFan2001Statistical , the contradiction between the noise problem and sample size has not been well solved.

It has been know that rarefied gas dynamics are represented properly by the Boltzmann equation due to its kinetic nature. That is continuum, slip and transition regimes can be described by one equation. Unfortunately, the Boltzmann equation is a 6-dimensional problem and the computational cost is formidable to solve such an equation by numerical method directly. In order to alleviate the heavy computational burden of directly solving Boltzmann equation, a variety of Boltzmann equation-based methods, such as the unified gas-kinetic scheme (UGKS)Xu2004Microchannel ; Xu2016A , the discrete velocity method (DVM)Yang2016Numerical , the discrete unified gas-kinetic scheme (DUGKS)Guo2013Discrete ; Guo2015Discrete , the lattice Boltzmann method (LBM)Zhang2005Gas ; Guo2007An ; Watari2009 ; Shan2006Kinetic ; Sofonea2005Boundary ; Tang2008Lattice ; Li2011Lattice ; li2016lattice ; Ansumali2002Kinetic ; Succi2002Mesoscopic ; Niu2004A ; Tang2005Lattice ; Guo2007Discrete , have been presented and well developed. Recently, Discrete Boltzmann Method (DBM) has also been developed and widely used in various complex flow simulationsXu2015Progess ; Xu2016Progess ; Xu2016Complex , such as multiphase flowsGan2015discrete , flow instabilitylai2016nonequilibrium ; chen2016viscosity , combustion and detonationlin2016double ; zhang2016kinetic , etc. From the viewpoint of numerical algorithm, similar to finite-different LBM, the velocity space is substituted by a limited number of particle velocities in DBM. However, the discrete distribution function in DBM satisfy more moment relations which make it fully compatible with the macroscopic hydrodynamic equations including energy equation. The macroscopic quantities, including density, momentum, and energy are calculated from the same set of discrete distribution functions. From the viewpoint of physical modeling, beyond the traditional macroscopic description, the DBM presents two sets of physical quantities so that the nonequilibrium behaviors can have a more complete and precise description. One set includes the dynamical comparisons of nonconserved kinetic moments of distribution function and those of its corresponding equilibrium distribution function. The other includes the viscous stress and heat flux. The former describes the specific nonequilibrium flow state, the latter describes the influence of current state on system evolution. The study on the former helps understanding the latter and the nonlinear constitutive relationsXu2014Progess . The new observations brought by DBM have been used to understand the mechanisms for formation and effects of shock wave, phase transition, energy transformation and entropy increase in various complex flowsLin2014Polar ; Gan2015discrete ; zhang2016kinetic , to study the influence of compressibility on Rayleigh-Taylor instabilitylai2016nonequilibrium ; chen2016viscosity . In a recent study, it is interesting to find that the maximum value point of the thermodynamic nonequilibrium strength can be used to divide the two stages, spinnodal decomposition and domain growth, of liquid-vapor separation.

Some of the new observations brought by DBM, for example, the nonequilibrium fine structures of shock waves, have been confirmed and supplemented by the results of molecular dynamicsLiu2016Molecular ; Liu2017Molecular ; Liu2017Recent . It should be pointed out that the molecular dynamics simulations can also gives microscopic view of points to the origin of the slip near boundary, such as the non-isotropic strong molecular evaporation flux from the liquidKang2008Thermal , which might help to develop more physically reasonable mesoscopic models for slip-flow regime.

In order to extend DBM to the micro-fluid, it is critical to develop a physically reasonable kinetic boundary condition. Many efforts have been made to devise mesoscopic boundary condition for LBM to capture the slip phenomenonAnsumali2002Kinetic ; Succi2002Mesoscopic ; Niu2004A ; Tang2005Lattice ; Guo2007Discrete ; Zhang2005Gas . However, the previous works are most suitable for two-dimensional (2D) models with a very small number of particle velocities and can not directly applicable to the DBM. On the other hand, those boundary conditions fail to capture flow characteristics in the Knudsen layer so the effective viscosity or effective relaxation time approach needs to be adoptedGuo2007An ; Tang2008Lattice . Besides, the results of LBM and DBM should be verified by the results of continuous Botlzmann equation. In 2009, WatariWatari2009 gave a general diffuse reflection boundary for his thermal LB modelWatari2006 and investigated the velocity slip and temperature jump in the slip-flow regime. Then, in his sequent workWatari2010 , he compared the relationship between accuracy and number of particle velocities in velocity slip. Two types of 2D models, octagon family and D2Q9 model, are used. It was found that D2Q9 model fails to represent a relaxation process in the Knudsen layer and the accuracy of the octagon family is improved with the increase in the number of particle velocities. However, all the boundary conditions were set as diffuse reflection wall and the tangential momentum accommodation coefficient (TMAC) was not taken into account.

Because of the dependence of the mean free path on microscopic details of molecular interaction, especially the collision frequency, the Knudsen number may have different values in various interaction models for the same macroscopic properties.

In this paper, we first clarify the definitions of Knudsen number and the connection between the hard sphere model and BGK model for three-dimensional (3D) condition so that the results obtained from various models can be compared dynamically. Then a general Maxwell-type boundary condition for DBM is represented and accommodation coefficient is introduced to implement a gas-surface interaction model. Two kinds of gas flows, Couette flow and Poiseuille flow, in a microchannel are simulated. In the section of Couette flow, the relation between the analysis solutions based on hard sphere and BGK model are verified. The simulation results with various Knudsen numbers and accommodation coefficients are compare with analytical ones based on linear Boltzmann equation in the literature not only on the velocity slip but on the Knudsen profiles. While in the section of Poiseuille flow, the simulation results are compare with analytical solution based on Navier-Stokes equation and the second order slip boundary condition.

II Models and Methods

II.1 Definition of Knudsen number

The Knudsen number is defined as the ratio of the free path of molecules ($\lambda$) to the characteristic length ($L$),

[TABLE]

Throughout the paper, we consider the characteristic length $L$ as unit, so the Knudsen number $Kn$ is equal to the value of $\lambda$.

For the hard sphere collision model, the molecules are considered as hard spheres with diameter $d$, the mean free path of molecules $\lambda_{HS}$ can be calculated by

[TABLE]

where $n$ is the number density of moleculesChingShen2005Rarefied . According to Chapman and EnskogChapman1953The , the viscosity coefficient $\mu$ of hard sphere molecules can be expressed by

[TABLE]

where $k$ is the Boltzmann constant, $m$ is the molecular mass, and $T$ is the temperature. It should be note that gas constant $R$ can be expressed by $R=k/m$. Combining the state equation of ideal gas ($p=\rho RT$), we have the following relationship between $\lambda_{HS}$ and macroscopic quantities:

[TABLE]

For the BGK model, the mean free path of molecules $\lambda_{BGK}$ is defined as

[TABLE]

where $\tau$ is the reciprocal of collision frequency and called relaxation time, $\overline{c}$ is the average thermal speed ChingShen2005Rarefied . The definition of $Kn$ in DBM is in accordance with the definition here.

According to the kinetic theory of gas molecules, $\overline{c}$ is expressed by

[TABLE]

in 3D case. From the Chapman-Enskog expansion, we know that $\tau$ has the following relation with macroscopic quantities:

[TABLE]

Consequently, it has

[TABLE]

The comparison of Eq.(4) and Eq.(8) yields the relationship between viscosity-based mean free path $\lambda_{HS}$ and $\lambda_{BGK}$,

[TABLE]

II.2 Discrete Boltzmann Model

The 3D discrete Boltzmann model taking into account the effect of the external force was presented based on the thermal model represented by WatariWatari2006 . The evolution of the discrete distribution function $f_{ki}$ for the velocity particle $\mathbf{v}_{ki}$ is given as

[TABLE]

where the variable $t$ is the time, $r_{\alpha}$ is the spatial coordinate and $\tau$ is the relaxation-time constant. $a_{\alpha}$ and $u_{\alpha}$ denote the macroscopic acceleration and velocity, respectively, in the $r_{\alpha}$ direction. $T$ denotes the temperature. $f_{ki}^{eq}$ is the local equilibrium distribution function. The subscript $k$ indicates a group of the velocity particles whose speed is $c_{k}$ and $i$ indicates the direction of the particles. The subscript $\alpha$ indicates an $x$, $y$, or $z$ component.

To recover the NS equations, the local equilibrium distribution function should retain up to the fourth order terms of flow velocity. The discrete local equilibrium distribution $f_{ki}^{eq}$ containing the fourth rank tensor is written as

[TABLE]

The velocity particles $\mathbf{v}_{ki}$ consist of a rest particle and 32 moving particles. Each moving particle has four speeds and can be obtained from the unit vectors in Table 1 multiplied by difference $c_{k}$. The speeds $c_{k}$ of moving particle is determined according to the method presented by Watari in RefWatari2009 .

The $F_{k}$ in Eq.(11) is the weighting coefficient for the particle velocity $v_{ki}$ and is determined by $c_{k}$ using the following equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

II.3 Boundary condition models

To solve the evolution equation (10), finite-difference method is adopted. The spatial derivative is solved by the second-order upwind scheme and time derivative is solved by the first-order forward scheme. Then the evolution equation (10) can be rewritten as

[TABLE]

The derivation at position $I$ (see Fig.1)is calculated by

[TABLE]

For ${v_{ki\alpha}}\geq 0$, the evolution equation (17) with Eq.(18) is applied from $I=3$ up to the right wall. However, at the node $I=2$, the second-order upwind scheme in Eq.(18) is not applicable. The first-order upwind scheme,

[TABLE]

is applied there. For ${v_{ki\alpha}}<0$, the evolution equation (17) with Eq.(18) is applied from the left wall to the node $I=N-2$. Likewise, at the node $I=N-1$, the first-order upwind scheme,

[TABLE]

is applied.

To solve the value of the distribution function on the left wall for ${v_{ki\alpha}}>0$ and on the right wall for ${v_{ki\alpha}}<0$, boundary condition models are required.

II.3.1 Diffuse reflection boundary

The complete diffuse reflection model assumes that the molecules leaving the surface with a local equilibrium Maxwellian distribution irrespective of the shape of the distribution of the incident velocity. It can be expressed as

[TABLE]

[TABLE]

The equilibrium distribution functions, ${f_{ki}}^{eq}(\rho_{w}^{R},u_{w}^{R},T_{w}^{R})$ and ${f_{ki}}^{eq}(\rho_{w}^{L},u_{w}^{L},T_{w}^{L})$, are determined from the wall conditions including the velocities and the surface temperatures. Using the zero-mass flow normal to the wallWatari2009 , the density $\rho_{w}^{R}$ and $\rho_{w}^{L}$ can be respectively calculated by the following two equations:

[TABLE]

[TABLE]

As a result, the distribution function on the left wall (${f_{ki,1}}$) for ${v_{ki\alpha}}>0$ and on the right wall (${f_{ki,N}}$) for ${v_{ki\alpha}}<0$ are solved under the diffuse reflection boundary condition.

II.3.2 Specular reflection boundary

The specular reflection model assumes that the incident molecules reflect on the wall surface as the elastic spheres reflect on the entirely elastic surface. The component of the relative velocity normal to the surfaces reverses its direction while the components parallel to the surface remain unchanged. As an example, the direction normal to the wall surface parallels to the $x$ axis, then the molecules leave the surface with a distribution as

[TABLE]

[TABLE]

Since the distribution function ${f_{ki,N}}(-{v_{kix}},{v_{kiy}},{v_{kiz}})$ for ${v_{kix}}<0$ and ${f_{ki,1}}(-{v_{kix}},{v_{kiy}},{v_{kiz}})$ for ${v_{kix}}>0$ can be solved by Eq.(17) with Eq.(18), the distribution function on the right wall (${f_{ki,N}}$) for ${v_{ki\alpha}}<0$ and on the left wall (${f_{ki,1}}$) for ${v_{ki\alpha}}>0$ are easy calculated from Eqs.(25) and (27).

II.3.3 Maxwell-type boundary

In practice, the real reflection of molecules on the body surfaces cannot be described properly by complete diffuse reflection or pure specular reflection. So the Maxwell-type reflection model which is composed of the two reflection modes is needed. The TMAC, $\alpha$ is introduced to measure the proportion of complete diffuse reflectionChingShen2005Rarefied . The $\alpha$ portion of the incident molecules reflect diffusely and the other ($1-\alpha$) portion reflect specularly. The value of TMAC is used to characterize the degree to which the reflected molecules has adjusted to the tangential momentum of the surface,

[TABLE]

where $\tau_{i}$ and $\tau_{r}$ are the tangential components of the momentum fluxes of the incident and reflected molecules, respectively. $\tau_{w}$ is the tangential momentum fluxes of the molecules in the wall. $\alpha=1$ corresponds to the case of complete tangential momentum accommodation and the molecules reflect with the Maxwellian distribution under wall condtion, $u_{w}$ and $T_{w}$. $\alpha=0$ corresponds to the the case when the incident molecules are entirely not adjusted to the conditions of the surface, $\tau_{r}=\tau_{i}$.

Under this boundary condition, the distribution function on the right wall, (${f_{ki,N}}$), for ${v_{ki\alpha}}<0$ and on the left wall, (${f_{ki,1}}$), for ${v_{ki\alpha}}>0$ are solved by the following equations, respectively,

[TABLE]

[TABLE]

III Simulation restlts

III.1 Couette flow

Consider a gas flow between two parallel walls, one at $x=-L$ and the other at $x=L$. The two plates are kept at uniform temperature $T_{0}$ and move with velocity $(0,-v,0)$ and velocity $(0,v,0)$, respectively. Velocity slip becomes more significant with the decrease of the distance between the two plates or with the increase of the mean free path of the molecules, more exactly, with the increase of Knudsen number.

The typical velocity profile between parallel plates in the slip-flow regime is depicted in Fig.2. Only right half of the profile is shown because of its antisymmetry. The gas flow away from the wall can be described by NS equations, and the corresponding flow area is referred to as the NS flow area. The flow near the wall possesses pronounced non-equilibrium characteristics, and the corresponding flow layer is known as the Knudsen layer whose thickness is of the order of the mean free path. In Fig.2, the linear portion $A$-$B$, whose gradient is $dv/dx$, corresponds to the NS flow area and the portion $B$-$D$ corresponds to the Knudsen layer. The line $B$-$C$ is extended from the line $A$-$B$ and the point $C$ is the cross point of the extended line with the right wall.

The slip velocity $v_{slip}$ is defined as the difference between the velocity of the right wall (the value the point $W$) and the velocity value of the point $C$. Considering the complete diffuse reflection, Sone gave the relation between $v_{slip}$ and the mean free path $\lambda$ in Ref.sone2007molecular .

For the hard sphere model,

[TABLE]

For the BGK model,

[TABLE]

The relationship between $\lambda_{HS}$ and $\lambda_{BGK}$ deduced in Sec.II.1 is verified by Eq.(30) and Eq.(31) since $1.2540\lambda_{HS}\approx 1.0162\lambda_{BGK}$.

Knudsen profile $\Delta v$ is defined as the difference between the curves, $B$-$D$ and $B$-$C$. Sonesone2007molecular gave also the relation between $\Delta v$ and $\lambda$,

[TABLE]

by introducing the so-called Knudsen layer function, $Y_{0}(\eta)$, where $\eta$ is a coordinate transformed from $x$ through the following conversion:

[TABLE]

The Knudsen layer function for hard sphere model, $Y_{0}^{HS}(\eta)$, and for BGK model, $Y_{0}^{BGK}(\eta)$, are both shown in Fig.3. The correction of the function $Y_{0}^{HS}(\eta)$ according to the relation in Eq.(9) is also plotted. It can be seen that, the profile of the corrected function is in excellent agreement with the profile $Y_{0}^{BGK}(\eta)$. As a consequence, the Eq.(9) is revalidated. In addition, the results based on hard sphere model can be compared with those from BGK model under same macro conditions by using the relation of Eq.(9).

Considering the Maxwell-type boundary condition, OnishiOnishi1974A gave the expression of slip velocity and Knudsen layer function $Y_{0}^{(\alpha)}(\eta)$ under various TMACs as

[TABLE]

[TABLE]

where

[TABLE]

$A_{i}$ and $k_{s}$ is the coefficient calculated by refined moment methodssone1973kinetic ; Onishi1974A . According to OnishiOnishi1974A , the solutions of $N=7$ are good approximations with high and sufficient accuracy to the exact ones. The coefficients for partial values of $\alpha$ are listed in Table 2. It should be noted that $Y_{0}^{(\alpha)}(\eta)$ is in complete agreement with ${Y_{0}}(\eta)$ shown in Fig.3 when $\alpha=1$ .

The DBM simulation results for the Couette flow with different Knudsen numbers under complete diffusion boundary condition are shown in Fig.4. Results for different Knudsen numbers are obtained by changing the relaxation-time constant $\tau$ according to Eq.(5).

Figure 4(a) shows that the phenomena of velocity slip become more significant with the increase of $Kn$. Comparison of the values of slip velocity normalized by $dv/dx$ between the DBM results and Sone’s results is shown in Fig.4(b). The two kinds of results are in excellent agreement with each other. The DBM accurately capture the velocity slip. Besides, the Knudsen layer is also well described by DBM. As shown in Fig.5, comparison of normalized Knudsen profiles calculated from DBM are also in excellent agreement with Sone s results. The Knudsen profiles $\Delta v$ in Fig.5 are normalized by $\frac{\sqrt{\pi}}{2}{\lambda}\frac{dv}{dx}$.

Taking the TMAC ($\alpha$) into consideration, the Maxwell-type boundary condition is adopted in the following simulation. The DBM simulation results with several different values of $\alpha$ are shown in Fig.6. From Fig.6(a), we can see that the phenomena of velocity slip are more significant with the decrease of $\alpha$. The values of velocity slip for different values of $\alpha$ are compared with those given by Eq.(34). Figure 6(b) shows good agreement of DBM simulation results with those of OnishiOnishi1974A .

The results of $Y_{0}^{\alpha}(\eta)$ for various values of $\alpha$ are compared in Fig.7. The DBM results also show good agreement with those of Eq.(35).

III.2 Poiseuille flow

Pressure driven gas flow, known as Poiseuille flow, in a microchannel is also very common in MEMS. In the slip-flow regime, the Navier-Stokes equations with slip boundary condition are applicable. Accurate second-order slip coefficients are of significantly since they directly determines the accuracy of the results given by Navier-Stokes equationsHadjiconstantinou2003Comment . The first second-order slip model was presented by Cercignani. Using the BGK approximation he obtained

[TABLE]

where $\theta=\mu/p\sqrt{2RT}$. It can be seen, the first-order coefficient is same with Eq.(31). Subsequently, HadjiconstantinouHadjiconstantinou2003Comment improved the model for a hard sphere gas by considering Knudsen layer effects. In his article, the viscosity-based mean free path, $\lambda=\mu/p\sqrt{\frac{\pi RT}{2}}$, was used. Then he obtained the following slip velocity

[TABLE]

However, in our DBM model, viscosity-based mean free path is defined as $\lambda=\mu/p\sqrt{\frac{8RT}{\pi}}$, so the first-order and second-order coefficients should be rescaled by $\pi/4$. The slip velocity formulate should be

[TABLE]

Considering the Maxwell-type boundary, the fully-developed velocity profile can be expressed by

[TABLE]

where $dp/dx$ is the pressure gradient in the streamwise direction. $H$ is the width of the micro-channel. Nondimensionalize the two sides of Eq.(40) by the mean channel velocity $\overline{u}$ gives

[TABLE]

It is clear that the pressure gradient and the viscosity coefficient vanish. Firstly, complete diffuse boundary condition is adopted.

The simulation results by DBM for different Knudsen numbers are shown in Fig.8. It can be found that the velocity slip is significant with the increase of Knudsen number. The nondimensional velocity has a higher maximum value for a smaller Knudsen number. The velocity profile described by Eq.(41) with $\alpha=1$ is also plotted for comparison. The simulation results show good agreement with the expression of Eq.(41).

Considering the behaviours of velocity slip under various TMACs. Then the Maxwell-type boundary condition is used. It can be seen from Fig.9 that the velocity slip is more significant with the decrease of TMAC and the nondimensional velocity has the highest maximum value when the complete diffuse reflection occurs. It concluded that the effect of Knudsen and TMAC on velocity is in the opposite direction. The numerical results are in well agreement with Eq.(41) for different values of $\alpha$ which verify the accuracy of the Maxwell-type boundary condition.

IV Conclusion

A discrete Boltzmann method with Maxwell-type boundary condition for slip flow is presented. The definition of Knudsen number is clarified for DBM. The relation between the Knudsen number based on hard sphere model and that based on BGK model is given. Two kinds of gas flows, including Couette flow and Poiseuille flow, are simulated to verify and validate the new model. The results show that the DBM with Maxwell-type can reasonably capture both the velocity slip and the flow characteristics in Kundsen layer under various Knudsen numbers and tangential momentum accommodation coefficients.

V acknowledge

The authors would like to thank Drs. Wei Jiang, Hongwei Zhu, Wei Wang and Ge Zhang for fruitful discussions on discrete Boltzmann modeling of slip flows. We acknowledge support of National Natural Science Foundation of China [under grant nos. 11475028 and 11772064], Science Challenge Project (under Grant No. JCKY2016212A501 and TZ2016002), and Science Foundation of Laboratory of Computational Physics.