Numerical investigation of curved channel Knudsen pump performance

International Journal of Heat and Mass Transfer 76 (2014) 1–15

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Numerical investigation of curved channel Knudsen pump performance D.M. Bond a,⇑, V. Wheatley a, M. Goldsworthy b a b

School of Mechanical and Mining Engineering, University of Queensland, St Lucia 4072, Australia CSIRO Energy Technology, Newcastle 2304, Australia

a r t i c l e

i n f o

Article history: Received 28 August 2013 Received in revised form 9 April 2014 Accepted 9 April 2014

Keywords: Knudsen pump Thermal creep Thermal transpiration Micro-channel Rarefied gas Unified gas kinetic scheme Validation

a b s t r a c t Thermal creep driven flows in curved micro-channels, commonly referred to as Knudsen pumps, are investigated across a range of rarefaction levels with particular focus on the effects of realistic gas coefficients and geometric configuration on performance. Two base geometries are investigated consisting of a previously proposed curved-straight channel and a newly developed double-curved channel with no straight sections. Use of the S-model kinetic equations enables investigation with realistic values of the Prandtl number and viscosity index for argon and nitrogen as well as for Maxwell molecules. For each gas, the pumping performance and flow structure of each base geometry are investigated for three channel aspect ratios and an array of Knudsen numbers finely spaced between 0.1 and 2.0 to allow performance extrema to be identified. The influence of Prandtl number is found to be significant with increased maximum mass flow rates for argon and nitrogen. The impact of viscosity index is found to be comparatively minor. The double-curved geometry is found to generate over twice the mass flow rate. Both geometries are also tested in pumping array configurations involving multiple sections with the newly introduced geometry again providing superior performance. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Fluidic micro scale devices are currently employed in many areas ranging from sensing duties to fluid metering and delivery. Many of these devices involve the flow of gases (e.g. micro gas chromatography) and so require a means of manipulating very small fluid volumes by some pumping apparatus. These micropumps have typically involved variation in some volume (displacement pumps) or motivated the fluid by dynamic means (centrifugal, electro-osmotic, etc.) [1]. An alternative to these traditional means of pumping was shown experimentally by Knudsen [2] and relies on the gas being in a rarefied state. The level of rarefaction of the gas may be characterized by the ratio of mean free path, k, to some characteristic length and is denoted the Knudsen number, Kn ¼ k=L. The Knudsen number provides an indication of the gas flow regime as shown in Table 1 [3]. Given a gas that is sufficiently rarefied (slip to free molecular) and contained in a pipe with an axial temperature gradient then the gas will proceed to flow in the direction of the gradient (cold to hot) via the mechanism of thermal creep. This is somewhat counter intuitive and a kinetic theory explanation following Sone

⇑ Corresponding author. Tel.: +61 408958505. E-mail addresses: [email protected] (D.M. Bond), [email protected] (V. Wheatley), [email protected] (M. Goldsworthy). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.04.021 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved.

[4] is as follows. From a wall with temperature gradient we take a small segment and estimate the momentum transfer from gas to the wall, see Fig. 1. The molecules impacting the wall arrive at the wall carrying the property of their original location, approximately one mean free path length back along the inverse trajectory. As the molecules arriving from the hotter region have, on average, a higher speed than those from the colder region there is a net tangential momentum transfer to the wall in the opposite direction to the temperature gradient. As a consequence of conservation of momentum the gas experiences a force in the direction of the temperature gradient and so the gas is induced to flow in that direction. If the gas is already moving then some momentum is transferred to the wall in the direction of the flow. When the momentum transfer due to differential temperature and gas flow is equal then a steady flow is established. This effect has been investigated experimentally [2,5–7] and numerically [8,9]. If two reservoirs at different temperatures, T 2 > T 1 , are connected, as shown in Fig. 2 then it can be shown [5] that the pressure P2 will increase to some equilibrium value such that P2 > P1 . To achieve a greater pressure rise than is available with the two reservoir system a series of stages, with the two reservoir system as it’s basic unit, may be employed as shown in Fig. 3. This configuration is the classic Knudsen pump. This method of pumping has the benefit of having no moving parts nor any fluids present that may contaminate the flow.

2

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

Table 1 Gas flow regimes. Continuum flow

Kn < 103

Slip flow

103 6 Kn < 101

Transition flow

101 6 Kn < 10 10 6 Kn

Free molecular flow

gas flow COLD

HOT

diffuse reflection momentum transfer to surface Fig. 1. Schematic showing thermal creep mechanism of a rarefied gas.

Cold P1 , T1

micro-channel

Hot P2 , T2

Fig. 2. Schematic showing thermal creep of a rarefied gas. Initial flow direction is shown through the micro-channel.

1.1. Curved channel geometry In the classic pump the opposing thermal creep in the wider sections of the series is supported by recirculating flow and so a positive net mass flow rate is maintained. A simpler geometry has been proposed by Aoki et al. [10,11] consisting of a channel with alternating curved and straight sections with periodic temperature gradient as shown in Fig. 4. Due to the curved sections generating a stronger thermal creep effect than the straight sections, a consequence of higher temperature gradient on the inside curve, a positive net mass flow rate has been numerically demonstrated [11,12]. The results shown by Aoki et al. are produced using the Bhattnagar–Gross–Krook (BGK) model equation for a monatomic gas, argon, with Prandtl number of unity, as a consequence of the BGK model, and viscosity index also set to unity to reflect the assumption of Maxwellian molecules. These assumptions are physically incorrect for argon which has a theoretical Prandtl number of two-thirds and viscosity index of around 0.81 [13]. Here we perform a numerical investigation of the performance of the curved-straight channel pump with realistic values of the ideal gas variables, namely the Prandtl number and viscosity index, and compare to results obtained with the gas properties used by Aoki

TL

TH

TL

et al. This investigation is carried out for a range of geometric configurations, involving variation in the ratio of channel width to length, as well as for Knudsen numbers ranging from 0.1 to 2. The performance of the pump is characterized by the mass flow rate through a periodic section and the variation in flow structure due to changes in gas properties are also investigated and visualized. The geometry proposed by Aoki et al. generates mass flow due to the imbalance of thermal creep effects in the curved and straight sections. If the geometry were modified even further to increase this imbalance then superior mass flow and pressure rise could be generated. To achieve this objective a design involving no straight sections is herein proposed. The general design of the new pump is shown in Fig. 5. The straight sections of the pump shown in Fig. 4 have been removed and the curved sections altered to introduce width variation along the flow path. This width variation is required to allow mass flow as without it the sections of opposing temperature gradient act in perfect opposition resulting in zero net mass flux. This geometry has a varying Knudsen number, based on the changing channel width, with the greater rarefaction occurring in those sections that drive the flow in the ‘forward’ direction. The zones that support opposing thermal creep are also ‘smooth’ in nature allowing for recirculation with minimal losses due to geometric irregularities. The geometry introduced herein is investigated in a similar manner to the geometry of Aoki et al. with simulation of mass flow rate and flow structure for varying gas properties and geometric configuration. Further to the closed-loop simulations an investigation of pressure rise through the formation of an array of the fundamental unit of each style of pump is carried out. These simulations cover closed arrays of four and eight sections with the average pressure across the channel recorded along the full length of the array. These traces and the total pressure rise obtained are then compared for the two designs presented. This study has been carried out using a two dimensional (2D) method that solves the S-model kinetic equations [14] using the unified gas kinetic scheme [15] which allows the Prandtl number to be set as desired. 1.2. Overview The remainder of the paper now follows with Section 2 giving a brief description of the kinetic equations that are numerically solved according to the method discussed in Section 2.4. The boundary conditions for the numerical method are then covered in Section 2.5. To validate the method three different tests are carried out in Section 3. These tests verify the ability of the method to capture the correct Prandtl number as well as the thermal creep effect. The results for the simulations described above for the curved-straight and double-curved channels are then presented in Section 4 before conclusions are presented in Section 5. 2. Numerical method To describe a monatomic dilute gas a velocity distribution func  tion, f ~ x; ~ n; t , can be used where the total mass per unit volume of

TH

TL

TH

Fig. 3. Schematic showing a typical Knudsen pump. Exterior arrows show thermal creep direction. Interior arrows show gas flow. Recirculating flow is shown in the thicker sections of the series. Dashed lines indicate possible extra units in the series.

3

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15 Table 2 Geometry conditions. Condition

j

Nc

Nt

I II III I–III

0.2 0.5 1 b ¼ 1=2

1600 400 400

24 48 96

Fig. 4. Schematic showing a curved-straight channel Knudsen pump. Interior arrows show direction of thermal creep.

Fig. 7. Diagram of thermal creep cell. General flow pattern shown. Channel widths of w ¼ 1 lm, 100 nm and 20 nm are simulated.

Fig. 5. Schematic showing a Knudsen pump with no straight sections. Exterior arrows show direction of thermal creep. Interior arrows show bulk flow.

n, is given by f d~ n. molecules at a point ~ x and time t, having velocity ~ The Boltzmann equation,

@f ~ @f þn ¼ Xðf Þ @t @~ x

ð1Þ

describes the time evolution of f at point ~ x, in terms of the convection of molecules, ~ n  @@f~x, and the effect of collisions, Xðf Þ. The collision operator can be replaced by the BGK approximation,

Xðf Þ ¼

fT  f

s

;

ð2Þ

which relaxes the known distribution, f, towards a target distribution, f T . Note that f T is a function of the moments of f. The relaxation time, s, determines the rate of relaxation of f towards f T . The standard BGK model [16] uses the Maxwellian equilibrium distribution, Eq. (4), as the target. Whenever the term ‘non-equilibrium’ is used it refers to the deviation of a distribution from the equilibrium Maxwell distribution. In this implementation the S-model [14], Eq. (3), is used as the target distribution, f T ¼ f S . Henceforth we will refer to the S-model distribution as f S which is defined below:

       ~ c ~ c q fS ~ 5 ð5pRT Þ ; x; ~ n; t ¼ f M ~ x; ~ n; t 1 þ ð1  PrÞ~ c ~ RT     ~ ~ q c  c fM ~ exp  ; x; ~ n; t ¼ 2RT ð2pRT Þ3=2

ð3Þ ð4Þ

where,

~ n~ c ¼~ N: The Prandtl number is given by,

Pr ¼

lc p k

;

ð5Þ

where l is the dynamic viscosity, cp is the specific heat at constant pressure and k is the thermal conductivity of the gas. In these equations the mean, or macroscopic, velocity vector is given by ~ q is the heat flux N ¼ ½U; V; W, while ~ c is the thermal velocity and ~ vector. All of the macroscopic properties such as density, q, temperature, T, and pressure, p, are moments of the distribution f. The specific gas constant is given as R. For a monatomic gas Pr is equal to two-thirds. The Chapman–Enskog solution of the BGK model [17] gives a relationship between viscosity and relaxation time of

l ¼ ps;

Fig. 6. Non-dimensional temperature profile for planar Couette flow, Kn = 0.01.

ð6Þ

4

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

Fig. 8. Pressure profiles for the current method and DSMC solutions [23] for a 5lm in length micro-channel.

Fig. 11. Schematic of flow disturbances observed in velocity flow field. Maxwellargon: R ¼ 208 J=kg K; c ¼ 5=3; Pr ¼ 1; x ¼ 1.

Fig. 9. Schematic of curved-straight channel ring section.

where p ¼ qRT is the pressure. The viscosity can be assumed [17] to follow a power law dependence on temperature according to,



l T ¼ T ref lref

x ð7Þ

;

where lref is the reference viscosity at a temperature of T ref . The constant, x, depends on the gas and may be found in reference texts [17,13].

as wall separation or by the length scale of flow property variations [13]. The non-dimensional variables can then be introduced:

^ ~ n ¼~ n=C ref ;

~¼~ N N=C ref ;

^t ¼ t=t ref ;

q^ ¼ q=qref ;

h ¼ T=T ref ;

l^ ¼ l=qref C ref L; ^f ¼ f

2.1. Nondimensionalization pffiffiffiffiffiffiffiffiffiffiffiffiffi We define a characteristic length, L, and speed, C ref ¼ 2RT ref , L and the characteristic time, t ref ¼ C , follows accordingly. A characref

teristic density, qref , and temperature, T ref , are also required. The characteristic length can be defined by a geometric feature such

^

~ x^ ¼ ~ x=L;

.

qref =C 3ref



~ ^ ¼ p=qref C 2ref ; q^ ¼ ~ q=qref C 3ref ; p .  ^ ¼ h=q ; ; g^ ¼ g qref =C 2ref ; h ref

^ give the reduced distribution functions that are The g^ and h introduced in Section 2.3. All of the above quantities from this point on, unless otherwise specified, will be in non-dimensional form and the carats (^ ) will be omitted.

Fig. 10. Non-dimensional mass flow rate for Maxwell-argon (c ¼ 5=3; Pr ¼ 1; x ¼ 1). Error bars on Aoki et al. data shows standard deviation.

5

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

If a nominal reference mean free path is defined according to

Table 3 Gases [13,24].

C ref lref ; pref

ð9Þ

kref Knref ¼ : 2L 2

ð10Þ

Gas

R, J/kg K

c

Pr

x

kref ¼

Maxwell-argon Argon Nitrogen

208 208 297

5=3 5=3 7=5

1 2=3 0.72

1 0.81 0.74

then

w¼

A local channel Knudsen number can also be defined according to this expression for mean free path,

pffiffiffiffiffi 2w hL KnL ¼ ; wqL

ð11Þ

where w is the non-dimensional width of the channel and the subscript L denotes a local variable. 2.3. Reduced distribution functions The numerical method is applied to 2D simulations in this paper. To reduce the demand on computer resources, both storage and computation, the single distribution function is split according to the procedure outlined by Chu [19]. By integrating out the w velocity dependence, Fig. 12. Schematic of curved-straight channel pump section.

g ðx; u; v ; tÞ ¼ hðx; u; v ; t Þ ¼

2.2. Relaxation time

1

Z1 1

  f x; ~ n; t dw;   w2 f x; ~ n; t dw;

1

To characterize the condition at which a simulation is carried out a non-dimensional constant, w ¼ lref =qref C ref L, is introduced. From Eqs. (1), (2), (6) and (7) we have,

the three dimensional equation (1) reduces to two simultaneous 2D equations and the equilibrium distributions become;

"

 # ~ q 2~ 4 c ~ c ~ c g ¼ g 1 þ ð1  PrÞ 2 4 ; 5 h qh "  # ~ q 2~ gM h 4 c ~ c ~ c S 1 þ ð1  PrÞ 2 2 ; h ¼ 2 5 h qh   ~ ~ c  c q : exp  gM ¼ h ph S

 @f ~ @f qh1x T þn ¼ f f : @t 2w @~ x

Z

ð8Þ

This w gives the reference non-dimensional viscosity and is a useful method of relating the relaxation time of the BGK model to a desired viscosity [18]. In this way a gas with arbitrary viscosity coefficient may be simulated.

M

ð12Þ ð13Þ ð14Þ

Fig. 13. Non-dimensional mass flow rate in the curved-straight channel for Maxwell-argon (c ¼ 5=3; Pr ¼ 1; x ¼ 1), argon (c ¼ 5=3; Pr ¼ 2=3; x ¼ 0:81) and nitrogen (c ¼ 7=5; Pr ¼ 0:72; x ¼ 0:74). Fitted curves follow the form ða exp ðbKnÞ þ c exp ðdKnÞÞ, see Table A.5 for coefficient values.

6

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

Fig. 14. Non-dimensional velocity with streamlines for Kn ¼ 1:0; j ¼ 1. For direction of flow see Fig. 9. Note: additional truncated streamlines added in areas of re-circulation to aid in flow visualization.

Fig. 15. Non-dimensional velocity with streamlines for Kn ¼ 0:5; j ¼ 1. For direction of flow see Fig. 9. Note: additional truncated streamlines added in areas of re-circulation to aid in flow visualization.

2.4. Unified gas kinetic scheme From this point on, it should be noted that all vector quantities refer to only two spatial variables i.e. ~ n ¼ ½u; v  and ~ x ¼ ½x; y. To efficiently approximate gases with internal degrees of freedom the method of Nie et al. [20] is implemented. This approach approximates variable internal degrees of freedom through applying a multiplier, K, to the energy contribution from h. The value of K may be calculated according to

K¼

4  2c ; c1

where c is the ratio of specific heats. Further details on this approach are outlined by Nie et al. [20] and Bond et al. [21]. At this point we have a means of calculating the flow variables as functions of spatial variation. To evolve these equations numerically the unified gas kinetic scheme [15] is implemented.

The unified gas kinetic scheme, henceforth referred to as UGKS, is a cell centered finite volume method with velocity space discretised according to the choice of numerical quadrature. Therefore, in each cell, the average conserved macroscopic properties are tracked along with values of g and h sampled at discrete velocity abscissa, ~ na . The number, d, and location of these abscissa being dependent on the quadrature scheme used. The macroscopic prop  erties in the cell can be calculated from these g a ¼ wa g ~ na and   ~ na ha ¼ wa h na where wa is the quadrature weight for abscissa ~ as follows;

q¼

d X ga ;

ð15Þ

a¼1

qNi ¼

d X na;i g a ; a¼1

ð16Þ

7

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

j ¼ 1; b ¼ 1=2. Note: This mesh is for illustrative purposes

Fig. 18. Mesh layout. only.

The time variation of the discrete distribution values, and by extension the macroscopic properties, in each cell are determined by the transport and interaction of the particles described by the discretised distribution functions. The updated discrete distributions are given by

g nþ1 ¼ g na þ a

1

"Z

U

t nþ1

tn

Z X um g^m ðtÞDSm dt þ

t nþ1

tn

m

#

Z

Xðg Þdt ;

ð19Þ

U

where m denotes the mth piecewise linear boundary of control volume with area U. The normal molecular velocity to the boundary is given by um and g^m is the distribution function found at the boundary. Boundary length is given by DSm and n denotes the current time step. This equation essentially describes the processes of transport and molecular collisions according to the Boltzmann equation, Eq. (1), and also applies to ha . The typical approach to solving Eq. (19) is to solve for the contribution to the cell average value of g a assuming free streaming of the interface distribution and then calculate the effect of the collision operator on the cell average value. The UGKS on the other hand calculates the form of the interface distribution, g^m , over the period of the time step taking into account both the free streaming and particle collisions. The interface distribution is then given by,

g^a ¼

    1 g Sa x0 ; t0 ; ~ na exp  ðt  t0 Þ dt 0 s tn s     1 na ðt  tn Þ; tn ; ~ þ exp  ðt  t 0 Þ g 0a x  ~ na ; 1

Z

t

ð20Þ

s

Fig. 16. Non-dimensional velocity with streamlines for Kn ¼ 0:1; j ¼ 1. For direction of flow see Fig. 9. Note: additional truncated streamlines added in areas of re-circulation to aid in flow visualization.

na ðt  t 0 Þ is the where x is a location on the cell interface, x0 ¼ x  ~ particle trajectory and g 0a is the initial gas distribution of g a at time tn . ^ a , multiplied by the ^a and h By integration (summation) of g correct terms according to Eqs. (15)–(18), over the time step Dt ¼ t nþ1  t n it is possible to calculate the flux of conserved properties across the interface and so update the cell average macroscopic quantities. The updated discrete distributions, g nþ1 , can be a calculated in a like manner. This allows the calculation of an Sðnþ1Þ updated equilibrium distribution g a and so the final updated discrete distribution is given by,

g nþ1 ¼ g na þ a

Dt þ 2

1

U



"Z

t nþ1

tn

1

snþ1

X um g^m ðtÞDSm dt m

g aSðnþ1Þ

 1 SðnÞ   g nþ1 þ n g a  g na a

s

 ;

ð21Þ

where snþ1 is calculated according to the updated cell average macnþ1 roscopic properties. This equation is also applied to calculate ha . Fig. 17. Schematic of double-curved channel.

h¼

 P    ~~ na  ~ na g a þ Jha  qN ðc  1Þ 12 da¼1 ~ N

qi ¼

q d

 1X na;i  Na;i ð~ ca  ~ ca g a þ Jha Þ: 2 a¼1

;

2.5. Boundary conditions

ð17Þ ð18Þ

The physics present in gas–surface interactions, under rarefied conditions, provides the mechanism for the Knudsen pump to operate. The approach used to model these interactions is therefore of critical importance to the effectiveness of the numerical simulation. In this work the fully diffuse wall model is used for consis-

8

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

Fig. 19. Non-dimensional mass flow rate in the double-curved channel for Maxwell-argon (c ¼ 5=3; Pr ¼ 1; x ¼ 1), argon (c ¼ 5=3; Pr ¼ 2=3; x ¼ 0:81) and nitrogen (c ¼ 7=5; Pr ¼ 0:72; x ¼ 0:74). Fitted curves follow the form ða exp ðbKnÞ þ c exp ðdKnÞÞ, see Table B.6 for coefficient values.

tency with the work of Aoki et al. In this model the surface is assigned a temperature, T w and a velocity, ~ Nw . Any molecule that interacts with the surface is then assumed to be fully accommodated and re-emitted into the flow with a random velocity chosen from the Maxwellian equilibrium distribution, Eq. (14), at the wall temperature and velocity. Assuming the Maxwellian density is set at unity, the discrete distributions being re-introduced to the flow (~ na  ~ n > 0) are given by,

when the slope of the residual on a log–log plot closely approached zero. The residual is calculated according to,

  ~ g a ¼ bg M a 1; Nw ; T w ;  h  ~ ha ¼ b g M a 1; Nw ; T w ; 2

~ gives the difference in M ~ over one time step. where DM

ð22Þ ð23Þ

where ~ n is the normal vector to the surface (directed into the flow). The value of b is set to ensure no mass flux across the boundary and so can be calculated according to,

P

~ ~ n   : ~ M ~ ~ ~ n>0 g a 1; Nw ; T w na  n na ~

b ¼ P

~ na ~ n<0 g a na

ð24Þ

3. Validation To validate the UGKS method three tests are carried out. The first test, planar Couette flow with comparison to an analytical solution of the Navier–Stokes equations with velocity slip and temperature jump walls by Sofonea et al. [22], is used. This test is intended to demonstrate the validity of the approximation of internal degrees of freedom. The second test is a comparison with DSMC of the pressure gradient in an enclosed chamber due to thermal creep effects. The third test is a series of simulations of curved channel thermal creep flow using Pr = 1 and l / T to allow a direct comparison with the numerical results of Aoki et al. [11].

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 P ~ ~ DM  DM Nc Nt R¼ ; P ~ M ~ ¼ ½q; U; V; 1=h; M

ð25Þ ð26Þ

3.2. Planar Couette flow To ensure that the numerical method is able to correctly simulate gas flows, with varying coefficients of specific heat ratio and Prandtl number, planar Couette flow is simulated and validated against an analytic solution [22]. The test case consisted of a two-dimensional domain constrained by two infinite walls moving in opposite directions, at a combined differential speed of U ¼ 0:1  C ref , with the walls held at the same initial temperature as the gas. The Knudsen number is given as Kn = 0.01 with the reference length being the channel width. The relatively low Knudsen number used in this test is due to the use of the Navier–Stokes equations in the analytical solution which are limited to low rarefaction flows even when implemented with specialized velocity slip and temperature jump boundary conditions. The numerical

 domain consisted of a grid with cell count of N x ¼ 2; N y ¼ 128 , fully diffuse walls at ymax and ymin , and periodic conditions in the x direction. The temperature profiles from these simulations are shown in Fig. 6. The simulation closely follows the analytical solution obtained via the Navier–Stokes equations with modified boundary conditions and thereby provides confidence that the numerical method is accurate. 3.3. Thermal creep cell

3.1. Convergence to steady state For all simulations carried out in this paper steady state solutions are presented. Steady state is assumed when all components of the global residual vector are reduced by a factor of 105 or

To verify that the thermal creep effect is recovered by the given system of equations a sealed 2D micro-channel with rectangular cross section, as shown in Fig. 7, is simulated and compared to results by the DSMC method according to Masters and Ye [23].

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

Fig. 20. Non-dimensional velocity with streamlines for Kn ¼ 1:0; direction of flow see Fig. 9.

j ¼ 1. For

The two ends of the channel are maintained at constant temperatures of T L ¼ 273 K and T H ¼ 573 K. The side wall temperatures varied linearly according to T ðx; yÞ ¼ ðT H  T L Þx=L þ T L . The pressure is initialised uniformly throughout the cell at P ¼ 101:325 kPa. The walls of the channel are modeled as fully diffusely reflecting. The gas being simulated is argon with a reference viscosity of lref ¼ 2:117  105 N s=m2 and viscosity index of x ¼ 0:81 [13]. Three channel widths of w ¼ 1 lm, 100 nm and 20 nm are simulated to steady state with Knudsen numbers at the low temperature end, according to Eq. (11), of 0.07, 0.7 and 3.5 respectively. The pressure variation along the center-line of the channel is then calculated and compared to a DSMC solution as shown in Fig. 8. From Fig. 8 it is clear that the current method recovers the pressure variation in the micro-channel closely with the trace lying almost entirely within the scatter of the DSMC solution for all cases simulated. This result indicates that the numerical method is able to accurately capture the thermal creep effect in micro-channels. 3.4. Curved-straight channel validation A significantly simplified version of the curved channel Knudsen pump geometry was presented by Aoki et al. [11] with accompanying mass flow rate data. This geometry consists of a

Fig. 21. Non-dimensional velocity with streamlines for Kn ¼ 0:3; direction of flow see Fig. 9.

9

j ¼ 1. For

half-segment of the repeated element used in the pump which is then mirrored to form a ring shaped channel as shown in Fig. 9(a). To evaluate this ring only the half-segment, seen in Fig. 9(b), needs to be simulated with periodic boundary conditions connecting E1 and E2 . To accommodate the rotation that this periodicity requires the fluxes undergo a reflection about the x–y line according to,

FðE1 ; u; v Þ ¼ FðE2 ; u; v Þ: The same geometric parameters are used as in [11] and these are given as follows;

D R

j¼ ; b¼

LS : LS þ pR

ð27Þ

The physical space flow domain is therefore fully defined according to D; j; b and the spatial discretization. According to Aoki et al. the value of b is a constant of b ¼ 1=2. The grid used to define the physical space is defined according to the number of cells along the center line of the domain, N c , and across the domain, N t . The spacing of the grid is such that the spacing along the walls is constant for each section of curvature (straight, small radius, large radius) and constant normal to the walls. Three geometry conditions are simulated for the curved-straight channel and

10

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

_ is calculated for Maxwell-argon and can be seen The value of m _ is found to be below in Fig. 10 where the standard deviation of m 2  104 for all cases. _ with respect to the reference soluThe percentage change in m tion for the data points calculated in Fig. 10 at Kn = 1.0 is 4.7%, 4.8% and 7.1% for j from 1 to 0.2. Overall the results calculated using the current method match closely with those generated for the previous work and so there is confidence in the current approach. As Knudsen number increases the flow structure is disturbed with perturbations in the streamlines observed along lines sketched in Fig. 11. These disturbances are also present in the visualizations presented by Aoki et al. [11] and are believed to provide a mechanism _ at high Kn. This divergence for the divergence of the computed m can be attributed to the inability of the discrete ordinate method of solution (that both methods employ) to fully represent the velocity space as Knudsen number increases. This results in the propagation of discontinuities throughout the gas that are directly related to the spacing of the velocity space discretization. As the Knudsen number increases the effect of relaxation is outweighed by that of free transport and so the discrete velocities propagate in semi-discrete waves away from any source of discontinuity in velocity space, such as the inner radius of the channel. For low Knudsen numbers this is not an issue and the effect is not observed in the flow for Kn 6 0:5. However, regardless of numerical artifacts, the results shown in Fig. 10 show reasonable agreement with the results of Aoki et al. and indicate that we are able to reliably reproduce similar simulation results in the curved-straight channel.

4. Results

Fig. 22. Non-dimensional velocity with streamlines for Kn ¼ 0:1; direction of flow see Fig. 5.

j ¼ 1. For

are listed in Table 2. These are chosen to correspond with those used by Aoki et al. [11]. Note that the width of the channel, D, is defined according to Eq. (11). Also following Aoki et al. the wall temperatures are set to T L ¼ 300 K and T H ¼ 900 K while the Prandtl number and viscosity index are both set to unity to emulate Maxwellian molecules. The ratio of specific heats is set to reflect a monatomic gas such that c ¼ 5=3. To characterize the flow around the channel the non-dimensional mass flow rate across section Xi is given by

_ Xi ¼ m

~ ~ qN ni Xi pffiffiffiffiffiffiffiffiffiffi ffi qav 2RT L D

ð28Þ

and the average around the channel by,

_ ¼ m

Nc X _X m i

Nc

i

;

ð29Þ

ni is the normal to the i-th section and qav is the average denwhere ~ sity in the channel. Following Aoki et al. numerical fluctuation along the channel is accounted for by computing this quantity for each column in the computational grid, which spans the channel in real space, and the average is taken as the characterizing value.

The results presented in this section cover the effect of realistic gas coefficients and the comparison of the two different geometric configuration of Knudsen pumps as discussed in Section 1. The effects of gas properties and comparisons of pump performance are primarily presented in terms of the mass flow rate through a single periodic pumping section. The secondary presentation covers the pressure rise capability of the different geometries for the peak mass flow rate condition found in the preceding discussion. The pertinent details of geometry are presented along with the results while the operating conditions used are presented in the next section. 4.1. Operating conditions The gases simulated are Maxwell-argon, argon and nitrogen with reference viscosity determined by the relation given for channel Knudsen number in Eq. (11) such that w ¼ Kn=2. The specific gas constant, R, ratio of specific heats, c, Prandtl number [24] and viscosity index x [25] are given in Table 3. The reference length is set to L ¼ 1 lm and as in Section 3.4 the low wall temperature is set to T L ¼ 300 K and T H ¼ 900 K. The velocity discretization used in all cases is a product Gauss– Hermite type quadrature scheme, symmetric about zero, with a total of 784 discrete velocities. This givespaffiffiffi square velocity space with nmax 6 ~ na 6 nmax where nmax  132 R and R is determined by the gas. This velocity lattice is chosen as it is the maximum allowable given the computational resources. 4.2. Curved-straight channel: mass flow The full channel of Aoki et al., visualized in Fig. 9(a), can be split into identical sections and then ‘chained’ to reproduce a full pump with characteristics as desired. One of these sections is shown in

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

11

Fig. 23. Dimensional mass flow rates for the curved-straight and double-curved geometries.

Fig. 12. For periodic flow the boundaries labeled E1 and E2 are connected giving a continuous loop. This geometry is closely related to that shown in Section 3.4 and so the same values for j and b, shown in Table 2, also apply. The computational mesh used to discretize the simulation domain also retains the same number of cells as previously defined for the ring configuration due to computational resource limitations. The non-dimensional mass flow rate of the curved-straight channel is calculated for the three gas types and geometry configurations according to Eq. (29). The results are shown in Fig. 13 with flow visualization shown in Figs. 14–16. Coefficients for the fitted lines shown in Fig. 13 are provided in Table A.5 to allow simple comparison of alternative results to the current data. Note that the flow visualization displays only a portion of the full simulation domain, according to Fig. 9(b), as the flow structure in the non-displayed section is essentially identical to that shown. The variation in flow structure due to changes in geometry from Section 3.4 can be observed when comparing Fig. 11 to the equivalent flow in Fig. 14(a). _ given in Fig. 13 demonstrate that real gas The profiles for m coefficients have a noticeable impact on the simulated perfor_ mance of the channel mass flow rate. The change in profile of m is, however, proportional to channel width. Comparing argon and nitrogen to Maxwell-argon for the thickest channel (j ¼ 1), there _ of 6% and 5% for argon and nitrois a difference in the fitted peak m gen respectively. The Kn values at which these maxima occur vary with Knudsen numbers of 0.43, 0.38 and 0.42 for Maxwell-argon, argon and nitrogen respectively. Using the fitted curves it is possible to also calculate approximate maxima values and locations for j ¼ 0:5. In this case all gases are within 2% of each other with peak location of 1:25  0:05 Kn. The fact that there exists a maximum in the mass flow rate, as opposed to a monotonic increase as suggested for plain channels by Sharipov [26] and Rojas-Cardenas et al. [7], can be explained through the presence of complicated flow structures that are not found in the plain channel experiments. The recirculation zones

that are present, see Figs. 14–16, vary in strength and size with Knudsen number and so influence the mass flow rate in a non-linear fashion. The influence of Prandtl number and specific heat ratio, c, can also be investigated. Maxwell-argon (c ¼ 5=3; Pr ¼ 1; x ¼ 1) and argon (c ¼ 5=3; Pr ¼ 2=3; x ¼ 0:81) both have the same specific heat ratio while they vary in peak mass flow rate by 6% for the j ¼ 1 case. In comparison the difference between argon and nitrogen (c ¼ 7=5; Pr ¼ 0:72; x ¼ 0:74) is marginal at 0.5%. Comparing the coefficients for these three gases there is a significant difference in peak mass flow rate when the difference in Prandtl number is (relatively) large and little impact with changes in specific heat ratio. This significance of Prandtl number is to be expected as the flow is entirely due to thermal effects in a viscous medium and the Prandtl number characterizes the relative impact of these two factors on the flow field. Therefore, as the Prandtl number is decreased, it is expected that thermal effects will dominate and mass flow rate will correspondingly increase, as demonstrated in the simulation results. A comparison of the flow structures visualized in Figs. 14–16 show a number of interesting features. The Maxwell-argon gas typically shows a greater degree of recirculation which can especially be seen in Fig. 15(a) in the lower curved section of channel. The presence of numerical disturbances in the flow, discussed earlier in Section 3.4 and Fig. 11, are also present in Fig. 14. Overall, however, the flow structure variation between the different gas types is minimal. Across Knudsen number the flow structure is seen to progress from highly recirculatory for low rarefaction to highly uniform for high rarefaction. The flow velocities also vary with Knudsen number as indicated by the mass flow rates shown in Fig. 13. 4.3. Double-curved channel: mass flow In a similar manner to the simplification performed in Section 4.2, to reduce the full cascade shown in Fig. 5 to a single section,

12

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

to equal the area of the full loop curved-straight geometry as shown in Fig. 9(a). These constraints lead to the following expressions for a and g,

g¼ a¼

Fig. 24. Contours of absolute velocity for the two geometries considered at Kn ¼ 0:35 and j ¼ 1 with four sections. Note that each pump is not shown to the same scale.

the pump is simplified to the configuration shown in Fig. 17. The widest and narrowest channel widths are given as Da and Db respectively while the intermediate radius is given by R. To fully define the geometry the narrowest channel width, Da , and the ratios a and g are specified according to

a¼

Da ; R

g¼

Da : Db

ð30Þ

The wall temperature is set to vary linearly from T L to T H with lines of highest and lowest temperature located at the center of the smaller radius curve as shown in Fig. 17. This configuration allows a more physically reasonable temperature distribution when compared to the channel described in Section 4.2. To aid in comparison with the geometry described in Section 4.2 the parameter Da is constrained to be the same value as D, such that Da ¼ D. This means the maximum Knudsen number of the channels are also equal. The minimum radius of the ring channel is also set equal to the inner radius of the curved-straight channel. The planar area of the double-curved geometry is also constrained

2ðb  1Þ þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 b2 ðj  2Þ þ ðj þ 2Þ  2bð4 þ j2 Þ bðj  4Þ  j  4

4gj : 4g þ jð1  gÞ

;

ð31Þ ð32Þ

The computational grid is defined in a like manner to the curved-straight channel with the cell distribution defined according to Table 2 and shown schematically in Fig. 18. The double-curved channel is simulated using the same condi_ is tions described previously in Sections 4.2 and 3.4. The value of m calculated for a range of Knudsen numbers and can be seen in Fig. 19 with curve fits according to Table B.6 to allow easy comparison of alternative results with the provided data. Note that the value of D in Eq. (28) is taken as the minimum channel width i.e. D ¼ Da . The effect of Prandtl number is clearly shown in Fig. 19 with peak mass flow rates for argon and nitrogen being at least 10% higher than for Maxwell-argon when j ¼ 0:5 and by at least 16.5% greater when j ¼ 1. The effect of x and c is seen to be negligible with differences in peak mass flow rate of less than 2.4% when comparing argon and nitrogen. The peak mass flow rates are found to occur at lower Knudsen number values as the channel widens from j ¼ 0:2 to 1. The peak values for j ¼ 0:2 are Kn ¼ 0:74  0:01 for argon and nitrogen and Kn = 1.19 for Maxwell-argon. For j ¼ 0:5 and 1 the peak values are at Kn ¼ 0:48  0:11 and Kn ¼ 0:3  0:003 respectively for all gases. Note that these values are the expected peak locations calculated from the curve fits given in Table B.6. The flow is also visualized and can be seen in Figs. 20–22. The flow patterns shown in these figures do not display the clear differences that are seen for the curved-straight channel. The key features of the flow include two large recirculation zones occupying the outer regions of the curved channel along with relatively uniform flow through the neck of the channel. The primary flow in the channel for all gases is observed to concentrate as the Knudsen number increases. The tortuosity of this primary flow is also observed to increase with decreasing Knudsen number. Recirculation zones attached to the walls near the transition from small to large radius reduce in strength as Knudsen number decreases until _ The recirculation is they are entirely absent at the maximum m. then found to relocate and strengthen as the Knudsen number decreases further. It can be seen that at the maximum mass flow rate the flow is highly uniform through the narrow section of channel with only the primary recirculation zones in evidence. 4.4. Flow rate comparison Comparing the non-dimensional mass flow rate of the doublecurved channel to the curved-straight channel the first noticeable feature is the significant difference in peak mass flow rate of all configurations. Taking the maximum, fitted, value of argon we find that the double-curved channel has a value 2.6 times that of the curved-straight channel. This marked increase in mass flow may be attributed primarily to the change in channel layout. Similar increases in flow rate are observed across the full range of gas and geometry configurations. While significant effort has been expended on matching geometric features to allow direct comparison of the non-dimensional mass flow rates (see Section 4.3) the dimensional flow rates have also been calculated to improve confidence in the analysis. The average dimensional mass flow rate through the pump is calculated according to

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

Fig. 25. Pressure rise for argon in a four and eight segment pump arrangement at Kn ¼ 0:35 and

No. Sec.

Geometry

^1 ^¼p ^2  p Dp

^1 ^ 2 =p p

4

CS DC

0.128 0.165

1.077 1.088

8

CS DC

0.135 0.218

1.081 1.118

Nc ~ ~ X qN ni Xi i

Nc

;

j ¼ 1. Note that bL p is the center-line length of a single pump segment.

The average pressure across the channel is measured for each span-wise column of cells according to,

Table 4 Pressure rise (CS = curved-straight, DC = double-curved).

_ ¼ m

13

ð33Þ

where the quantities are all dimensional and Xi is the width of section i along a column of cells as described in Section 3.4. To relate the dimensional and non-dimensional quantities the relation

pffiffiffiffiffiffiffiffiffiffiffi  ^_ ¼ m= _ qav 2RT L D is given with the average density given by m

qav ¼ p=RT L , a pressure of p = 1 atm and temperature of T L ¼ 300 K. Following Eq. (11) the channel width, D, is defined 3 according to D=L ¼ ðT L =T ref Þ2 with reference length, L ¼ 106 m, and temperature, T ref ¼ 273 K. The dimensional mass flow rates, in 2D, are therefore shown in Fig. 23. When comparing the location of the maximum flow rates for j ¼ 1 the maxima for the double-curved channel are also somewhat closer to the continuum regime, Kn  0:3 as opposed to Kn  0:4, however the maxima for lower values of j are found at markedly lower values of Kn. Overall, the doublecurved channel is shown to possess much higher mass flow capacity. This capacity is also found at lower Knudsen numbers, especially as the average channel width decreases, although the sensitivity to Knudsen number is increased with peak values found within a much narrower band than found with the curved-straight geometry. 4.5. Pressure rise The pressure rise that could be obtained by joining multiple instances of the modeled geometries is also investigated. This is performed for a condition close to that for peak mass flow rate for both geometries calculated above. Four and eight full sections are joined together with the ends blocked off as shown for the four section configuration in Fig. 24.

PNc pi j p ¼ ; pref Nc qav RT L

ð34Þ

where Nc is the number of cells spanning the channel. The average pressure for each span-wise column is then plotted against the corresponding normalized center-line distance as shown in Fig. 25. The performance characteristics are noted explicitly in Table 4. The traces shown in Fig. 25 shows clear evidence that the double-curved layout has greater pumping performance, at the conditions specified, than the curved-straight channel. In both four and eight segment configurations the slope of the least squares fit for the double-curved channel is at least twice that of the curvedstraight channel. The percentage pressure rise is also greater by one percent for the four section configuration (CS 8%, DC 9%) and by four percent for the eight section configuration (CS 8%, DC 12%). These figures show that the pressure rise induced by adding further sections to the pump is greater for the doublecurved channel although both geometries do not show the linear gains originally shown by Aoki et al. These non-linearities can be partly assigned to the change in effective Knudsen number along the channel as the pressure increases. The difference to the results of Aoki et al. is believed to be due to the use of argon, not Maxwellian argon, as the test gas which introduces significant changes in the balance of thermal and momentum effects. This non-linearity is evident when viewing the results in Fig. 25 where the general trend in pressure rise, for both geometries, appears to ‘flatten out’ slightly in the mid sections of the pressure trace with the main pressure gains at either end of the pump. This effect is most apparent for the curved-straight geometry for eight sections seen in Fig. 25. The double-curved geometry does, however, show somewhat less of this effect and so demonstrates less sensitivity to the variation in operating conditions along the length of the channel. 4.6. Discussion The results show the curved-straight geometry proposed by Aoki et al. induces flow in a rarefied gas. The double-curved geometry proposed here improves on that design by providing two and a half times the peak non-dimensional mass flow rate of the curved-

14

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15

straight channel under similar operating conditions (argon, j ¼ 1; Kn  0:3). When simulating an actual pumping configuration the pressure rise is also found to be greater for the doublecurved pump for both four and eight section pumps under the simulated conditions. The double-curved pump also shows less sensitivity to variation in operating conditions along the length of the pump. With the objective of applying these geometries to real world applications the double-curved geometry uses a temperature distribution that is straightforward to apply. This is in contrast to the curved-straight geometry which requires a complex temperature gradient to be applied along each wall. Overall the double-curved geometry appears to be the superior pump in terms of mass flow rate and pressure rise capability for the conditions presented in this paper. To implement this technology in a real world situation further optimization for the particular use case would have to be carried out. As noted by Aoki et al. the temperature gradients used in this study (T H =T L ¼ 3 over distances of order 106 m) are quite severe and would most likely need to be reduced leading to lower mass flow rates. Increasing the pressure rise capability of the pumping array through decreasing channel width to maintain a constant Knudsen number throughout the array is one possibility of ameliorating this effect. Further modifications of the geometry as well as the relative positioning of the spatially varying temperature field could also be investigated. Axisymmetric, or fully three dimensional, simulations would also be required for pipe flow cases were the two dimensional case no longer holds. As the flow is driven by the thermal creep effect at the interface between gas and solid it is this area that should also be investigated further. The boundary interaction is of vital importance to the performance of these devices and more sophisticated numerical models should be implemented. Implementations such as the Cercignani–Lampis reflection kernel [27] with varied accommodation coefficients as well as mass adsorption and diffusion at the surface could be used. Multi-species gas effects are also of critical importance if air is to be used in situations requiring micro-scale pumping under atmospheric conditions. As always the final verdict on application of this technology must come from experimental validation and so practical application in an experimental setting is called for.

Knudsen number at which this maxima occurs. The viscosity index and specific heat ratio, on the other hand, does not have as great an impact. Variation in the aspect ratio of the channel shows the maximum flow rate is found to occur in a wide channel with Knudsen number in the slip to transitional regime. For thinner channels a lower peak mass flow is observed at later Knudsen numbers. A pumping configuration is also tested using four and eight fundamental elements with argon as the test gas and pressure rises of 7.7% and 8.8% are generated respectively. The geometry introduced in this paper, consisting of a width varying, continuously curved, micro-channel is also simulated for varying gas type and aspect ratio. The same conclusions about the influence of gas properties are made for the double-curved channel as for the curved-straight channel. The mass flow rates for comparable geometric configuration are found to exceed the mass flow rates for the curved-straight channel by two and a half times. The percentage pressure increase when operating as a pump is also found to be greater with increases of 8.8% and 11.8% for four and eight section pumps respectively. Overall these simulations demonstrated the ability of a curved channel to support a non-negligible net mass flow rate and pressure rise with no moving parts through the application of a tangential temperature gradient. The geometry introduced is shown to be suitable for applications requiring mass transport and pressure increase with a favorable comparison to a previously proposed geometric configuration. The significant changes in characterizing mass flow rate due to the addition of accurate Prandtl number also indicates that simulations allowing physical insight into the processes occurring in micro-channels should be carried out with methods that allow flexible Prandtl number. It is also noted that this geometry with no open boundaries, and associated conditions, can be considered a highly useful benchmark for the comparison of numerical methods in the slip and transitional regime.

5. Concluding remarks

Appendix A. Curved-straight channel mass flow fit coefficients

The geometry introduced by Aoki et al., consisting of a microchannel made up of a curved and straight section, is revisited with simulations carried out using real gas values for Prandtl number and viscosity index. The introduction of variable Prandtl number is found to be a significant factor in increasing the maximum mass flow rate of the curved-straight channel and also affecting the

Acknowledgement The authors would like to acknowledge Dr. M. Macrossan for helpful discussions related to the content of this paper.

See Table A.5.

Appendix B. Double-curved channel mass flow fit coefficients See Table B.6.

Table A.5 Curved-straight channel fit coefficients ða exp ðbKnÞ þ c exp ðdKnÞÞ. Gas

j

a

b

c

d

Maxwell

0.2 0.5 1

705=919 3=772 3=557

162=649 977=529 5518=731

293=382 2=561 1=226

112=451 71=253 211=495

Argon

0.2 0.5 1

649=876 3=794 3=496

118=509 943=481 7769=893

20=27 1=289 1=212

211=915 113=397 139=279

Nitrogen

0.2 0.5 1

717=995 4=965 4=685

123=638 87=52 3241=426

116=161 1=260 3=634

78=407 185=603 326=709

15

D.M. Bond et al. / International Journal of Heat and Mass Transfer 76 (2014) 1–15 Table B.6 Double-curved channel fit coefficients. ða exp ðbKnÞ þ c exp ðdKnÞÞ. Gas

j

a

b

c

d

Maxwell

0.2 0.5 1.0

1=234 3=358 1=105

2107=768 154=659 291=712

2=399 3=391 9=803

97=954 4149=751 11929=983

Argon

0.2 0.5 1.0

3=614 4=415 9=790

51=505 352=995 577=999

4=899 9=926 2498=399

5408=991 596=75 40727=539

Nitrogen

0.2 0.5 1.0

4=803 8=845 11=991

92=731 149=478 359=702

4=855 4=431 3=122

3853=783 5459=812 11164=589

References [1] D.J. Laser, J.G. Santiago, A review of micropumps, J. Micromech. Microeng. 14 (2004) R35. [2] M. Knudsen, Eine revision der gleichgewichtsbedingung der gase. thermische molekularstrmung, Ann. Phys. 336 (1909) 205–229. [3] A.A. Rostami, A.S. Mujumdar, N. Saniei, Flow and heat transfer for gas flowing in microchannels: a review, Heat Mass Transfer 38 (2002) 359–367, http:// dx.doi.org/10.1007/s002310100247. [4] Y. Sone, Kinetic Theory and Fluid Dynamics, Birkhäuser Boston Inc., 2002. [5] O. Reynolds, On certain dimensional properties of matter in the gaseous state, Proc. R. Soc. London 28 (1878) 303–321. [6] S. McNamara, Y. Gianchandani, On-chip vacuum generated by a micromachined knudsen pump, J. Microelectromech. Syst. 14 (2005) 741–746. [7] M. Rojas-Cardenas, I. Graur, P. Perrier, J.G. Meolans, Time-dependent experimental analysis of a thermal transpiration rarefied gas flow, Phys. Fluids 25 (2013) 072001. [8] F. Sharipov, G. Bertoldo, Poiseuille flow and thermal creep based on the Boltzmann equation with the Lennard–Jones potential over a wide range of the Knudsen number, Phys. Fluids 21 (2009) 067101. [9] V. Titarev, Rarefied gas flow in a planar channel caused by arbitrary pressure and temperature drops, Int. J. Heat Mass Transfer 55 (2012) 5916–5930. [10] K. Aoki, P. Degond, L. Mieussens, S. Takata, H. Yoshida, A diffusion model for rarefied flows in curved channels, Multiscale Model. Simul. 6 (2008) 1281– 1316. [11] K. Aoki, P. Degond, L. Mieussens, Numerical simulations of rarefied gases in curved channels: thermal creep, circulating flow, and pumping effect, Commun. Comput. Phys. 6 (2009) 919–954. [12] V. Leontidis, J.J. Brandner, L. Baldas, S. Colin, Numerical analysis of thermal creep flow in curved channels for designing a prototype of Knudsen micropump, J. Phys.: Conf. Ser. 362 (2012) 012004. [13] G.A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford University Press, 1994.

[14] E.M. Shakhov, Generalization of the Krook kinetic relaxation equation, Fluid Dyn. 3 (1968) 95–96. [15] J.-C. Huang, K. Xu, P. Yu, A unified gas-kinetic scheme for continuum and rarefied flows. II: Multi-dimensional cases, Commun. Comput. Phys. 12 (2012) 662–690. [16] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I: Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [17] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-Uniform Gases, Cambridge University Press, 1970. [18] W.G. Vincenti, C.H. Kruger, Introduction to Physical Gas Dynamics, Rober E. Krieger Publishing Company, 1975. [19] C.K. Chu, Kinetic-theoretic description of the formation of a shock wave, Phys. Fluids 8 (1965) 12–22. [20] X. Nie, X. Shan, H. Chen, Thermal lattice Boltzmann model for gases with internal degrees of freedom, Phys. Rev. E 77 (2008) 035701. [21] D. Bond, V. Wheatley, M. Macrossan, M. Goldsworthy, Solving the discrete Smodel kinetic equations with arbitrary order polynomial approximations, J. Comput. Phys. 259 (2014) 175–198. [22] V. Sofonea, R.F. Sekerka, Diffuse-reflection boundary conditions for a thermal lattice Boltzmann model in two dimensions: evidence of temperature jump and slip velocity in microchannels, Phys. Rev. E 71 (2005) 066709. [23] N.D. Masters, W. Ye, Octant flux splitting information preservation {DSMC} method for thermally driven flows, J. Comput. Phys. 226 (2007) 2044–2062. [24] F.P. Incropera, D.P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley & Sons, Inc., 2002. [25] F.M. White, Fluid Mechanics, McGraw Hill, 2003. [26] F. Sharipov, Non-isothermal gas flow through rectangular microchannels, J. Micromech. Microeng. 9 (1999) 394. [27] C. Cercignani, M. Lampis, Kinetic models for gas–surface interactions, Transp. Theory Stat. Phys. 1 (1971) 101–114.