Application of a hybrid kinetic-continuum solver to the near wall modelling

A hybrid method dynamically coupling the direct numerical solution of the S-model kinetic equation and Navier-Stokes equations is applied to a numerical simulation of the flow through the channel of a finite length due to arbitrarily pressure ratios and for a wide range of Knudsen number. The decomposition of the physical domain into kinetic and hydrodynamic sub-domains is updated at each time step. The solution is advanced in time simultaneously in both kinetic and hydrodynamic domains: the coupling is achieved by matching half fluxes at the interface of the kinetic and Navier-Stokes domains, thus taking care of the conservation of momentum, energy and mass through the interface. Solver efficiency is increased via MPI (Message Passing Interface) parallelization. Accuracy and reliability of the method, for different decomposition criteria, are assessed via comparison with a pure kinetic solution.


Introduction
At microscale, the analysis of gaseous flow may require to take into account rarefaction effects. The first area where such effect becomes significant is the near-wall region. This is due to the well known fact that the flow near a solid surface can be divided into a thin boundary layer, of the order of a few mean free path (so thin to be negligible for macro configuration), which is a rarefied regime area and the internal core flow, which is a continuum one. The rarefied domain is naturally described by kinetic equation for the velocity distribution function, which involves a considerable effort in terms of CPU time and memory requirements, due to the discretization in both physical and velocity spaces. On the other hand, the continuum domain is well described by the fluid Navier-Stokes (NS) equations in terms of mere average flow velocity, gas density and temperature. These equations are more efficient, but less accurate in critical rarefied areas. The inaccuracy of NS equations in boundary layer can be partially overcome by introducing slip boundary and temperature jump conditions on the solid surface. However, as indicated in the literature, the slip conditions are valid only for the local Knudsen number Kn ≤ 0.1, and any attempt to increase their range, resorting to higher order slip boundary conditions, is not trivial and highly geometrical dependent. Therefore, the development of hybrid solvers applying kinetic model in the region with high rarefaction, while keeping continuum model in the rest of the flow allow us to simulate flows with an accuracy close to the full kinetic solution and computational expenses close to the continuum one. Thus, this became an important area of research over the last decade. Potential applications of such solvers range from gas flows in micro systems to aerospace applications, such as high altitude flights. Major challenges in the development of hybrid code are the identification of kinetic and continuum domains and the choice of the coupling technique.
The most popular methods, in the open literature, are based on domain decomposition in physical space: the computational domain is decomposed into kinetic and continuum sub-domains using appropriate criteria [1][2][3][4][5][6][7]. Typically, particle methods such as DSMC or Molecular Dynamics are used in regions with strong deviations from equilibrium, and a continuum fluid (Euler or NS, depending on problem features) solver is used elsewhere [4]. Nonetheless, the Direct Numerical Solution (DNS) of kinetic equations is a viable alternative to DSMC and may even be preferable to it, for coupling purposes, since both kinetic and continuum models use similar numerical techniques. Recent efforts to combine the DNS of kinetic equation with a NS solver used a priori domain decomposition in [1,2] and later evolved in dynamically updated decomposition in [3]. In [3] a wall boundary independent problem (the gas flow through the slit) was considered and the decomposition of the physical domain into kinetic and continuum sub-domains was done by computing gradient-length Knudsen number Kn GL , based on the local Knudsen number and macroparameters gradients proposed in [4,6].
In the present paper the same hybrid solver is applied to the simulation of the rarefied gas flow through a channel of finite length caused by an arbitrary pressure ratio, where the flow is dominated by the wall boundary effect. The hybrid results are compared with pure kinetic and NS solutions, for different pressure ratios and a wide range of Knudsen number. The results are discussed in terms of both accuracy and computational efficiency.

Statement of the problem
The two-dimensional pressure-driven monoatomic gas flow through a planar microchannel of width H and length l = 10H connecting two reservoirs of L x × L y = 30H × 15H is considered. The gas in reservoirs far from the channel is in equilibrium at constant pressures p 0 and p e , (p 0 > p e ), and temperature T 0 . The temperature of the wall T w is equal to the temperature in reservoirs T 0 . Since the flow is symmetric about y = 0 only a half of the domain is simulated. The interface I c between kinetic and NS sub-domains is dynamically updated during the computation. The gas flow is determined by the pressure ratio p 0 /p e and the rarefaction parameter δ: v 0 = (2RT 0 ) 0.5 is the thermal speed, R is the specific gas constant, μ 0 is the dynamic viscosity at T 0 . The dimensionless flow rate W is introduced as a global characteristic of the flow: m fm is the analytically deduced in the limit of free molecular regime mass flow rate. The mass flow rate through the channel is constant and computed at channel mid-section as: where ρ(x, y) is the gas density, u(x, y) is the gas bulk velocity.
The friction factor f r is defined as an average value over a channel section of ∆L ≈ 9 (i.e., skipping small regions ∆l ≈ 0.5 both at the inlet and exit of channel to avoid entrance and exit effects): where D h is the hydraulic diameter, D h = 2H, overbar means average values over cross section, ∆ and av are the difference and the average between inlet and outlet values. The Poiseuille number Po is then written in terms of the local Re as: Due to the conservation of mass flow rate, Re will change only slightly along the channel, because of the dynamic viscosity µ variations with temperature.

Coupling algorithm and numerical method
The coupling strategy between kinetic and continuum solvers is applied to the S-model of Boltzmann equation [8] and NS equations. The S-model can be written in dimensionless form as: where f = f(t, x, ξ) is the velocity distribution function, i.e. the probability of finding a molecule with velocity ξ = (ξ x , ξ y , ξ z ) in the position x = (x, y) at the time moment t. S is the standard Shakhov function [8], ρ is the gas density, T is the temperature, c = ξ -V is the relative speed of a molecule against a background gas with bulk velocity V = (u, v). The problem was recast in terms of dimensionless variables using the inlet reservoir equilibrium values and the height of the channel H.
Assuming the hard-sphere molecular model the dimensionless viscosity coefficient is μ = T 0.5 . The macroscopic density, momentum and internal energy e per unit mass, respectively, are computed as: Since the flow considered is two-dimensional it is possible to eliminate the third velocity component [1][2][3]. The S-model equation is discretized in both velocity and physical spaces and solved using the explicit-implicit numerical scheme [1][2][3]. In particular, the transport term in equation (6) is treated explicitly and approximated by a standard finite volume scheme. Numerical fluxes are determined by the TVD scheme with minmod limiter.
According to the Chapman-Enskog (CE) theory [9] the NS equations are the first order approximation of the Boltzmann equation. Thus, introducing dimensionless CE function f CE as: where M(c) is the standard Maxwellian, f 1 is the correction term, τ ij is the shear stress tensor, δ ij is the Kronecker delta, κ = 15/4µ is the thermal conductivity for the hard sphere molecular model. Substituting f CE for f in equation (6), multiplying by the collision vector φ(ξ) = (1, ξ, ξ 2 /2) and integrating over the whole velocity domain [-∞;∞] we recover the conservative form of NS equations: where U is the vector of macroscopic values, F(U) is the flux. The Navier-Stokes solver is based on a hybrid finite difference-finite volume method and has formal second order accuracy in space and time [1][2][3]. Fluxes are defined via neighbouring values averaging, and an artificial dissipation term is added to prevent checker boarding and numerical instabilities. Artificial dissipation terms are given by a blend of second and fourth order differences, scaled by the maximum eigenvalue of jacobian matrix of the convective terms in vector F, as suggested in [10]. Second order terms are switched on near discontinuities. Viscous flux vectors are evaluated with second order finite differences centred at mesh midpoint i+1/2. The solution is advanced in time via Crank Nicolson integration scheme. The use of the spatially factored ADI scheme proposed by Beam and Warming leads, at each time step, to the resolution of two series of cheap block tridiagonal algebraic systems, rather than the original pentadiagonal block system arising from flux discretization.

Breakdown parameter
The choice of the breakdown criterion defining the size of Ω K region is important, since a wrong domain decomposition could even lead to a non-positive velocity distribution function. One of the options proposed in the literature is the gradient-length Knudsen number Kn GL (x, y) [4,6]: where Φ(x, y) represents the parameter of interest: the local density ρ, the magnitude |V| of local bulk velocity and the local temperature T. The kinetic solution is performed where Kn GL (x, y) > ε = 0.05. It should be noticed that here Kn GL , based on the gradient of the velocity magnitude |V|, is different from what proposed in [4,6]. The new definition takes into account share stresses which are essential for the considered problem. Moreover, the velocity magnitude gradient in Kn GL |V| is normalized via the maximum between |V| and a minimum value v min , here chosen as 5% of the theoretical isentropic exit velocity, in order to avoid singularities where the velocity approaches to zero. Another breakdown parameter is based on the CE distribution function [5,6]. When f 1 (equation (8)) is small Ψ(c) slightly deviates from 1 and NS equations are valid; when Ψ(c) is sufficiently far from unity then NS equations may fail, and a kinetic approach is required. The direct evaluation of f 1 is difficult and numerically expensive because it is a function of not only flow field gradients but also the random velocity c and so either an average or maximum value of f 1 must be evaluated over the full distribution function. Thus, in [5,6] an approximate breakdown parameter based on the normalized shear stress τ * ij and heat flux q * i was proposed: The kinetic solution is activated when B qτ > 0.1. Tiwari [7] has proposed a similar criterion based on the norm of additional term in CE distribution function ||f 1 || which can be written in dimensionless form as: For the continuum approximation to hold this parameter should be much less than unity.

Coupling procedure
To solve a kinetic equation in a physical cell x i+1 we need the distribution function of incoming particles from the neighbouring cell x i (see Fig. 1). On the coupling interface I c this information has to be provided by the continuum solution: therefore, we impose f CE distribution function: where η(x) is the inward normal vector to the boundary of Ω K . Macroscopic values ρ, V, T appearing in f CE (x i ) are computed at the grid point x i ∈ Ω NS and the evaluation of parameters gradients involves also values in the neighbor to x i points: x i-1 and x i+1 . The interface I c position is updated at each time step and depends on the chosen breakdown criterion and threshold value ε; thus, some nodes considered as continuum at previous time step may become kinetic ones at current time step. In this case, the kinetic distribution function f in "new" kinetic nodes will be initialized as f CE (x) computed using current macroscopic values and their gradients. To solve NS equations at the coupling interface I c the following boundary condition is imposed:

F i (U)·η(x) and F o (U)·η(x)
are half fluxes associated with incoming (from Ω K to Ω NS ) and outcoming (from Ω NS to Ω K ) molecules: where f(t, x, ξ) ≡ f K is the current solution of kinetic equation for molecules exiting from Ω K . If a boundary node at the inlet/outlet falls in continuum domain the conditions in left and right reservoirs specify the inlet total temperature and pressure at the inlet, or the static pressure at the outlet. If an inlet/outlet boundary node falls in the kinetic domain a standard Maxwell distribution function, based on inlet/outlet total temperature and pressure, is assumed for incoming particles.
At solid wall in kinetic domain Ω K Maxwell diffuse reflecting boundary condition with the full accommodation is applied. In Ω NS domain the Maxwell first order slip boundary condition and a Dirichlet temperature boundary condition with Smoluchowski temperature jump are imposed.
The hybrid code is parallelized in order to improve its efficiency using MPI message passing protocol. The solution of kinetic equation (5) is local in velocity node, and therefore completely parallelizable in the velocity space. For the solution of NS equation a set of block tridiagonal algebraic systems is solved sequentially: the solution of each system may thus be considered as independent and can easily be parallelized. The software code was written in C++ (for the kinetic and coupling coding) and Fortran (for the NS related procedures) with the use of MPI. Computations were carried out on Multi Core system consisting of 2 quad core processors, Intel(R) Xeon(R) E5520 CPU, 2.27 (2.93) GHz, 8 MB Cache, for a maximum of 8 concurrent cores.

Parameters of modeling
The size of uniform two-dimensional velocity grid should be selected large enough to capture all of the important features of the problem: here, the velocity space boundary satisfies the condition v max ≥ max(|u|, |v|) + 3. direction (with minimum grid spacing 0.008H), is used. Grid independence test has been done using a coarse grid with 240 × 40 nodes. The mass flow difference between finer and coarse meshes is less than 1%. The time step is unique for both solvers and it should satisfy the stricter stability (or accuracy) constraint Δt = min(Δt K , Δt NS ). The solution is considered to be converged when the criterion ||U n+1 -U n || L2 < Δ is fulfilled with L 2 norm and Δ = 10 -7 .

Breakdown criteria comparison
In [3], for the flow through a slit, the breakdown criterion Kn GL > 0.05 guaranteed a difference between hybrid and kinetic solutions within 1%. Here, we evaluate the effectiveness of modified Kn GL (using gradient of velocity magnitude |V|) and the two more complex criteria B qτ and ||f 1 || applied to the boundary value problem, via comparison with pure kinetic Boltzmann computations [11]. In table 1 mass flow rates W h and global Po h using different criteria at different pressure ratios and rarefaction levels are presented. At p 0 /p e = 2 and δ = 10 the value of ε, for each criterion, was chosen such as to provide hybrid solution close to kinetic one while using the same number of kinetic points N K . This 'optimal' value is retained for further computations with the same criterion, i.e.: ε = 0.1 for Kn GL , ε = 0.015 for B qτ and ε = 0.05 for ||f 1 ||. For both pressure ratios the maximum deviation of mass flow rate W h from kinetic value is around 1% and 2% for Po number. At the same time the difference between hybrid solutions computed with different criteria is less than 1%. It should be noticed that for low rarefaction level (δ = 50) the use of criterion ||f 1 || and B qτ generates larger number of the kinetic points than the use of Kn GL with little or no gain in accuracy. As can be seen in figure 2 (top), showing contour line of local Knudsen number, all criteria create almost identical kinetic regions at p 0 /p e = 2 and δ = 10. On the other hand, for high rarefaction δ = 50 at the same pressure ratio (figure 2, bottom) criteria ||f 1 || and B qτ create significantly larger kinetic region: 6800 points for ||f 1 ||, 5300 for B qτ , 3030 for Kn GL . Thus, for further hybrid computations breakdown parameter Kn GL with ε = 0.1 will be used.

Hybrid results discussion
Assuming the kinetic solution as the reference one, we may define the relative error induced by the use of hybrid and pure NS solvers as: for mass flow rate W and Poiseuille number Po, respectively. These relative differences ΔW NS , ΔPo NS , ΔW h and ΔPo h at pressure ratios p 0 /p e of 10, 3, 2 and 1.1 and rarefaction parameter δ ranging from 10 to 100 are shown in figure 3. In the slip regime (δ ≥ 50) mass flow rates and Po obtained by hybrid, kinetic and NS solvers are close to each other: maximum difference between results is less than 2%. Starting from δ = 20 the difference becomes noticeable, around 4-6%. When rarefaction increases the difference growth up to 10-14% at δ = 10. When gradients of macroparameters become higher the local Knudsen number increases: thus, at the same rarefaction parameter δ, we can have higher local Knudsen for higher pressure ratio; e.g., for rarefaction δ = 10 and for pressure ratio p 0 /p e = 2 the difference between NS and kinetic results is around 10%; for p 0 /p e = 3 the difference is around 12% and 14% for p 0 /p e = 10. On the other hand, the hybrid code successfully reproduces kinetic values for both W and Po under any condition (difference is within 1% for W and 2% for Po).
Profiles of density and temperature, normalized by inlet values, and velocity along the symmetry axis y = 0 are shown in figure 4 for selected cases, chosen in order to demonstrate noticeable difference between NS and S-model solutions.
Density (or pressure) variations are qualitatively similar in all cases. Before and after the channel density tends to upstream and downstream conditions, while decreasing along the channel. At high pressure ratio the density profile is nonlinear. The axial velocity, close to zero in the large reservoirs, grows along the channel. Accordingly, the temperature decreases through the channel, while in the low velocity reservoirs it approaches to the reference temperature. The maximum variations of macroparameters appear for the highest pressure ratio p 0 /p e = 10 ( figure 4, top). The hybrid solutions are close enough to the kinetic ones for all values of rarefaction and pressure ratios.  and the coupling procedures. Since the sum of CPU times relative to NS solution and coupling computation is quite small in comparison to the kinetic solution requirements, time t h is essentially dictated by the number of N K kinetic points should be solved during computation. Figure 8 shows the hybrid code speedup compared to the pure kinetic S-model code as a function of rarefaction parameter. The hybrid code offers a CPU time speedup ranging from 2, at a rarefaction level where NS solution would be completely inadequate, up to 5 in the slip flow regime. For example, pure kinetic computation at p 0 /p e = 10 and δ = 20 on physical mesh of 360×40 and velocity mesh of 40×40 requires 2.51 s, while hybrid time step t h is 1.3 s (number of kinetic points N K = 7050) and the NS step is t NS = 0.14 s. The total CPU advantage over the whole computation, with the dynamic coupling, is however appreciably greater, since at the beginning of the computation the kinetic area is typically much smaller than at converged state. It can be seen in figure 8 that with decrease in δ the speedup decreases, since the number of kinetic points in Ω K increases.

Conclusion
A hybrid algorithm based on the direct numerical solution of the kinetic S-model equation coupled with a Navier-Stokes model was applied to the numerical investigation of a gas flow through the channel. It was shown that hybrid code gives results close to full kinetic solutions for flow regimes where the Navier-Stokes solution clearly fails. It was found that the breakdown criterion based on the modified gradient length Knudsen number is more suitable for the boundary value problem. The obtained speed up, with respect to full S -model solutions, is between 2 and 5, although strongly depends on the size of the kinetic region. Moreover, the accuracy of the NS solution provided with slip boundary conditions is reasonable up to Kn ≤ 0.05 (deviation from kinetic solution does not exceed 5%). Thus, for higher Kn the use of first order slip boundary condition is not recommended.