When the flow is sufficiently rarefied, a temperature gradient, for example, between two walls separated by a few mean free paths, induces a gas flow—an observation attributed to the thermostress convection effects at the microscale. The dynamics of the overall thermostress convection process is governed by the Boltzmann equation—an integrodifferential equation describing the evolution of the molecular distribution function in six-dimensional phase space—which models dilute gas behavior at the molecular level to accurately describe a wide range of flow phenomena. Approaches for solving the full Boltzmann equation with general intermolecular interactions rely on two perspectives: one stochastic in nature often delegated to the direct simulation Monte Carlo (DSMC) method and the others deterministic by virtue. Among the deterministic approaches, the discontinuous Galerkin fast spectral (DGFS) method has been recently introduced for solving the full Boltzmann equation with general collision kernels, including the variable hard/soft sphere models—necessary for simulating flows involving diffusive transport. In this work, the deterministic DGFS method, Bhatnagar-Gross-Krook (BGK), Ellipsoidal statistical BGK (ESBGK), and Shakhov kinetic models, and the widely used stochastic DSMC method, are utilized to assess the thermostress convection process in micro in-plane Knudsen radiometric actuator—a microscale compact low-power pressure sensor utilizing the Knudsen forces. The BGK model underpredicts the heat-flux, shear-stress, and flow speed; the S-model overpredicts; whereas, ESBGK comes close to the DSMC results. On the other hand, both the statistical/DSMC and deterministic/DGFS methods, segregated in perspectives, yet, yield inextricable results, bespeaking the ingenuity of Graeme Bird who laid down the foundation of practical rarefied gas dynamics for microsystems.

## I. INTRODUCTION

In microscale flows, the length scale dictates the type of forces governing the physical phenomena. The surface to the volume ratio is high and hence the surface forces dominate. The Reynolds number is low and the viscous shear stresses are significantly increased.^{1} Under sufficiently rarefied flow conditions, an application of temperature gradient, say, between two parallel plates separated by few mean free paths, induces a low velocity gas flow commonly identified as thermostress convection effects.^{2} A necessary condition to induce a sufficiently *useful* gaseous velocity requires the characteristic length scale of the thermal gradients *T*/|∇_{x}*T*| to be comparable to the molecular mean free path *λ*. At the macroscale, such magnitudes are prohibitive, necessitating thermal gradients on the order of 10^{6} K/m. However, at the microscale, such conditions are readily achieved allowing the thermostress effects to overcome the classically dominant viscous forces.^{3}

From a historical and experimental viewpoint, Knudsen, in 1910, explored the possibility of gas actuation under the influence of temperature gradients using evacuated glass bulbs separated by a long narrow tube, wherein heating one of the bulbs resulted in a pumping action creating a high pressure at the hot end and low pressure at the cold end.^{4,5} In 1950s,^{6} Knudsen carried out various experiments using Crooke’s radiometer,^{7} wherein a device consisting of a long thin and narrow platinum band with dark (hot) and bright (cold) sides, in a rarefied environment, exhibits a net force due to momentum imbalance of particles reflecting from the dark and bright sides. Without being exhaustive, we refer to Ref. 8 for a comprehensive review of the radiometric phenomenon. From a theoretical viewpoint, Maxwell hypothesized that one of the possible causes of radiometric effects are temperature stresses. However, based on linearized kinetic theory and corresponding reduced macroscopic equations of motion (see Sec. 15 in Ref. 9), the author concluded that no motion can be produced by temperature stresses,^{2,9} which, in general, is incorrect. Later, Kogan, in 1976, introduced the theory of thermostress convection, wherein the bulk velocity is attributed to the presence of higher order terms of temperature stresses [see Eq. (2.6) in Ref. 2], arrived in part by the second order Chapman-Enskog expansion commonly identified as the Burnett approximations. In the multispecies context, however, the phenomenon and the effect of thermostress convection on the flow concentration (and the subsequent induced velocity) is more apparent.

Chapman,^{10} as early as 1953, developed the theory of diffusion processes [see Eq. (8.4, 7) in Ref. 10 again derived using Chapman-Enskog expansion] wherein the difference in concentrations of two species is proportional to the thermal gradient term *k*_{T}∇ ln *T*, where *k*_{T} is the thermal diffusion factor. At normal conditions, this coefficient is very low and is therefore not accounted for in practice. For instance, as a classical example, Bird^{11} devised a self-diffusion test case (see Sec. 12.6) where the diffusion coefficient was measured by ignoring the thermal gradient term *k*_{T}∇ ln *T* of Eq. (8.4, 7) in Ref. 10. Note, however, that there is considerable thermal gradient in self-diffusion cases; see for instance Ref. 12, where we presented the results for temperature variation for self-diffusion cases. Although the temperature gradient is unaccounted for, the diffusion coefficient, which is measured by a self-diffusion simulation, matches well with the experimentally^{10} observed diffusion coefficient. This suggests that *k*_{T} is potentially low—which is indeed the case, for instance, see Ref. 13, wherein the authors noted thermal diffusion coefficient on order of 10^{−3}. In microscale flows where the per unit temperature drop can easily reach 10^{6} K/m, as noted earlier, *k*_{T}∇ ln *T* can have appreciable contributions. This type of process has been interpreted in terms of thermostress convection due to concentration inhomogeneities by Kogan.^{2} The overall thermostress convection phenomenon/effect is highly coupled and exhibits highly rich flow structures (as will be shown in Sec. IV C), and an in-depth understanding can prove to be very useful for development of next generation of microsystems.

To summarize, Sone^{14} identified three broad groups of the temperature driven flow based on its application in microsystems: (a) thermal creep flow^{15–17} which is an induced flow around a body with nonuniform temperature; (b) thermal stress slip flow, which is induced by nonuniform temperature gradient over the boundary;^{2,18–23} and (c) nonlinear thermal stress flow,^{2} which is important only when the temperature gradient in the gas is high, and nonlinear terms of temperature variations in stress tensor should be taken into account. The present study is delegated to the third, i.e., nonlinear thermal stress flow.

From a practical engineering viewpoint, thermostress convection has been applied for microstructure actuation. Passian,^{24,25} in 2003, demonstrated a microcantilever suspended over a substrate, which when heated via a pulsed laser generated deflections at the cantilever tip as a consequence of the Knudsen forces in the gap between the substrate and microcantilever. Foroutan,^{26} in 2014, demonstrated untethered levitation in concave microflying robots relying on the Knudsen force. The phenomenon has been further explored in small satellite and spacecraft attitude control devices^{27} and high-altitude propulsion systems.^{28}

The dynamics of the overall thermostress convection process is governed by the Boltzmann equation—an integrodifferential equation describing the evolution of the distribution function in six-dimensional phase space—which models the dilute gas behavior at the molecular level to accurately describe a wide range of noncontinuum flow phenomena. In the present work, we assess the thermostress convection process using the fundamental microscopic full Boltzmann equation. The approaches for numerical solution of the Boltzmann equation date back to as early as the 1940s.^{29} However, it was not until the 1960s that the numerical simulations were feasible. In practice, the numerical simulations of the Boltzmann equation were made possible by the introduction of direct simulation Monte Carlo (DSMC) method.^{11,30} Over sufficient small intervals, by decoupling the molecular motion and interaction processes, DSMC first advects the particles deterministically according to their velocities, also termed *free transport*, and then describes the collisions by statistical models with a specified interaction potential.

The choice of the interaction potential substantially affects the simulation fidelity and computational complexity. Early implementations of the DSMC method relied on the purely repulsive hard sphere (HS) interaction model.^{30} The HS model, however, deviates from experimental observations for common gases^{31} due to a square-root viscosity variation with temperature. The variable hard sphere (VHS) model proposed by Bird^{11} results in a more general power-law viscosity variation with temperature and has been widely used for DSMC simulations of *single-species* gas flows due to its computational efficiency and ease of implementation. The VHS model, however, deviates from experimental observations for common *multispecies* flows^{32,33} involving diffusive transport. Later, several variations of the VHS model were proposed, including, the variable soft sphere (VSS),^{33} M-1,^{34} generalized soft sphere (GSS),^{35} all of which belong to a class of repulsive interactions. The VSS model modifies the scattering law of the VHS model by using a scattering parameter (*α*) that allows reproduction of measured diffusion coefficients in addition to the viscosity coefficient. The M-1 model is a modification of the VHS model to have a linear distribution of scattering angles in terms of the impact parameter. This modification allows M-1 to reproduce correct viscosity and diffusivity without the need of an additional parameter (*α*).^{36} The GSS model, although general, needs additional parameters for reproducing the viscosity and diffusion coefficients (see Refs. 11 and 37 for additional details/equations for these models). In particular, the complexity of the DSMC algorithm is independent of the number of species in the mixture, as well as the mass of the individual species. This makes the method highly useful and efficient for modeling sufficiently *fast*^{2} nonequilibrium flows. From a usage perspective, there is a growing number of applications requiring DSMC simulations. Statistics (see Fig. 1) show that over 1000 papers on DSMC are now published every year. Many of these papers are based on the codes made available to the research community in the 1990s by Bird. Over the years—in accordance with the predictions of Moore’s law—the number of collisions performed per hour have increased exponentially.

However, it is the stochastic nature of the DSMC that introduces high statistical noise in low-speed flows. In the present work, we study the thermostress convection process using the recently developed deterministic discontinuous Galerkin fast spectral (DGFS) method^{12,38} as well as DSMC: the primary tool for rarefied flow simulations. DGFS allows arbitrary unstructured geometries; high order accuracy in physical space time, and velocity space; arbitrary collision kernels, including, the well-known VSS model^{12} and provides excellent nearly linear scaling characteristics on massively parallel architectures.^{39,40} DGFS produces noise-free solutions and can simulate low-speed flows encountered in thermostress convection dominated devices.

From a flow modelling viewpoint, Loyalka,^{41} using a linearized Boltzmann equation, calculated the longitudinal and transversal Knudsen forces on the cylindrical surfaces of a hanging wire of a vacuum microbalance. The authors noted the Knudsen force maximum in the transitional regime for helium–an observation attributed to the bimodal nature of radiometric forces.^{42} Fierro^{43} studied the problem using a Bhatnagar-Gross-Krook (BGK) model for the range of Knudsen numbers and different molecular species noting an inverted parabolic profile for variation of Knudsen force with pressure (which can be reinterpreted in terms of Knudsen number since a fixed size geometry was used for all cases). The authors observed a peak Knudsen force in the 10–100 N/m^{2} pressure range for helium, krypton, hydrogen, oxygen, and carbon dioxide. Alexeenko^{44} carried out numerical simulations around heated microbeams using the conventional Navier-Stokes incorporating the first order Maxwell slip and Smoluchowski temperature jump boundary conditions, DSMC, and primarily using a deterministic kinetic ellipsoidal statistical Bhatnagar-Gross-Krook (ESBGK) model employing a finite-difference-discrete-velocity scheme. The gas-damping coefficients on a moving microbeam for quasistatic isothermal conditions were estimated by the three numerical methods for Knudsen numbers from 0.1 to 1.0. It was concluded that the Navier-Stokes simulations overestimate the gas-damping force for Knudsen numbers larger than 0.1, while the ESBGK and DSMC methods are in good agreement for the slip and transitional flow regimes. Moreover, the Knudsen force peaks in the transitional regime at Kn ≈ 2, and the numerically predicted variation of the force is consistent with experimental observations of the displacement of a heated micro-beam. Zhu^{45} analyzed the problem specifically using DSMC in the slip, transition, and free molecular regimes noting *qualitative* agreements between DSMC and experimental results of Passian.^{24,25} Nabeth^{46} analyzed the problem using the ESBGK model within a finite volume framework. Notably, the authors devised a semiempirical relation between the force and the Knudsen number based on dynamic similarity. Anikin^{47} studied the radiometric forces via a direct solution of the Boltzmann equation on 2-D velocity grids via a discrete ordinate projection method.^{48} More recently, Lotfian^{49} analyzed the various arrangements for radiometric pumps featuring vane and ratchet structures, including, zigzag triangular fins, using DSMC and the finite volume based BGK-Shakhov model.

In more complex scenarios, one can stack an array of microheaters to significantly enhance the Knudsen force output.^{3,50,51} Strongrich^{51} demonstrated the possibility of amplifying the Knudsen forces as well as reversing its direction by combining thermal gradients between several solid bodies. The idea was further explored, resulting in the development of a Microscale In-Plane Knudsen Radiometric Actuator (MIKRA) sensor for flow actuation and measurement.^{3,52,53} MIKRA consists of an array of hot and cold microbeams termed the heater and shuttle arm. When the heater arm is heated under the application of an electric current, the Knudsen force is generated in the gap between the shuttle and heater arm. The displacement of the shuttle arm is then measured using a capacitor (specific details to follow in Sec. IV). MIKRA presents an interesting problem for analyzing thermostress convection due to temperature gradients as well as concentration inhomogeneties; see Ref. 54 where authors observed species separation in MIKRA which *might* be, in part, due to be the effect of the *k*_{T}∇ ln *T* term. We believe it is too early to make a definite conclusion on the topic.

A key question, and a subject of ongoing research is the following: How well can the kinetic equations/methods/models, for instance, the McCormack model,^{55} Lattice Boltzmann method (LBM),^{56} Bhatnagar-Gross-Krook (BGK),^{57–61} Ellipsoidal statistical Bhatnagar-Gross-Krook (ESBGK),^{62,63} BGK-Shakhov (S-model),^{64} Unified Gas Kinetic Scheme (UGKS),^{65,66} Discontinuous Galerkin Fast Spectral (DGFS),^{38} and direct simulation Monte Carlo (DSMC),^{11} describe the thermostress convection process, including, their applicability regimes at a wide range of rarefaction levels and temperature gradients, and required computational cost for reproducing the correct induced low speed velocity profile on *common standard* benchmark problems such as MIKRA where the experimental results are readily available. As noted by Kogan,^{2} the overall thermostress convection process is *complicated* function of concentration, of mass-ratio, molecule-collision cross section, etc. An in-depth understanding of the overall thermostress convection process at the microscale may potentially prove useful for development of a series of new MEMS devices without any moving parts (see, for instance, Refs. 52 and 67). This paper, in part, focuses on quantifying the fidelity of results recovered from BGK, ESBGK, S-model, DGFS and DSMC for the Knudsen radiometric actuator MIKRA. To the best of our knowledge, analysis of single/multispecies complex flows such as MIKRA has not been carried out using deterministic full Boltzmann.

Finally, we mention that there is another class of low-speed gas flow problems which appears within the context of gas damping:^{68,69} an effect generally observed in radio frequency (RF) switches^{70} which are used in radar systems, wireless communication systems, and other instrumentation fields. Gas damping sharply influences the dynamic behavior of MEMS devices, including, mechanical quality factors of microfabricated resonators, switching time, impact velocity, and bounceback of contacting MEMS.^{71} A commonly identified trade-off in gas damping is as follows: the gas damping must be minimized to achieve high sensitivity of MEMS, such as resonators, and it must be maximized to mitigate the shock response and transient performance of MEMS. This is yet another class of problems where the fidelity of different modelling approaches, in particular, the deterministic ones, can be tested.

The rest of this paper is organized as follows. In Sec. II, we give an overview of the multispecies Boltzmann equation, the self/cross collision integrals, and the phenomenological VHS/VSS collision kernels used in practical engineering applications. Extensive numerical verification for the BGK, ESBGK, S-model, and DGFS against DSMC are performed and discussed in Sec. III. Section IV provides the description, problem statement, and results for the thermostress convection enabled MIKRA sensor. Section V presents the analysis of *multispecies* thermostress convection in the MIKRA sensor. Concluding remarks are given in Sec. VI.

## II. BOLTZMANN EQUATION

In this section, we give a brief overview of the multispecies Boltzmann equation. Readers are referred to Ref. 12 for more details.

Suppose we consider a gas mixture of *s* species (*s* ≥ 2), each represented by a distribution function *f*^{(i)}(*t*, *x*, $v$), where *t* ≥ 0 is the time, $x\u2208\Omega \u2282R3$ is the position, and $v\u2208R3$ is the particle velocity [*f*^{(i)}d*x*d$v$ gives the number of particles of species *i* to be found in an infinitesimal volume d*x*d$v$ centered at the point (*x*, $v$) of the phase space]. The time evolution of *f*^{(i)} is described by the multispecies Boltzmann equation written as^{72,73}

Here $Q(ij)$ is the collision operator that models the binary collisions between species *i* and *j*, and acts only in the velocity space

where ($v$, $v*$) and ($v$′, $v*\u2032$) denote the precollision and postcollision velocity pairs. During collisions, the momentum and energy are conserved

where *m*_{i}, *m*_{j} denote the mass of particles of species *i* and *j*, respectively. Hence, one can parameterize $v$′ and $v*\u2032$ as follows:

with *σ* being a vector varying on the unit sphere *S*^{2}. $Bij=Bji(\u22650)$ is the collision kernel characterizing the interaction mechanism between particles. It can be shown that

where *χ* is the deviation angle between $v$ − $v*$ and $v$′ − $v*\u2032$.

Given the interaction potential between particles, the specific form of *B*_{ij} can be determined using the classical scattering theory

where Σ_{ij} is the differential cross section given by

with *b*_{ij} being the impact parameter.

With a few exceptions, the explicit form of Σ_{ij} can be hard to obtain since *b*_{ij} is related to *χ* implicitly. However, as stated in the Introduction, the choice of the interaction potential substantially affects the simulation fidelity and computational complexity. Proposed as a modification of the Bird’s VHS model, Koura and Matsumoto^{33} introduced the so-called VSS model by assuming

where *α*_{ij} is the scattering parameter and *d*_{ij} is the diameter borrowed from the VHS model [Eq. (4.79) in Ref. 11]

with Γ being the Gamma function, $\mu ij=mimjmi+mj$ the reduced mass, *d*_{ref,ij}, *T*_{ref,ij}, and *ω*_{ij}, respectively, the reference diameter, reference temperature, and viscosity index. Substituting Eqs. (7)–(9) into (6), one can obtain *B*_{ij} as

where $b\omega ij,\u2009\alpha ij$ is a constant given by

In particular, the VHS kernel is obtained when *α*_{ij} = 1 and 0.5 ≤ *ω*_{ij} ≤ 1 (*ω*_{ij} = 1: Maxwell molecules; *ω*_{ij} = 0.5: HS) and the VSS kernel is obtained when 1 < *α*_{ij} ≤ 2 and 0.5 ≤ *ω*_{ij} ≤ 1.

Given the distribution function *f*^{(i)}, the number density, mass density, velocity, and temperature of species *i* are defined as

The total number density, mass density, and velocity are given by

Furthermore, the diffusion velocity, stress tensor, and heat flux vector of species *i* are defined as

where *c* = $v$ − *u* is the peculiar velocity. Finally, the total stress, heat flux, pressure, and temperature are given by

### A. Stochastic modelling

From a stochastic viewpoint, DSMC, as introduced by Bird,^{11,30} incorporates four principal steps: (a) index, (b) move, (c) collide, and (d) sample. The flowchart of a standard DSMC algorithm is illustrated in Fig. 2. DSMC uses a spatial grid to contain the simulated molecules and perform sampling. The algorithm starts with the distribution of molecules in the spatial domain according to the prespecified initial condition: bulk velocity *u*_{ini}, temperature *T*_{ini}, and number densities $nini(i)$. The individual molecules must be tracked, and therefore an indexing mechanism is used to track which molecules are in which cell of the spatial domain. Repeated calling of the index subroutine is necessitated by molecular movement.

The *move* subroutine moves each molecule according to their velocity a distance appropriate for the specified time step. This velocity is assumed constant over each time step. The velocity of a molecule is changed either by external forces, such as electrostatic forces, or by the scattering resulting from a collision. In the absence of external forces, the velocity will only change as a result of a collision. In the present work, we consider elastic collisions, i.e., collisions in which the total kinetic energy is unchanged.

The *collide* subroutine randomly selects a pair of molecules to collide using the acceptance-rejection method.^{11,37} The collision is accepted with a probability

where $(|v\u2212v*|\u2009\Sigma ij)max$ is the maximum, effective volume swept out by a molecule. This maximum value is recorded in each cell such that each cell may have a separate collision frequency. Taking this into account, the goal of the *collide* is to determine scattering angles and postcollision energies, as well as to obtain correct collision frequencies and microscopic properties.

Finally, *sample* performs sampling over all cells to determine macroscopic properties. The microscopic properties from each simulated molecule, for instance, molecular velocities and translational energy, are averaged in each cell to compute the macroscopic properties such as density, bulk velocity, pressure, translational temperature, etc. For example, in a given *cell*, consider $Mcell(i)$ particles of species *i*, the mass density of the species is simply the total mass per unit cell volume, i.e.,

Details about sampling of other macroscopic properties, specifically in the stochastic context, can be found in Ref. 11.

### B. Deterministic modelling

From a deterministic viewpoint, we use the recently introduced discontinuous Galerkin fast spectral (DGFS) method.^{12,38} DGFS directly approximates the Boltzmann equation (1), where the transport term (spatial derivative) is discretized by the classical DG method and the collision term (integral in $v$) is discretized by the fast Fourier spectral method.^{12,74} The discretized system is then advanced in time using the Runge-Kutta method.

The coupling of two kinds of methods (DG in the physical space and spectral method in the velocity space) is possible due to the special structure of the Boltzmann equation—the collision operator acts only in $v$ wherein *t* and *x* can be treated as parameters. Simply speaking, given the distribution functions *f*^{(i)} and *f*^{(j)} of species *i* and *j* at *N*^{3} velocity grid, the fast Fourier spectral method produces $Q(ij)(f(i),f(j))$ at the same grid with *O*(*MN*_{ρ}*N*^{3} log *N*) complexity, where *M* ≪ *N*^{2} is the number of discretization points on the sphere and *N*_{ρ} ∼ *O*(*N*) is the number of Gauss-Legendre quadrature/discretization points in the radial direction needed for low-rank decomposition. Further details can be found in Refs. 12 and 38.

The overall DGFS method is simple from mathematical and implementation perspective; highly accurate in both physical and velocity spaces as well as time, robust, i.e., applicable for general geometry and spatial mesh, exhibits nearly linear parallel scaling, and directly applies to general collision kernels needed for high fidelity modelling. Due to these features, we use DGFS for deterministic modelling of flows considered in this work.

## III. BGK/ESBGK/S-MODEL/DGFS: VERIFICATIONS

Due to the nonlinearity and the complexity of Boltzmann collision term $Q$, the collision operator is often simplified for practical reasons—a major motivation behind the development of kinetic models.^{57,59,62,64,65,75} In this section, we shall restrict our discussion to the single-species system, i.e., *s* = 1, *i* = {1}. We will drop superscripts ^{(i)} and ^{(ij)} for simplicity.

While devising kinetic models for single species system, the collision term $S$—in this work, we denote kinetic models by symbol $S$ to differentiate it from the full Boltzmann collision integral $Q$—is expected to have the following four properties:^{72,75,76}

It guarantees the conservation of mass, momentum, and energy, i.e.,

The entropy production is always positive, i.e.,

Due to specific form of $S$, the phase density in equilibrium is a Maxwellian, i.e.,

where

The Prandtl number is close to 2/3 for monoatomic gases, i.e.,

where *μ* and *κ*, respectively, refer to the dynamic viscosity and thermal conductivity.

Among popular kinetic models, BGK/ESBGK collision operators are relaxation type kernels given as

where *f*_{γ} is the local equilibrium function, and *ν* = Pr *p*/*μ* is the collision frequency. Here, *p* denotes pressure. For BGK, *f*_{γ} is a local Maxwellian given as

whereas for ESBGK, *f*_{γ} is anisotropic Gaussian given as

where Id is an identity matrix. For the S-model, *f*_{γ} is given as

It can be easily shown that the BGK, ESBGK, and S-model satisfy conditions (17)–(19). The S-model and ESBGK satisfy (21), whereas BGK does not.

To put things more concretely, we consider five Fourier-Couette flow cases and a flow over a microelectronic chip to verify the BGK/ESBGK/S-model/DGFS method against DSMC.

### A. Verification: Fourier-Couette flows

In the current test case, we consider the effect of velocity and temperature gradients on the solution. The coordinates are chosen such that the walls are parallel to the *y* direction and *x* is the direction perpendicular to the walls. The geometry as well as boundary conditions are shown in Fig. 3. Specific case details have been provided in Tables I and II. Figure 4 illustrates the velocity and temperature along the domain length, wherein we note an excellent agreement between DGFS and DSMC. The velocity profiles from BGK/ESBGK are in good agreement with DGFS and DSMC, whereas the temperature profiles from ESBGK are in good agreement with DGFS and DSMC. The deviation in BGK temperature profiles is due to its Prandtl number defect.

Common parameters . | . |
---|---|

Molecular mass: m_{1} (×10^{27} kg) | 66.3 |

Nondimentional physical space | [0, 1] |

Spatial elements | 2 |

DG order | 3 |

Time stepping | Euler |

Viscosity index: ω | 0.81 |

Scattering parameter: α | 1.4 |

Ref. diameter: d_{ref} (×10^{10} m) | 4.17 |

Ref. temperature: T_{ref} (K) | 273 |

Ref. viscosity: μ_{ref} (Pa s) | 2.117 × 10^{−5} |

Characteristic mass: m_{0} (×10^{27} kg) | 66.3 |

Characteristic length: H_{0} (mm) | 1 |

Characteristic velocity: u_{0} (m/s) | 337.2 |

Characteristic temperature: T_{0} (K) | 273 |

Characteristic no. density: n_{0} (m^{−3}) | 3.538 × 10^{22} |

Initial conditions | |

Velocity: u (m/s) | 0 |

Temperature: T (K) | 273 |

Number density: n (m^{−3}) | 3.538 × 10^{22} |

Knudsen number:^{a} (Kn) | 0.036 |

Common parameters . | . |
---|---|

Molecular mass: m_{1} (×10^{27} kg) | 66.3 |

Nondimentional physical space | [0, 1] |

Spatial elements | 2 |

DG order | 3 |

Time stepping | Euler |

Viscosity index: ω | 0.81 |

Scattering parameter: α | 1.4 |

Ref. diameter: d_{ref} (×10^{10} m) | 4.17 |

Ref. temperature: T_{ref} (K) | 273 |

Ref. viscosity: μ_{ref} (Pa s) | 2.117 × 10^{−5} |

Characteristic mass: m_{0} (×10^{27} kg) | 66.3 |

Characteristic length: H_{0} (mm) | 1 |

Characteristic velocity: u_{0} (m/s) | 337.2 |

Characteristic temperature: T_{0} (K) | 273 |

Characteristic no. density: n_{0} (m^{−3}) | 3.538 × 10^{22} |

Initial conditions | |

Velocity: u (m/s) | 0 |

Temperature: T (K) | 273 |

Number density: n (m^{−3}) | 3.538 × 10^{22} |

Knudsen number:^{a} (Kn) | 0.036 |

Parameter . | Case FC-01 . | Case FC-02 . | Case FC-03 . | Case FC-04 . | Case FC-05 . |
---|---|---|---|---|---|

Nondimentional velocity space^{a} | [−5, 5]^{3} | [−5, 5]^{3} | [−5, 5]^{3} | [−5, 5]^{3} | [−8, 8]^{3} |

{N^{3}, N_{ρ}, M}^{b} | {24^{3}, 6, 6} | {24^{3}, 6, 6} | {24^{3}, 6, 6} | {24^{3}, 6, 6} | {32^{3}, 16, 6} |

Left wall (purely diffuse) conditions | |||||

Velocity: u_{l} (m/s) | (0, −50, 0) | (0, −50, 0) | (0, −250, 0) | (0, −250, 0) | (0, −250, 0) |

Temperature: T_{l} (K) | 273 | 223 | 273 | 223 | 173 |

Right wall (purely diffuse) conditions | |||||

Velocity: u_{r} (m/s) | (0, 50, 0) | (0, 50, 0) | (0, 250, 0) | (0, 250, 0) | (0, 250, 0) |

Temperature: T_{r} (K) | 273 | 323 | 273 | 323 | 373 |

Parameter . | Case FC-01 . | Case FC-02 . | Case FC-03 . | Case FC-04 . | Case FC-05 . |
---|---|---|---|---|---|

Nondimentional velocity space^{a} | [−5, 5]^{3} | [−5, 5]^{3} | [−5, 5]^{3} | [−5, 5]^{3} | [−8, 8]^{3} |

{N^{3}, N_{ρ}, M}^{b} | {24^{3}, 6, 6} | {24^{3}, 6, 6} | {24^{3}, 6, 6} | {24^{3}, 6, 6} | {32^{3}, 16, 6} |

Left wall (purely diffuse) conditions | |||||

Velocity: u_{l} (m/s) | (0, −50, 0) | (0, −50, 0) | (0, −250, 0) | (0, −250, 0) | (0, −250, 0) |

Temperature: T_{l} (K) | 273 | 223 | 273 | 223 | 173 |

Right wall (purely diffuse) conditions | |||||

Velocity: u_{r} (m/s) | (0, 50, 0) | (0, 50, 0) | (0, 250, 0) | (0, 250, 0) | (0, 250, 0) |

Temperature: T_{r} (K) | 273 | 323 | 273 | 323 | 373 |

### B. Verification: Flow around a micro-electronic chip

In the current test case, we consider the effect of temperature gradients on a solid substrate placed in a rarefied environment. The problem schematic, geometry, as well as boundary conditions are shown in Fig. 5. Case details have been provided in Table IV. We mention that, to the best of our knowledge, an analysis of such a low-temperature flow has not been carried out previously within a deterministic full Boltzmann framework.

#### 1. Numerical details

We employ DSMC and DGFS to carry out simulation of flow around a micro-electronic chip. The simulation specific numerical parameters as well as differences between stochastic (DSMC) and deterministic (DGFS) modelling is described next.

**DSMC**: SPARTA^{77}has been employed for carrying out DSMC simulations in the present work. It implements the DSMC method as proposed by Bird.^{11}The solver takes into account the translational/rotational/vibrational kinetic energies associated with the molecular motion. The solver has been benchmarked^{77}and widely used for studying hypersonic, subsonic and thermal^{12,36,38,39,78–81}gas flow problems. In this work, cell size less than*λ*/3 has been ensured in all the test cases. The no-time collision (NTC) algorithm is used in conjunction with Bird’s VHS scattering model. The simulations are first run for 200 000 unsteady steps wherein the particles move, collide, and allowed to equilibrate. No sampling is performed at this stage. Next, the simulation is run for another 4 000 000 steady steps wherein the samples of flow properties namely number density, flow velocity, temperature, stress, and heat-flux, are taken for sufficiently long time so as to produce a meaningful bulk properties as well as minimize the statistical noise therein. In the present case, the DSMC domain is discretized with a uniform cell size of 0.2*μ*m, with 300 particles per cell on average during initialization. A time step of 10^{−9}s is used during*move*step of DSMC algorithm throughout the course of simulation. N_{2}is used as the working gas in simulations. The properties of the working gas is given in Table III. We want to emphasize that for DSMC simulations, we take rotational/vibrational degrees of freedom into account, i.e., N_{2}is treated as a diaotomic species. DSMC simulations on 30 cores of Intel(R) Xeon(R) CPU E5-2670 v2 2.50GHz, took ∼73 h.**DGFS**: We use the DGFS implementation described in Ref. 38. The spatial domain consists of 281 uniform square cells of 1*μ*m each. Since we are seeking a steady state solution, the time-step is selected based on the Courant-Friedrichs-Lewy (CFL) constraints of the forward Euler scheme. Other case specific DGFS parameters have been provided in Table IV. Note that, we employ N_{2}as the working gas in simulations, since MIKRA experiments^{3}were performed in N_{2}medium. N_{2}is diatomic, however, DGFS, as of now, is applicable for monoatomic gases only. Since the working temperature range is low, we anticipate the effects of vibrational degrees of freedom to be negligible. DGFS simulations on 2 Nvidia-P100 GPUs took ∼9 h.

Mass: m (kg) | 46.5 × 10^{−27} |

Viscosity index: ω (−) | 0.74 |

Scattering index: α (−) | 1.0 |

Ref. diameter: d_{ref} (m) | 4.17 × 10^{−10} |

Ref. temperature: T_{ref} (K) | 273 |

Ref. viscosity: μ_{ref} (Pa s) | 1.656 × 10^{−5} |

DSMC specific parameters | |

Rotational degrees of freedom: ζ_{R} (−) | 2 |

Rotational relaxation: Z_{R} (−) | 2 |

Vibrational degrees of freedom: ζ_{V} (−) | 2 |

Vibrational relaxation Z_{V} (−) | 1.901 14 × 10^{−5} |

Vibrational temperature T_{V} (K) | 3371 |

Mass: m (kg) | 46.5 × 10^{−27} |

Viscosity index: ω (−) | 0.74 |

Scattering index: α (−) | 1.0 |

Ref. diameter: d_{ref} (m) | 4.17 × 10^{−10} |

Ref. temperature: T_{ref} (K) | 273 |

Ref. viscosity: μ_{ref} (Pa s) | 1.656 × 10^{−5} |

DSMC specific parameters | |

Rotational degrees of freedom: ζ_{R} (−) | 2 |

Rotational relaxation: Z_{R} (−) | 2 |

Vibrational degrees of freedom: ζ_{V} (−) | 2 |

Vibrational relaxation Z_{V} (−) | 1.901 14 × 10^{−5} |

Vibrational temperature T_{V} (K) | 3371 |

Parameters . | MEC-01 . |
---|---|

Spatial elements | 190 quadrilaterals |

DG order | 3 |

Time stepping | Euler |

Points in velocity mesh: N^{3} | 24^{3} |

Points in radial direction:^{a} N_{ρ} | 6 |

Points on half sphere:^{a} M | 6 |

Size of velocity mesh^{b} | [−5, 5]^{3} |

Characteristic length: H_{0} (μm) | 3 |

Characteristic velocity: u_{0} (m/s) | 402.54 |

Characteristic temperature: T_{0} (K) | 273 |

Characteristic no. density: n_{0} (m^{−3}) | 4.894 × 10^{23} |

Initial conditions | |

Velocity: u (m/s) | 0 |

Temperature: T (K) | 273 |

Number density: n (m^{−3}) | 4.894 × 10^{23} |

Knudsen number:^{c} (Kn) | 0.881 58 |

Parameters . | MEC-01 . |
---|---|

Spatial elements | 190 quadrilaterals |

DG order | 3 |

Time stepping | Euler |

Points in velocity mesh: N^{3} | 24^{3} |

Points in radial direction:^{a} N_{ρ} | 6 |

Points on half sphere:^{a} M | 6 |

Size of velocity mesh^{b} | [−5, 5]^{3} |

Characteristic length: H_{0} (μm) | 3 |

Characteristic velocity: u_{0} (m/s) | 402.54 |

Characteristic temperature: T_{0} (K) | 273 |

Characteristic no. density: n_{0} (m^{−3}) | 4.894 × 10^{23} |

Initial conditions | |

Velocity: u (m/s) | 0 |

Temperature: T (K) | 273 |

Number density: n (m^{−3}) | 4.894 × 10^{23} |

Knudsen number:^{c} (Kn) | 0.881 58 |

#### 2. Results and discussion

Figure 6 illustrate the contours of various flow properties for the flow around the solid chip/substrate. Ignoring the statistical noise, we observe excellent agreement between DSMC and DGFS. In particular, DGFS reproduces noise-free smooth results.

Next, we compute the force acting on the substrate as a result of the temperature gradients initially present in the flow. In general, the pressure force on a surface is given by

where *n* is the unit surface normal, *p* is the pressure on the surface, and *A* is the area of the surface.

Table V presents the *x*-component of force on the substrate for the micro-electronic chip verification case. Again, we note reasonable agreement between the values recovered from DSMC and DGFS simulations.

### C. Verification: Flow in short microchannels

The present test case closely follows case-I(a) from Ref. 82. In the current test case, two reservoirs filled with N_{2} gas, at different temperatures, are connected by a two-dimensional capillary tube, both with a finite length *L* and height *H*/2, are considered. The problem schematic, geometry, as well as boundary conditions are shown in Fig. 7. Case details have been provided in Table VI. Note in particular, we introduce a linearly decreasing temperature profile at the top wall. Previous studies for flow in short microchannels have been restricted to model kinetic equations^{82–84} or moment methods.^{85} To the best of our knowledge, an analysis of such a low-temperature flow has not been carried out previously within a deterministic full Boltzmann framework.

Parameters . | SM-01 . |
---|---|

Spatial elements | 127 quadrilaterals |

DG order | 3 |

Time stepping | Euler |

Points in velocity mesh: N^{3} | 32^{3} |

Points in radial direction:^{a} N_{ρ} | 8 |

Points on half sphere:^{a} M | 6 |

Size of velocity mesh^{b} | [−5.72, 5.72]^{3} |

Characteristic length: H_{0} (μm) | 1 |

Characteristic velocity: u_{0} (m/s) | 421.98 |

Characteristic temperature: T_{0} (K) | 300 |

Characteristic no. density: n_{0} (m^{−3}) | 6.62 × 10^{24} |

Initial conditions | |

Velocity: u (m/s) | 0 |

Temperature: T (K) | 300 |

Number density: n (m^{−3}) | 6.62 × 10^{24} |

Knudsen number:^{c} (Kn) | 0.2 |

Inlet condition | |

Velocity: u_{in} (m/s) | 0 |

Temperature: T_{in} (K) | 600 |

Number density: n_{in} (m^{−3}) | 3.31 × 10^{24} |

Pressure: p_{in} (N/m) | 27 420 |

Outlet condition | |

Velocity: u_{out} (m/s) | 0 |

Temperature: T_{out} (K) | 300 |

Number density: n_{out} (m^{−3}) | 6.62 × 10^{24} |

Pressure: p_{out} (N/m) | 27 420 |

Parameters . | SM-01 . |
---|---|

Spatial elements | 127 quadrilaterals |

DG order | 3 |

Time stepping | Euler |

Points in velocity mesh: N^{3} | 32^{3} |

Points in radial direction:^{a} N_{ρ} | 8 |

Points on half sphere:^{a} M | 6 |

Size of velocity mesh^{b} | [−5.72, 5.72]^{3} |

Characteristic length: H_{0} (μm) | 1 |

Characteristic velocity: u_{0} (m/s) | 421.98 |

Characteristic temperature: T_{0} (K) | 300 |

Characteristic no. density: n_{0} (m^{−3}) | 6.62 × 10^{24} |

Initial conditions | |

Velocity: u (m/s) | 0 |

Temperature: T (K) | 300 |

Number density: n (m^{−3}) | 6.62 × 10^{24} |

Knudsen number:^{c} (Kn) | 0.2 |

Inlet condition | |

Velocity: u_{in} (m/s) | 0 |

Temperature: T_{in} (K) | 600 |

Number density: n_{in} (m^{−3}) | 3.31 × 10^{24} |

Pressure: p_{in} (N/m) | 27 420 |

Outlet condition | |

Velocity: u_{out} (m/s) | 0 |

Temperature: T_{out} (K) | 300 |

Number density: n_{out} (m^{−3}) | 6.62 × 10^{24} |

Pressure: p_{out} (N/m) | 27 420 |

#### 1. Numerical details

**DSMC**: The no-time collision (NTC) algorithm is used in conjunction with Bird’s VHS scattering model. The simulations are first run for 500 000 unsteady steps wherein the particles move, collide, and allowed to equilibrate. Next, the simulation is run for another 100 000 steady steps wherein the samples of flow properties are taken. In the present case, the DSMC domain is discretized with a uniform cell size of 0.01*μ*m, with 30 particles per cell on average during initialization (Note that SPARTA uses hierarchical Cartesian grid over the simulation domain: used to track particles and to co-locate particles in the same grid cell for performing collision and chemistry operations. At the junction, where the walls join the inlet and outlet regions, one can identify two boundary cells. We further refine, specifically, these two boundary cells into 10 × 10 subcells. These two cells are unique, i.e., for each of these cells, the top face is marked as inlet, and the left face is marked as solid wall. The cell-size has been made smaller to avoid any potential leakage). A time step of 10^{−10}s is used. N_{2}is (Table III) used as the working gas in simulations.**DGFS**: The spatial domain consists of 127 nonuniform quadrilateral elements as shown in Fig. 7(b). Case specific DGFS parameters have been provided in Table VI.

#### 2. Results and discussion

Figure 8 illustrate the contours of various flow properties for the flow around the solid chip/substrate. Ignoring the statistical noise, we gain note excellent agreement between DSMC and DGFS. In particular, minor differences in *x*-component of heat-flux i.e., *Q*_{x} can be attributed to the fact that DSMC simulations consider rotational degrees-of-freedom of N_{2} into account, whereas DGFS does not.

Figure 9 shows the variation of flow properties over the vertical centerline, wherein we again observe an excellent agreement.

## IV. MIKRA: MICRO IN-PLANE KNUDSEN RADIOMETRIC ACTUATOR

MIKRA, acronym for Micro In-Plane Knudsen Radiometric Actuator, is a microscale compact low-power pressure sensor. A computer-aided design (CAD) representation of the device has been illustrated in Fig. 10. Simply speaking, the device consists of an array of (twelve) microbeams labelled the *Shuttle Arm* and *Heater Arm* in Fig. 10. The heater arm is heated, and a thermal motion is induced in the gap between the heater and the shuttle. Subsequently, the shuttle arm experiences forces on order of few micro-newtons. This force is commonly identified as Knudsen force. Depending on the temperature of the heater, the shuttle gets displaced, and this displacement is measured capacitively. The magnitude of displacement is then used to estimate the ambient pressure. Specific details on MIKRA can be found in Refs. 3, 52, and 54.

### A. Problem statement

The flow configuration is shown in Fig. 11. Consider the 2D uniform flow of N_{2} with freestream velocity *U*_{∞}, freestream temperature *T*_{∞}, and freestream pressure *p*_{∞} over two two-dimensional square vanes, each with side lengths of 50 *μ*m, separated by a gap of 20 *μ*m (also used as the nondimensionalizing length scale). The vanes are modeled as purely diffuse solid walls. The left vane, indicated in blue, is kept at a lower/cold temperature which we denote by *T*_{C}. The right vane, indicated in red, is kept at a higher/hot temperature which we denote by *T*_{H}. The substrate, indicated in green, forms the lower boundary of the domain, and is modelled as a purely diffuse solid wall. The end goal is to simulate the motion of gas flows in the gap between the two vanes, subject to different initial pressures *p*_{∞}, hot (*T*_{H}) and cold (*T*_{C}) vane temperatures as listed in Table VII, in order to identify the correct circulation, induced low velocity, temperature gradient, and Knudsen forces from the vanes. The results are to be obtained from both stochastic (DSMC) and deterministic (DGFS) simulations.

. | Cases . | ||
---|---|---|---|

Parameter . | M-01 . | M-02 . | M-03 . |

Pressure: p (Torr) | 1.163 | 2.903 | 7.246 |

Number density: n (×10^{21} m^{−3}) | 37.8609 | 94.5058 | 235.8901 |

Knudsen number:^{a} Kn | 1.85 | 0.74 | 0.30 |

Cold vane temperature: T_{C} (K) | 306 | 306 | 304 |

Hot vane temperature: T_{H} (K) | 363 | 356 | 331 |

DGFS parameters | |||

Points in velocity mesh: N^{3} | 24^{3} | 24^{3} | 24^{3} |

Points in radial direction:^{b} N_{ρ} | 6 | 6 | 6 |

Points on half sphere:^{b} M | 6 | 6 | 6 |

Size of velocity mesh^{c} | [−5, 5]^{3} | [−5, 5]^{3} | [−5, 5]^{3} |

BGK/ESBGK/S-model parameters | |||

Points in velocity mesh: N^{3} | 48^{3} | 24^{3} | 24^{3} |

Size of velocity mesh^{c} | [−7, 7]^{3} | [−5, 5]^{3} | [−5, 5]^{3} |

. | Cases . | ||
---|---|---|---|

Parameter . | M-01 . | M-02 . | M-03 . |

Pressure: p (Torr) | 1.163 | 2.903 | 7.246 |

Number density: n (×10^{21} m^{−3}) | 37.8609 | 94.5058 | 235.8901 |

Knudsen number:^{a} Kn | 1.85 | 0.74 | 0.30 |

Cold vane temperature: T_{C} (K) | 306 | 306 | 304 |

Hot vane temperature: T_{H} (K) | 363 | 356 | 331 |

DGFS parameters | |||

Points in velocity mesh: N^{3} | 24^{3} | 24^{3} | 24^{3} |

Points in radial direction:^{b} N_{ρ} | 6 | 6 | 6 |

Points on half sphere:^{b} M | 6 | 6 | 6 |

Size of velocity mesh^{c} | [−5, 5]^{3} | [−5, 5]^{3} | [−5, 5]^{3} |

BGK/ESBGK/S-model parameters | |||

Points in velocity mesh: N^{3} | 48^{3} | 24^{3} | 24^{3} |

Size of velocity mesh^{c} | [−7, 7]^{3} | [−5, 5]^{3} | [−5, 5]^{3} |

^{a}

Based on hard-sphere definition (see Ref. 11).

### B. Numerical details

The simulation is carried out at wide range of Knudsen number for flows in early slip to early free molecular regime. The simulation specific numerical parameters as well as differences between stochastic (DSMC) and deterministic (DGFS) modelling is described next.

**DSMC**: SPARTA^{77}has been employed for carrying out DSMC simulations in the present work. The simulations are first run for 200 000 unsteady steps wherein the particles move, collide, and allowed to equilibrate. No sampling is performed at this stage. Next, the simulation is run for another 5 000 000 steady steps wherein the samples of flow properties namely number density, flow velocity, temperature, stress, and heat-flux, are taken for sufficiently long time so as to produce a meaningful bulk properties as well as minimize the statistical noise therein. In the present case, the DSMC domain is discretized into 300 × 150 cells, resulting in a uniform cell size of 2*μ*m, with 50 particles per cell on average during initialization. A time step of 10^{−9}s is used during*move*step of DSMC algorithm throughout the course of simulation. Note that these DSMC parameters have been taken from Ref. 3 wherein the authors performed multiple verification cases with different time-steps, grid-size, domain length, particles per cell, etc. N_{2}is used as the working gas in simulations, since MIKRA experiments^{3}were performed in N_{2}medium. The properties of the working gas is given in Table III. DSMC simulations treat N_{2}as diatomic species, and takes rotational degrees of freedom into account.**DGFS/BGK/ESBGK/S-model**: We use the DGFS implementation described in Ref. 38. The spatial domain consists of 849 elements [39 × 23 (total) - 2 × 6 × 4 (remove the vane regions)]. We use a linearly refined structured grid as illustrated in Fig. 12. While structured grids might seem inflexible compared to unstructured grids, they are known to produce more stable scheme with superior convergence rates,^{86,87}are amenable to highly efficient adaptive*h*/*p*mesh refinement via recursive element splitting^{88}(Nevertheless, DGFS is more general, and test cases on general grids will be reported in future works). Since we are seeking a steady state solution, the time-step is selected based on the CFL constraints of the forward Euler scheme. Other case specific DGFS parameters have been provided in Table VII.

It is worth noting that both the methods have different cell size requirements. In DSMC method, the contribution of particle collision to the transport properties is affected by strict spatial cell size requirements. In DGFS, however, the transport properties are strongly affected by local 3-D velocity space resolution rather than spatial resolution. As we show later, one can resolve the flow properties with fewer cells using DGFS.

### C. Results and discussion

#### 1. Flow pattern

Figures 13 and 14 illustrate the contour plot of various flow properties for the highest pressure case Kn = 0.3 (left column) and Kn = 1.85 (right column). For each of these plots, the DSMC and DGFS contours have been overlaid, wherein DSMC results have been indicated by thin black lines, and DGFS results have been indicated with thick red lines. Since the flow is strictly driven by temperature gradients, we expect very small deviation in the number density from the equilibrium value of 235.8901 × 10^{21} m^{−3}, as is also evident from Fig. 13(a). In terms of temperature, in Fig. 13(c), we observe a rather familiar flow expansion, in the sense that, the hot vane dissipates heat to the surrounding acting as a source, thereby giving rise to a spiral with spiral’s origin at the hot vane. Observe the interaction of contour lines (isotherms at 305 K and 310 K) with the cold vane in the region (250 *μ*m ≤ *X* ≤ 300 *μ*m, 25 *μ*m ≤ *Y* ≤ 60 *μ*m). We notice sharply curved isotherms near the top and right sides of the cold vane (see Fig. 15). Taking into account the Knudsen number of 0.3 and the characteristic length scale of system of 20 *μ*m, the Knudsen layer should extend few mean free paths from the solid surfaces i.e., *O*(*λ*) ≈ *O*(6 *μ*m). Therefore, one should expect some temperature jump, and therefore nonlinearity in the temperature in the near-wall region. More interestingly, we note an inflection in the isotherms at the top surface of the cold vane. This is essentially because the cold vane surface temperature is 304 K, while the free-stream is at 296 K. Hence, near to the heating source, say top-right end of the cold vane, the surface temperature is lower than the temperature of a layer of molecules just above the surface; and far away from the heating source, say top-left end of the cold vane, the surface temperature is higher than the temperature of a layer of molecules just above the surface. Therefore, an inflection in isotherms is expected somewhere between the top-left and top-right corner of the cold vane.

The origin of Knudsen force can be appreciated as follows. Consider a differential area *dS* over the cold vane as shown Fig. 15. The molecules impinging on the area *dS* can be thought as made up of two types of molecules: molecules coming from colder point *A* and molecules coming from hotter point *B*, both separated by few mean free paths. Near to the top right end of the cold vane, nearer to the hot vane, one should expect larger concentration of molecules of type *B*, and smaller concentration of molecules of type *A*. Conversely, near to the top left end of the cold vane, which is (relatively) far away from the hot vane, one should expect a smaller concentration of molecules of type *B*, and larger concentration of molecules of type *A*. Specifically, at the top left end of the cold vane, due to this imbalance of particles hitting the surface area, the momentum transferred to the surface element *dS* is in the opposite direction to the temperature gradient; however the gas flow is induced in the direction of the temperature gradient.^{14,23} This overall momentum imbalance contributes to the Knudsen force.

Figure 13(e) illustrates the variation of off-diagonal (*xy*) component of stress tensor at Kn = 0.3. First, we note the development of four ovals/ellipses originating at the four corners/edges of the hot vane. The effect is more pronounced at the right end (top-right and bottom-right corners) of the hot vane i.e., the length of the semi-major axis is larger for the ellipses on the right. At the top-left corner of the hot vane, in particular, we observe interaction of ovals with the top-right edge of the cold vane (note the distorted shape of the oval/ellipse at the top-left boundaries of the hot vane). Since the Knudsen number is in the *slip/early-transition* regime (Kn = 0.3), consider the expression for the stress-tensor, arrived in part by second order Chapman-Enskog expansion^{2}

where *P*_{ij}, *p*, *u*, *μ*, *ρ* are stress tensor, pressure, velocity, dynamic viscosity, and density respectively. *δ*_{ij} is the Kronecker delta function, *τ*_{ij} is the off-diagonal term of the stress tensor, and *K*_{i} ≈ 1, *i* = {1, 2, 3} are species/molecular-interaction specific constants.^{2} This yields

Let us consider four points in the flow: *A* (top-right corner of cold vane), *B* (top-left corner of hot vane), *C* (third vertex of equilateral triangle Δ*ABC* s.t. $BC\u20d7\xd7CA\u20d7/\Vert BC\u20d7\xd7CA\u20d7\Vert =k^$), and *D* (mid point of *A* and *B*) as shown in Fig. 16. Based on isotherms in Fig. 15, it can be inferred that the temperature difference between points *A* and *B* is ≈10 K, whereas the temperature difference between points *C* and *D* is ≈5 K. Consistent with Eq. (28), theoretically, we expect the thermal stresses (and therefore *P*_{xy}) to be larger between points *A* and *B* since *∂T*/*∂x*|_{AB} ≫ *∂T*/*∂y*|_{CD} (more formally: ∥∇*T*∥_{AB} ≫ ∥∇*T*∥_{CD}, $\Vert \u22072T\Vert AB\u226b\Vert \u22072T\Vert CD$). Hence, the distorted ellipse. A more subtle observation is as follows: Why, *precisely*, should an isocontour line of *xy* component of stress, start from top-left corner of the hot-vane (i.e., point *B*) and end at the top-right corner of cold-vane (i.e., point *A*). What happens to the entire flow field if we introduce roughness on the walls, or smooth the vane corners–few questions that we delegate to a future study.

Next, Figs. 14(a) and 14(c) depict the variation of *x* and *y* components of heat flux. We want to reemphasize that DSMC simulations consider the rotational degrees of freedom of N_{2} into account, whereas DGFS, being in very early stages of research, does not. Nevertheless, we observe a fair agreement between DSMC and DGFS. In Fig. 14(a), in the region (250 *μ*m ≤ *X* ≤ 320 *μ*m, 50 *μ*m ≤ *Y* ≤ 90 *μ*m), we again note presence of isocontour lines between the top-left and top-right corners of the cold and hot vanes. A more subtle observation is as follows: Multiple isocontours, for instance −2000 W/m,^{2} −4000 W/m,^{2} −6000 W/m^{2} (the unlabeled contour just below −4000 W/m^{2} isocontour), differing by large magnitudes, start at *approximately* the top-left corner of the hot vane, and end at the top-right corner of the cold vane, resulting in sharply curved isocontours. A partial explanation of such effects appears in Ref. 20, wherein the author attributed the observation to simply *edge effects*, basing the argument on the imbalance of particles of type A (cold) and type B (hot) near to the edges, as was mentioned earlier in the discussion.

Figure 14(e) illustrates the flow speed in the domain. We notice significant statistical fluctuations in DSMC (thin black lines), to an extent that removing DGFS contour lines in red, would make it difficult, if not impossible, to decipher the overall flow structure. A more complete picture of the flow is presented through transient DGFS streamlines in Fig. 17. First, we note the streamlines pointing in the upward direction. This is essentially due to the heating of the molecules (and therefore the thermal energy imparted to them) in the lower portions of the domain. In the process, four characteristic vortexes appear at the four corners of the heated vane, relatively early during the course of the simulation, for instance, see Fig. 17(a) at 1.25 ms. Over the time, secondary vortexes appear in the flow, most notably, a larger vortex at the top of the cold vane, and a smaller vortex near the top-right corner of hot vane.

Figure 18 shows the *steady state* speed contours at different Knudsen numbers with the corresponding flow streamlines overlaid. With increase in Knudsen number from Kn = 0.3 to Kn = 0.74, we note sharp increase in flow velocity, approximately by a factor of two. Consequently, the vortexes grow in size. The change in flow speed, however, from Kn = 0.74 to Kn = 1.85, although appreciable, is relatively mild.

Finally, we compare the variation of flow properties along the vertical centerline (*x* = 300 *μ*m, 0 ≤ *y* ≤ 300 *μ*m) in Figs. 19–22 for various models. We observe a fair agreement between DSMC and DGFS results ignoring the statistical noise (see Figs. 13 and 14). In particular, in Fig. 19(b), we observe peak temperatures near the *edges* of hot and cold vanes i.e., in the region *x* = 300 *μ*m, 30 ≤ *y* ≤ 60 *μ*m. Through Figs. 19(c) and 19(d), we infer that the thermal gradients are stronger in the *x*-direction. More notably, we observe the highest thermal-stress in the *edge* region (note the valley in the region *x* = 300 *μ*m, 40 ≤ *y* ≤ 60 *μ*m). We conjecture the trough of the valley to be shallower if the vane edges ought to be made smoother. A slightly peculiar observation is as follows: the trough of the valley is deeper at Kn = 0.74 compared to Kn = 0.30, and shallower at Kn = 1.85 compared to Kn = 0.74. This could be explained as follows: at Kn = 0.30 the temperature difference, *T*_{H} − *T*_{C}, is lower than the one correspoding to the Kn = 0.74 case and therefore the thermal stress increases in the latter case. For the Kn = 1.85 and Kn = 0.74 cases, wherein the temperature difference is approximately same, the peak thermal-stress decreases owing to the bimodal nature of the Knudsen forces.

## V. MULTISPECIES MIKRA

In the present section, we carry out the MIKRA simulations for binary mixture consisting of N_{2} and H_{2}O using the variable soft sphere model.

### A. Problem statement

The flow configuration remains the same as shown in Fig. 11. We consider the 2D uniform flow of binary mixture of N_{2} and H_{2}O. The end goal is to simulate the motion of gas flows in the gap between the two vanes, subject to initial pressure *p*_{∞}, hot (*T*_{H}) and cold (*T*_{C}) vane temperature as listed in Table VIII, in order to identify the correct circulation, induced low velocity, temperature gradient, Knudsen forces, and heat transfer rate from the vanes. The results are to be obtained from both stochastic (DSMC) and deterministic (DGFS) simulations.

. | Cases . |
---|---|

Parameter . | MSM-01 . |

Pressure: p (Torr) | 1.163 |

Total number density: n (×10^{21} m^{−3}) | 37.860 91 |

Concentration: ($n(N2)/n,\u2009n(H2O)/n$) | (0.05, 0.95) |

Knudsen number:^{a} Kn | 1.85 |

Cold vane temperature: T_{C} (K) | 306 |

Hot vane temperature: T_{H} (K) | 363 |

DGFS parameters | |

Points in velocity mesh: N^{3} | 32^{3} |

Points in radial direction:^{b} N_{ρ} | 8 |

Points on full sphere:^{b} M | 12 |

Size of velocity mesh^{c} | [−6, 6]^{3} |

. | Cases . |
---|---|

Parameter . | MSM-01 . |

Pressure: p (Torr) | 1.163 |

Total number density: n (×10^{21} m^{−3}) | 37.860 91 |

Concentration: ($n(N2)/n,\u2009n(H2O)/n$) | (0.05, 0.95) |

Knudsen number:^{a} Kn | 1.85 |

Cold vane temperature: T_{C} (K) | 306 |

Hot vane temperature: T_{H} (K) | 363 |

DGFS parameters | |

Points in velocity mesh: N^{3} | 32^{3} |

Points in radial direction:^{b} N_{ρ} | 8 |

Points on full sphere:^{b} M | 12 |

Size of velocity mesh^{c} | [−6, 6]^{3} |

### B. Numerical details

The multispecies simulations are carried out for flows in transition regime. The specific differences between stochastic (DSMC) and deterministic (DGFS) modelling is described next.

**DSMC**: SPARTA^{77}has been employed for carrying out DSMC simulations in the present work. The geometric parameters remain the same as described in Sec. IV B. A minimum of 300 DSMC simulator particles per cell is used in conjunction with the no-time collision (NTC) algorithm and VSS scattering model. The simulations are first run for 200 000 unsteady steps, and subsequently another 5 000 000 steady steps wherein the flow sampling is performed. Similar to the previous single-species MIKRA case, the DSMC domain is discretized into 300 × 150 cells, resulting in a uniform cell size of 2*μ*m, with 285 particles of H_{2}O and 15 particles of N_{2}, per cell on average during initialization. A time step of 10^{−9}s is used during*move*step of DSMC algorithm throughout the course of simulation. N_{2}and H_{2}O are used as the working gases in simulations. The properties of the working gas is given in Table IX. Note that for N_{2}, we consider*ζ*_{R}= 2 rotational degrees of freedom, rotational relaxation*Z*_{R}= 0.2,*ζ*_{V}= 2 vibrational degrees of freedom, vibrational relaxation*Z*_{V}= 1.901 14 × 10^{−5}, and vibrational temperature $Tv$ = 3371 K; and for H_{2}O, we consider*ζ*_{R}= 3 rotational degrees of freedom, rotational relaxation*Z*_{R}= 0.2,*ζ*_{V}= 3 vibrational degrees of freedom, vibrational relaxation*Z*_{V}= 1.901 14 × 10^{−5}, and vibrational temperature $Tv$ = 5261 K.**DGFS**: We use the DGFS implementation described in Ref. 12. The geometrical parameters remain the same as described in Sec. IV B. Multispecies case specific DGFS parameters have been provided in Table VIII. Note that, we employ N_{2}and H_{2}O as the working gas in simulations. N_{2}is diatomic, and H_{2}O is triatomic, however, DGFS, as of now, is applicable for monoatomic gases only. Since the working temperature range is low, we anticipate the effects of vibrational degrees of freedom to be negligible.

. | N_{2}
. | H_{2}O
. |
---|---|---|

Mass: m (kg) | 46.5 × 10^{−27} | 29.9 × 10^{−27} |

Viscosity index:^{a} ω_{i}, (−) | 0.74 | 1.00 |

Scattering index: α_{i}, (−) | 1.36 | 1.00 |

Ref. diameter: d_{ref,i} (m) | 4.07 × 10^{−10} | 5.78 × 10^{−10} |

Ref. temperature: T_{ref,i} (K) | 273 | 273 |

. | N_{2}
. | H_{2}O
. |
---|---|---|

Mass: m (kg) | 46.5 × 10^{−27} | 29.9 × 10^{−27} |

Viscosity index:^{a} ω_{i}, (−) | 0.74 | 1.00 |

Scattering index: α_{i}, (−) | 1.36 | 1.00 |

Ref. diameter: d_{ref,i} (m) | 4.07 × 10^{−10} | 5.78 × 10^{−10} |

Ref. temperature: T_{ref,i} (K) | 273 | 273 |

### C. Results and discussion

#### 1. Flow pattern

Figures 23 and 24 illustrate the contour plot of various flow properties for the MSM-01 case in transition regime, wherein the N_{2} and H_{2}O are in 0.05:0.95 concentration ratio. Similar to the single species case, for each of these plots, the DSMC and DGFS contours have been overlaid, wherein DSMC results have been indicated by thin black lines, and DGFS results have been indicated with thick red lines. Since the flow is strictly driven by temperature gradients, we expect very small deviation in the number density from the equilibrium values of 35.967 864 5 × 10^{21} m^{−3} for H_{2}O and 1.893 045 5 × 10^{21} m^{−3} for N_{2}, as is also evident from Figs. 23(a) and 23(b). In terms of temperature, in Figs. 23(c) and 23(d), we again observe a rather familiar flow expansion, in the sense that, the hot vane dissipates heat to the surrounding acting as a source, thereby giving rise to a spiral with spiral’s origin at the hot vane. From the fundamental mass/momentum conservation principles, one can infer that, in the presence of temperature gradients, the heavier species, here N_{2}, moves slower and the lighter species, here H_{2}O, moves faster giving rise to the well-known thermal diffusion. This explain why the isotherms for H_{2}O spread farther apart compared to those of N_{2}.

Figures 24(c) and 24(d) illustrate the variation of off-diagonal (*xy*) component of stress tensor. Again, we observe the development of four ovals/ellipses originating at the four corners/edges of the hot vane, wherein the effects are more pronounced at the right end (top-right and bottom-right corners) of the hot vane. The stress is higher for H_{2}O compared to N_{2}. Figures 24(a) and 24(b) illustrate the flow speed in the domain. We notice significant statistical fluctuations in DSMC (thin black lines) contour lines for N_{2} due to lower number of DSMC simulator particles. In particular, we observe that DGFS results/contours are insusceptible to the concentration of the individual species, thereby opening the possibility of its application for simulating flows involving species in trace concentrations.

Finally, we compare the variation of individual species flow properties along the vertical centerline (*x* = 300 *μ*m, 0 ≤ *y* ≤ 300 *μ*m). We observe a fair agreement between DSMC and DGFS results ignoring the statistical noise for the bulk properties. In Fig. 19(b), we note peak temperatures near the *edges* of hot and cold vanes i.e., in the region *x* = 300 *μ*m, 30 ≤ *y* ≤ 60 *μ*m. More specifically, the temperature is higher for H_{2}O compared to N_{2}, an observation consistent with fundamental conservation principles. One can infer that the magnitude of the thermal gradients are stronger in the *x*-direction. Notably, in Fig. 25(c), we observe the highest thermal-stress in the *edge* region (note the valley in the region *x* = 300 *μ*m, 40 ≤ *y* ≤ 60 *μ*m). Finally, consistent with aforementioned observations, we observe higher velocity at *x* = 300 *μ*m, *y* ≈ 54 *μ*m—the location of the top edges of the two vanes.

## VI. CONCLUSIONS

We have presented an application of the recently introduced deterministic discontinuous Galerkin fast spectral (DGFS) method for assessing the flow phenomenon in the thermostress convection enabled microscale device MIKRA—a compact low-power pressure sensor. We carried out MIKRA simulations in the slip-to-transition regime gas flows at different Knudsen numbers. The single-species cases are run with a variable hard sphere scattering model. We conclude that the results obtained with DGFS and DSMC are inextricable ignoring the statistical noise. The DSMC provides verification benchmark for the solution of the Boltzmann equation with real gas effects. The overall DGFS method is simple from a mathematical and implementation perspective, highly accurate in both physical and velocity spaces as well as time, robust, i.e., applicable for general geometry and spatial mesh, exhibits nearly linear parallel scaling, and directly applies to general collision kernels, for instance, Bird’s variable hard/soft sphere models, needed for high fidelity modelling. DGFS presents a viable alternative for simulation of highly information rich thermostress convection processes at the microscale.

## ACKNOWLEDGMENTS

S.J. and J.H.’s research is partially supported by NSF Grant No. DMS-1620250 and NSF CAREER Grant No. DMS-1654152. Support from NSF CDS&E Grant No. CBET-1854829 is also gratefully acknowledged.

## REFERENCES

_{2}H

_{2}O gas mixtures

_{2}-O

_{2}deflagration waves