Edinburgh Explorer Hybrid continuum-molecular modelling of multiscale internal gas flows

We develop and apply an eﬃcient multiscale method for simulating a large class of low- speed internal rareﬁed gas ﬂows. The method is an extension of the hybrid atomistic– continuum approach proposed by Borg et al. (2013) [28] for the simulation of micro/nano ﬂows of high-aspect ratio. The major new extensions are: (1) incorporation of ﬂuid compressibility; (2) implementation using the direct simulation Monte Carlo (DSMC) method for dilute rareﬁed gas ﬂows, and (3) application to a broader range of geometries, including periodic, non-periodic, pressure-driven, gravity-driven and shear-driven internal ﬂows. The multiscale method is applied to micro-scale gas ﬂows through a periodic converging–diverging channel (driven by an external acceleration) and a non-periodic channel with a bend (driven by a pressure difference), as well as the ﬂow between two eccentric cylinders (with the inner rotating relative to the outer). In all these cases there exists a wide variation of Knudsen number within the geometries, as well as substantial compressibility despite the Mach number being very low. For validation purposes, our multiscale simulation results are compared to those obtained from full-scale DSMC simulations: very close agreement is obtained in all cases for all ﬂow variables considered. Our multiscale simulation is an order of magnitude more computationally eﬃcient than the full-scale DSMC for the ﬁrst and second test cases, and two orders of magnitude more eﬃcient for the third case.

for internal-flow applications, which typically have a geometry of high-aspect ratio (i.e. are extremely long, relative to their cross-section). The high-aspect ratio creates a formidable multiscale problem: processes need to be resolved occurring over the smallest characteristic scale of the geometry (e.g. a channel's height), as well as over the largest characteristic scale of the geometry (e.g. the length of a long channel network), simultaneously.
Domain decomposition (DD) is a methodology that applies a molecular-based model where non-equilibrium/noncontinuum effects are located (statically or dynamically); for example, in the vicinity of boundaries, interfaces and moving shock waves. This localised molecular treatment is coupled to a continuum model applied over the rest of the domain (often an Euler solver for rarefied gas flows); the coupling is typically performed by matching hydrodynamic flux or state properties in an overlapping region (see Fig. 1(a), (b)). Despite their demonstrated success in a number of applications, DD methods have some disadvantages in certain classes of flow [26]. For example, a major disadvantage of using domain decomposition to solve the kind of highly-confined micro-/nano-channel problems we investigate in this paper is that the majority (if not all) of the flow can be considered 'near-wall', and thus highly non-equilibrial. As illustrated in Fig. 1(b) for a generic flow geometry, the molecular component (e.g. DSMC) would need to describe the entire length of the channel walls in order to capture the slip and other near-wall phenomena. As a consequence this severely restricts the maximum channel lengths (and the minimum channel heights) that can be usefully tackled by DD-type methods.
The heterogeneous multiscale method, known as HMM [27], is the main hybrid alternative to DD (see Fig. 1(c)). In HMM a continuum description is regarded as applicable in the whole domain and defined everywhere (although no assumptions regarding constitutive models or boundary conditions are made). Small 'micro elements' (individual molecular solvers) are dispersed on the nodes of the continuum computational grid, with the sole objective to provide the microscopic information that the macro model requires for closure (i.e. constitutive and boundary information). For a flow in a channel or pipe, for example, molecular micro elements would be placed on the bounding surfaces and in the fluid region, as illustrated in Fig. 1(c). Molecular simulations are performed using the local continuum strain-rate as a constraint, and from which only the necessary information close to the continuum governing equations is measured and supplied to the macro model. The shear stress measured in the micro elements located in the bulk of the fluid replaces any form of constitutive model or requirement for transport properties, while the wall micro elements produce accurate slip velocities and shear stresses, replacing any need for phenomenological slip boundary conditions [18].
The main drawback of HMM that makes it unsuitable for highly-confined micro channels (where the molecular scale is comparable to the transverse scale) is that micro elements, which have a minimum size of the order of the molecular scale (e.g. the mean free path in a gas flow), will be forced to overlap (illustrated in Fig. 1(c)). This is less efficient and less accurate than a full molecular treatment. It is less accurate for gases, for example, because micro elements in the bulk (whose state is dictated by local continuum variables) would not be able to capture non-equilibrium Knudsen-layer effects, which typically extend one or two mean free paths into the channel from the walls.
To tackle this problem the internal-flow multiscale method (IMM) has been developed [28] (see Fig. 1(d)). The IMM is a hybrid method tailored to flows in high-aspect-ratio channels. The exploitable multiscale feature of this type of flow configuration is the existence of length-scale separation between hydrodynamic processes occurring along the direction of the flow, and molecular processes occurring on scales transverse to the flow direction. As such, in this method every micro element covers the entire cross-section of the channel/tube, enabling the use of simple parallel-wall molecular simulations. In IMM, micro elements do not communicate directly with each other, instead they interact indirectly via constraints applied by the macroscopic conservation laws. No constitutive or slip models are required in the continuum formulation because this is provided (indirectly) by the micro elements. The IMM algorithm consists of a series of iterations, involving a two-way exchange of information between the macro domain and micro sub-domains, until convergence is obtained.
In this paper we extend the IMM to compressible fluids and dilute gas flows. A substantial modification to the methodology is needed in order to account for fluid compressibility, and to capture streamwise variations in density in order to predict mass flow rates and velocity profiles accurately. Even in low Reynolds number and low Mach number micro-scale flows, to which we restrict our attention, significant compressibility can occur [29], particularly in dilute gases. In fact, in the original IMM paper [28] a significant error due to the incompressibility assumption was identified, even though a dense Lennard-Jones fluid was considered. A simple explanation for this very low-Mach compressibility is that in highaspect-ratio geometries, significant variations in density are generated from a combination of high viscous losses and long channel lengths. Examples of such compressibility effects are provided below.
This paper is structured as follows. In Section 2, the modified IMM is introduced alongside details of a new numerical solution procedure. In Section 3 results obtained from the multiscale method for three examples cases are presented. These results are compared to equivalent results obtained using a full DSMC simulation, serving as a validation of the new IMM method. In Section 4, recommendations are made and conclusions drawn.

Multiscale and numerical methodology
In the internal-flow multiscale (IMM) approach [28], individual micro sub-domains, covering the full cross-section of the channel or tube, are distributed in the streamwise direction with a spacing sufficient to resolve any streamwise variation in the channel/tube geometry. The key assumption in the IMM -and in fact what it is designed to exploit -is that there is a degree of scale separation between hydrodynamic variation along streamlines (i.e. in the s-direction, see Fig. 2(a)) and microscopic processes transverse to the flow direction (i.e. in the y and z directions). For this scale separation to exist, the section geometry must vary slowly in the streamwise direction. A dimensionless number indicating the degree of such scale separation can be developed: where L y and L z are length scales in the transverse directions that characterise the cross-section geometry, and R y and R z are the radii of curvature of the centre streamline. Provided S 1, it can safely be assumed that, in small streamwise sections of the tube/channel, the walls are approximately parallel to the centre streamline (and consequently all other local streamlines). This is referred to in [28] as a local parallel-flow assumption, and it enables the representation of small sections of the channel/tube geometry by micro sub-domains with exactly parallel walls, as illustrated in Fig. 2(b). These sub-domains are otherwise identical to the full geometry in terms of the cross-section dimensions of the channel at that location. In terms of boundary conditions for the sub-domains, again, these are the same as for the full geometry, except that periodicity is imposed in the streamwise direction (the consequences of this are discussed later).
The objective of the method is for the flow conditions in these otherwise independent micro sub-domains to represent the conditions that would be occurring at the same point in the whole flow geometry, i.e. in the macro domain. To do this a coupling needs to occur between the micro sub-domains and an appropriate macroscopic representation of the flowfield.

The hydrodynamic conservation equations
Attention is restricted here to low-speed (low Reynolds number, low Mach number), steady-state flows, as these are of primary importance to the majority of applications outlined in Section 1. In such conditions the equation for conservation where τ sy and τ sz are perpendicular components of shear stress, p is the hydrostatic pressure, ρ is the mass density, and a ext is an external acceleration (ρa ext is a body force per unit volume in the flow direction). For individual decoupled micro sub-domains (i.e. which do not involve, as yet, any macro connection) a different momentum balance is formed: The difference between Eqs. (2) and (3) is the existence of a local pressure-gradient term, which in the macro domain is generated by a combination of the applied pressure difference (p B -p A , in Fig. 2(a)), and the requirement for mass continuity through the entire geometry -this term does not feature in Eq. (3) because of the streamwise periodicity of the micro sub-domains. To facilitate a coupling mechanism between the micro and macro domains, the missing pressure gradient can be emulated within the micro sub-domain by a fictitious body force, Φ, which is referred to as a pressure-gradient correction: Provided the difference in pressure over the streamwise length of a micro domain is relatively small, the pressure-gradient correction force, Φ, is hydrodynamically equivalent to the effects of the missing pressure gradient, i.e.
Note, because of the locally-parallel-flow assumption, the pressure gradient (and consequently Φ) is constant across the cross-section.
The aim of the incompressible IMM presented in [28] is to find (by iteration) values of Φ for each micro sub-domain that collectively ensure: (a) the boundary conditions on pressure are correctly imposed (p B -p A ), and (b) mass conservation is satisfied, i.e. the same mass flow rate is generated in each of the micro simulations. Where the method in this present paper differs is that it must also obtain, as part of the iterative procedure, values of gas density in each micro sub-domain (ρ i ). Gas density directly affects viscosity, rarefaction effects, and mass flow rate, and so obtaining these values accurately is essential for the prediction of compressible flows, even at very low Mach numbers.

Compressible IMM
The aim of the IMM is to find individual sub-domain states that collectively satisfy macroscopic mass conservation, i.e. all sub-domains should have the same mass flow rate: whereṁ i is the mass flow rate in the ith of Π sub-domains andṁ is a target value (an output of the multiscale simulation).
The body forces (Φ i ) and density states (ρ i ) in each sub-domain are found via an iterative algorithm, with the target of satisfying Eq. (6). In order for this iterative process to be efficient, we need to make some prediction of the general response of the mass flow rate to changes in Φ and ρ. To this end we assume that the difference between the total mass flow rate and the mass flow rate due directly to wall motion is proportional to the net momentum flux introduced by external acceleration (a ext , e.g. gravity) and the pressure-gradient correction, i.e.
where Q wall is the volumetric flow rate that is due to wall motion alone (i.e. the volumetric flow rate in the absence of any pressure gradient or body force). The constants of proportionality k i can be estimated at the beginning of the algorithm by performing initial simulations for which ρ i and Φ i are known (or which can be chosen arbitrarily) and from which mass flow rate is measured; k i are then calculated from Eq. (7).
Note, the prediction in Eq. (7) is exact for near-equilibrium flows at the macro scale, but only approximate at the micro scale for rarefied gas flows. However, the accuracy of this prediction only determines the convergence characteristics of the method and does not affect the accuracy of the final solution (provided one can be found).
Eq. (7) is then used to estimate the change in mass flow rate between successive iterations: where l is the iteration index and angular brackets denote a value measured from the sub-domains. To move the mass flow rates in each sub-domain towards the same single value, we substituteṁ l+1 i forṁ, to obtain, after rearrangement: As Φ and ρ converge (ρ l+1 → ρ l and Φ l+1 → Φ l ), Eq. (9) ensures mass conservation (6) is satisfied, i.e. ṁ i l →ṁ.
For simplicity, we focus on ideal gases, for which an equation of state exists, allowing the re-expression of Eq. (5) as: where R is the specific gas constant and T is the temperature (here assumed to be constant in the streamwise direction). At each iteration, Eqs. (9) and (10) are solved simultaneously. We therefore need a discrete representation of the derivative of density in Eq. (10) (which is not required in the incompressible IMM formulation of [28]). Sub-domain simulations are very expensive, and the discrete representation of Eq. (10) should be as spatially accurate as is practical to implement. To this end we adopt a spectral collocation method. Spectral methods have extremely good error characteristics, exhibiting so-called spectral convergence (the fastest possible). Spectral convergence means that increasing resolution (i.e. increasing then number of sub-domains per streamwise length, Π/L s ) decreases the error exponentially (i.e. ∝ (L s /Π ) Π ) as opposed to algebraically as for finite-difference methods (i.e. ∝ (L s /Π ) r , where r depends on the order of the scheme). In fact, even though it is more accurate, the implementation of the spectral method is not significantly more complex than a standard finite-difference scheme for the IMM.

Spectral collocation scheme
The most suitable polynomial type for the collocation scheme for Eq. (10) depends on the configuration of the particular problem. 1 In this paper we consider periodic and non-periodic problems, but here we present the scheme for periodic configurations (used in Section 3.1 and Section 3.3), for which the most natural choice is Fourier polynomials. An equivalent scheme for non-periodic configurations (used in Section 3.2) is presented in Appendix A.
The density variation in a periodic domain is represented by the following Fourier polynomial: whereρ 1,..,Π are the Fourier coefficients (note, Π must be odd) and L s is the total streamwise length of the geometry. From Eq. (10) we can directly obtain a Fourier polynomial description of Φ: Eqs. (11) and (12) are evaluated at the sub-domain locations (collocation points): 1 A mapped spatial variable can also be used to redistribute resolution of the chosen polynomial.
to give variables in physical space at the new iteration: where iteration indices for the Fourier coefficients are omitted for clarity. These expressions are substituted into Eq. (9) to obtain:ṁ − ρ 1 + Eq. (16) represents Π equations for Π + 1 unknowns (ρ 1,..,Π andṁ). The additional required equation comes from knowledge of the total mass within the system, M (for periodic systems, this is fixed from the outset): where A i are the local sectional areas of the geometry at the Π sub-domain positions. The integration here is using the trapezoidal rule (i.e. the mean of the variables at the collocation points multiplied by L s ), which is very accurate for periodic systems. The set of linear equations (16) and (17) can be solved at each iteration using LU decomposition, or similar. After each iteration, the variables at collocation points (in physical rather than Fourier space) can be obtained from Eqs. (14) and (15).

Compressible IMM algorithm
The solution to the multiscale problem is an iterative one. Once converged, it enforces macroscopic mass conservation (momentum is automatically conserved at every iteration). The algorithm stops once a convergence criterion (e.g. on standard deviation of the mass flow rates in the micro simulations) is satisfied. The iterative procedure has the following steps: (0) Initial setup. Set the density states in all Π sub-domains to a constant value that satisfies the total mass constraint (17) for periodic cases, or satisfies boundary conditions on pressure/density for non-periodic cases. Since the first iteration is primarily used to obtain an estimate of k i for each micro sub-domain, the initial values of Φ i are arbitrary, provided Φ i = ρ i a ext (see Eq. (7)).
(1) Initialise all micro (DSMC) simulations with the density and pressure-gradient correction values for the current iteration: (2) Simulate all Π micro simulations in two stages: a transient period and an averaging period. The transient period must be long enough for the flow in the micro sub-domains to reach steady state, and the averaging period must be long enough so that a satisfactorily noise-free mass flow rate measurement can be obtained: ṁ i . As a convergence criterion, we take: where ζ tol is a predetermined tolerance. For the simulations in this paper we choose ζ tol = 0.025.

Micro simulation model
For simulating the rarefied gas flow in each of the micro sub-domains ( Fig. 2(b)) we use the direct simulation Monte Carlo (DSMC) method [2]. The DSMC code used in the present work has been implemented in OpenFOAM, an open-source C++ tool box for computational fluid dynamics which is free to download from [30]. This solver, called dsmcFoam, has been validated for a variety of benchmark cases including hypersonic and micro-channel flows [31][32][33].

Results and discussion
In this section we apply the compressible IMM described in Section 2 to three rarefied, low-speed, high-aspectratio, micro-channel flow configurations: streamwise-periodic converging-diverging channel flow (Section 3.1); non-periodic pressure-driven channel flow around a corner (Section 3.2); and eccentric cylindrical Couette flow (Section 3.3). For validation purposes, in each case the multiscale results are compared to results obtained from a full DSMC simulation of the complete geometry using the same DSMC code. From this perspective, the full-scale DSMC simulations are considered the 'accurate' results, but are, as will be shown later, far more computationally expensive than our compressible IMM. Note, DSMC has itself been extensively validated against experiment for a range of applications.
The DSMC simulations presented in this paper are two-dimensional for both the sub-domain and the full-scale simulations.

Test case: converging-diverging channel
Our first test case is a converging-diverging channel flow, with periodicity in the streamwise direction, driven by a uniform and constant acceleration; the configuration is shown in Fig. 3. The length of the geometry is L s = 24 μm, and the thickest and narrowest points of the channel are 1 μm and 0.4 μm, respectively. The temperature of the walls is kept constant at T = 300 K, and the average density in the domain is ρ = 0.442 kg/m 3 . The simulated gas is argon.
For reference, a full DSMC simulation of the flow in the configuration is performed with ∼6 × 10 5 DSMC particles ( P full ), chosen to give ∼100 particles per cell on average, with each cell having dimensions considerably smaller than the mean free path, λ. The flow has a maximum Reynolds number Re = 0.44 (based on local centreline velocity and local channel thickness, h), a maximum Mach number Ma = 0.13, and Knudsen numbers ranging from Kn = 0.4 to 1.1 (Kn = λ/h). The applied acceleration is a ext = 10 10 m/s 2 ; this large value is used so that a manageable signal-to-noise ratio from the expensive full DSMC simulation can be obtained. The measured mass flow rate isṁ = 60.05 ± 0.2 × 10 −15 kg/s. , with slip velocity increasing concomitantly. What is also obvious from Fig. 5 is that the density within the domain varies substantially (the highest density is almost 2.5× greater than the lowest value), demonstrating the need for our new compressible extension to the IMM. Note, this compressibility occurs despite the low Mach number conditions (maximum Ma = 0.13) due to a combination of high viscous forces and long channel lengths, as discussed in [29].  To accurately estimate the computational saving that the IMM has afforded us, we can compare the total particletimesteps (G) computed in each method. For the full-scale reference DSMC this is straightforward: G full = P full (T su,full + T av,full ), (19) where P full is the total number of DSMC particles simulated, T su is the number of timesteps required for the flow to become developed (a 'start-up' time), and T av is the number of timesteps used in the averaging phase of the simulation. For the compressible IMM: G IMM = P IMM (T su,IMM + T av,IMM )I , (20) where I is the number of iterations of the multiscale algorithm performed (here I = 3). The ratio of Eqs. (19) and (20) provides the total computational savings factor: For consistency, we set our averaging times to be equal: T av,full = T av,IMM = T av (the timestep is also kept the same). In the converging-diverging test case of this section, the development times for both the multiscale simulations and the full-scale simulation are much smaller than the chosen averaging time, and thus: C ≈ P full /(P IMM I). From this we can calculate that the overall computational cost of the compressible IMM is approximately 6× less than the full-scale DSMC.
To investigate the influence of the number of sub-domains on the accuracy of the compressible IMM, the same converging-diverging channel simulation has been performed for Π = 5, 7 and 9 sub-domains. The calculated mass flow rates are:ṁ = 60.54 × 10 −15 kg/s (Π = 5);ṁ = 60.48 × 10 −15 kg/s (Π = 7); andṁ = 60.3 × 10 −15 kg/s (Π = 9) -all these are within 1% of the full DSMC value. Streamwise variation of density through the geometry is plotted in Fig. 6 for each of the multiscale simulations. It is clear that the accuracy of the method is still very good when as few as Π = 5 sub-domains are employed, for which the computation is over 10× cheaper than the full DSMC.

Test case: pressure-driven channel flow around a corner
The second test case is also a converging-diverging channel flow, but unlike in Section 3.1 the flow is driven by a pressure difference and the geometry turns smoothly through 90 degrees. This non-periodic configuration is shown in Fig. 7. The length of the geometry is L s = 24 μm, and the thickest and narrowest points of the channel are 1 μm and 0.4 μm, respectively. The inlet pressure is p 1 = 46.3 kPa (Kn = 0.12) and the outlet pressure p Π = 8.9 kPa (Kn = 0.6). The temperature of the walls, and the inlet/outlet, are kept constant at T = 300 K. The simulated gas is argon.
A full DSMC reference simulation of the configuration is performed with ∼18 × 10 5 DSMC particles ( P full ), chosen to give ∼100 particles per cell on average, with each cell having dimensions considerably smaller than the mean free path, λ. The flow in the geometry has a maximum Mach number Ma = 0.056, and Knudsen numbers ranging from Kn = 0.12 to 0.72. The measured mass flow rate isṁ = 18.9 ± 0.2 × 10 −15 kg/s.

Test case: eccentric-cylinder Couette flow
Now we consider the flow generated in a gas contained between two eccentric cylinders (the outer stationary, the inner rotating at a constant rate), as in Fig. 10; this is representative of a plain bearing configuration. The inner and outer cylinders have diameters of 40 μm and 38 μm, respectively, and the inner cylinder's axis is offset by 0.5 μm, making the maximum and minimum gap between the two cylinders 1.5 μm and 0.5 μm, respectively. The constant angular velocity of the inner cylinder creates a tangential wall velocity of 19 m/s. The contained gas is argon, and the average density in the domain is ρ = 0.354 kg/m 3 . Both walls are held at a constant temperature of T = 300 K.
A full reference DSMC simulation of this configuration is performed using 2.88 × 10 6 DSMC particles. The flow has a maximum Reynolds number Re = 0.16 (based, again, on local centreline velocity and local channel thickness), a maximum Mach number of Ma = 0.042, and Knudsen numbers ranging from Kn = 0.35 to 0.6. The measured mass flow rate iṡ m = 27.27 ± 0.54 × 10 −15 kg/s. A compressible IMM simulation of the geometry is also performed with Π = 11 sub-domains, whose circumferential positions (s i ) are given by Eq. (13) where L s = 122.5 μm (the circumference of the centre streamline); sub-domain positions and scale are illustrated in Fig. 10. The minimum value of the scale separation number in the geometry, as defined in Eq. (1),  is S = 13. The total number of DSMC particles used in the multiscale simulation (i.e. summed over all 11 sub-domains) is 44 × 10 3 , which is roughly 65× fewer than the full DSMC simulation. Here, again, the number of DSMC particles has been chosen to give comparable resolution.
Based on the convergence criterion given in Eq. (18), the IMM reaches a solution after just two iterations (i.e. after two separate simulations of the Π = 11 sub-domains). The predicted mass flow rate from the compressible IMM iṡ m = 27.16 × 10 15 , which is within 0.5% of the prediction of the full-scale DSMC. Fig. 11 shows velocity profiles at a number of circumferential positions and Fig. 12 shows the circumferential variation of density. In the narrower portion of the configuration there is a negative curvature (at around i = 2) due to a favourable pressure gradient; this curvature becomes positive in the wider regions as the pressure gradient becomes adverse. As in the previous cases, very close agreement is found between the compressible IMM and the full-scale DSMC calculation. The noise in the full DSMC velocity profiles appears larger than that in the multiscale simulations because a smaller bin size is used to create the velocity profile; however, the DSMC cell sizes are comparable.  Eq. (21) is used to evaluate the computational saving of the IMM. In this test case the development time (i.e. the time to reach steady state) for the full-scale geometry is far longer than the development time of the slowest micro element (T su,full ∼ 200T su,IMM ), and greater than the chosen averaging time (T su,full ∼ 8T av ). This relatively long development time of the full-scale configuration is due to the high scale separation and the weak gradients in the streamwise direction as compared to gradients across the channel. It is the streamwise gradients in mass flow rate that promote temporal changes in density, but because these gradients are weak this occurs very slowly relative to viscous-driven development across the channel. Not having to simulate this lengthy development period is a major computational advantage of the compressible IMM for steady-state flows, and results from having coupled the micro sub-domains to a steady-state form of the continuum conservation equations. We find the total computational cost of the IMM, as calculated from Eq. (21), is approximately 300× less than the full-scale DSMC.

Conclusions
In this paper we have presented a major extension to the internal-flow multiscale method (IMM) [28] to take into account fluid compressibility (which can be very significant in micro-channel geometries at low Mach numbers). We have applied the method to dilute gas flows (using DSMC as a micro model) for periodic and non-periodic internal-flow geometries.
In all three test cases we considered, the compressible IMM gives flow predictions that are in very close agreement with those of a full-scale DSMC of the same geometry. Most importantly, the compressible IMM is approximately 6×, 50× and 300× more efficient than full-scale DSMC simulations of these cases.
The major savings afforded by our multiscale treatment of the final case (the eccentric cylindrical Couette flow of Section 3.3) is in part due to the extremely long development time of the flow in the full-scale simulation. The IMM avoids this lengthy development period because the macro conservation equations used in the coupling algorithm are all steady-state. The IMM is currently not appropriate for transient problems. However, in principle a time-accurate macro/micro coupling could be developed, with the time-stepping in each model controlled to exploit time-scale separation where it exists, as described in [25]. This is for future work.