Multiscale simulation of heat transfer in a rarefied gas

We present a new hybrid method for dilute gas ﬂows that couples a continuum-ﬂuid description to the direct simulation Monte Carlo (DSMC) technique. Instead of using a domain-decomposition framework, we adopt a heterogeneous approach with micro resolution that can capture non-equilibrium or non-continuum ﬂuid behaviour both close to bounding walls and in the bulk. A continuum-ﬂuid model is applied across the entire domain, while DSMC is applied in spatially-distributed micro regions. Using a ﬁeld-wise coupling approach, each micro element provides a local correction to a continuum sub-region, the dimensions of which are identical to the micro element itself. Interpolating this local correction between the micro elements then produces a correction that can be applied over the entire continuum domain. Key advantages of this method include its suitability for ﬂow problems with varying degrees of scale separation, and that the location of the micro elements is not restricted to the nodes of the computational mesh. Also, the size of the micro elements adapts dynamically with the local molecular mean free path. We demonstrate the method on heat transfer problems in dilute gas ﬂows, where the coupling is performed through the computed heat ﬂuxes. Our test case is micro Fourier ﬂow over a range of rarefaction and temperature conditions: this case is simple enough to enable validation against a pure DSMC simulation, and our results show that the hybrid method can deal with both missing boundary and constitutive information.


Introduction
While the conventional hydrodynamic equations are generally excellent for modelling the majority of fluid flow problems, the presence of localised regions of non-continuum or non-equilibrium flow can result in some degree of inaccuracy. Such regions appear when the flow is far from local thermodynamic equilibrium, for example, when there are large gradients in fluid properties, or when surface e↵ects become dominant. Although molecular simulation tools can provide an accurate modelling alternative in these cases, they are usually much too computationally expensive for resolving engineering spatial and temporal scales. Multiscale methodologies that exploit 'scale separation' have therefore been developed over the past decade. Scale separation occurs when the variation of hydrodynamic properties across small regions of space or periods of time is only very loosely coupled with the flow behaviour on a much larger spatial or temporal scale.
Often referred to as 'hybrids', these multiscale methods combine continuum and molecular descriptions of the flow. A traditional continuum description is employed in macro flow regions, and a molecular treatment is applied in small-scale micro or nano regions. Essentially, the aim is to combine the best of both solvers: the computational e ciency associated with continuum methods, and the detail and accuracy of molecular techniques.
In the literature, two di↵erent hybrid frameworks have emerged for fluid flows: a) the domain-decomposition technique, and b) the Heterogeneous Multiscale Method (HMM). For liquids, molecular dynamics (MD) is the appropriate molecular simulation tool. This deterministic method is, however, ine cient for dilute gas flows, and the direct simulation Monte Carlo (DSMC) method [1] can instead provide a coarse-grained molecular description. Founded on the kinetic theory of dilute gases, DSMC reduces computational expense by adopting a stochastic approximation for the molecular collision process.
Typically, thermodynamic non-equilibrium e↵ects occur in the vicinity of bounding surfaces or other interfaces. Recognising this behaviour, domain-decomposition has become the most popular hybrid framework for both liquids [2,3,4,5,6] and dilute gases [7,8,9,10,11,12,13,14,15,16,17]. In this method, the simulation domain is partitioned -a molecular solver is applied in the regions closest to the surfaces, while a conventional continuum fluid solver is implemented in the remainder. These micro and macro sub-domains are independent but communicate through an overlap region that enables mutual coupling. This coupling is typically established by matching fluxes of fluid mass, momentum, and energy, or by matching state properties. However, despite its popularity, a fundamental disadvantage of this micro-macro decomposition approach is that computational e ciency can be increased above that of a full molecular simulation only when non-continuum flow is confined to 'near-wall' regions. In highly non-equilibrium flows, or when the temperature dependence of fluid properties such as dynamic viscosity or thermal conductivity is not known a priori (for example, in unusual chemically-reacting gas mixtures), the traditional fluid constitutive relations may be inaccurate in the bulk of the domain.
Domain-decomposition techniques are also inappropriate for simulating flow through micro or nanoscale geometries that have a high aspect ratio, e.g. one dimension of the geometry is orders of magnitude larger than another. This class of flow presents a challenge as it requires simultaneous solution of the microscopic processes occurring over the smallest dimension and the macroscopic processes occurring over the largest dimension, and is often too computationally intensive for a full molecular approach. Such flows are generally beyond the reach of domain-decomposition as the majority (or perhaps all) of the flowfield can be considered 'near-wall'.
The less-common HMM framework overcomes the limitations of domain-decomposition by adopting a micro-resolution approach that can be employed anywhere in the domain -near bounding surfaces, or in the bulk flow. In this case, a continuum model is applied across the entire flow field, and the molecular solver is applied in spatially-distributed micro regions. These micro regions provide the missing data that is required for closure 2 of the local continuum model, either in the form of unknown boundary conditions, or unknown constitutive information. Existing HMM studies in the literature are mainly for liquid flows, using MD as the molecular solver [18,19,20,21]. Generally, these studies consider flow problems where momentum transport is dominant and the transfer of heat is negligible. Coupling is therefore based on momentum: velocity fields or strain-rates are prescribed in each micro region, and the resultant stress is used to apply a correction back into the hydrodynamic momentum equation [18,19]. Each HMM micro region supplies information to a computational node on the continuum mesh. This point-wise coupling approach [20] is ideal when there is a large degree of spatial scale separation in the system, providing significant computational savings over a pure molecular simulation. However, in flow problems with smaller, or mixed, degrees of spatial scale separation, the molecular resolution required can result in the micro elements overlapping, making HMM more expensive and less accurate than a pure molecular treatment. Despite its advantages over domain-decomposition techniques, there has been little development of HMM-type hybrids that use DSMC as the molecular treatment. In 2010, Kessler et al. [22] proposed the Coupled Multiscale Multiphysics Method (CM 3 ) that couples both momentum and heat transfer, with DSMC providing corrections to both the hydrodynamic momentum and energy equations. This method was, however, developed to simulate transient flows where a time-accurate solution is sought, and so any computational advantage over a full DSMC simulation is achieved only through decoupling of the time scales; the length scales remain fully coupled, with both the continuum description and DSMC employed over the same region of space. More recently, Patronis et al. [23] adapted the Internal-flow Multiscale Method (IMM) to simulate dilute gas flows with DSMC. Originally developed by Borg et al. [21] for liquid flows, IMM adopts a framework similar to HMM but is tailored to model flows in high-aspect-ratio channels. Although large computational savings are presented, this method is tailored specifically to deal with cases where the length scale in the direction of flow is significantly larger than the length scales transverse to the flow direction.
In this paper we propose a new form of the HMM technique, with DSMC providing the molecular description. We adopt the field wise coupling (HMM-FWC) approach developed by Borg et al. [24] for liquid flows: rather than supplying a correction to a node on the continuum mesh, each micro element instead corrects a continuum subregion, the spatial dimensions of which are identical to those of the micro element itself. This means that, unlike point wise coupling, HMM-FWC is suitable for dealing with flow problems with varying degrees of spatial scale separation. Also, the location of each micro element is not restricted to the nodes of the computational mesh: both the position and size of the micro elements can be optimised for each problem, independently of the continuum mesh.
Our method is designed to cope with inaccuracy in the traditional flow boundary conditions and/or constitutive relations, and is therefore able to deal with problems which are beyond the reach of domain-decomposition. This includes problems where the fluid behaviour is unknown in the bulk, i.e. the traditional constitutive relations fail due to non-equilibrium e↵ects (for example, in the wake of a re-entry vehicle), or the transport properties are unknown (for example, in unusual gas mixtures). The method is also suitable for simulating high aspect ratio geometries. While the IMM is designed to simulate problems where the largest length scale is in the flow direction, our new method has no such restriction and so provides a more general approach. It could therefore be 3 useful when the largest length scale is transverse to the flow direction; for example, the flow through microscale cracks in valves. The form of the method we present in this paper is tailored to model heat transfer problems in dilute gas flows. As a starting point we consider problems in which the gas is essentially motionless, with large applied temperature gradients placing the focus on heat transfer. With negligible transport of momentum, our coupling is performed through the heat flux: we impose the local temperature fields on the micro elements and measure the consequent heat flux from the relaxed DSMC particle ensembles. A suitable correction is then applied to the hydrodynamic conservation of energy equation. (Full coupling of mass, momentum, and heat transfer is a subject for future work.) The level of translational non-equilibrium in a rarefied gas is generally characterized by the Knudsen number Kn, defined as the ratio of the gas molecular mean free path to a characteristic system dimension L. Typically, the traditional 'no-temperaturejump' boundary condition at a bounding surface (wall) is only valid when Kn < 0.001; the flow is then in thermodynamic equilibrium as the frequency of both intermolecular and molecule-wall collisions is very high. As Kn increases above 0.001, this collision frequency decreases, resulting in a temperature discontinuity between the wall and its adjacent gas. For low Kn, the conventional conservation equations can be extended to account for this by employing von Smoluchowski temperature-jump boundary conditions [25]. However, as Kn increases, molecule-wall collisions become more frequent than intermolecular collisions and a thermal Knudsen layer develops. This is essentially a region of non-equilibrium that extends from the wall into the domain, with its thickness determined by the degree of rarefaction. In this layer, the gas behaviour deviates from the conventional linear heat flux/temperature-gradient constitutive relation. Even with the use of temperature-jump boundary conditions, the conventional hydrodynamic equations cannot model this phenomenon. Our hybrid method, however, aims to capture both the temperature-jump and the thermal Knudsen layer at a lower computational cost than a full-domain DSMC treatment.
This paper is organised as follows. In Section 2, we discuss the macro and micro descriptions of the flow. We then present the general three-dimensional form of our multiscale coupling framework and its iterative algorithm. For simplicity, one-dimensional heat transfer is considered when validating the method in Section 3, with results compared with pure DSMC simulations at equivalent conditions. In Section 4 we draw our conclusions.

Continuum description
For simplicity in this paper, we restrict our attention here to steady-state problems where the gas remains stationary, i.e. it has no streaming velocity. With negligible transport of mass and momentum, the continuum description is based on the conservation of energy, where q is the heat flux vector. To close this equation, a suitable constitutive relation is required. For conventional gas flows, the traditional Navier-Stokes-Fourier (NSF) constitutive relations are accurate and so Fourier's law can be used, i.e.
where  is the thermal conductivity of the gas. The energy equation then reads, Assuming no-temperature-jump boundary conditions, solution of this conventional energy equation will produce a continuum NSF temperature field T NSF . This field is taken as an initial condition for our hybrid method, and provides a starting point from which the coupling framework iterates towards the correct temperature field. While Eq. (2) is valid for typical flow problems, it fails in certain flow conditions, for instance, in flows of complex fluids or in conditions of thermodynamic non-equilibrium. Inaccuracy in the conventional constitutive model can therefore be quantified by a heatflux-correction field , so that The flux-correction field incorporates not only the departure of the gas state from equilibrium, but also any additional inaccuracy due to the assumed thermal conductivity model. Using this relation to close Eq. (1) results in a 'flux-corrected' energy equation, The general strategy of our heterogeneous hybrid approach is therefore as follows. Across each individual micro element, the heat flux and temperature fields are measured. The flux-correction field across each element is then computed using Eq. (4). By interpolating between all micro elements, the full flux-correction field across the entire flowfield is approximated. With this, and the boundary information obtained from micro simulations located at the bounding walls, an appropriate continuum method (e.g. finite di↵erence, finite element, or finite volume) can be used to solve Eq. (5). This produces a flux-corrected continuum temperature field T across the domain. With continuing iterations, this temperature field should converge towards that which would be obtained from a full molecular simulation of the problem.

DSMC technique
As we focus on heat transfer in dilute gases, we use DSMC as our molecular model in the micro elements. DSMC has become the dominant method for simulation of dilute gas flows that lie in the continuum-transition regime. The fundamental concept is to track a large number of numerical particles, storing their position, velocity, and internal state as they move through a computational mesh. During a simulation, the particles collide with each other and with bounding surfaces while maintaining conservation of mass, momentum, and energy. A major advantage of the method is that two key approximations significantly reduce computational expense: (a) each simulated particle typically represents a large number of real gas molecules, and (b) molecular motion and intermolecular collisions can be decoupled over small time intervals. Particle movements are computed deterministically, while interparticle collisions are treated statistically within numerical mesh cells. 5 Using expressions provided by Bird [1], local hydrodynamic properties (including the temperature and the heat flux) can be recovered in DSMC by averaging microscopic data over all of the particles in each cell. However, the inherent statistical scatter associated with the method means that a large number of independent samples are usually required to capture smooth fields, particularly for low speed flows. Therefore, time averaging is typically used for steady state problems, while ensemble averaging is used for transient flow problems. For continuum-molecular hybrids, the smoothing of statistical fluctuations is particularly important as the transfer of noisy data to the continuum fluid description may result in instability. Consequently, each DSMC simulation in this paper is performed in two stages: a transient period enables the simulation to reach steadystate, then a longer averaging period reduces the statistical scatter in the measured hydrodynamic properties.
The open-source C++ toolbox OpenFOAM [26] incorporates a DSMC solver, dsmc-Foam, which has been validated for various benchmark cases, including hypersonic and microchannel flows [27,28,29]. This is used to perform all the micro and full-scale DSMC simulations in this paper.

Coupling framework
In our framework, the continuum description is applied across the entire macro domain, while DSMC is performed in dispersed micro elements (the arrangement of which depends on the flow problem being investigated). As discussed in the Introduction, we propose a new strategy to achieve two-way coupling between these domains. The underlying methodology presented here is general and may be applied to 1D, 2D, or 3D problems.

Macro-to-micro coupling: constraint of the micro elements
It is crucial that the flux-correction information is extracted from DSMC elements in which the gas state is properly representative of the local conditions in the macro domain. There is, however, a fundamental problem in creating such elements: the particle distribution required at the boundaries of the element cannot be extracted directly from the macroscopic continuum flowfield. For this reason, macro-to-micro coupling remains a challenge in multiscale hybrids. In our method, we circumvent this problem by not imposing an exact particle distribution at these boundaries; instead, we allow the distribution to be dictated by the macroscopic continuum fields alone.
A natural evolution to the correct particle distribution and macroscopic state at the boundaries of an element is possible by establishing an artificial relaxation region (zone) around the element. It is not essential that the gas state in these relaxation zones accurately represent the conditions in the corresponding macro domain; their sole purpose is to develop the boundary conditions for the core of the element. Sampling of property fields is performed only in the core region, which we refer to as the 'sampling zone'.
Suitable boundary conditions to the sampling zone are generated by enforcing the continuum variation of macroscopic properties (i.e. in our case, the continuum temperature field) across the surrounding relaxation zone. This is done by implementing particle controllers [30]. A particle distribution does, however, need to be applied at the outer boundaries of the relaxation zone. This imposed distribution introduces error in the local distribution close to these boundaries, so it is important that the particle state within the sampling zone is su ciently independent of this introduced error. If the relaxation zone is large enough (i.e. several molecular mean free paths), the form of the imposed distribution 1 is irrelevant -its e↵ect will decay through the zone as the particle state gradually relaxes through interparticle collisions. Complete relaxation needs to occur over the extent of the relaxation zone such that the particle distribution in the sampling zone is dictated solely by the applied continuum state. Once the hybrid method has converged, the artificiality of the outer relaxation zone dissolves seamlessly into the true particle distribution and fluid state in the inner sampling zone. Figure 1(a) shows a typical 'bulk' micro element in 2D. Similarly, Fig. 1(b) shows an example computational domain set-up for a 2D problem, including both bulk and 'nearwall' micro elements. In summary, the sampling zone in a micro element is surrounded by the relaxation zone over which the local continuum fields are imposed, and a chosen particle distribution is applied at the outer boundaries. However, in order to capture regions of non-equilibrium that often appear at solid bounding surfaces or walls, the sampling zone in a near-wall micro element must be adjacent to the wall itself. The arrangement of the micro elements is clearly crucial to the resulting accuracy of our method, and the appropriate arrangement will vary from problem to problem.
To impose the local continuum temperature field, we divide the relaxation zone in each micro element into a grid of control cells as shown in Fig. 1. Using particle controllers [30], the appropriate temperature is then set in each control cell to create the desired temperature variation. Similarly, the sampling zone of each micro element is divided into a grid of measurement cells. We then extract the averaged hydrodynamic properties from each measurement cell, enabling us to capture the property fields across the entire sampling zone. Boundary information is also obtained from these measurement cells: we extract the temperature of the gas in contact with the wall from the wall-adjacent cell faces.
Essentially, these measurement and control cells form a measurement/control mesh that is completely independent of the computational mesh used by DSMC itself. This has the benefit of enabling us to define the resolution of the macroscopic fields that we both extract and impose, and to control the noise in our measurements, without a↵ecting the accuracy of the DSMC calculations (i.e. the particle collision rate). In Fig. 1, the measurement and control cells are shown to have the same dimensions, but this does not have to be the case. Also, in Fig. 1(b), the measurement and control cells are shown to be collocated with the computational cells of the continuum mesh; although this simplifies the transfer of data between the micro elements and the continuum description, it is only for convenience and is not essential. If these cells are not collocated, we simply interpolate between the meshes.

Micro-to-macro coupling: correcting the continuum description
The ease in converting microscopic particle information into macroscopic fields means that, generally, micro-to-macro coupling is less problematic than macro-to-micro coupling. Based on the property fields extracted from the DSMC solver, a suitable correction can be applied to the continuum description. In our coupling strategy, this correction is in the form of a heat-flux-correction field , as discussed in section 2.1. In our DSMC simulations, the average translational temperature (for a single species gas) in each measurement cell is computed from, where k B is the Boltzmann constant, m is the molecular mass and c 0 x , c 0 y , and c 0 z are the x, y, and z components of the thermal velocity vector c 0 . Assuming the gas has no vibrational energy, the average heat flux q in each measurement cell is obtained from, where j represents the x, y, and z components, n is the number density of the gas, and " rot is the rotational energy of a single molecule. Note that we consider monatomic gas flows in this paper, and so " rot = 0 and the macroscopic temperature T = T tr . By computing the average macroscopic temperature and heat flux in each measurement cell, we capture the temperature and heat flux fields across the sampling zone of each micro element. Substituted into Eq. (4), these then provide the flux-correction field across this zone. However, the continuum description is applied across the full simulation domain, and so Eq. (5) requires the flux-correction field across the full domain. This field can be approximated using appropriate interpolation between all the sampling zones.
This hybrid methodology not only corrects for inaccurate constitutive information, but also provides missing boundary information, i.e. the temperature jump. In each near-wall sampling zone, we measure the temperature of the gas in contact with the wall by summing over all particles that strike the surface [31], i.e.
where c n is the particle velocity normal to the wall, c t is the particle velocity tangential to the wall, kck is the velocity magnitude, and U slip is the slip velocity which, with zero wall velocity, is given by We then interpolate this boundary gas temperature between the near-wall sampling zones to produce an estimate of the gas temperature at all bounding surfaces. Using this boundary information and the full flux-correction field , solution of Eq. (5) then produces a new (flux-corrected) temperature field T across the whole domain. The process is then repeated until T relaxes to a solution that is close to that obtained from a full DSMC simulation.

Iterative algorithm
The general iterative coupling procedure of this method is: (0) Assuming no-temperature-jump at bounding surfaces, solve the conventional energy equation (3) to obtain an initial estimate for the temperature field across the entire simulation domain, T NSF .
(1) Constrain each micro element by applying boundary conditions: (a) Using particle control in each control cell, enforce the local continuum temperature variation throughout the relaxation zone.
(b) At the outer boundaries of the relaxation zone, impose Maxwellian particle distributions at the local continuum temperature.
(2) Execute DSMC in each micro element as described in section 2.2. When steadystate is reached, perform averaging of properties in all measurement cells across the sampling zone of each element. Extract the temperature and the heat flux values from each of these cells. Also, from each wall-adjacent measurement cell face, extract the temperature of the gas at the wall surface.
(3) Compute the flux-correction in each measurement cell using the temperature and heat flux values extracted from that cell in the flux-corrected constitutive relation (4). From this, obtain the flux-correction field across each sampling zone.
(4) Carry out appropriate interpolations between sampling zones to approximate the flux-correction field across the full simulation domain. Similarly, perform appropriate interpolations between near-wall sampling zones to obtain a gas temperature at all bounding walls.
(5) Using the boundary gas temperature information and the full flux-correction field, solve the flux-corrected energy equation (5) across the full domain to obtain a new flux-corrected temperature field T .
(6) Repeat from Step (1) until T does not change between iterations to within a user-defined tolerance.

Validation: one-dimensional Fourier flow
A simple validation test case is chosen in order to assess the performance of the coupling method, and to easily validate it against a full DSMC simulation. We therefore demonstrate the method on the case of one-dimensional heat transfer, i.e. the classical Fourier flow problem. This problem has a motionless gas (in this case argon) confined between two infinite parallel planar walls at di↵erent temperatures, T cold and T hot .

Numerical implementation
We assume that the variation of the gas thermal conductivity  with temperature is unknown. We therefore take a reasonable reference value  r that is independent of the temperature field across the system domain; as discussed in section 2.1, the fluxcorrection field will automatically adjust for any error resulting from this assumption.
Discretizing the domain in one-dimensional space, the macroscopic computational mesh consists of N macro nodes, including a node at each wall. An example computational domain set-up is shown in Fig. 2, where we have a micro element at each wall, and one in the bulk. This element arrangement is an example only -the appropriate arrangement will in practise depend on the case itself, as will be discussed further in section 3.2.1. For this 1D problem, each near-wall element comprises a single sampling zone and a single relaxation zone, while each bulk element consists of a single sampling zone with a relaxation zone on either side, as indicated in Fig. 2.
Each sampling zone is divided into a number of measurement bins. Similarly, each relaxation zone is divided into a number of control bins. In order to keep the transfer of data between the macroscopic mesh and the micro elements as simple as possible here, the bin arrangement in each element is set such that the centre of each bin coincides exactly with a macro node, i = 1, 2, ..., N . Again this is merely for convenience, and does not generally need to be imposed. The length of each 1D bin is then equal to the spacing between each macro node, x.
Application of the general coupling algorithm of section 2.4 to this 1D system is as follows: Sampling zone (0) An initial estimate for the temperature field across the simulation domain is computed by solving the 1D conservation of energy equation, where q x is the streamwise component of the heat flux vector. This expression can be closed using the 1D form of Fourier's law, Using the reference thermal conductivity  r , the 1D energy equation becomes, and can be approximated using a second order central finite di↵erence scheme. With no temperature-jump at the solid walls (i.e. T 1 = T cold and T N = T hot ), solution results in the initial continuum temperature field T NSF .
(1) As the micro element size can change at each iteration (as will be discussed in section 3.2.1), each element is initialised at equilibrium before boundary conditions are applied. This is done by sampling particle velocities from a Maxwellian distribution. For consistency between the micro and macro domains, the temperature and density for initialisation are obtained from the local continuum solution: continuum values of these properties are averaged over a sub-region of the macro domain that corresponds to the particular micro element.
Each micro element is then constrained: (a) The local continuum temperature field is imposed throughout each relaxation zone by implementing a thermostat in each control bin.
(b) Maxwellian particle distributions are imposed at the outer boundaries of each relaxation zone via a di↵use solid wall condition at the local continuum temperature.
(2) DSMC is executed in each micro element. When steady-state conditions are reached, time averaging of properties is performed in each measurement bin. From each bin, time-averaged values of the temperature and the heat flux are extracted. Also, from both near-wall sampling zones, the temperature of the gas in contact with the bounding wall is extracted.
(3) Using the temperature and the heat flux extracted from each measurement bin, the flux-correction in each bin can be computed by applying the 1D flux-corrected constitutive relation with constant conductivity  r , For this calculation, the temperature gradient in each measurement bin is approximated using a central finite di↵erence scheme (or a forward/backward finite di↵erence scheme for bins at the edges of the sampling zone). By computing the flux-correction in each measurement bin, the flux-correction field across each sampling zone is captured.
(4) To obtain x at every point in the domain, interpolation between sampling zones is required. A study was performed showing that a simple linear interpolation provides an adequate representation of the full flux-correction field for this 1D case. Note that, for this 1D geometry, interpolation of the boundary gas temperature is not required.
which can be approximated using a second order central finite di↵erence scheme. As discussed in section 2.3.2, temperature-jump boundary conditions are applied: the temperature of the gas at each bounding wall is set to be that measured in the corresponding near-wall element during Step (2). Solving this equation then results in the flux-corrected temperature field T across the entire domain.
(6) Replacing the initial temperature field T NSF with T , the process is repeated from Step (1). This continues until T converges to within a user defined tolerance, i.e.
where N is the number of macro nodes, l is the iteration index, and ⇣ tol is the tolerance parameter.

Results
The overall level of non-equilibrium in the gas is characterized here by a global Knudsen number Kn, where the characteristic dimension is the separation between the heated walls L, i.e.
We use a Variable Hard Sphere (VHS) collision model in our DSMC simulations, and so the gas mean free path is given by [1], where d is the VHS molecular diameter, and n is the number density of the gas. Localised regions of non-equilibrium can also result from large gradients in the macroscopic fluid properties, for example, temperature gradients. Varying both Kn and the temperature gradient across the system, we simulate a range of test cases here. These explore the ability of the hybrid method to deal with both missing constitutive and boundary information, while remaining simple enough for validation against an equivalent full-scale DSMC simulation. For simplicity, we simulate monatomic argon gas with the following VHS parameters: a reference temperature T ref = 273 K, a reference diameter d ref = 4.17⇥10 10 m, and viscosity exponent ! = 0.81. For all test cases, the separation between the walls L is 1 µm. In order to accurately capture the variation of the property fields, we set the number of macro nodes N across the domain to be 201 for all cases, making the constant spacing x between each macro node 5 nm. For each test case, we then set the gas density and the temperature of the walls to obtain the desired Kn and temperature gradient.
To ensure a fair comparison between each hybrid simulation and the equivalent fullscale DSMC simulation, we use the same cell-size and time-step in both the hybrid DSMC elements and the full-scale simulation. The cell size is set as a fraction of the gas mean free path . The DSMC time-step must be a fraction of the mean collision time t mc : a time-step t = 1⇥10 12 s is su ciently small for all cases simulated in this paper. In our testing, we found that an initial start-up run of 3 million time-steps allowed all DSMC simulations to reach steady state. To minimise the statistical scatter associated with the averaging of the field properties, all simulations were then run for a further 50 million time-steps.
As discussed in section 3.1, we monitor convergence of the hybrid procedure using Eq. (15), which quantifies the di↵erence in the temperature solution between the current iteration and the previous. The tolerance parameter ⇣ tol depends on the case itself, with typical values of O(10 2 ) and below. We then require a measure of the overall accuracy for the converged hybrid solution when comparing with the corresponding fullscale DSMC solution T Full . Here we consider the mean percentage error✏ in the hybrid temperature profile, i.e.✏ = 1 The value of this error depends on the hybrid's ability to capture the 'true' flux-correction field across the simulation domain. The true flux-correction is that which can be computed from the full-scale DSMC solution x, Full , using the measured temperature field T Full and the measured heat flux field q x, Full in Eq. (13).

Micro resolution
The temperature-jump and associated thermal Knudsen layer are modelled within the micro element at each wall. To obtain a high level of accuracy for this Fourier flow problem, the sampling zone in each of these elements should extend to capture the entire Knudsen layer. This means that the required size of our near-wall elements increases with the level of rarefaction, i.e. Kn. Note that, for high values of Kn, the required size of the near-wall micro elements may result in the hybrid approach (over a number of iterations) becoming less e cient than a full-scale DSMC treatment. For low values of Kn (where the Knudsen layers remain close to the walls), micro elements could be required in the bulk to capture non-equilibrium behaviour caused by strong temperature gradients. Even when the bulk of the domain is near-equilibrium, bulk micro elements may still be needed to provide a correction for any error in the assumed thermal conductivity  r . The micro resolution (i.e the number of micro elements ⇧ and their size) required for a particular problem is therefore dependent on both the level of rarefaction, and the temperature conditions.
To explore the e↵ect of the micro element arrangement, we consider an example test case with Kn = 0.01. For a separation L = 1µm, we require a global mean free path = 0.01µm, corresponding to a gas number density n = 1.295⇥10 26 m 3 from Eq. (17). Setting an average gas temperature T av = 273 K, we assume a value of  r = 0.0164 W/mK [32]. The temperature di↵erence between the walls T is set to 50 K (i.e. T cold = 248 K and T hot = 298 K), resulting in a temperature gradient of 50⇥10 6 K/m across the domain.
For this initial test case, with this Knudsen number and temperature gradient, the bulk of the domain will be in equilibrium. However, the flux-correction field x in the bulk will be non-zero to correct for the assumption of a constant thermal conductivity. Over a temperature range of 200 K to 350 K, the thermal conductivity of argon varies approximately linearly with temperature, so x will also be approximately linear in the bulk. Direct linear interpolation between the sampling zones of the near-wall elements should therefore provide an adequate representation of the full flux-correction field, and bulk micro elements are not needed for this initial case. Two near-wall elements should be su cient, i.e. ⇧ = 2.
Ideally, the sampling zones of these near-wall micro elements should capture the thermal Knudsen layers fully. Also, the relaxation zones must be large enough to allow full relaxation of the particle distribution. The size of both of these zones is therefore key to the accuracy of the hybrid method, and it is important that we are able to make an estimate of the appropriate size depending on the flow conditions. Typically, thermal Knudsen layers extend several mean free paths from a wall surface. Similarly, the relaxation of a particle state typically occurs over several mean free paths. We therefore express the extents of both the sampling and relaxation zones as some number of local mean free paths l . Two sensitivity studies have been performed to find the appropriate extent for each zone. In the first study, we keep the length of the sampling zone L SZ at 5 l and consider the e↵ect of the relaxation zone length L RZ by incrementally increasing its value from 1 l to 7 l . Then, in the second study, we consider the impact of L SZ by increasing it from 1 l to 7 l , keeping L RZ constant at 5 l . Note that, for the first iteration, l is assumed equal to the global mean free path, = KnL. In subsequent iterations, l is updated for each micro element using Eq. (17), where n is the average number density measured in the sampling zone of the element in the previous iteration.
With elements in this size range, we need to see convergence of the hybrid method within 3 to 4 iterations for there to be any computational advantage over a full DSMC 14 simulation. For this test case, we set the tolerance parameter ⇣ tol to be 0.001. Figure 3(a) shows that, for all relaxation zone lengths, convergence of the method is reached within 4 iterations. However, convergence is not reached within 4 iterations when the sampling zone lengths are 1 l or 3 l , as shown in Fig. 3(b). We therefore require larger sampling zones of either 5 l or 7 l , both of which provide convergence in only 2 to 3 iterations. Note that, for all simulations, the level of convergence fluctuates slightly due to noise. The converged temperature profiles from the hybrid are shown in Fig. 4; also plotted is the initial temperature field T NSF , and the full DSMC temperature field T Full . The accuracy of the hybrid method is, however, more apparent from the mean percentage error at each iteration, presented in Fig. 5. The error at iteration l = 0 is the error in the initial hydrodynamic NSF solution. It is seen that the final mean error is less than 0.1% when relaxation and sampling zone extents are both 5 l or 7 l . This accuracy is due to the hybrid method's ability to capture the flux-correction field across the domain, as shown in Fig. 6.
Essentially, the element size depends on the desired balance between accuracy and computational savings. Although sampling and relaxation zone extents of 7 l provide slightly higher accuracy, we select an extent of 5 l for both zones. This provides a su cient level of accuracy while also enabling greater computational saving over a fullscale DSMC simulation. For the problems we study in this paper, each near-wall element therefore has a total extent of 10 l , while each bulk element has a total extent of 15 l .

Various rarefaction and temperature conditions
We now demonstrate the method's ability to capture temperature jump and the associated thermal Knudsen layer under various rarefaction conditions (study A) and various temperature conditions (study B). In each test case, we maintain a wall separation L = 1µm, and an average gas temperature T av = 273 K.
In study A, we consider the global Knudsen numbers: Kn = 0.01, 0.02, 0.03. This is achieved by varying the gas density for each case according to Eqs. (17) and (16). For all three cases, the temperature di↵erence T between the walls is set to 50 K. In study B, we simulate a range of temperature conditions by varying the temperature di↵erence between the walls: T = 50 K, 100 K, 150 K. For all three cases, the gas density is set to maintain Kn = 0.01. For every test case considered, both Kn and the temperature gradient are small enough that the bulk of the domain will be in equilibrium. As in section 3.2.1, two near-wall elements should therefore be su cient to capture the flux-correction field.
For all the cases, convergence of the hybrid solution occurs within 3 iterations to a tolerance parameter of ⇣ tol = 0.001. This is shown in Figs B, respectively. In Fig. 8(a), we present the converged temperature profiles from all three hybrid simulations of study A, along with the corresponding full DSMC temperature profiles. Also shown is the initial NSF temperature profile: as this is independent of Kn, it is the same for all test cases in study A. Similarly, the converged temperature profiles from all three hybrid simulations of study B are shown in Fig. 8(b). The corresponding full-scale DSMC and initial NSF temperature profiles are also shown, with the NSF solution di↵erent for each value of T . Again, a clearer measure of accuracy is to consider the mean percentage error of the hybrid solution (compared with the corresponding full DSMC solution) at each iteration. This is presented in Fig. 9 and, for each case, the error at iteration l = 0 represents the error in the initial NSF solution. For all three values of Kn in study A, the hybrid technique reduces✏ to approximately 0.07% as shown in Fig. 9(a). Figure 10(a) shows that, for all Kn, the final hybrid flux-correction representation is in fairly good agreement with the corresponding full DSMC flux-correction. However, as T is increased in study B, Fig. 9(b) shows that the overall accuracy of the hybrid method decreases. This is a reflection of the lower quality of the hybrid flux-correction representation, as seen in Figure 10(b). However it should be noted that, for all values of T , the hybrid solution provides a considerable improvement over the corresponding NSF solution.
For higher temperature gradients, the accuracy of the hybrid could be improved by slightly increasing the extent of the near-wall elements (i.e. larger L RZ and L SZ ), or by adding a element in the bulk. This would, however, reduce the computational savings obtained. As discussed in section 3.2.1, there must be a balance exercised between the accuracy required from, and the computational speed-up o↵ered by, the hybrid method.

Extreme temperature conditions
For the cases in sections 3.2.1 and 3.2.2, we used only two near-wall elements in our hybrid approach. However, even if the bulk of the domain is near-equilibrium, higher temperature conditions might result in a non-linear variation of conductivity with temperature. If so, the use of linear interpolation of the flux-correction field between opposite near-wall elements is not likely to produce the most accurate solution. Accuracy should be improved by the addition of micro elements within the bulk, i.e. increasing the micro resolution.
To verify the method's ability to deal with more extreme temperature conditions, we consider a case where the global Knudsen number is set to Kn = 0.01, but T av is increased to 500 K. The assumed value of  r is therefore increased to 0.027 W/mK [32]. The temperature di↵erence T between the walls is then set equal to 600 K (i.e. T cold = 200 K and T hot = 800 K). We consider two hybrid configurations: in the first, we have only two near-wall elements, i.e. ⇧ = 2; in the second, we add a bulk element in the centre of the system, i.e. ⇧ = 3.
Setting ⇣ tol = 0.01 for this test case, convergence occurs within 3 iterations for both configurations, as seen in Fig. 11. The converged temperature profiles from both hybrid simulations, along with the initial NSF and the full-scale DSMC temperature profiles, are presented in Fig. 12.
Once again, the accuracy of the hybrid is more clearly visible from the mean per-18 centage error, shown for each iteration in Fig. 13. As we would expect, the addition of a bulk micro element results in a more accurate temperature solution. In fact, in comparison with the initial NSF solution, this configuration provides an order of magnitude reduction in the mean error. This increased accuracy is the result of the higher quality flux-correction representation, as shown in Fig. 14.
The accuracy of the hybrid could be further improved by increasing the molecular resolution either by increasing the element lengths, or by adding more bulk elements, or both. It should be noted however that even with only two near-wall elements the hybrid solution still provides a considerable improvement over the NSF solution.

Computational savings
In this paper we have not exploited time-scale separation -the full-scale and the hybrid element simulations are run for the same number of DSMC time-steps. Computational savings therefore come only from length-scale separation. The dimensions of the Fourier flow test cases considered here have been restricted by the need for a full-scale DSMC simulation in order to validate the hybrid method. With the degree of lengthscale separation fairly small, the required micro elements occupy a significant portion of the domain. The computational savings provided by the hybrid on this particular heat transfer problem are therefore modest.
Note that M is the total number of time-steps over both the transient and steady-state periods and, for all cases in this paper, M Hybrid = M Full = 53⇥10 6 . The computational cost of the hybrid approach depends on both the micro resolution, and the number of iterations required. As the physical extent of each element is updated after each iteration based on the local mean free path, the average clock time per time step t c for a particular element will di↵er for each iteration. The total average clock time per time-step for the hybrid approach t c,Hybrid is therefore calculated as the sum of t c for all micro elements ⇧, over all iterations I, i.e.
For all test cases in section 3.2.2, our hybrid requires only two near-wall elements and convergence is reached inside 3 iterations. With ⇧ = 2 and I = 3, the computational speed-up S is presented in Table 1 for each case. These speed-ups are determined by the extents of the near-wall subdomains, (L SZ + L RZ ) left and (L SZ + L RZ ) right which are each 10 l . These extents during the final iteration (I = 3) of each case are also presented in Table 1. If Kn = 0.01, (L SZ + L RZ ) left and (L SZ + L RZ ) right are small enough that we see a modest computational speed-up (S > 1). However, when Kn is increased to 0.02 and 0.03, these extents become larger so that we see no speed-up at all (S < 1). 21    With I = 3, the computational speed-ups for this case for both ⇧ = 2 and ⇧ = 3 are shown in Table 2. For the configuration with only near-wall elements, (L SZ + L RZ ) left and (L SZ +L RZ ) right are again small enough for us to see a modest computational speedup. However, the increase in accuracy that comes with the addition of the bulk element increases the computational cost. In fact, with a bulk element length (L SZ + 2L RZ ) bulk corresponding to 15 l , there is no computational speed-up at all, and S <   In summary, while several cases in this paper demonstrate moderate computational savings by using the hybrid method, a number of cases are more computationally expensive than the equivalent full-scale DSMC simulation. It is important to note, however, that these Fourier flow cases have been chosen simply to test and validate the hybrid methodology itself. Future simulations of larger, more realistic, 2D and 3D problems will highlight any computational advantages of our multiscale approach. In transient flows, when temporal variations in the macroscopic property fields are much slower than the molecular time scales, computational savings can also be obtained by exploiting time scale separation, i.e. by controlling the time-steps in the macro and micro models as described by Lockerby et al. [33].

Conclusions
Based on a HMM-FWC approach, we have proposed a hybrid multiscale method that couples a continuum-fluid solver with the DSMC particle technique. The key advantage of this hybrid over existing continuum-DSMC hybrids is that DSMC micro elements of any size can be placed at any location in the flow, i.e. close to walls or in the bulk of the domain, independent of the continuum description. Therefore, unlike traditional HMM techniques, our method is able to simulate problems with any degree of scale separation. The micro resolution can be adjusted for each problem to obtain the desired balance 24 between accuracy and computational cost. Another useful feature of our approach is that the size of the micro elements adapts dynamically depending on the local mean free path of the gas. For heat-transfer problems, the coupling in our method is performed through the computed heat fluxes. As a simple validation test, the method was demonstrated on micro Fourier flow. For a range of rarefaction and temperature conditions, we have shown the method's ability to compensate for inaccurate boundary and constitutive information. The hybrid procedure was found to converge very quickly, inside 3 iterations for all test cases considered. Generally, good agreement with the equivalent full-scale DSMC simulations was observed, with the exact level of accuracy depending on the micro element arrangement.
Due to the small degree of scale separation in the Fourier flow cases considered here, the computational speed-ups observed (over the full DSMC simulations) were modest. This problem was, however, selected only to provide a simple example on which to test and validate the underlying methodology. Much greater speed-ups can be expected when simulating larger, more realistic engineering problems.
Before tackling 2D and 3D problems, the method needs to be developed to enable full coupling of mass and momentum, as well as heat transfer. Extension to 2D/3D then presents some additional challenges, including the application of mass, momentum, and heat transfer boundary conditions within 2D/3D relaxation zones, and 2D/3D interpolation of the flux-correction field between sampling zones. However, with this development, the method has the potential to tackle a variety of problems, such as the flow of complex gas mixtures, or the flow through micro heat exchangers. The modelling of thin bow shockwaves at the front of a planetary re-entry vehicle is another example which could benefit from the development of this hybrid as it encompasses both complex fluid behaviour (including chemical reactions of many gas species) and a high aspect ratio geometry. Our method also has the potential to be used in an inverse manner to obtain the transport properties of a gas.