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.

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/|∇xT| to be comparable to the molecular mean free path λ. At the macroscale, such magnitudes are prohibitive, necessitating thermal gradients on the order of 106 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 kT∇ ln T, where kT 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, Bird11 devised a self-diffusion test case (see Sec. 12.6) where the diffusion coefficient was measured by ignoring the thermal gradient term kT∇ 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 experimentally10 observed diffusion coefficient. This suggests that kT 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 106 K/m, as noted earlier, kT∇ 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, Sone14 identified three broad groups of the temperature driven flow based on its application in microsystems: (a) thermal creep flow15–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 devices27 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 gases31 due to a square-root viscosity variation with temperature. The variable hard sphere (VHS) model proposed by Bird11 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 flows32,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 fast2 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.

FIG. 1.

A holistic view of growth in DSMC publications over years; and the exponential increase in number of simulated collision events per hour.

FIG. 1.

A holistic view of growth in DSMC publications over years; and the exponential increase in number of simulated collision events per hour.

Close modal

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) method12,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 model12 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 Fierro43 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/m2 pressure range for helium, krypton, hydrogen, oxygen, and carbon dioxide. Alexeenko44 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. Zhu45 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 Nabeth46 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. Anikin47 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, Lotfian49 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 Strongrich51 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 kT∇ 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) switches70 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.

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ΩR3 is the position, and vR3 is the particle velocity [f(i)dxdv gives the number of particles of species i to be found in an infinitesimal volume dxdv centered at the point (x, v) of the phase space]. The time evolution of f(i) is described by the multispecies Boltzmann equation written as72,73

tf(i)+vxf(i)=j=1sQ(ij)(f(i),f(j)),i=1,,s.
(1)

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

Q(ij)(f(i),f(j))(v)=R3S2Bij(vv*,σ)f(i)(v)f(j)(v*)f(i)(v)f(j)(v*)dσdv*,
(2)

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

miv+mjv*=miv+mjv*,mi|v|2+mj|v*|2=mi|v|2+mj|v*|2,
(3)

where mi, mj denote the mass of particles of species i and j, respectively. Hence, one can parameterize v′ and v* as follows:

v=v+v*2+(mimj)2(mi+mj)(vv*)+mj(mi+mj)|vv*|σ,v*=v+v*2+(mimj)2(mi+mj)(vv*)mi(mi+mj)|vv*|σ,
(4)

with σ being a vector varying on the unit sphere S2. Bij=Bji(0) is the collision kernel characterizing the interaction mechanism between particles. It can be shown that

Bij=Bij(|vv*|,cosχ),cosχ=σ(vv*)|vv*|,
(5)

where χ is the deviation angle between vv* and v′ − v*.

Given the interaction potential between particles, the specific form of Bij can be determined using the classical scattering theory

Bij(|vv*|,cosχ)=|vv*|Σij(|vv*|,χ),
(6)

where Σij is the differential cross section given by

Σij(|vv*|,χ)=bijsinχdbijdχ,
(7)

with bij being the impact parameter.

With a few exceptions, the explicit form of Σij can be hard to obtain since bij 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 Matsumoto33 introduced the so-called VSS model by assuming

χ=2cos1{(bij/dij)1/αij},
(8)

where αij is the scattering parameter and dij is the diameter borrowed from the VHS model [Eq. (4.79) in Ref. 11]

dij=dref,ij2kBTref,ijμij|vv*|2ωij0.51Γ(2.5ωij)1/2,
(9)

with Γ being the Gamma function, μij=mimjmi+mj the reduced mass, dref,ij, Tref,ij, and ωij, respectively, the reference diameter, reference temperature, and viscosity index. Substituting Eqs. (7)–(9) into (6), one can obtain Bij as

Bij=bωij,αij|vv*|2(1ωij)(1+cosχ)αij1,
(10)

where bωij,αij is a constant given by

bωij,αij=dref,ij242kBTref,ijμijωij0.5αijΓ(2.5ωij)2αij1.
(11)

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

n(i)=R3f(i)dv,ρ(i)=min(i),u(i)=1n(i)R3vf(i)dv,T(i)=mi3n(i)kBR3(vu(i))2f(i)dv.
(12)

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

n=i=1sn(i),ρ=i=1sρ(i),u=1ρi=1sρ(i)u(i).
(13)

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

vD(i)=1n(i)R3cf(i)dv=u(i)u,P(i)=R3miccf(i)dv,q(i)=R312mic|c|2f(i)dv,

where c = vu is the peculiar velocity. Finally, the total stress, heat flux, pressure, and temperature are given by

P=i=1sP(i),q=i=1sq(i),p=nkBT=13tr(P).
(14)

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 uini, temperature Tini, 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.

FIG. 2.

Flowchart of the DSMC process.

FIG. 2.

Flowchart of the DSMC process.

Close modal

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

P=|vv*|Σij(|vv*|Σij)max,
(15)

where (|vv*|Σ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.,

ρcell(i)=1Vcellj=1Mcell(i)mj,cell(i).
(16)

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

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 N3 velocity grid, the fast Fourier spectral method produces Q(ij)(f(i),f(j)) at the same grid with O(MNρN3 log N) complexity, where MN2 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.

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

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

R3Sdv=R3vSdv=R3v2Sdv=0.
(17)
  1. The entropy production is always positive, i.e.,

R3ln(f)Sdv0.
(18)
  1. Due to specific form of S, the phase density in equilibrium is a Maxwellian, i.e.,

S=0R3ln(f)Sdv=0f=M,
(19)

where

M=n(2πRT)3/2exp(vu)22RT.
(20)
  1. The Prandtl number is close to 2/3 for monoatomic gases, i.e.,

Pr=52kBmμκ,
(21)

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

S=ν(fγf),
(22)

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

fγBGK=M,
(23)

whereas for ESBGK, fγ is anisotropic Gaussian given as

fγESBGK=ndet(2πT)exp12(vu)TT1(vu),ρT=1PrρRTId+11Prρ,ρ=R3ccfdv,ρRT=R3ccfγBGKdv,
(24)

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

fγS-model=fγBGK1+1Pr5ScpRTc22RT52,S=R3cc2fdv.
(25)

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.

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.

FIG. 3.

Numerical setup for 1D Fourier-Couette flow.

FIG. 3.

Numerical setup for 1D Fourier-Couette flow.

Close modal
TABLE I.

Common numerical parameters for Fourier-Couette flow.

Common parameters
Molecular mass: m1 (×1027 kg) 66.3 
Nondimentional physical space [0, 1] 
Spatial elements 
DG order 
Time stepping Euler 
Viscosity index: ω 0.81 
Scattering parameter: α 1.4 
Ref. diameter: dref (×1010 m) 4.17 
Ref. temperature: Tref (K) 273 
Ref. viscosity: μref (Pa s) 2.117 × 10−5 
Characteristic mass: m0 (×1027 kg) 66.3 
Characteristic length: H0 (mm) 
Characteristic velocity: u0 (m/s) 337.2 
Characteristic temperature: T0 (K) 273 
Characteristic no. density: n0 (m−33.538 × 1022 
Initial conditions  
Velocity: u (m/s) 
Temperature: T (K) 273 
Number density: n (m−33.538 × 1022 
Knudsen number:a (Kn) 0.036 
Common parameters
Molecular mass: m1 (×1027 kg) 66.3 
Nondimentional physical space [0, 1] 
Spatial elements 
DG order 
Time stepping Euler 
Viscosity index: ω 0.81 
Scattering parameter: α 1.4 
Ref. diameter: dref (×1010 m) 4.17 
Ref. temperature: Tref (K) 273 
Ref. viscosity: μref (Pa s) 2.117 × 10−5 
Characteristic mass: m0 (×1027 kg) 66.3 
Characteristic length: H0 (mm) 
Characteristic velocity: u0 (m/s) 337.2 
Characteristic temperature: T0 (K) 273 
Characteristic no. density: n0 (m−33.538 × 1022 
Initial conditions  
Velocity: u (m/s) 
Temperature: T (K) 273 
Number density: n (m−33.538 × 1022 
Knudsen number:a (Kn) 0.036 
a

Based on variable hard-sphere definition (see Refs. 11 and 38).

TABLE II.

Numerical parameters for Fourier-Couette cases.

ParameterCase FC-01Case FC-02Case FC-03Case FC-04Case FC-05
Nondimentional velocity spacea [−5, 5]3 [−5, 5]3 [−5, 5]3 [−5, 5]3 [−8, 8]3 
{N3, Nρ, M}b {243, 6, 6} {243, 6, 6} {243, 6, 6} {243, 6, 6} {323, 16, 6} 
Left wall (purely diffuse) conditions 
Velocity: ul (m/s) (0, −50, 0) (0, −50, 0) (0, −250, 0) (0, −250, 0) (0, −250, 0) 
Temperature: Tl (K) 273 223 273 223 173 
Right wall (purely diffuse) conditions 
Velocity: ur (m/s) (0, 50, 0) (0, 50, 0) (0, 250, 0) (0, 250, 0) (0, 250, 0) 
Temperature: Tr (K) 273 323 273 323 373 
ParameterCase FC-01Case FC-02Case FC-03Case FC-04Case FC-05
Nondimentional velocity spacea [−5, 5]3 [−5, 5]3 [−5, 5]3 [−5, 5]3 [−8, 8]3 
{N3, Nρ, M}b {243, 6, 6} {243, 6, 6} {243, 6, 6} {243, 6, 6} {323, 16, 6} 
Left wall (purely diffuse) conditions 
Velocity: ul (m/s) (0, −50, 0) (0, −50, 0) (0, −250, 0) (0, −250, 0) (0, −250, 0) 
Temperature: Tl (K) 273 223 273 223 173 
Right wall (purely diffuse) conditions 
Velocity: ur (m/s) (0, 50, 0) (0, 50, 0) (0, 250, 0) (0, 250, 0) (0, 250, 0) 
Temperature: Tr (K) 273 323 273 323 373 
a

Nondimensional (see Refs. 12 and 38 for details on nondimensionalization).

b

Required only in the fast Fourier spectral low-rank decomposition for DGFS method (see Refs. 38 and 74).

FIG. 4.

Variation of flow properties along the domain for Fourier-Couette flow cases obtained with BGK, DGFS, and DSMC for Argon: (a) normalized temperature, cases FC-0{1, 2, 3, 4} and (b) normalized y-velocity, and temperature for case FC-05.

FIG. 4.

Variation of flow properties along the domain for Fourier-Couette flow cases obtained with BGK, DGFS, and DSMC for Argon: (a) normalized temperature, cases FC-0{1, 2, 3, 4} and (b) normalized y-velocity, and temperature for case FC-05.

Close modal

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.

FIG. 5.

Numerical setup for the flow around a micro-electronic chip: (a) Schematic, (b) mesh for DGFS simulations. For DSMC simulations, we subdivide each cell of the mesh above into 5 × 5 subcells.

FIG. 5.

Numerical setup for the flow around a micro-electronic chip: (a) Schematic, (b) mesh for DGFS simulations. For DSMC simulations, we subdivide each cell of the mesh above into 5 × 5 subcells.

Close modal

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: SPARTA77 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 benchmarked77 and widely used for studying hypersonic, subsonic and thermal12,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. N2 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., N2 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 N2 as the working gas in simulations, since MIKRA experiments3 were performed in N2 medium. N2 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.

TABLE III.

N2 gas VHS parameters used in 2-D single-species DSMC and DGFS simulations. Note that DGFS, being in very early stage of research, treats N2 as a monoatomic species.

Mass: m (kg) 46.5 × 10−27 
Viscosity index: ω (−) 0.74 
Scattering index: α (−) 1.0 
Ref. diameter: dref (m) 4.17 × 10−10 
Ref. temperature: Tref (K) 273 
Ref. viscosity: μref (Pa s) 1.656 × 10−5 
DSMC specific parameters  
Rotational degrees of freedom: ζR (−) 
Rotational relaxation: ZR (−) 
Vibrational degrees of freedom: ζV (−) 
Vibrational relaxation ZV (−) 1.901 14 × 10−5 
Vibrational temperature TV (K) 3371 
Mass: m (kg) 46.5 × 10−27 
Viscosity index: ω (−) 0.74 
Scattering index: α (−) 1.0 
Ref. diameter: dref (m) 4.17 × 10−10 
Ref. temperature: Tref (K) 273 
Ref. viscosity: μref (Pa s) 1.656 × 10−5 
DSMC specific parameters  
Rotational degrees of freedom: ζR (−) 
Rotational relaxation: ZR (−) 
Vibrational degrees of freedom: ζV (−) 
Vibrational relaxation ZV (−) 1.901 14 × 10−5 
Vibrational temperature TV (K) 3371 
TABLE IV.

Numerical parameters for flow around microelectronic chip.

ParametersMEC-01
Spatial elements 190 quadrilaterals 
DG order 
Time stepping Euler 
Points in velocity mesh: N3 243 
Points in radial direction:aNρ 
Points on half sphere:aM 
Size of velocity meshb [−5, 5]3 
Characteristic length: H0 (μm) 
Characteristic velocity: u0 (m/s) 402.54 
Characteristic temperature: T0 (K) 273 
Characteristic no. density: n0 (m−34.894 × 1023 
Initial conditions  
Velocity: u (m/s) 
Temperature: T (K) 273 
Number density: n (m−34.894 × 1023 
Knudsen number:c (Kn) 0.881 58 
ParametersMEC-01
Spatial elements 190 quadrilaterals 
DG order 
Time stepping Euler 
Points in velocity mesh: N3 243 
Points in radial direction:aNρ 
Points on half sphere:aM 
Size of velocity meshb [−5, 5]3 
Characteristic length: H0 (μm) 
Characteristic velocity: u0 (m/s) 402.54 
Characteristic temperature: T0 (K) 273 
Characteristic no. density: n0 (m−34.894 × 1023 
Initial conditions  
Velocity: u (m/s) 
Temperature: T (K) 273 
Number density: n (m−34.894 × 1023 
Knudsen number:c (Kn) 0.881 58 
a

Nondimensional (see Refs. 12 and 38 for details on nondimensionalization).

b

Required only in the fast Fourier spectral low-rank decomposition for DGFS method (see Refs. 38 and 74).

c

Based on variable hard-sphere definition (see Refs. 11 and 38).

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.

FIG. 6.

Flow properties at steady state for micro-electronic chip obtained from DSMC and DGFS using VHS collision model. For each of these figures, DSMC results (mirrored along y-axis) have been shown in the second quadrant (−17 μm ≤ x < 0 μm), whereas DGFS results have been illustrated in the first quadrant (0 μm ≤ x < 17 μm). Observe the legend for number-density: (a) Number density (m−3); (b) Temperature (K); (c) Speed (m/s); (d) xy-component of stress (N/m2); (e) x-component of heat-flux (W/m2); (f) y-component of heat-flux (W/m2).

FIG. 6.

Flow properties at steady state for micro-electronic chip obtained from DSMC and DGFS using VHS collision model. For each of these figures, DSMC results (mirrored along y-axis) have been shown in the second quadrant (−17 μm ≤ x < 0 μm), whereas DGFS results have been illustrated in the first quadrant (0 μm ≤ x < 17 μm). Observe the legend for number-density: (a) Number density (m−3); (b) Temperature (K); (c) Speed (m/s); (d) xy-component of stress (N/m2); (e) x-component of heat-flux (W/m2); (f) y-component of heat-flux (W/m2).

Close modal

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

F=dApndA,
(26)

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.

TABLE V.

x-component of force on the substrate for MEC-01 case, obtained using DSMC and DGFS simulations.

Force (μN/μm)
Pressure (Pa)KnDSMCDGFS
2000 0.881 58 −0.040 008 843 −0.040 010 413 
Force (μN/μm)
Pressure (Pa)KnDSMCDGFS
2000 0.881 58 −0.040 008 843 −0.040 010 413 

The present test case closely follows case-I(a) from Ref. 82. In the current test case, two reservoirs filled with N2 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 equations82–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.

FIG. 7.

Numerical setup for the flow in short microchannels. On the horizontal channel walls, we impose a linearly decreasing temperature profile similar to case I(a) in Ref. 82: (a) Schematic, (b) Mesh for DGFS simulations.

FIG. 7.

Numerical setup for the flow in short microchannels. On the horizontal channel walls, we impose a linearly decreasing temperature profile similar to case I(a) in Ref. 82: (a) Schematic, (b) Mesh for DGFS simulations.

Close modal
TABLE VI.

Numerical parameters for flow in short microchannels.

ParametersSM-01
Spatial elements 127 quadrilaterals 
DG order 
Time stepping Euler 
Points in velocity mesh: N3 323 
Points in radial direction:aNρ 
Points on half sphere:aM 
Size of velocity meshb [−5.72, 5.72]3 
Characteristic length: H0 (μm) 
Characteristic velocity: u0 (m/s) 421.98 
Characteristic temperature: T0 (K) 300 
Characteristic no. density: n0 (m−36.62 × 1024 
Initial conditions  
Velocity: u (m/s) 
Temperature: T (K) 300 
Number density: n (m−36.62 × 1024 
Knudsen number:c (Kn) 0.2 
Inlet condition  
Velocity: uin (m/s) 
Temperature: Tin (K) 600 
Number density: nin (m−33.31 × 1024 
Pressure: pin (N/m) 27 420 
Outlet condition  
Velocity: uout (m/s) 
Temperature: Tout (K) 300 
Number density: nout (m−36.62 × 1024 
Pressure: pout (N/m) 27 420 
ParametersSM-01
Spatial elements 127 quadrilaterals 
DG order 
Time stepping Euler 
Points in velocity mesh: N3 323 
Points in radial direction:aNρ 
Points on half sphere:aM 
Size of velocity meshb [−5.72, 5.72]3 
Characteristic length: H0 (μm) 
Characteristic velocity: u0 (m/s) 421.98 
Characteristic temperature: T0 (K) 300 
Characteristic no. density: n0 (m−36.62 × 1024 
Initial conditions  
Velocity: u (m/s) 
Temperature: T (K) 300 
Number density: n (m−36.62 × 1024 
Knudsen number:c (Kn) 0.2 
Inlet condition  
Velocity: uin (m/s) 
Temperature: Tin (K) 600 
Number density: nin (m−33.31 × 1024 
Pressure: pin (N/m) 27 420 
Outlet condition  
Velocity: uout (m/s) 
Temperature: Tout (K) 300 
Number density: nout (m−36.62 × 1024 
Pressure: pout (N/m) 27 420 
a

Nondimensional (see Refs. 12 and 38 for details on non-dimensionalization).

b

Required only in the fast Fourier spectral low-rank decomposition for DGFS method (see Refs. 38 and 74).

c

Based on variable hard-sphere definition (see Refs. 11 and 38).

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. N2 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., Qx can be attributed to the fact that DSMC simulations consider rotational degrees-of-freedom of N2 into account, whereas DGFS does not.

FIG. 8.

Flow properties for short microchannel test-case obtained from DSMC and DGFS using VHS collision model. For each of these figures, DSMC results have been shown in the first quadrant (0 μm ≤ y < 2 μm), whereas DGFS results (mirrored along x-axis) have been illustrated in the fourth quadrant (−2 μm ≤ y < 0 μm). Differences in Qx can be attributed to the fact that DSMC simulations consider rotational degrees-of-freedom of N2 into account, whereas DGFS does not: (a) Number density (m−3); (b) Temperature (K); (c) Speed (m/s); (d) xy-component of stress (N/m2); (e) x-component of heat-flux (W/m2); (f) y-component of heat-flux (W/m2).

FIG. 8.

Flow properties for short microchannel test-case obtained from DSMC and DGFS using VHS collision model. For each of these figures, DSMC results have been shown in the first quadrant (0 μm ≤ y < 2 μm), whereas DGFS results (mirrored along x-axis) have been illustrated in the fourth quadrant (−2 μm ≤ y < 0 μm). Differences in Qx can be attributed to the fact that DSMC simulations consider rotational degrees-of-freedom of N2 into account, whereas DGFS does not: (a) Number density (m−3); (b) Temperature (K); (c) Speed (m/s); (d) xy-component of stress (N/m2); (e) x-component of heat-flux (W/m2); (f) y-component of heat-flux (W/m2).

Close modal

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

FIG. 9.

Flow properties on the horizontal centerline (y = 0 μm) for short microchannel test-case obtained from DSMC and DGFS using VHS collision model: (a) Number density (m−3), and Temperature (K); (b) Speed (m/s), and the xy-component of stress (N/m2).

FIG. 9.

Flow properties on the horizontal centerline (y = 0 μm) for short microchannel test-case obtained from DSMC and DGFS using VHS collision model: (a) Number density (m−3), and Temperature (K); (b) Speed (m/s), and the xy-component of stress (N/m2).

Close modal

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.

FIG. 10.

The CAD model for Gen1 Micro In-Plane Knudsen Radiometric Actuator (MIKRA).3 

FIG. 10.

The CAD model for Gen1 Micro In-Plane Knudsen Radiometric Actuator (MIKRA).3 

Close modal

The flow configuration is shown in Fig. 11. Consider the 2D uniform flow of N2 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 TC. The right vane, indicated in red, is kept at a higher/hot temperature which we denote by TH. 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 (TH) and cold (TC) 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.

FIG. 11.

Schematic for numerical simulation of thermostress convection in MIKRA Gen1.3 The interior dashed thin black lines indicate the blocks used for structured mesh generation. Specifically for deterministic DGFS simulations, a linear gradient is applied within blocks such that the cells are finer in the near-vane region.

FIG. 11.

Schematic for numerical simulation of thermostress convection in MIKRA Gen1.3 The interior dashed thin black lines indicate the blocks used for structured mesh generation. Specifically for deterministic DGFS simulations, a linear gradient is applied within blocks such that the cells are finer in the near-vane region.

Close modal
TABLE VII.

Numerical parameters for thermostress convection in MIKRA Gen1 simulations for DSMC and DGFS using VHS collision model for N2 molecules.

Cases
ParameterM-01M-02M-03
Pressure: p (Torr) 1.163 2.903 7.246 
Number density: n (×1021 m−337.8609 94.5058 235.8901 
Knudsen number:a Kn 1.85 0.74 0.30 
Cold vane temperature: TC (K) 306 306 304 
Hot vane temperature: TH (K) 363 356 331 
DGFS parameters 
Points in velocity mesh: N3 243 243 243 
Points in radial direction:bNρ 
Points on half sphere:bM 
Size of velocity meshc [−5, 5]3 [−5, 5]3 [−5, 5]3 
BGK/ESBGK/S-model parameters 
Points in velocity mesh: N3 483 243 243 
Size of velocity meshc [−7, 7]3 [−5, 5]3 [−5, 5]3 
Cases
ParameterM-01M-02M-03
Pressure: p (Torr) 1.163 2.903 7.246 
Number density: n (×1021 m−337.8609 94.5058 235.8901 
Knudsen number:a Kn 1.85 0.74 0.30 
Cold vane temperature: TC (K) 306 306 304 
Hot vane temperature: TH (K) 363 356 331 
DGFS parameters 
Points in velocity mesh: N3 243 243 243 
Points in radial direction:bNρ 
Points on half sphere:bM 
Size of velocity meshc [−5, 5]3 [−5, 5]3 [−5, 5]3 
BGK/ESBGK/S-model parameters 
Points in velocity mesh: N3 483 243 243 
Size of velocity meshc [−7, 7]3 [−5, 5]3 [−5, 5]3 
a

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

b

Required only in the fast Fourier spectral low-rank decomposition (see Refs. 38 and 74).

c

Nondimensional (see Refs. 12 and 38 for details on nondimensionalization).

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: SPARTA77 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. N2 is used as the working gas in simulations, since MIKRA experiments3 were performed in N2 medium. The properties of the working gas is given in Table III. DSMC simulations treat N2 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 splitting88 (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.

FIG. 12.

Spatial mesh used for carrying out DGFS simulations for MIKRA Gen1 device. A linear gradient is applied within blocks such that the cells are finer in the near-vane region. A 3rd order nodal/sem DG scheme has been used.

FIG. 12.

Spatial mesh used for carrying out DGFS simulations for MIKRA Gen1 device. A linear gradient is applied within blocks such that the cells are finer in the near-vane region. A 3rd order nodal/sem DG scheme has been used.

Close modal

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.

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 × 1021 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.

FIG. 13.

Variation of flow properties along the domain for MIKRA Gen1 cases (M-01: Kn = 1.85, and M-03: Kn = 0.3) obtained from DSMC (thin black lines), DGFS (thick red lines), BGK (thick blue lines), ESBGK (thick green lines), and S-model (thick orange lines): (a) Kn = 0.3, Number density (m−3); (b) Kn = 1.85, Number density (m−3); (c) Kn = 0.3, Temperature (K); (d) Kn = 1.85, Temperature (K); (e) Kn = 0.3, xy-component of stress (N/m2); (f) Kn = 1.85, xy-component of stress (N/m2).

FIG. 13.

Variation of flow properties along the domain for MIKRA Gen1 cases (M-01: Kn = 1.85, and M-03: Kn = 0.3) obtained from DSMC (thin black lines), DGFS (thick red lines), BGK (thick blue lines), ESBGK (thick green lines), and S-model (thick orange lines): (a) Kn = 0.3, Number density (m−3); (b) Kn = 1.85, Number density (m−3); (c) Kn = 0.3, Temperature (K); (d) Kn = 1.85, Temperature (K); (e) Kn = 0.3, xy-component of stress (N/m2); (f) Kn = 1.85, xy-component of stress (N/m2).

Close modal
FIG. 14.

Continuation of Fig. 13: (a) Kn = 0.3, x-component of heat-flux (W/m2); (b) Kn = 1.85, x-component of heat-flux (W/m2); (c) Kn = 0.3, y-component of heat-flux (W/m2); (d) Kn = 1.85, y-component of heat-flux (W/m2); (e) Kn = 0.3, Speed (m/s); (f) Kn = 1.85, Speed (m/s).

FIG. 14.

Continuation of Fig. 13: (a) Kn = 0.3, x-component of heat-flux (W/m2); (b) Kn = 1.85, x-component of heat-flux (W/m2); (c) Kn = 0.3, y-component of heat-flux (W/m2); (d) Kn = 1.85, y-component of heat-flux (W/m2); (e) Kn = 0.3, Speed (m/s); (f) Kn = 1.85, Speed (m/s).

Close modal
FIG. 15.

Sharp curvature in isotherms near the surface of cold vane at Kn = 0.3. This can be interpreted in terms of imbalance of molecules of type A (cold) and type B (hot) at the top-left/top-right ends of the cold vane.

FIG. 15.

Sharp curvature in isotherms near the surface of cold vane at Kn = 0.3. This can be interpreted in terms of imbalance of molecules of type A (cold) and type B (hot) at the top-left/top-right ends of the cold vane.

Close modal

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 expansion2 

Pij=pδij+τij(1)+τij(2)+τij(1)=2μuixjτij(2)=K2μ2ρT2Txixj+K3μ2ρT2TxiTxjτijT:Thermalstresstensor+K1μ2ρukxkuixj,i,j,k{1,2},
(27)

where Pij, 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 Ki ≈ 1, i = {1, 2, 3} are species/molecular-interaction specific constants.2 This yields

P12=τ12(1)+τ12(2)+=Pxy,τ12T=K2μ2ρT2Tx1x2+K3μ2ρT2Tx1Tx2=K2μ2ρT2Txy+K3μ2ρT2TxTy.
(28)

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×CA/BC×CA=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 Pxy) to be larger between points A and B since ∂T/∂x|AB∂T/∂y|CD (more formally: ∥∇TAB ≫ ∥∇TCD, 2TAB2TCD). 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.

FIG. 16.

xy component of stress tensor at Kn = 0.3: origin of oval/ellipses at the edges of the vanes. Note the distorted shape of the ellipse between the top-right corner of cold vane and top-left corner of hot vane. Since the temperature gradient is stronger between point A and B, compared to point C and D, we expect the semi-major axis of the ellipse to be larger than the semi-minor axis, and hence the distorted ellipse/oval – an observation consistent with Eqs. (27) and (28) since ∂T/∂x∂T/∂y.

FIG. 16.

xy component of stress tensor at Kn = 0.3: origin of oval/ellipses at the edges of the vanes. Note the distorted shape of the ellipse between the top-right corner of cold vane and top-left corner of hot vane. Since the temperature gradient is stronger between point A and B, compared to point C and D, we expect the semi-major axis of the ellipse to be larger than the semi-minor axis, and hence the distorted ellipse/oval – an observation consistent with Eqs. (27) and (28) since ∂T/∂x∂T/∂y.

Close modal

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 N2 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/m2 (the unlabeled contour just below −4000 W/m2 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.

FIG. 17.

Instantaneous streamlines near the vanes of MIKRA Gen1 device at Kn = 0.30 obtained from DGFS using VHS collision model. Observe the vortex formation above the cold vane, and top right corner of hot vane: (a) t = 1.25 × 10−3 s; (b) t = 2.5 × 10−3 s; (c) t = 3.75 × 10−3 s; (d) t = 4.375 × 10−3 s.

FIG. 17.

Instantaneous streamlines near the vanes of MIKRA Gen1 device at Kn = 0.30 obtained from DGFS using VHS collision model. Observe the vortex formation above the cold vane, and top right corner of hot vane: (a) t = 1.25 × 10−3 s; (b) t = 2.5 × 10−3 s; (c) t = 3.75 × 10−3 s; (d) t = 4.375 × 10−3 s.

Close modal

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.

FIG. 18.

Variation of flow speed (m/s) at steady state for MIKRA Gen1 cases obtained from DSMC and DGFS using VHS collision model: (a) Kn = 0.30, DSMC; (b) Kn = 0.30, DGFS; (c) Kn = 0.74, DSMC; (d) Kn = 0.74, DGFS; (e) Kn = 1.85, DSMC; (f) Kn = 1.85, DGFS.

FIG. 18.

Variation of flow speed (m/s) at steady state for MIKRA Gen1 cases obtained from DSMC and DGFS using VHS collision model: (a) Kn = 0.30, DSMC; (b) Kn = 0.30, DGFS; (c) Kn = 0.74, DSMC; (d) Kn = 0.74, DGFS; (e) Kn = 1.85, DSMC; (f) Kn = 1.85, DGFS.

Close modal

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, THTC, 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.

FIG. 19.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and DGFS (lines) using VHS collision model: (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

FIG. 19.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and DGFS (lines) using VHS collision model: (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

Close modal
FIG. 20.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and BGK (lines): (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

FIG. 20.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and BGK (lines): (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

Close modal
FIG. 21.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and ESBGK (lines): (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

FIG. 21.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and ESBGK (lines): (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

Close modal
FIG. 22.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and Shakhov (lines): (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

FIG. 22.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for MIKRA Gen1 cases obtained from DSMC (symbols) and Shakhov (lines): (a) Number density (m−3); (b) Temperature (K); (c) x-component of heat-flux (W/m2); (d) y-component of heat-flux (W/m2); (e) xy-component of stress (N/m2); (f) Speed (m/s).

Close modal

In the present section, we carry out the MIKRA simulations for binary mixture consisting of N2 and H2O using the variable soft sphere model.

The flow configuration remains the same as shown in Fig. 11. We consider the 2D uniform flow of binary mixture of N2 and H2O. The end goal is to simulate the motion of gas flows in the gap between the two vanes, subject to initial pressure p, hot (TH) and cold (TC) 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.

TABLE VIII.

Numerical parameters for thermostress convection in MIKRA Gen1 simulations for DSMC and DGFS using VSS collision model for N2/H2O binary mixture.

Cases
ParameterMSM-01
Pressure: p (Torr) 1.163 
Total number density: n (×1021 m−337.860 91 
Concentration: (n(N2)/n,n(H2O)/n(0.05, 0.95) 
Knudsen number:a Kn 1.85 
Cold vane temperature: TC (K) 306 
Hot vane temperature: TH (K) 363 
DGFS parameters 
Points in velocity mesh: N3 323 
Points in radial direction:bNρ 
Points on full sphere:bM 12 
Size of velocity meshc [−6, 6]3 
Cases
ParameterMSM-01
Pressure: p (Torr) 1.163 
Total number density: n (×1021 m−337.860 91 
Concentration: (n(N2)/n,n(H2O)/n(0.05, 0.95) 
Knudsen number:a Kn 1.85 
Cold vane temperature: TC (K) 306 
Hot vane temperature: TH (K) 363 
DGFS parameters 
Points in velocity mesh: N3 323 
Points in radial direction:bNρ 
Points on full sphere:bM 12 
Size of velocity meshc [−6, 6]3 
a

Based on hard-sphere definition using total number density (see Ref. 11).

b

Required in the fast Fourier spectral low-rank decomposition (see Ref. 39).

c

Nondimensional (see Refs. 12 and 38 for details on nondimensionalization).

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: SPARTA77 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 H2O and 15 particles of N2, 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. N2 and H2O are used as the working gases in simulations. The properties of the working gas is given in Table IX. Note that for N2, we consider ζR = 2 rotational degrees of freedom, rotational relaxation ZR = 0.2, ζV = 2 vibrational degrees of freedom, vibrational relaxation ZV = 1.901 14 × 10−5, and vibrational temperature Tv = 3371 K; and for H2O, we consider ζR = 3 rotational degrees of freedom, rotational relaxation ZR = 0.2, ζV = 3 vibrational degrees of freedom, vibrational relaxation ZV = 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 N2 and H2O as the working gas in simulations. N2 is diatomic, and H2O 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.

TABLE IX.

N2 and H2O gas VSS parameters used in MIKRA Gen1 DSMC and DGFS simulations.

N2H2O
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: dref,i (m) 4.07 × 10−10 5.78 × 10−10 
Ref. temperature: Tref,i (K) 273 273 
N2H2O
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: dref,i (m) 4.07 × 10−10 5.78 × 10−10 
Ref. temperature: Tref,i (K) 273 273 
a

For cross-collision (see Refs. 11 and 12): Ψij = 0.5(Ψi + Ψj), ij, where Ψ = {ω, α, dref, Tref}.

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 N2 and H2O 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 × 1021 m−3 for H2O and 1.893 045 5 × 1021 m−3 for N2, 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 N2, moves slower and the lighter species, here H2O, moves faster giving rise to the well-known thermal diffusion. This explain why the isotherms for H2O spread farther apart compared to those of N2.

FIG. 23.

Variation of flow properties along the domain for multispecies MIKRA Gen1 case (MSM-01: Kn = 1.85) obtained with DSMC (thin black lines) and DGFS (thick red lines) using VSS collision model. We want to reemphasize that DSMC simulations consider the rotational degrees of freedom of N2 and H2O into account, whereas DGFS, being in very early stages of research, does not; and therefore we expect some differences between DSMC and DGFS results: (a) N2, Number density (m−3); (b) H2O, Number density (m−3); (c) N2, Temperature (K); (d) H2O, Temperature (K).

FIG. 23.

Variation of flow properties along the domain for multispecies MIKRA Gen1 case (MSM-01: Kn = 1.85) obtained with DSMC (thin black lines) and DGFS (thick red lines) using VSS collision model. We want to reemphasize that DSMC simulations consider the rotational degrees of freedom of N2 and H2O into account, whereas DGFS, being in very early stages of research, does not; and therefore we expect some differences between DSMC and DGFS results: (a) N2, Number density (m−3); (b) H2O, Number density (m−3); (c) N2, Temperature (K); (d) H2O, Temperature (K).

Close modal
FIG. 24.

Continuation of Fig. 23: (a) N2, Speed (m/s); (b) H2O, Speed (m/s); (c) N2, xy-component of stress (N/m2); (d) H2O, xy-component of stress (N/m2).

FIG. 24.

Continuation of Fig. 23: (a) N2, Speed (m/s); (b) H2O, Speed (m/s); (c) N2, xy-component of stress (N/m2); (d) H2O, xy-component of stress (N/m2).

Close modal

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 H2O compared to N2. 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 N2 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 H2O compared to N2, 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.

FIG. 25.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for multispecies MIKRA Gen1 case obtained with DSMC (symbols) and DGFS (lines) using VSS collision model: (a) Number density (m−3); (b) Temperature (K); (c) xy-component of stress (N/m2); (d) Speed (m/s).

FIG. 25.

Variation of flow properties along the domain vertical centerline (X = 300 μm) for multispecies MIKRA Gen1 case obtained with DSMC (symbols) and DGFS (lines) using VSS collision model: (a) Number density (m−3); (b) Temperature (K); (c) xy-component of stress (N/m2); (d) Speed (m/s).

Close modal

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.

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.

1.
C.-M.
Ho
and
Y.-C.
Tai
, “
Micro-electro-mechanical-systems (MEMS) and fluid flows
,”
Annu. Rev. Fluid. Mech.
30
,
579
612
(
1998
).
2.
M.
Kogan
,
V.
Galkin
, and
O.
Fridlender
, “
Stresses produced in gases by temperature and concentration inhomogeneities. New types of free convection
,”
Usp. Fiz. Nauk
119
,
111
125
(
1976
).
3.
A.
Strongrich
,
A.
Pikus
,
I. B.
Sebastião
, and
A.
Alexeenko
, “
Microscale in-plane knudsen radiometric actuator: Design, characterization, and performance modeling
,”
J. Microelectromech. Syst.
26
,
528
538
(
2017
).
4.
M.
Knudsen
, “
Thermischer molekulardruck der gase in röhren
,”
Ann. Phys.
338
,
1435
1448
(
1910
).
5.
G.
Karniadakis
,
A.
Beskok
, and
N.
Aluru
,
Microflows and Nanoflows: Fundamentals and Simulation
(
Springer Science & Business Media
,
2006
), Vol. 29.
6.
M.
Knudsen
,
The Kinetic Theory of Gases: Some Modern Aspects
(
Methuen
,
1950
).
7.
W.
Crookes
 et al, “
On attraction and repulsion resulting from radiation
,”
Philos. Trans. R. Soc. London
164
,
501
527
(
1874
).
8.
A.
Ketsdever
,
N.
Gimelshein
,
S.
Gimelshein
, and
N.
Selden
, “
Radiometric phenomena: From the 19th to the 21st century
,”
Vacuum
86
,
1644
1662
(
2012
).
9.
J. C.
Maxwell
, “
On stresses in rarified gases arising from inequalities of temperature
,”
Philos. Trans. R. Soc. London
170
,
231
256
(
1879
).
10.
S.
Chapman
and
T.
Cowling
,
The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases
(
Cambridge Mathematical Library, Cambridge University Press
,
1970
), Vol. 1, pp.
27
52
.
11.
G. A.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows
(
Clarendon Press
,
Oxford
,
1994
).
12.
S.
Jaiswal
,
A. A.
Alexeenko
, and
J.
Hu
, “
A discontinuous Galerkin fast spectral method for the multi-species full Boltzmann equation
,”
Comput. Methods Appl. Mech. Eng.
352
,
56
84
(
2019
).
13.
A. V.
Pavlov
, “
Diffusion and thermodiffusion of atmospheric neutral gases: A review
,”
Surv. Geophys.
40
,
247
276
(
2019
).
14.
Y.
Sone
,
Kinetic Theory and Fluid Dynamics
(
Springer Science & Business Media
,
2012
).
15.
E. H.
Kennard
 et al,
Kinetic Theory of Gases, with an Introduction to Statistical Mechanics
(
McGraw-Hill
,
1938
).
16.
Y.
Sone
, “
Thermal creep in rarefied gas
,”
J. Phys. Soc. Jpn.
21
,
1836
1837
(
1966
).
17.
F.
Sharipov
and
V.
Seleznev
, “
Data on internal rarefied gas flows
,”
J. Phys. Chem. Ref. Data
27
,
657
706
(
1998
).
18.
K.
Aoki
,
Y.
Sone
, and
Y.
Waniguchi
, “
A rarefied gas flow induced by a temperature field: Numerical analysis of the flow between two coaxial elliptic cylinders with different uniform temperatures
,”
Comput. Math. Appl.
35
,
15
28
(
1998
).
19.
Y.
Sone
, “
Flow induced by thermal stress in rarefied gas
,”
Phys. Fluids
15
,
1418
1423
(
1972
).
20.
Y.
Sone
and
M.
Yoshimoto
, “
Demonstration of a rarefied gas flow induced near the edge of a uniformly heated plate
,”
Phys. Fluids
9
,
3530
3534
(
1997
).
21.
N.
Selden
,
C.
Ngalande
,
S.
Gimelshein
,
E.
Muntz
,
A.
Alexeenko
, and
A.
Ketsdever
, “
Area and edge effects in radiometric forces
,”
Phys. Rev. E
79
,
041201
(
2009
).
22.
K.
Fowee
,
A.
Ibrayeva
,
A.
Strongrich
, and
A.
Alexeenko
, “
Experimental measurements and numerical modeling of a thermostress convection-based actuator
,”
AIP Conf. Proc.
1786
,
200004
(
2016
).
23.
A.
Ibrayeva
, “
Numerical modeling of thermal edge flow
,” M.S. thesis,
Purdue University
,
2017
.
24.
A.
Passian
,
A.
Wig
,
F.
Meriaudeau
,
T.
Ferrell
, and
T.
Thundat
, “
Knudsen forces on microcantilevers
,”
J. Appl. Phys.
92
,
6326
6333
(
2002
).
25.
A.
Passian
,
R.
Warmack
,
T.
Ferrell
, and
T.
Thundat
, “
Thermal transpiration at the microscale: A Crookes cantilever
,”
Phys. Rev. Lett.
90
,
124503
(
2003
).
26.
V.
Foroutan
,
R.
Majumdar
,
O.
Mahdavipour
,
S.
Ward
, and
I.
Paprotny
, “
Levitation of untethered stress-engineered microflyers using thermophoretic (Knudsen) force
,” in
Technical Digest of the Hilton Head Workshop
(
Transducer Research Foundation
,
2014
), pp.
105
106
.
27.
R. T.
Nallapu
,
A.
Tallapragada
, and
J.
Thangavelautham
, “
Radiometric actuators for spacecraft attitude control
,” in
2017 IEEE Aerospace Conference
(
IEEE
,
2017
).
28.
B. M.
Cornella
,
A. D.
Ketsdever
,
N. E.
Gimelshein
, and
S. F.
Gimelshein
, “
Analysis of multivane radiometer arrays in high-altitude propulsion
,”
J. Propul. Power
28
,
831
839
(
2012
).
29.
H.
Grad
, “
On the kinetic theory of rarefied gases
,”
Commun. Pure Appl. Math.
2
,
331
407
(
1949
).
30.
G.
Bird
, “
Approach to translational equilibrium in a rigid sphere gas
,”
Phys. Fluids
6
,
1518
1519
(
1963
).
31.
G. C.
Maitland
and
E. B.
Smith
, “
Critical reassessment of viscosities of 11 common gases
,”
J. Chem. Eng. Data
17
,
150
156
(
1972
).
32.
K.
Koura
,
H.
Matsumoto
, and
T.
Shimada
, “
A test of equivalence of the variable-hard-sphere and inverse-power-law models in the direct-simulation Monte Carlo method
,”
Phys. Fluids A
3
,
1835
1837
(
1991
).
33.
K.
Koura
and
H.
Matsumoto
, “
Variable soft sphere molecular model for inverse-power-law or Lennard-Jones potential
,”
Phys. Fluids A
3
,
2459
2465
(
1991
).
34.
A.
Kersch
,
W.
Morokoff
, and
C.
Werner
, “
Selfconsistent simulation of sputter deposition with the Monte Carlo method
,”
J. Appl. Phys.
75
,
2278
2285
(
1994
).
35.
J.
Fan
, “
A generalized soft-sphere model for Monte Carlo simulation
,”
Phys. Fluids
14
,
4399
4405
(
2002
).
36.
S.
Jaiswal
,
I.
Sebastião
, and
A. A.
Alexeenko
, “
DSMC-SPARTA implementation of M-1 scattering model
,” in
Proceedings of 31st Rarefied Gas Dynamics Symposium
(
AIP
,
2018
), to appear: http://goo.gl/N7qFao.
37.
A. B.
Weaver
, “
Assessment of high-fidelity collision models in the direct simulation Monte Carlo method
,” Ph.D. thesis,
Purdue University
,
West Lafayette
,
2015
.
38.
S.
Jaiswal
,
A.
Alexeenko
, and
J.
Hu
, “
A discontinuous Galerkin fast spectral method for the full Boltzmann equation with general collision kernels
,”
J. Comput. Phys.
378
,
178
208
(
2019
).
39.
S.
Jaiswal
,
A. A.
Alexeenko
, and
J.
Hu
, “
Fast deterministic solution of the full Boltzmann equation on graphics processing units
,” in
Proceedings of 31st Rarefied Gas Dynamics Symposium
(
AIP
,
2018
), to appear: http://goo.gl/x4A7sy.
40.
S.
Jaiswal
,
J.
Hu
,
J. K.
Brillon
, and
A. A.
Alexeenko
, “
A discontinuous Galerkin fast spectral method for multi-species full Boltzmann on streaming multi-processors
,” in
Proceedings of the Platform for Advanced Scientific Computing Conference, PASC’19
(
ACM
,
New York, NY, USA
,
2019
), pp.
4-1
4-9
.
41.
S.
Loyalka
, “
Knudsen forces in vacuum microbalance
,”
J. Chem. Phys.
66
,
4935
4940
(
1977
).
42.
N.
Selden
,
C.
Ngalande
,
N.
Gimelshein
,
S.
Gimelshein
, and
A.
Ketsdever
, “
Origins of radiometric forces on a circular vane with a temperature gradient
,”
J. Fluid Mech.
634
,
419
431
(
2009
).
43.
J. G.
Fierro
and
A. A.
Garcia
, “
Gas dynamics at low pressures in a vacuum microbalance
,”
Vacuum
31
,
79
84
(
1981
).
44.
A.
Alexeenko
,
E. P.
Muntz
,
M.
Gallis
, and
J.
Torczynski
, “
Comparison of kinetic models for gas damping of moving microbeams
,” in
36th AIAA Fluid Dynamics Conference and Exhibit
(
AIAA
,
2006
), p.
3715
.
45.
T.
Zhu
and
W.
Ye
, “
Origin of Knudsen forces on heated microbeams
,”
Phys. Rev. E
82
,
036308
(
2010
).
46.
J.
Nabeth
,
S.
Chigullapalli
, and
A. A.
Alexeenko
, “
Quantifying the Knudsen force on heated microbeams: A compact model and direct comparison with measurements
,”
Phys. Rev. E
83
,
066306
(
2011
).
47.
Y. A.
Anikin
, “
Numerical study of radiometric forces via the direct solution of the Boltzmann kinetic equation
,”
Comput. Math. Math. Phys.
51
,
1251
1266
(
2011
).
48.
F.
Tcheremissine
, “
Conservative evaluation of Boltzmann collision integral in discrete ordinates approximation
,”
Comput. Math. Appl.
35
,
215
221
(
1998
).
49.
A.
Lotfian
and
E.
Roohi
, “
Radiometric flow in periodically patterned channels: Fluid physics and improved configurations
,”
J. Fluid Mech.
860
,
544
576
(
2019
).
50.
N.
Gimelshein
,
S.
Gimelshein
,
A.
Ketsdever
, and
N.
Selden
, “
Impact of vane size and separation on radiometric forces for microactuation
,”
J. Appl. Phys.
109
,
074506
(
2011
).
51.
A. D.
Strongrich
,
W. J.
O’Neill
,
A. G.
Cofer
, and
A. A.
Alexeenko
, “
Experimental measurements and numerical simulations of the Knudsen force on a non-uniformly heated beam
,”
Vacuum
109
,
405
416
(
2014
).
52.
A.
Alexeenko
and
A.
Strongrich
, “
Microelectromechanical gas sensor based on Knudsen thermal force
,” U.S. patent 15/183,259 (
25 December 2018
).
53.
A.
Strongrich
and
A.
Alexeenko
, “
Microstructure actuation and gas sensing by the Knudsen thermal force
,”
Appl. Phys. Lett.
107
,
193508
(
2015
).
54.
A.
Pikus
,
I. B.
Sebastião
,
A.
Strongrich
, and
A.
Alexeenko
, “
Characterization of a Knudsen force based vacuum sensor for N2H2O gas mixtures
,”
Vacuum
161
,
130
137
(
2019
).
55.
F. J.
McCormack
, “
Construction of linearized kinetic models for gaseous mixtures and molecular gases
,”
Phys. Fluids
16
,
2095
2105
(
1973
).
56.
L.-S.
Luo
and
S. S.
Girimaji
, “
Theory of the lattice Boltzmann method: Two-fluid model for binary mixtures
,”
Phys. Rev. E
67
,
036302
(
2003
).
57.
P. L.
Bhatnagar
,
E. P.
Gross
, and
M.
Krook
, “
A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems
,”
Phys. Rev.
94
,
511
(
1954
).
58.
L.
Sirovich
, “
Kinetic modeling of gas mixtures
,”
Phys. Fluids
5
,
908
918
(
1962
).
59.
P.
Andries
,
K.
Aoki
, and
B.
Perthame
, “
A consistent BGK-type model for gas mixtures
,”
J. Stat. Phys.
106
,
993
1018
(
2002
).
60.
J. R.
Haack
,
C. D.
Hauck
, and
M. S.
Murillo
, “
A conservative, entropic multispecies BGK model
,”
J. Stat. Phys.
168
,
826
856
(
2017
).
61.
A. V.
Bobylev
,
M.
Bisi
,
M.
Groppi
,
G.
Spiga
, and
I. F.
Potapenko
, “
A general consistent BGK model for gas mixtures
,”
Kinet. Relat. Models
11
,
1377
(
2018
).
62.
L. H.
Holway
, Jr.
, “
New statistical models for kinetic theory: Methods of construction
,”
Phys. Fluids
9
,
1658
1673
(
1966
).
63.
S.
Brull
, “
An ellipsoidal statistical model for gas mixtures
,”
Commun. Math. Sci.
13
,
1
13
(
2015
).
64.
E.
Shakhov
, “
Generalization of the Krook kinetic relaxation equation
,”
Fluid Dyn.
3
,
95
96
(
1968
).
65.
K.
Xu
and
J.-C.
Huang
, “
A unified gas-kinetic scheme for continuum and rarefied flows
,”
J. Comput. Phys.
229
,
7747
7764
(
2010
).
66.
Z.
Guo
,
K.
Xu
, and
R.
Wang
, “
Discrete unified gas kinetic scheme for all Knudsen number flows: Low-speed isothermal case
,”
Phys. Rev. E
88
,
033305
(
2013
).
67.
A.
Alexeenko
,
A. G.
Cofer
, and
S. D.
Heister
, “
Microelectronic thermal valve
,” U.S. patent application 15/370,633 (
8 June 2017
).
68.
J. W.
Lee
,
R.
Tung
,
A.
Raman
,
H.
Sumali
, and
J. P.
Sullivan
, “
Squeeze-film damping of flexible microcantilevers at low ambient pressures: Theory and experiment
,”
J. Micromech. Microeng.
19
,
105029
(
2009
).
69.
R. A.
Bidkar
,
R. C.
Tung
,
A. A.
Alexeenko
,
H.
Sumali
, and
A.
Raman
, “
Unified theory of gas damping of flexible microcantilevers at low ambient pressures
,”
Appl. Phys. Lett.
94
,
163117
(
2009
).
70.
M. A.
Gallis
and
J. R.
Torczynski
, “
An improved Reynolds-equation model for gas damping of microbeam motion
,”
J. Microelectromech. Syst.
13
,
653
659
(
2004
).
71.
A.
Granaldi
and
P.
Decuzzi
, “
The dynamic response of resistive microswitches: Switching time and bouncing
,”
J. Micromech. Microeng.
16
,
1108
(
2006
).
72.
C.
Cercignani
,
The Boltzmann Equation and its Applications
(
Springer-Verlag
,
New York
,
1988
).
73.
S.
Harris
,
An Introduction to the Theory of the Boltzmann Equation
(
Dover Publications
,
2004
).
74.
I.
Gamba
,
J.
Haack
,
C.
Hauck
, and
J.
Hu
, “
A fast spectral method for the Boltzmann collision operator with general collision kernels
,”
SIAM J. Sci. Comput.
39
,
B658
B674
(
2017
).
75.
L.
Mieussens
and
H.
Struchtrup
, “
Numerical comparison of Bhatnagar–Gross–Krook models with proper Prandtl number
,”
Phys. Fluids
16
,
2797
2813
(
2004
).
76.
L.
Mieussens
, “
Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries
,”
J. Comput. Phys.
162
,
429
466
(
2000
).
77.
M. A.
Gallis
,
J. R.
Torczynski
,
S. J.
Plimpton
,
D. J.
Rader
,
T.
Koehler
, and
J.
Fan
, “
Direct simulation Monte Carlo: The quest for speed
,”
AIP Conf. Proc.
1628
,
27
36
(
2014
).
78.
M. A.
Gallis
,
N. P.
Bitter
,
T. P.
Koehler
,
J. R.
Torczynski
,
S. J.
Plimpton
, and
G.
Papadakis
, “
Molecular-level simulations of turbulence and its decay
,”
Phys. Rev. Lett.
118
,
064501
(
2017
).
79.
M. A.
Gallis
,
T. P.
Koehler
,
J. R.
Torczynski
, and
S. J.
Plimpton
, “
Direct simulation Monte Carlo investigation of the Rayleigh-Taylor instability
,”
Phys. Rev. Fluids
1
,
043403
(
2016
).
80.
I. B.
Sebastiao
,
L.
Qiao
, and
A. A.
Alexeenko
, “
Direct simulation Monte Carlo modeling of H2-O2 deflagration waves
,”
Combust. Flame
198
,
40
53
(
2018
).
81.
S.
Jaiswal
,
I.
Sebastião
,
A.
Strongrich
, and
A. A.
Alexeenko
, “
FEMTA micropropulsion system characterization by DSMC
,” in
Proceedings of 31st Rarefied Gas Dynamics Symposium (RGD-31)
,
Glasgow, UK
,
23–27 July 2018
(
AIP
, to be published).
82.
A. A.
Alexeenko
,
S. F.
Gimelshein
,
E. P.
Muntz
, and
A. D.
Ketsdever
, “
Kinetic modeling of temperature driven flows in short microchannels
,”
Int. J. Therm. Sci.
45
,
1045
1051
(
2006
).
83.
H.
Yamaguchi
,
P.
Perrier
,
M. T.
Ho
,
J. G.
Méolans
,
T.
Niimi
, and
I.
Graur
, “
Mass flow rate measurement of thermal creep flow from transitional to slip flow regime
,”
J. Fluid Mech.
795
,
690
707
(
2016
).
84.
H.
Akhlaghi
and
E.
Roohi
, “
Mass flow rate prediction of pressure–temperature-driven gas flows through micro/nanoscale channels
,”
Continuum Mech. Thermodyn.
26
,
67
78
(
2014
).
85.
Q.
Sheng
,
G.-H.
Tang
,
X.-J.
Gu
,
D. R.
Emerson
, and
Y.-H.
Zhang
, “
Simulation of thermal transpiration flow using a high-order moment method
,”
Int. J. Mod. Phys. C
25
,
1450061
(
2014
).
86.
R.
Biswas
,
K. D.
Devine
, and
J. E.
Flaherty
, “
Parallel, adaptive finite element methods for conservation laws
,”
Appl. Numer. Math.
14
,
255
283
(
1994
).
87.
B.
Cockburn
,
G. E.
Karniadakis
, and
C.-W.
Shu
, “
The development of discontinuous Galerkin methods
,” in
Discontinuous Galerkin Methods
(
Springer
,
2000
), pp.
3
50
.
88.
A.
Bakhtiari
,
D.
Malhotra
,
A.
Raoofy
,
M.
Mehl
,
H.-J.
Bungartz
, and
G.
Biros
, “
A parallel arbitrary-order accurate AMR algorithm for the scalar advection-diffusion equation
,” in
Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC’16
(
IEEE Press
,
Piscataway, NJ, USA
,
2016
), pp.
44-1
44-12
.