## Abstract

Switching criteria for hybrid hydrodynamic/molecular gas flow solvers are developed, and are demonstrated to be more appropriate than conventional criteria for identifying thermodynamic non-equilibrium. For switching from a molecular/kinetic solver to a hydrodynamic (continuum-fluid) solver, the criterion is based on the difference between the hydrodynamic near-equilibrium fluxes (i.e. the Navier–Stokes stress and Fourier heat flux) and the actual values of stress and heat flux as computed from the molecular solver. For switching from hydrodynamics to molecular/kinetic, a similar criterion is used but the values of stress and heat flux are approximated through higher order constitutive relations; in this case, we use the R13 equations. The efficacy of our proposed switching criteria is tested within an illustrative hybrid kinetic/Navier–Stokes solver. For the test cases investigated, the results from the hybrid procedure compare very well with the full kinetic solution, and are obtained at a fraction of the computational cost.

## 1. Introduction

In a gas flow, the Knudsen number—the ratio between the mean free path *λ* of a gas molecule and a meaningful characteristic length of the process—is a dimensionless parameter that represents the degree of rarefaction. More specifically, it describes the extent to which the gas flow departs from local thermodynamic equilibrium. The classical hydrodynamic constitutive equations, i.e. the laws of Navier–Stokes and Fourier, follow from the Boltzmann kinetic equation in the limit of small Knudsen number by means of the Chapman–Enskog method (Chapman & Cowling 1970). That the Knudsen number can be used to characterize different rarefied regimes is therefore implicit in this method, and the modelling methods appropriate for these regimes follow from the analysis. Gas processes with small Knudsen number can be described with the equations of classical hydrodynamics, while processes with larger Knudsen number must be described by more elaborate means, e.g. the Boltzmann equation itself, or higher order continuum models (Reese *et al*. 2003; Struchtrup 2005).

The Knudsen number, in its global form, depends on a somewhat arbitrary choice of the macroscopic length scale. In a two-dimensional channel flow, for example, it is not clear whether this length should be half the channel height or the full channel height—this choice affects the Knudsen number by a factor of 2. More complicated flow processes will contain a multitude of relevant length scales, and assessment of the Knudsen number becomes correspondingly more difficult.

Additionally, some parts of a flow domain might be in the hydrodynamic regime, while others are in the Boltzmann kinetic regime. Since hydrodynamic codes are far faster than molecular/kinetic solvers, the computational effort can be minimized by using hybrid codes, which apply the appropriate solvers in the respective flow regions (Wang & Boyd 2003). However, these codes require switching criteria, also called ‘breakdown parameters’ (Kolobov *et al*. 2007), to identify these regions and for switching from hydrodynamic (continuum-fluid) to molecular/kinetic solvers and vice versa.

An often-used switching criterion is a local Knudsen number, where the length scale is a formulation based on the local spatial gradients of hydrodynamic variables, *viz*.(1.1)where *λ* is the molecular mean free path and *ϕ* is a significant flow quantity, typically density, temperature or pressure (Wang & Boyd 2003).

There are a number of ways to define the local Knudsen number, and each can give significantly different values. So the choice of what to use as a switching criterion has proven a problem in itself (Wang & Boyd 2003). This is especially apparent when considering microflows. For example, based on the definition given in equation (1.1), a low-speed micro gas flow will have a uniformly negligible local Knudsen number because the gradients of the flow variables are negligibly small. Nevertheless, non-equilibrium effects are far from negligible in these cases (Lockerby *et al*. 2005*a*). For hypersonic aerodynamic flows it has been suggested (Macrossan 2006) that *Ma*.*Kn*_{L}, where Ma is the local Mach number, is a more appropriate breakdown parameter; however, this parameter would be inappropriate for low-speed gas flows, which have a very small Mach number but can still be quite rarefied.

These issues prompt the question: is there a local Knudsen number definition, and hence a switching criterion, which is appropriate both for hypersonic and for micro gas flows?

In this paper, we address this question by introducing new local Knudsen number definitions as switching criteria measuring the departure from hydrodynamic behaviour. For switching from a molecular/kinetic solver to a hydrodynamic (continuum-fluid) solver, the switching criterion is based on the difference between the hydrodynamic near-equilibrium fluxes, i.e. the Navier–Stokes stress and Fourier heat flux, and the actual values of stress and heat flux as computed from the molecular solver. For switching from a hydrodynamic to a molecular solver, a similar criterion is used but the values of stress and heat flux are approximated through higher order constitutive relations. For the latter, we adopt the regularized 13 moment equations, which are an extension of hydrodynamics to third-order accuracy in the Knudsen number (i.e. super Burnett order) (Struchtrup & Torrilhon 2003; Struchtrup 2004, 2005).

## 2. Departure from hydrodynamics

Our arguments are based on the philosophy of the Chapman–Enskog method, and we recall some basic elements of kinetic theory and the Chapman–Enskog expansion before we proceed (Chapman & Cowling 1970; Struchtrup 2005). In this paper, we consider monatomic ideal gases exclusively.

The central quantity in kinetic theory is the velocity distribution function *f*(** x**,

*t*,

**), where**

*c**f*d

**d**

*c***is the number of molecules with microscopic velocities in the interval (**

*x***,**

*c***+d**

*c***) in the space element (**

*c***,**

*x***+d**

*x***) at time**

*x**t*. To find the distribution function we need to solve the Boltzmann equation.

The hydrodynamic quantities are velocity moments of the distribution function, with the following definitions for mass density *ρ*, velocity *v*_{i}, temperature *T*, pressure deviator (or stress) *σ*_{ij} and heat flux *q*_{i}:Here, *R* denotes the gas constant; *m* is the molecular mass of the gas; and *C*_{i}=*c*_{i}−*v*_{i} is the peculiar velocity; indices in angular brackets denote trace-free symmetric tensors. The pressure deviator is related to the stress tensor *t*_{ij} normally used in hydrodynamics by *σ*_{ij}=−(*t*_{ij}+*pδ*_{ij}), where *p*=*ρRT* is the ideal gas pressure.

The Chapman–Enskog method aims to find an approximate velocity distribution function from the Boltzmann equation. To this end, the distribution function is formally expanded in the Knudsen number,(2.1)where *f*_{E} denotes the Maxwellian equilibrium distribution and the additional terms obey the compatibility conditions(2.2)These imply that the Knudsen order corrections do not affect the basic hydrodynamic fields, *ρ*, *v*_{i}, *T*, but only higher moments, in particular stress *σ*_{ij} and heat flux *q*_{i}.

The stress and heat flux have no contribution from the equilibrium distribution *f*_{E}. The first-order correction, *f*^{(NSF)}, yields the laws of Navier–Stokes and Fourier,(2.3)with viscosity *μ* and heat conductivity *κ*. Deviation from conventional hydrodynamic behaviour is therefore described by the higher order terms *f*^{(2)}, *f*^{(3)}, … . Accordingly, we can write the stress and heat flux as(2.4)where and are the non-hydrodynamic contributions to stress and heat flux that result from *f*^{(2)}, *f*^{(3)}, … . From the Chapman–Enskog theory, we conclude that these are of second order or higher in the Knudsen number.

To be useful as a switching criterion, a local Knudsen number should indicate the degree of departure from equilibrium at a given point in the flow field. More precisely, if it is to be used to establish when the Navier–Stokes–Fourier equations can or cannot be employed, it should indicate the degree of departure from near-local equilibrium (i.e. the thermodynamic state that must exist for the Navier–Stokes–Fourier equations to be valid). Furthermore, it is not departure from near equilibrium in an absolute sense that we are interested in here; for example, low-speed flows of different flow magnitude will have the same rarefaction characteristics but depart from near equilibrium at different absolute levels—however, we would require the same modelling for each. So it is departure relative to the level of near equilibrium that we are interested in, i.e. relative to Navier–Stokes–Fourier hydrodynamics.

This therefore suggests, as measures for relative deviation from hydrodynamic behaviour, the following local Knudsen numbers:(2.5)where ‖.‖ denotes suitable norms to be discussed below. Since the hydrodynamic expressions and are of first order in the Knudsen number, and the non-hydrodynamic corrections and are at least of second order in the Knudsen number, their ratio should be of first order (or higher) in the Knudsen number. In general, mechanical and thermal effects may occur on different scales, so an overall local Knudsen number could be the maximum of the above values, i.e.

To complete our definition of the local Knudsen number,1 equation (2.5), it remains to choose suitable norms. For the definition of *Kn*_{q}, we propose using the length of the vectors, and for the definition of *Kn*_{σ}, we propose the square root of the second invariant of the (trace-free and symmetric) tensors, which gives a similar measure, i.e.(2.6)Both norms are invariant under linear transformations and thus independent of the choice of reference frame.

To test the performance of these local Knudsen number definitions as switching criteria, switching in both directions must be considered: molecular-to-continuum; and continuum-to-molecular. In some respects, determining appropriate switching from a molecular to a continuum solution is hypothetical, since the more accurate solution (the molecular one) has already been obtained. However, this may not be so clear cut in time-dependent simulations.

## 3. Local Knudsen numbers based on the R13 equations

In molecular-to-continuum switching, the actual values of stress, *σ*_{ij}, and heat flux, *q*_{i}, can be computed from the molecular/kinetic solution (the velocity distribution *f*) so the calculation of the local Knudsen numbers, equation (2.5), is straightforward. We may use the molecular/kinetic results for the velocity and temperature fields to estimate and as(3.1)

On the other hand, when switching from a hydrodynamic solver to a molecular-based simulation, the actual values of stress and heat flux are unknown: hydrodynamic solvers only produce the first-order contributions , , see equation (2.3).

However, higher order constitutive relations can be used to estimate the higher order contributions and from the hydrodynamic result. For this, we require a higher order theory that performs well in *both* low-speed and high-speed problems.

There are many competing high-order equation sets in the literature; space precludes a detailed discussion here, instead see Reese *et al*. (2003), Struchtrup & Torrilhon (2003), Lockerby *et al*. (2005*b*), Struchtrup (2005) and Lockerby & Reese (2008). Best known, perhaps, are the Burnett equations (Burnett 1936; Chapman & Cowling 1970; Struchtrup 2005), variants of which have been shown to accurately reproduce the viscous structure of one-dimensional shock waves (Reese *et al*. 1995). Lockerby & Reese (2008) tested a number of different high-order continuum-type equations against a simple low-speed benchmark case with no bounding surfaces. They concluded that the regularized 13 moment (R13) equations, proposed by Torrilhon and Struchtrup (Struchtrup & Torrilhon 2003; Struchtrup 2004, 2005) as a development of Grad's original 13 moment technique, provided the best model among the several tested. The R13 equations have also shown good predictive capabilities in high-speed flows (Torrilhon & Struchtrup 2004). Only recently, a set of boundary conditions for the R13 equations has been provided (Struchtrup & Torrilhon 2007; Torrilhon & Struchtrup 2008), and it was shown that the equations can also predict Knudsen layer dominated problems to some extent (Struchtrup & Torrilhon 2008; Taheri *et al*. 2009).

We therefore propose the following local Knudsen number definitions based on the R13 equations:(3.2)If the R13 equations are being solved within the hydrodynamic solver, then these Knudsen numbers are straightforward to calculate: and form part of the result, while and can be calculated from the R13 results for velocity and temperature, i.e. , in the same way as described above for the molecular solver.

However, the local Knudsen numbers may also be estimated cheaply within a conventional Navier–Stokes–Fourier hydrodynamic solver, as(3.3)where and are calculated from the R13 constitutive relations, but using the hydrodynamic variables provided from the solver. The error this introduces is acceptable for switching purposes, as an estimate of the local Knudsen number is sufficient. The full R13 equations are given in appendix A, where we also state which terms will be computed from the hydrodynamic fields.

The idea of using the R13 equations with data from hydrodynamics recalls the testing method for higher order models proposed in Zheng *et al*. (2006), where results from a microscopic solver were inserted into the higher order models.

## 4. Examples

We consider several benchmark flow examples to explore the usefulness of our definitions for the local Knudsen number.

### (a) Example I. Shock structure using the Burnett equations

As an illustrative introduction to our new approach, we first examine the classical shock structure problem and draw conclusions about the local Knudsen number derived using the classical Burnett equations. Normal shock waves are essentially one-dimensional structures. The (trace-free) stress tensor is diagonal, with and the Navier–Stokes normal stress iswhere *x* is the direction longitudinally through the shock and *v* is the *x*-component of velocity. The dominant nonlinear term featuring in the Burnett expression for the stress is (Chapman & Cowling 1970)(4.1)where *A* is a dimensionless constant for the particular gas. If this Burnett expression is used to estimate the non-equilibrium stress , we find from equation (2.5) that the local Knudsen number is(4.2)Here, we have used the mass balance d(*ρv*)/d*x*=0, the mean free path (for convenience, without the factor , which is often found in literature) and the local Mach number . The coefficient *α* is a number of order unity.

The form given in equation (4.2) is equivalent to the usual local Knudsen number, equation (1.1), multiplied by the Mach number. This is the same local breakdown parameter identified by Bird (1970) for high-speed expanding flows, and closely related to Tsien's (1946) parameter, which was also identified by Macrossan (2006) as a better indicator for high-speed flows than the local Knudsen number, equation (1.1), alone.

So our present proposal can be related directly to accepted approaches for nonlinear high-speed flows, but has the additional advantage that it should also be relevant for low-speed micro flows, as will become clear in the following sections.

### (b) Example II. Nonlinear shear flow with second-order hydrodynamics

For nonlinear shear flow between two parallel plates, both the Burnett and R13 equations (to second order in *Kn*) yield (Struchtrup 2005; Struchtrup & Thatcher 2007),(4.3)where *y* is the direction normal to the plates. Note that *q*_{1} is a heat flux perpendicular to the temperature gradient. In classical hydrodynamics, *σ*_{12} and *q*_{2} have the values given above, while *σ*_{11}=*σ*_{22}=*q*_{1}=0. From equation (2.5) the local Knudsen numbers are(4.4)which are in essence agreeing with equation (4.2). The major difference lies not in the numerical factors but in the fact that for shear flow, the velocity gradient cannot be expressed through the density gradient.

### (c) Example III. Linear Poiseuille flow with the R13 equations

We now show that the local Knudsen number, equation (2.5), gives meaningful results also for linear flows.

In Struchtrup & Torrilhon (2007), linear force-driven Poiseuille flow was considered using the R13 equations with jump and slip boundary conditions. In dimensionless form, the analytical result reads(4.5)where *F* is the dimensionless driving force and is the global Knudsen number based on the channel height *H*. All other non-equilibrium quantities vanish in the linear regime, and the temperature is constant; the Fourier heat flux vanishes, making *Kn*_{q} irrelevant.

The Navier–Stokes stress computed using the R13 result for velocity reads(4.6)where and *σ*_{12} differ only due to Knudsen layer effects. The local Knudsen number, equation (2.5), then becomes(4.7)and figure 1 shows this function of the dimensionless space coordinate and the global Knudsen number for global Knudsen numbers between 0.01 and 1. Note that while the R13 equations are accurate only for Knudsen numbers below 0.5 (Struchtrup & Torrilhon 2008), they give a valid qualitative description of flows with larger Knudsen numbers. So equation (4.7) gives a useful estimate for the local Knudsen number.

For very small global Knudsen number, *Kn*=0.01, the local Knudsen number *Kn*_{σ} is zero in the bulk but is one order of magnitude larger than *Kn* close to the walls, i.e. in the Knudsen layers. This indicates that hydrodynamics is a sufficient model of the bulk flow, but a more sophisticated approach is needed near the walls. Owing to the relatively small extent of the Knudsen layer in this case, it might not be necessary to resolve the Knudsen layer, and proper higher order slip boundary conditions might be sufficient (Cercignani 1990; Lockerby *et al*. 2004; Struchtrup & Torrilhon 2008).

As the global Knudsen number grows to *Kn*=0.05, the Knudsen layer extends further into the gas, with the local Knudsen number below *Kn*_{σ}=0.1 at the walls and zero towards the middle of the channel. For higher global Knudsen numbers, the local Knudsen numbers *Kn*_{σ} are larger as well and become almost constant across the channel. For this particular problem, the values of *Kn*_{σ} lie below those of *Kn* and converge towards 0.25.

### (d) Example IV. Linear Poiseuille flow with a BGK kinetic solver

We now consider a numerical solution of the Bhatnagar–Gross–Krook (BGK) kinetic model equation for the case of low-speed Poiseuille flow. The global Knudsen number, based on the full channel height, is *Kn*=0.08.

In figure 2, a hydrodynamic solution with second-order slip boundary conditions (Cercignani 1990; Struchtrup & Torrilhon 2008) is plotted with a solution to the BGK kinetic equation (see Chapman & Cowling 1970) obtained using a discrete velocity method. The reader is referred to Valougeorgis (1988) and Valougeorgis & Naris (2003) for a detailed description of the mathematical and numerical formulation of this method (and which also deal with the application of more sophisticated model equations than the BGK approximation used here). In both simulations, 500 grid points are used. The new local Knudsen number, from equation (2.5), is also plotted, where(4.8)with *σ*_{12} the shear stress and *v* the flow velocity, both computed from the BGK solver.

As in the previous example, the local Knudsen number *Kn*_{σ} is identifying the non-equilibrium in the Knudsen layers close to the confining surfaces clearly and distinctly from the near-equilibrium bulk flow (where switching to a continuum solution could occur without significant error). Conventional local Knudsen number definitions based on equation (1.1) could not provide this information since their values are negligibly small throughout the channel.

### (e) Example V. Switching from a linear NSF solver to a molecular solver

For a simple linear shear flow with, again, a dimensionless driving force *F*, classical hydrodynamics gives(4.9)in dimensionless variables, while the linearized R13 equations reduce to(4.10)To apply the definition of the local Knudsen number in equation (3.3), we use the solution of the Navier–Stokes equations, equation (4.9), in the right-hand side of the expression for in equation (4.10), giving(4.11)or, more instructively,(4.12)Clearly, when the force is constant in the *y* direction, as would be the case for the force-driven channel flow of the previous example, we have , and the computation of the local Knudsen number according to equation (3.3) simply yields *Kn*_{σ}=0. As was seen in example III, in the linear regime, the difference between the two hydrodynamic models is only due to Knudsen layers: in the linear regime, our local Knudsen numbers can identify when to switch from a higher order theory or molecular/kinetic solver to classical hydrodynamics, but not when to switch in the other direction.

While our proposed local Knudsen numbers cannot identify the inability of the Navier–Stokes–Fourier model to reproduce Knudsen layers, the next example will show that the method can identify rarefaction in the bulk when, for example, the force is a function of space.

### (f) Example VI. A hybrid flow calculation

Following the discussion in §4*e*, we finally investigate the steady-state channel flow response to the non-constant body force,where *x* is in the flow direction (and perpendicular to *y*), , *a*_{2}=10^{3}, and is the global Knudsen number based on the full channel height *H*. The variation of this normalized body force through the channel is shown in figure 3 for the global Knudsen number *Kn*=0.04. This simple forcing function, while somewhat artificial, has been chosen because it generates a shear flow, shown in figure 4, that exhibits both near-equilibrium and strong non-equilibrium behaviour in the bulk flow (away from the walls). This flow is reminiscent, in some respects, of the velocity variation through a stationary monopole vortex.

From equation (4.11), the local Knudsen number, equation (3.3), becomes(4.13)For this geometry, our definition of the local Knudsen number gives a quadratic relationship between the local Knudsen number *Kn*_{σ} and the global Knudsen number *Kn*, since there are no second-order contributions in the equation for .

To test the performance of equation (4.13) as a switching parameter, we construct a simple hybrid solution procedure. First, the Navier–Stokes equations are solved, shown by the dotted line in figures 5 and 6 for global *Kn*=0.02 and 0.08, respectively. The local Knudsen number is then calculated from equation (4.13) and, based on a local Knudsen number threshold of 0.2, the high-*Kn*_{σ} portion of the domain is handled by the BGK kinetic solver (at the same spatial resolution). As before, the solver of Valougeorgis (1988) and Valougeorgis & Naris (2003) is used with the standard BGK approximation to the collision integral. In this simple case, the near-equilibrium (Navier–Stokes) distribution is (Chapman & Cowling 1970, §6.6)(4.14)where *c*_{x} and *c*_{y} are components of molecular velocity in the streamwise and cross-stream directions, respectively. Enforcing a flow gradient in the BGK solver as a boundary only requires the odd part of the distribution (with respect to *c*_{x} and *c*_{y}) to be prescribed. So, if the Navier–Stokes velocity gradient is [d*v*/d*y*]_{NS} at the switching points, in the BGK solver this velocity gradient (and therefore near equilibrium) can be enforced using(4.15)where *f*(*c*_{y}) and *f*(−*c*_{y}) are the distributions of molecules entering and exiting the domain, respectively. The velocity profile, *v*, calculated by the BGK solver is then uniformly scaled (by a few per cent), so that the velocities at the switching points match the velocities of the Navier–Stokes solution at the same points. The solutions are then combined across the domain.

Figures 5 and 6 show the solutions from this illustrative hybrid method, compared with pure Navier–Stokes and BGK kinetic solutions, for global *Kn*=0.02 and 0.08. The hybrid simulation at *Kn*=0.02 is of comparable accuracy to the complete BGK kinetic solution, but has been obtained with a kinetic domain size that is five times smaller. The relative proportions of the deconstructed domain are 20%/80% (kinetic/Navier–Stokes) for *Kn*=0.02; 34%/66% for *Kn*=0.04; and 57%/43% for *Kn*=0.08. In all cases, the discrepancy between the complete BGK solution and the hybrid method could be reduced by adopting a lower switching threshold on the local Knudsen number, i.e. by solving a greater proportion of the domain using the BGK solver.

Note that had we solved the R13 equations fully, a smaller kinetic domain would have been required (and a higher local Knudsen number threshold would be permissible). This may offset, partially at least, the additional computational expense of their solution.

## 5. Modified switching parameters

There are certain situations in which a modified version of the local Knudsen number discussed in §2 might be warranted—for example, for flow configurations that have wide-ranging degrees of non-equilibrium (say, in stress). In the same simulation, there could exist regions of high local Knudsen number (calculated from equation (2.5)) at relatively low stress, and other regions with much greater stress but a negligible local Knudsen number. In such a case, it would not be efficient for a hybrid solver to switch to a molecular model to calculate stress in the regions of high local Knudsen number because values for the stress, although highly non-equilibrial, would be negligible in magnitude compared with those in other regions of the flowfield. In such circumstances, a refined set of switching criteria are(5.1)where the subscript ‘max’ denotes a maximum across the spatial domain (although quantities other than the maximum could reasonably be used). Figure 7 shows this modified local Knudsen number calculated for linear Poiseuille flow using BGK kinetic data (for a global *Kn*=0.08); these results are directly comparable to those shown earlier in figure 2. The Knudsen layers are again clearly identified, although a lower threshold would be required for switching with this definition. The illustrative hybrid simulations shown in figures 5 and 6 can be reproduced identically using our modified local Knudsen number of equation (5.1), but with a threshold of 0.01, rather than 0.2.

In circumstances where the spatial maximum does *not* reflect the range of non-equilibrium in a solution, and/or situations where the range of non-equilibrium is not reflective of the grid resolution required by the user, then an alternative modification to equation (5.1) is required.

Another situation where the use of equation (2.5) alone might be problematic is when and , or and , are zero in different spatial locations. In these situations, our definition of the local Knudsen number tends to infinity. However, this can be handled relatively easily by using the modified Knudsen number in equation (5.1) instead, or by the following purely local modification:(5.2)Effectively, this definition limits the local Knudsen number to values less than or equal to one. For switching purposes, we are interested in threshold values significantly lower than one, so this modification seems reasonable.

## 6. Conclusions

A set of local Knudsen numbers has been defined, which have been demonstrated to be more appropriate than the conventional ones for identifying micro *and* hypersonic gas flow non-equilibrium. The problem of choosing appropriate molecular/hydrodynamic switching criteria has been addressed by adopting a local Knudsen number definition based on higher order constitutive relations; here the regularized 13 moment (R13) equations were chosen.

Several benchmark tests showed that our proposed local Knudsen number gives meaningful results. For strongly nonlinear flows, the definition agrees with the usual definitions in the literature.

We also described a procedure to estimate efficiently the R13 local Knudsen number within a Navier–Stokes solver, and the efficacy of this as a switching criterion has been tested within an illustrative hybrid BGK/Navier–Stokes solver. For the test case investigated, the results from the hybrid solver compare very well with the full BGK solution, and are obtained at a fraction (that depends on the global *Kn*) of the computational cost.

Since the local Knudsen number measures the relative deviation from classical hydrodynamics, it can be considered as the error associated with the use of classical hydrodynamics (more precisely, the error in the determination of the local stress and heat flux). It depends on the flow application how much error the user considers to be acceptable, and so at what value of the local Knudsen number switching should occur.

## Acknowledgements

The authors would like to thank Prof. Dimitris Valougeorgis for providing a well-documented DVM code that generated the BGK results presented in this paper. J.M.R. gratefully acknowledges the support of the UK's Engineering and Physical Sciences Research Council under grant no. EP/D007488/1. H.S. gratefully acknowledges support by the Natural Sciences and Engineering Research Council (NSERC). We also thank the referees of this paper for their very helpful comments.

### Appendix A. The R13 equations

For convenience, we present the complete R13 equations in this appendix; see Struchtrup & Torrilhon (2003) and Struchtrup (2004, 2005) for the detailed derivation and deeper discussion.

The basic equations are the conservation laws for mass, momentum and energy, which can be written as(A1)where *θ*=*RT* is the temperature in energy units; is the convective time derivative; *F*_{i} is the external body force; andThe balance equations for stress and heat flux can be written as(A2)(A3)with the additional constitutive equations

When the R13 equations themselves are to be solved, and on the left-hand sides of (A 2) and (A 3) are just the stress and heat flux, i.e. and .

However, when the R13 equations are used to estimate the deviation from hydrodynamics for use in equation (3.3), then the values on the right-hand side of the above equations must be replaced by the hydrodynamic values, and .

## Footnotes

↵It may be helpful to refer to this as a local

*transport*Knudsen number in order to distinguish it from the classical local Knudsen number; in this paper, though, in the interest of brevity, we simply use ‘local Knudsen number’.- Received December 5, 2008.
- Accepted January 22, 2009.

- © 2009 The Royal Society