Coupling heterogeneous continuum-particle fields to simulate non-isothermal microscale gas flows

This paper extends the hybrid computational method proposed by Docherty et al. (2014) for simulating non-isothermal rarefied gas flows at the microscale. Coupling a continuum fluid description to a direct simulation Monte Carlo (DSMC) solver, the original methodology considered the transfer of heat only, with validation performed on 1D micro Fourier flow. Here, the coupling strategy is extended to consider the transport of mass, momentum, and heat, and validation in 1D is performed on the high-speed micro Couette flow problem. Sufficient micro resolution in the hybrid method enables good agreement with an equivalent pure DSMC simulation, but the method offers no computational speed-up for this 1D problem. However, considerable speed-up is achieved for a 2D problem: gas flowing through a microscale crack is modelled as a microchannel with a high-aspect-ratio cross-section. With a temperature difference imposed between the walls of the cross-section, the hybrid method predicts the velocity and temperature variation over the cross-section very accurately; an accurate mass flow rate prediction is also obtained. 2016 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
The behaviour of fluid flows through and around micro-and nano-scale devices is still not well understood; traditional continuum fluid mechanics is often inaccurate due to the importance of the fluid's molecular nature at these scales. Deterministic molecular dynamics (MD) may be employed to model liquid flows at the nano-scale, while the stochastic direct simulation Monte Carlo (DSMC) method is the most popular computational tool for dilute gases at the micro-or nano-scale. Both of these techniques are, however, too computationally intense to resolve the spatial and temporal scales in real engineering flow problems. The DSMC method is significantly cheaper than MD, but the expense of tracking and computing collisions between thousands/millions of DSMC particles can require months (or even years) of computing time. Therefore, continuum-molecular 'hybrid' methods are being developed to reduce this expense. Hybrid methods combine the efficiency of a continuum-fluid description with the resolution and accuracy of a molecular treatment: the molecular tool is applied over micro/nano scales to resolve the molecular flow behaviour, while the continuum-fluid description is employed over macro scales to resolve macroscopic flow variations.
The majority of current continuum-molecular hybrid methods are based on a domain decomposition (DD) framework [2][3][4][5][6], where the molecular solver is applied in a 'micro' sub-domain (which is typically close to a bounding wall) and the continuumfluid solver is applied in the remaining 'macro' sub-domain; an overlap region then enables coupling of the two solvers, as illustrated in Fig. 1(a). These DD methods are, however, limited to flow problems where microscopic resolution is required only in localised regions. The less common 'equation-free' [7] and 'heterogeneous' [8] multiscale approaches are able to simulate problems that require the molecular solver to provide information to the continuum solver everywhere in the flow, i.e. when the conventional linear fluid-constitutive relations fail in the bulk of a flowfield -this could be the case if non-equilibrium flow appears in the bulk, or if the transport properties are unknown in an unusual gas mixture. In this case, a continuum-fluid description is applied across the entire flowfield, and spatially-distributed micro elements (in which the molecular solver is implemented) are then deployed to 'correct' this continuum description.
The Heterogeneous Multiscale Method (HMM) [8] uses a 'point-wise coupling' approach, illustrated in Fig. 1(b). The micro elements supply information (consisting of updated boundary conditions and fluid-constitutive information) directly to the nodes of the macro grid, both at the bounding walls and in the bulk. At continuum solution at the collocated macro node. This point-wise coupling means that the position and size of the micro elements is restricted by the placement and density of the macro nodes. It also means that, while the HMM is effective when the spatial scales are highly separated (i.e. when the variation of the flow properties relative to the physical extent of a micro element is small), it is inefficient and potentially inaccurate for flow problems that exhibit mixed degrees of spatial scale separation.
To overcome the restrictions of point-wise coupling, the HMM with field-wise coupling (HMM-FWC) was proposed by Borg et al. [9], which builds on the equation-free [7] and heterogeneous [8] multiscale frameworks. A continuum description is again applied over the entire flowfield, but each micro element now represents a field that correlates directly with an identicallysized continuum sub-region. The coupling is then performed via fields rather than nodal points: the local continuum property fields are imposed over the micro elements, and local constitutive correction fields (i.e. corrections to properties derived from linear constitutive-fluid relations) are extracted 1 ; these local corrections are then interpolated to provide global corrections across the entire flowfield. Essentially, the HMM-FWC takes advantage of the fact that properties like stress and heat flux (and hence their corrections) often vary slowly in space, and so computational savings can be made by interpolating these fields between micro elements that are more sparsely located than the macro nodes. As in the HMM, near-wall micro elements also provide the continuum-fluid description with updated boundary information. The HMM-FWC can be considered a more general heterogeneous approach than the HMM -the position and size of the micro elements is not restricted by the macro nodes, as indicated in Fig. 1(c), and flows with mixed degrees of spatial scale separation can be simulated efficiently.
A specific class of the HMM-FWC has recently been developed to simulate efficiently the flow through long micro/nanochannels, which are a fairly common feature of emerging micro and nano devices. Presented by Borg et al. [10], the framework and coupling strategy of the Internal-flow Multiscale Method (IMM) is tailored to exploit the large length scale separation that exists in the streamwise direction of these flows. A continuum-fluid description is applied over the entire channel, and very short micro elements occupying the entire cross-section (i.e. the full channel height in 2D problems) are distributed along the channel length, as shown in Fig. 1(d). Compared with the other heterogeneous methods, the coupling is simplified -pressure gradients are imposed over the periodic micro elements via body-forcing, and the resulting mass flux is used to correct the continuum description.
Both the HMM-FWC and the IMM were originally implemented using MD as the micro solver for the simulation of liquids, and with assumptions of incompressible and isothermal flow. The IMM has seen significant development since: it has been implemented for a continuum-DSMC coupling [11,12], and extended to compressible [11], non-isothermal [12], and unsteady flows [13]. Note, however, that the non-isothermal coupling strategy in [12] is not general, and is applicable only to the long micro/nanochannel flows tackled by the IMM.
In 2014, Docherty et al. [1] adapted the HMM-FWC for a continuum-DSMC coupling, tailoring the constraint of the micro elements to suit the use of a DSMC solver. The focus was on heat transfer problems, and non-isothermal coupling was achieved via the conservation of energy. However, this coupling strategy limited the method to 'stationary' heat transfer problems, i.e. where the  Real gas flow problems often involve the transport of momentum and heat. This paper makes progress towards a general hybrid method that can be applied to a wide range of dilute gas flows, by extending the coupling strategy of the continuum-DSMC method of [1] to consider the transport of mass, momentum, and heat, simultaneously. This extended HMM-FWC method can therefore be used to simulate steady, non-stationary, compressible, and non-isothermal flows.
To date, the HMM-FWC has been implemented only for 1D flow problems. This paper demonstrates its implementation in 2D by considering a simplified 2D model of a gas flow through a microscale crack. In addition, this 2D flow problem highlights the method's ability to deal efficiently with mixed degrees of spatial scale separation.
The structure of this paper is as follows. The methodology of our extended HMM-FWC approach is discussed in Section 2 -this discussion applies to 1D, 2D, and 3D problems. In Section 3, 1D validation is performed on a high-speed micro Couette flow problem. Section 4 then applies the method to the larger and more realistic 2D crack flow problem mentioned above. Conclusions are drawn in Section 5.

Micro solver: the DSMC method
The DSMC method tracks a large number of gaseous particles that move through a computational mesh, storing their position, velocity, and internal state. Compared with deterministic MD, the DSMC method reduces computational expense through two assumptions. First, each simulated particle represents a large number of real gas molecules. Secondly, particle collisions are treated stochastically, rather than deterministically: particle movements and collisions are decoupled over small time intervals, with the movements computed deterministically based on the velocity vector and the DSMC time step, and the collisions are then modelled statistically.
The open source C++ toolbox OpenFOAM [14] includes a DSMC solver named dsmcFoam that offers the key features of any modern DSMC code. This code has been extended within the research collaboration of the authors [15], and this solver dsmcFoamStrath has been validated for various benchmark cases, including hypersonic and microchannel flows [16][17][18]. It is used to perform all of the hybrid and pure DSMC simulations in this paper. The Variable Hard Sphere (VHS) collision model is adopted for all of these simulations.
Note that the dsmcFoamStrath code includes new measurement and state controller tools that are useful when implementing the coupling algorithm of our hybrid method. Volume-based measurements are performed by averaging over all particles in the volume [19], while wall measurements are obtained by summing over all particles that strike the surface [20,21]. State controllers achieve a desired macroscopic state in a prescribed volume by manipulating the molecular behaviour; a simple feedback loop algorithm is used [22], where the measured property (for example, the velocity) is compared against its target value, and the difference provides the basis for the controlling action. The need for these tools is discussed in Section 2.3.1.

Macro solver: the continuum-fluid description
The steady-state conservation equations of fluid mechanics form the basis of our hybrid method, i.e.
and r Á ðqeuÞ ¼ Àr Á ðpI Á uÞ þ r Á ðs Á uÞ À r Á q; ð3Þ where q is the fluid mass density, u the velocity vector, s the viscous stress tensor, f an external force per unit volume, e the specific energy, and q the heat flux vector; extending the hybrid method to unsteady flow will be the subject of future work. For most conventional flow problems, the linear Navier-Stokes-Fourier (NSF) constitutive relations are used to close these continuum equations. The focus of this paper, however, is on flows where the linear NSF relations fail.
In our heterogeneous hybrid method, the continuum equations above are unmodified. The constitutive relations used for closure are, however, the NSF relations augmented by 'correction fields', i.e.
s ¼ lðruÞ þ lðruÞ T À 2 3 lðr Á uÞI þ X; where X is a stress correction tensor and U is a heat flux correction vector, l is the gas dynamic viscosity and j is the gas thermal conductivity. Closing Eqs. (2) and (3) with these 'corrected' constitutive relations gives, and r Á ðqeuÞ ¼ Àr Á ðpI Á uÞ þ r Á ½lðruÞ Á u þ r Á ½lðruÞ T Á u respectively. Along with the gas equation of state, Eqs. (1), (6), and (7) can then be used to describe the entire gas flowfield. 2 The correction fields X and U are supplied by distributed DSMC micro elements, which thereby incorporate in the flow model any effects due thermodynamic non-equilibrium and/or any inaccuracies in the transport properties l and j.

Coupling strategy
At the start of a hybrid simulation, the continuum-fluid description is applied over the entire flowfield; the correction fields are equal to zero and the conventional no-slip/no-temperature-jump boundary conditions are imposed. Micro elements are then dispersed across the flowfield (with the arrangement depending on the flow problem), and we begin the iterative coupling algorithm of the hybrid method.
The two-way coupling exchange performed in each iteration is summarised in Fig. 2. Each micro element is constrained by the local continuum solution, and DSMC is performed in each element. Constitutive correction fields and updated boundary conditions are then extracted from the micro elements and passed to the continuum-fluid description. Solving the continuum equations then produces 'corrected' continuum property fields -these are used to constrain the micro elements again, and the process is repeated. With continuing iterations, the corrected continuum solution converges towards a pure DSMC solution of the same problem, as will be demonstrated.

Macro-to-micro coupling: constraining the micro elements
The corrections to the constitutive relations must be extracted from micro elements that properly represent the local conditions in the macro domain. For dilute gases, this is achieved when the correct particle distribution is imposed at the boundaries of the micro element. However, the particle distribution required at the element boundaries cannot be extracted directly from the macroscopic continuum solution -such detail is not available from a continuum-fluid description. This problem is circumvented here by introducing an artificial 'relaxation' region around each micro element that facilitates a natural relaxation to the correct particle distribution (and macroscopic state) within the core of the element. It is not essential that the gas state in the relaxation region accurately represents the conditions in the corresponding macro domain; its sole purpose is to develop the boundary conditions needed for the core region of the element. Sampling of the property fields is then performed only in this core 'sampling region'. A schematic of a 3D bulk micro element is shown in Fig. 3.
In this paper, the imposition of the micro constraints is tailored to dilute gases (i.e. to a DSMC micro solver) and is achieved in two stages. First, the spatial variation of the local continuum properties is imposed throughout the relaxation region by implementing state controllers for the duration of the micro DSMC simulation. Secondly, a particle distribution is applied at the outer boundaries of the relaxation region; local Maxwellian distributions are applied here for simplicity. 3 The relaxation region must then be large enough that the particle state relaxes fully (via particle collisions) across its extent. If this is the case, the particle distribution in the sampling region will then be dictated solely by the applied continuum state with no effects from the imposed Maxwellian distribution. When the hybrid solution converges, the artificiality of the relaxation region dissolves seamlessly into the true particle distribution and fluid state in the core sampling region.
Force-driven flows are modelled in pure DSMC simulations by applying an external acceleration to all particles in the system; although the acceleration is applied continuously, particle-wall interactions mean that the macroscopic velocity will reach a constant steady-state. In an equivalent hybrid simulation, we apply the same acceleration to all of the particles in the micro elements (in both the relaxation and sampling regions), in addition to the micro constraints discussed above.
In summary, the sampling region of a micro element is surrounded by a relaxation region: the local continuum property fields  are imposed in this relaxation region, and a particle distribution is applied at its outer boundaries; if required, an external acceleration is also applied to all particles in the element. The sampling region of a near-wall micro element must, however, be adjacent to the wall itself in order to capture the non-equilibrium behaviour that exists at bounding walls. Fig. 4 shows bulk and near-wall micro elements in an example computational set-up for a 2D problem.

Sampling region
Our coupling strategy requires that the spatial variation of the flow properties is both imposed on and extracted from different regions of each DSMC micro element. We achieve this by using state controller and measurement tools based on the method of bins (MOB): a volume (the entire domain or a user-defined zone) is split into a number of 'bins', and the spatial average of the properties in each bin is controlled or measured. Each relaxation region is therefore divided into a grid of 'control bins', and state controllers impose the continuum properties in each bin to create the desired spatial variation across the region. Similarly, each sampling region is divided into a grid of 'measurement bins' and the property fields are extracted by averaging over the particles in each bin; in the near-wall micro elements, boundary information is also extracted from the wall-adjacent faces of the measurement bins by summing over all particles that strike the surface. The key advantage of this binning method is its independence from the computational mesh used by the DSMC algorithm, which enables us to define the resolution of the extracted/imposed macroscopic fields, and to control the level of noise, without affecting the DSMC computations. As well as being independent of the DSMC mesh, the measurement and control bins can also be independent of each other, and of the continuum mesh -if they are not collocated with the cells of the continuum mesh, then data can simply be interpolated between the meshes. Sufficient micro resolution is crucial to obtaining an accurate result from the hybrid method, and the appropriate micro element configuration (i.e. the number, location, and size of the micro elements) will vary from problem to problem; this is discussed in more detail in [1,24]. However, for a given number of micro elements at predefined locations, the extents of both the sampling and relaxation regions can be made to adapt dynamically depending on the local flow if these extents are defined as some number of local mean free paths k l -in this paper, the appropriate number for each region is determined from some initial testing. At each iteration, the dimensions of each element are then set depending on k l : in the first iteration, k l is assumed equal to the global mean free path k gl ¼ Kn=l, where Kn is the Knudsen number and l is some characteristic dimension of the system; in subsequent iterations, k l is updated for each micro element according to the VHS expression [19], where x is the viscosity exponent of the gas. The temperature T av and the density q av are spatial averages measured in the sampling region of the element in the previous iteration.

Micro-to-macro coupling: correcting the continuum-fluid description
The data transferred from the spatially-distributed micro elements to the continuum-fluid description must not contain too much statistical scatter, as this could lead to instability of the coupling algorithm. Each micro DSMC simulation is therefore performed in two stages: a transient period enables the flow to reach steady-state, and a longer cumulative averaging period then reduces scatter in the measured properties. 4 The correction fields are extracted from the micro elements (both near-wall and bulk) in two steps. In the first step, the local correction fields across each sampling region are computed, i.e. the corrections in each measurement bin and at each wall boundary are calculated from the time-averaged flow properties measured in that bin or at that wall boundary. The stress correction is calculated according to Eq. (4) based on the measured stress and the velocity gradient; similarly, the heat flux correction is calculated from Eq. (5) using the measured heat flux and the temperature gradient. Note that the property gradients in each bin and at each wall boundary are approximated by finite difference representations here (based on the properties in the adjacent bins or the adjacent wall boundary). The second step is then to obtain approximations of the global correction fields across the entire flowfield. In this paper, linear interpolations between adjacent local correction fields are used for simplicity.
In addition to constitutive correction fields, the micro elements also provide the continuum-fluid description with updated boundary information. The properties (i.e. the velocity, temperature, and pressure) of the gas at the wall are extracted from the walladjacent faces of the measurement bins in each near-wall sampling region. This local boundary information is then linearly interpolated (or extrapolated) to produce estimates at all bounding walls. The configuration of the near-wall micro elements is dictated by the need to construct a reasonable approximation of the boundary properties with this type of interpolation.

Iterative algorithm for steady flows
The general algorithm of the hybrid method therefore consists of the following steps: (0) Assuming no constitutive corrections (i.e. X ¼ U ¼ 0) and no-slip/no-temperature-jump at the bounding walls, solve the continuum Eqs. (1), (6), and (7), and the gas equation of state (8), to obtain an initial NSF velocity field u NSF , temperature field T NSF , pressure field p NSF , and density field q NSF across the entire system.
(1) Initialise and constrain each DSMC micro element: (a) Impose the local continuum property fields (i.e. the velocity, temperature, and density) across the relaxation region by implementing state controllers in each control bin. (b) Apply a Maxwellian particle distribution at the outer boundaries of the relaxation region based on the local continuum temperature and velocity.  (7), and the gas equation of state (8), to obtain a new, corrected velocity field u corr , temperature field T corr , pressure field p corr , and density field q corr across the system.
(7) Using this corrected solution to constrain the micro elements, repeat from Step (1) until the corrected fields do not change between iterations to within user-defined tolerances.

1D validation: high-speed micro Couette flow
We now validate our hybrid technique on a high-speed micro Couette flow problem. A gas is confined between two parallel planar walls that have the same temperature but are moving with different velocities: with a Mach number Ma > 0:2, viscous heating occurs in the gas.

Applying the hybrid method in 1D
Fig . 5 shows the computational set-up for this 1D flow problem: the continuum grid consists of M y macro nodes, and an example micro element arrangement is shown. The measurement and control bins are set to have the same height dy. For convenience, the centre of each bin b coincides exactly with a macro node j (meaning that the bin height dy is equal to the macro node spacing Dy) -while this collocation of the bins and macro nodes is not a requirement of the method, it simplifies the transfer of data between the continuum mesh and the micro elements.
The implementation of the hybrid algorithm for this high-speed Couette flow problem is described in A. Note that mass conservation is satisfied automatically, and so density coupling is required only in the macro-to-micro direction: the continuum density field is imposed on the micro elements, but no DSMC data is needed to satisfy mass conservation. Both momentum and heat flux coupling are, however, performed in both directions.
Reasonable and constant reference values for the gas dynamic viscosity and thermal conductivity are assumed here (l ref and j ref , respectively). Any variations of these transport properties are then modelled indirectly through the constitutive correction fields which, for this flow problem, reduce to X xy ; X yy , and U y .
Linear interpolations between the sampling regions approximate these corrections everywhere in the flowfield, while interpolations of the gas boundary information are not required in this simple flow case.
The hybrid method provides a corrected streamwise velocity field u corr , pressure field p corr , temperature field T corr , and density field q corr . For a property s (i.e. s ¼ u; T; p, and q), convergence occurs when, where k is the iteration index, j ¼ 1; 2; . . . M y , and f tols is the property tolerance value that depends on both the property and the case itself; the hybrid algorithm is converged when all properties satisfy

Results
Monatomic argon is chosen as the working gas for simplicity: a molecular mass m ¼ 66:3 Â 10 À27 kg, a viscosity exponent Mach number Ma ¼ ju wall jðcRT wall Þ À0:5 % 5; although such hypersonic conditions are unlikely in the types of problems that Couette flow might represent, the aim here is to test our hybrid coupling strategy under extreme conditions. From Eq. (9), the gas density is set to obtain a global Knudsen number Kn gl ¼ k gl =H ¼ 0:01. While this value of Kn gl is small, the high-strain rates that produce viscous heating in the bulk of the flowfield also produce local non-equilibrium effects. Viscous heating means that the average gas temperature T av will be considerably higher than T wall : with no velocity slip or temperature jump, the NSF solution predicts T av ¼ 1294 K. Based on this temperature, a reference viscosity l ref ¼ 6:25 Â 10 À5 kg/ms and a reference conductivity j ref ¼ 0:049 W/m K are assumed [25].
The same cell dimensions and DSMC time step dt ¼ 1 Â 10 À12 s are adopted in both the hybrid micro simulations and the pure DSMC simulation of the entire flowfield that is used for validation purposes. The start-up and averaging periods are also the same for both: a start-up run of 3 million DSMC time steps enables all simulations to reach steady-state, and an averaging period of 20 million DSMC time steps reduces scatter in the measured flow properties.
As mentioned above, adequate micro resolution is crucial to ensuring an accurate result from the hybrid method; in this paper we use a trial-and-error approach. Three micro element configurations are considered: in the first, there are only two near-wall elements (i.e. P ¼ 2); in the second, a bulk element centred at y ¼ 0:5H is added (i.e. P ¼ 3); in the third, there are two bulk elements centred at y ¼ 0:33H and y ¼ 0:67H (i.e. P ¼ 4). Based on initial testing [24], the extents of the near-wall and bulk sampling regions are set to 15k l and 10k l , respectively; extents of 5k l are sufficient for the relaxation regions.
With tolerance values of f tolu ¼ 0:05 and f tol T ¼ f tolp ¼ f tolq ¼ 0:002, Fig. 6 shows that all hybrid property fields converge inside 3 iterations, for all three configurations.
It is difficult to assess the accuracy of the hybrid method from observing the final converged property fields in Fig. 7. Instead, the mean error in these fields, presented in Fig. 8 (for each  iteration), gives a much clearer indication of the method's accuracy. The final correction fields constructed by the hybrid method are presented in Fig. 9.
When simulating this case with only two near-wall elements (P ¼ 2), the mean errors in the final hybrid velocity and pressure fields are u ¼ 1:53% and p ¼ 4:47%, respectively. Adding micro elements in the bulk reduces these errors to approximately 0:1%, for both P ¼ 3 and P ¼ 4 -this is due to improved approximations of X xy and X yy , as indicated in Fig. 9(a) and (b). The largest difference in accuracy between the P ¼ 3 and P ¼ 4 configurations lies in the temperature solution, largely due to the capture of U y : the use of one bulk element results in T ¼ 2:42%, while the use of two bulk elements gives T ¼ 1:59%. This increase in accuracy might not be considered to be worth the increase in the computational cost; the cost of all three configurations is discussed in the following section.

Computational savings
A measure of the computational speed-up S provided by the hybrid method is given by the ratio of the total processing time required for the pure DSMC simulation to the total processing time required for the hybrid approach. In the hybrid approach, the time taken to compute the continuum solution is negligible compared with the computational time required by the DSMC micro simulations. The total processing time of each approach can therefore be estimated as [the total number of DSMC time steps G tot ] Â [the average clock time per DSMC time step t cl ], i.e. S ¼ G tot;pure t cl;pure G tot;hyb t cl;hyb : ð12Þ Note that G tot is the total number of time steps for both the transient and averaging periods. The total average clock time per time step for the hybrid approach t cl;hyb depends on both the micro element configuration and the number of iterations required: as the physical extent of an element can differ at each iteration, t cl;hyb is calculated as the sum of t cl for all micro elements h ¼ 1; 2; . . . P, over all iterations k ¼ 1; 2; . . . I, i.e.
In this paper, we focus on obtaining computational savings by exploiting spatial scale separation, and both the full and micro element simulations run for the same number of DSMC time steps; so for this Couette case, G tot;hyb ¼ G tot;pure ¼ 23 Â 10 6 . Table 1 presents the speed-up from each of the three micro element configurations considered. With convergence inside I ¼ 3 iterations for all cases, S is determined by the spatial extents of the near-wall micro elements, H lower and H upper , and the spatial extents of any bulk micro elements, H bulk . The extents presented in Table 1 are those at convergence (i.e. at k ¼ 3); note the changes from their original extents of 0.2 lm due to variations in k l during the simulation.
For this Couette flow test case, low spatial scale separation and the need to capture multiple constitutive correction fields means that the micro elements must occupy a significant portion of the flowfield. This makes the hybrid method more expensive than the pure DSMC simulation, i.e. S < 1, for Pi = 2, 3, and 4. Nevertheless, the purpose of this test case was to test/validate our extended coupling strategy in 1D, and simulations of larger 2D and 3D flow problems will highlight the computational advantages of the method (as demonstrated in the following section).

2D validation: flow through a microscale crack-type geometry
Channels with high-aspect-ratio cross-sections can be found in MEMS devices, and are also representative of narrow cracks. Consider the channel geometry shown in Fig. 10: the width W in the x-direction is much greater than the height H in the y-direction, i.e. W=H ) 1; the streamwise length L is then much greater than the width.
To date, the study of gas flows in such crack-type geometries has focused mainly on the continuum and slip regimes [26], with conventional fluid mechanics used to predict the flow behaviour. However, if the mean free path of the gas is of the same order as the crack opening (i.e. the height H), the flow will be nonequilibrium and predictions from conventional fluid mechanics will be poor. While the DSMC method could provide accurate flow predictions for such conditions, pure DSMC simulations are likely to be prohibitively expensive for very high cross-sectional aspect-ratios. Previous DSMC studies [27,28] have reduced the expense by assuming negligible variation of the flow in the direction of the width (essentially considering a 2D model in height and length), but this assumption will not hold true if there is significant geometric variation in the direction of the width. In reality, the cross-section of a crack will not be perfectly rectangular, for example, it could be more open at one end than the other; surface roughness or even surface defects could also have substantial effects on the flow behaviour.
A continuum-DSMC hybrid approach could be useful for modelling more realistic crack-type flows that see flow variation in the direction of the width. As the entire flowfield is 'near-wall', the popular DD framework would not be suitable. The point-wise coupling approach of the original HMM would be inefficient and inaccurate as the micro resolution required over the channel or crack height would force the micro elements to overlap. While the IMM could exploit high scale separation in the streamwise z-direction, it could not exploit scale separation in the direction of the width (the x-direction), and would require each micro element to occupy the entire cross-sectional area.
However, our HMM-FWC method has the potential to deal effectively with the mixed degrees of spatial scale separation that exist over the cross-section of these flows: it is able to exploit high scale separation in the x-direction by distributing the micro elements across the width; simultaneously, it can cope with the low scale separation in the y-direction by setting the micro elements to occupy the entire crack height. The higher the scale separation (i.e. the more gradual the variation of the flow properties) in the direction of the width, the more efficient our method will be.
To maintain simplicity and limit the expense of the pure DSMC simulation that we use for validation, we demonstrate the hybrid method's ability to cope with this mixed spatial scale separation by applying it only to the 2D cross-section of the crack-type channel. Periodicity is assumed in the streamwise z-direction and an external acceleration a z is used to drive the flow. Also for simplicity, we assume the channel cross-section to be perfectly rectangular, i.e. there is no geometric variation in the x-direction. Instead, a temperature difference is imposed between the left and right walls of the channel, with an associated linear temperature gradient imposed along the lower and upper walls. Although this is not necessarily representative of real crack flows, this test case enables us to consider variation of the flow properties over the crack width, while also enabling validation of our simultaneous momentum and heat flux coupling in 2D. More realistic geometric configurations should be considered in future work.

Applying the hybrid method in 2D
Compared with its application to 1D flow problems, applying our hybrid method in 2D presents additional challenges: the continuum conservation equations become more complex; the micro elements must be constrained by applying boundary conditions within 2D relaxation regions; and the interpolations of the constitutive correction fields between the sampling regions must be 2D.
The computational set-up for this 2D flow problem is shown in Fig. 11: a uniform continuum-fluid grid is imposed over the entire cross-section, consisting of M x macro nodes across the width and M y macro nodes over the height (including nodes at the bounding walls); the horizontal and vertical node spacings are Dx and Dy, respectively. The micro elements span the full channel height H. The element arrangement shown is an example, consisting of a 'side-wall' element at each of the left and right walls, and one element in the 'bulk'. All sampling regions are discretised into a grid of 2D measurement bins, while all relaxation regions are split Iteration, k into a grid of 2D control bins. Here, the measurement and control bins are set to have the same width dx and height dy.
A full description of the hybrid algorithm for this nonisothermal test case is given in B 5 . With subsonic flow and periodicity in the streamwise z-direction, mass conservation is satisfied automatically, and the transport of momentum in the x-and y-directions is assumed negligible, i.e. the gas pressure and density are assumed constant over the flow cross-section. Once again, constant reference values for the transport properties are adopted (l ref and j ref ) with the constitutive corrections X zx ; X zy ; U x , and U y capturing the true variation of these properties as the iterations proceed.
In this flow case, constraining the micro elements involves applying the external acceleration a z to the particles, while imposing only the velocity and temperature fields across each relaxation region. After executing DSMC and extracting the required information from the local sampling regions, we approximate the global correction fields and the boundary information at all bounding walls by using simple linear interpolations across the channel width.
A corrected streamwise velocity field w corr and a corrected temperature field T corr are produced by the hybrid method. For this 2D problem, the convergence of each property s (i.e. s ¼ w; T) is monitored at each iteration k according to, Both w and T are expected to vary considerably over the crosssection, and so the denominator of this expression is the flow property range in the pure DSMC solution. As mentioned above, the magnitude of these errors is determined by the accuracy of the correction fields constructed by the method, which should ideally match the pure DSMC correction fields (which are calculated by using the pure DSMC property fields in the corrected constitutive relations).
When investigating fluid flows in cracks and general microchannels, it is often important to predict the mass flow rate _ m through the crack or along the channel. Here, the mass flow rate predicted by each approach (i.e. the NSF approach, the hybrid method, and the pure DSMC simulation) is approximated using a surface integral, where i=1; 2; . . . M x and j=1; 2; . . . M y . As well as using the different velocity fields from each solution (w NSF ; w hyb ¼ w corr , and w pure ), this mass flow rate calculation also requires the different density fields: q NSF and q hyb are simply equal to the average density q av everywhere, while q pure is the density field predicted in the pure DSMC simulation. For this low-speed test case, q pure will vary only slightly from the value of q av .

Results
We consider a geometry where the width W ¼ 200 lm and the height H ¼ 0:1 lm, giving a cross-sectional aspect-ratio of 2000.
Greater macroscopic resolution is needed in the y-direction (height) than in the x-direction (width); M y ¼ 26 macro nodes over the height (including nodes at the lower and upper walls) give a vertical node spacing Dy ¼ 4 nm, while M x ¼ 801 nodes across the width (including nodes at the left and right walls) give a horizontal spacing Dx ¼ 250 nm. The micro elements of the hybrid method capture an even higher resolution: the height of each measurement/control bin dy is equal to the vertical macro node spacing of 4 nm, but the bin width dx is set to 50 nm. In all the micro DSMC simulations in the hybrid method, and in the validating pure DSMC simulation, we apply an acceleration a z ¼ 10 Â 10 10 m/s 2 to each particle to generate flow in the z-direction -this corresponds to a body forcing f z ¼ a z q ¼ 8:5 Â 10 10 N/m 3 . Note that it is common practice to apply such large, unrealistic forcing in molecular-type simulations (up to the point that the response behaves nonlinearly) so that the resulting signal-to-noise ratio is manageable. A time step dt ¼ 5 Â 10 À12 s is used in all the DSMC simulations (micro and pure), and steady state is reached inside a start-up period of 0.6 million time steps, while a further 1.4 million time steps reduces scatter in the measured flow properties.
The no-slip/no-temperature-jump NSF description of this flow problem predicts a streamwise velocity field w NSF with a maximum of approximately 5 m/s (i.e. Ma % 0:016); the effects of the temperature difference are not incorporated and so this maximum is in the centre of the cross-section, as shown in Fig. 12(a). The pure DSMC simulation, on the other hand, predicts a velocity field w pure with a much larger magnitude and a skew: Fig. 12(b) shows a maximum of approximately 55 m/s (i.e. Ma % 0:18) at the colder left wall. The NSF temperature field T NSF assumes a linear variation in the x-direction, as shown in Fig. 12(c). Although it is difficult to observe from Fig. 12(d), the pure DSMC temperature field T pure features small thermal Knudsen layers at the left and right walls.
The micro element configuration in the hybrid method for this flow case consists of three micro elements (i.e. Pi = 3): two side-wall elements, and one bulk element at x ¼ 0:5W. Following some initial testing [24], horizontal extents of 5k l were deemed   sufficient for both the sampling and the relaxation regions, with both adapting dynamically in spatial extent at each iteration depending on k l . Adopting tolerance values of f tolw ¼ 3 Â 10 À3 and f tol T ¼ 3 Â 10 À4 , the velocity and temperature fields produced from this hybrid configuration reach convergence inside 6 iterations; see Fig. 13.
The final hybrid property fields are compared with the pure DSMC property fields by observing the profiles across the channel width at the centreline of its height (i.e. at y ¼ 0:5H), and over the height at the centreline of the width (i.e. at x ¼ 0:5W); the velocity profiles are displayed in Fig. 14(a) and (b), while the temperature profiles are shown in Fig. 14(c)   μm) x (μ m) y (10 -2 μm) x (μ m) Fig. 12. Non-isothermal crack-type flow: surface plot of (a) the NSF velocity field wNSF, (b) the pure DSMC velocity field wpure, (c) the NSF temperature field TNSF, and (d) the pure DSMC temperature field Tpure.
hybrid prediction agrees almost exactly with the pure DSMC solution -the method predicts the skew of the velocity in the xdirection, and the small thermal Knudsen layers at the left and right walls. The method's accuracy is more apparent from the mean percentage error (with respect to the pure DSMC solution) of the predicted property fields. The mean error w in the velocity field at each iteration is shown in Fig. 15(a): while the NSF velocity solution has an error of 97.59%, the hybrid method achieves an error of only 0.81% inside 5 iterations. Fig. 15(b) shows the mean error T in the temperature field at each iteration; although the NSF temperature field has an error of only 3.11%, the hybrid method is able to reduce this to 0.25% inside 6 iterations.
The hybrid method also provides an excellent approximation to the mass flow rate. Based on the pure DSMC velocity and density fields, Eq. (16) gives a prediction of _ m pure ¼ 73:97 Â 10 À11 kg/s. Compared with this value, the hybrid method's final prediction of _ m hyb ¼ 74:1 Â 10 À11 kg/s has an error of only 0.18%. On the other hand, the NSF prediction of _ m NSF ¼ 5:74 Â 10 À11 kg/s equates to an error of 92.24%.
The high accuracy of the hybrid solution reflects the method's ability to approximate all four constitutive correction fields, even with only simple linear interpolations between the sampling regions. In Fig. 16 the final correction fields constructed using the hybrid method are compared with the equivalent fields from the pure DSMC solution (with the figure showing profiles for all four corrections across the width at y ¼ 0:5H, and over the height at x ¼ 0:5W). Note that substantial scatter makes it difficult to observe the true variation of the heat flux correction fields. Scatter  in the bulk micro element results in a discontinuity in the hybrid method's approximation to U x near the centre of the channel width (see Fig. 16(c)) but the effects of this discontinuity on the other flow properties are not significant.

Computational savings
As with the Couette flow problem, the computational speed-up S provided by the hybrid method for this 2D case can be computed from Eqs. (12) and (13). Once again, the same number of DSMC time steps were performed in both the hybrid micro DSMC simulations and the pure DSMC validation simulation, i.e. G tot;hyb ¼ G tot;pure ¼ 2 Â 10 6 . Table 2 shows that, with convergence in I ¼ 6 iterations, exploitation of spatial scale separation enables the hybrid method to achieve a speed-up of nearly 10Â. Note that the element extents stated in this table are those at convergence, i.e. at iteration k ¼ 6.
While the hybrid configuration adopted here includes a bulk element, two near-wall elements with linear interpolations of the constitutive and boundary information may have been sufficient: the temperature gradient imposed on the upper and lower walls is linear and so (for the temperature range considered) the transport properties vary approximately linearly over the width. Bulk elements do, however, provide essential micro resolution if there are surface defects, or if the cross-section is non-rectangular; a more extreme temperature range might also require bulk micro resolution.

Conclusions and future work
We have extended the coupling strategy of the HMM-FWC (continuum-DSMC) method of [1] so that the transport of mass, momentum, and heat can now be modelled. Combined with the generality of the HMM-FWC approach itself, this extension makes the method suitable for a wide range of flows, including compressible and non-isothermal flows, and flows with mixed levels of spatial scale separation.
Validation of this extended coupling strategy has been demonstrated in 1D: for a high-speed micro Couette flow problem, the method (with sufficient micro resolution) achieved excellent agreement with a pure DSMC solution. Although low spatial scale separation meant that there was no computational speed-up, the  purpose of this 1D case was simply to test and validate the coupling strategy. The method has then been implemented in 2D, for a simplified model of gas flowing through a microscale crack. A temperature difference between the walls of the cross-section required simultaneous momentum and heat flux coupling in 2D; the resulting predictions of the property fields and the mass flow rate were within 1% of the pure DSMC benchmark solution. As well as demonstrating validation in 2D, this crack flow test case highlighted the method's ability to cope efficiently with mixed spatial scale separation, and a computational speed-up of almost 10Â has been achieved. This work opens up a number of directions for future work. For the crack-type problem, it would be interesting to consider a non-rectangular cross-section, or the presence of surface defects/ roughness. A 3D model could employ the HMM-FWC method along the streamwise length. Alternatively, if the streamwise scale separation was very high, then significant computational savings would be achieved by combining the HMM-FWC method with the IMM: the HMM-FWC would comprise the 'micro elements' spaced intermittently along the streamwise direction of an IMM implementation.
Development of the hybrid method itself could focus on the simulation of polyatomic gases and gas mixtures, and perhaps on the simulation of external high-speed, high-altitude gas flows (for example, around hypersonic space access vehicles); for very high temperatures, this would require models for the chemical reactions. An extension to unsteady flows could be achieved using a similar approach to that of Borg et al. [13].
From Eq. (A.5), the normal stress correction in each bin or at each wall is simply equal to the normal stress in that bin or at that wall. (5) Linear interpolations of the correction fields are performed between the sampling regions in order to approximate X xy ; X yy , and U y everywhere in the system. Interpolations of the gas boundary information are not required for this 1D problem. The boundary conditions are updated to the gas velocities, pressures, and temperatures measured at the bounding walls during Step (3). Solution of these three equations then provides the new, corrected streamwise velocity field u corr , pressure field p corr , and temperature field T corr , respectively.
The new density field q corr is computed from the equation of state, i.e. q corr ¼ p corr =RT corr .
(7) Replacing the NSF fields with these corrected fields, the process is repeated from Step (1); iterations continue until all fields converge within user-defined tolerances.