# Extended Macroscopic Study of Dilute Gas Flow within a Microcavity

Authors: Mohamed Hssikou; Jamal Baliti; Mohammed Alaoui
Journal: Modelling and Simulation in Engineering
 (2016)

---

## Abstract

The behaviour of monatomic and dilute gas is studied in the slip and early transition regimes using the extended macroscopic theory. The gas is confined within a two-dimensional microcavity where the longitudinal sides are in the opposite motion with constant velocity±Uw. The microcavity walls are kept at the uniform and reference temperature T0. Thus, the gas flow is transported only by the shear stress induced by the motion of upper and lower walls. From the macroscopic point of view, the regularized 13-moment equations of Grad, R13, are solved numerically. The macroscopic gas proprieties are studied for different values of the so-called Knudsen number (Kn), which gives the gas-rarefaction degree. The results are compared with those obtained using the classical continuum theory of Navier-Stokes and Fourier (NSF).


---

## Body

## 1. Introduction

Recently, the technology of the Microelectromechanical Systems (MEMS) has greatly developed and they have wide areas of application [1–3]. This fast growth of MEMS use is not followed enough by the physical understanding of rarefied gas flows in these microdevices. For this purpose, several studies have been recently focused on for more understanding of the physical phenomena involved in these small devices [4]. In fact, the performances of MEMS often defy the predictions made using the scaling laws developed for large systems. In fact, the gas flows inside the MEMS, under the standard conditions, are usually characterized by a mean free path λ comparable to the system characteristic-length L. Thus, the so-called Knudsen number Kn~λ/L of gas flow is in the slip-transition regimes range; that is, 0.001<Kn≤10. In this case, the conventional computational fluid dynamics (CFD) scheme, based on the classical Navier-Stokes and Fourier (NSF) equations, becomes inappropriate to describe the gas flow behaviour in MEMS devices. Therefore, the Knudsen number, in MEMS, is not sufficiently small to guarantee the validity of the NSF equations and the processes in MEMS need to be modelled with more accurate transport models. Similar rarefaction effects can be found in the problems of gas flows under low pressure and atmospheric conditions [5]. For gas flows outside the hydrodynamic regime (Kn>0.001) [6], many interesting rarefaction effects such as velocity-slip and temperature-jump at the walls [7–10], Knudsen paradox, Knudsen layers [11], transpiration flow [12, 13], thermal stress [14], and heat flux without temperature gradients can take place [15]. Hence, there is a pressing need to develop the more accurate methods allowing a good description of gas-dynamic processes into these microsystems. The direct simulation Monte Carlo (DSMC) is the largely kinetic method used to simulate a rarefied gas flow where the behaviour is mainly described by the Boltzmann equation [16]. The accuracy of this method is proved by many previous studies especially with the actual computers capabilities. But, the computational cost and fluctuations noises, especially in the low-signal flows, remain the major inconveniences of this kinetic method [17]. Indeed, many macroscopic approaches are proposed such as the Chapman-Enskog (CE) expansion and the Grad moments theory. At the first order of CE both approaches lead to the famous laws of Navier-Stokes and Fourier. However, on one hand, the instability of Burnett equations obtained at second order of CE expansion is the main problem of this approach. On the other hand, the Grad 13-moment equations are hyperbolic in nature, yielding finite wave speeds, and discontinuous subshock structures when the Mach number lies above Ma≥1.65 [18]. We also note that the Grad 13-moment equations, for nonlinear problems, lack suitable boundary conditions. But, when the rarefaction degree becomes more intense, the Knudsen number value is in the range of 1<Kn<10, and a purely kinetic approach is needed for describing the gas flows [19].Keeping the benefits of both approaches and to avoid their failures, Struchtrup has adopted recently the combination of the above approaches [19]. This leads to the set of regularized 13-moments equations (R13) used as higher order of continuum solution, that is, O(Kn3), to capture the rarefaction effects, described by the nonlinear terms. The main goal of this paper is to investigate the behaviour of a dilute gas flow inducing only the longitudinal shear stress using the classical theory of NSF, with slip and jump boundary conditions, and the regularized 13-moment equations of Grad approaches. In this study, the rarefaction effects are evaluated in the slip and early transition regimes range. We assume that the gas flow is induced with no synergetic contributions from external force fields.


## 2. Statement of Problem

A monatomic and Maxwell-molecules gas, where the collisions rate is independent of the collision-patterns velocities, is confined within a square microcavity. In this gas, the particles are interacting via a potentialV(r)~r-4 with r being the interparticles distance [16]. The orthogonal cross section and the origin of the coordinate system are shown in Figure 1. The upper and lower sides are in the opposite motion with a constant velocity Uw=±50m/s. The microcavity walls are kept at uniform and environmental temperature T0=273K. The macroscopic proprieties of the gas are evaluated for different values of the Knudsen number in the slip and early transition regime, 0.01≤Kn≤0.5.Figure 1
Microcavity geometry.In the kinetic point of view, the macroscopic gas proprieties can be estimated by solving the Boltzmann equation given by(1)∂∂t+c·∇rf=Jf,f′.In the steady state, the main unknown of this equation is the velocity distribution function f(r,c) in which r=(x,y) is the position vector and c=(cx,cy,cz) is the molecular velocity vector. The distribution function is defined such that the quantity f(r,c)drdc is the number of particles in the phase-volume element drdc near the point (r,c). The analytical solution of this kinetic equation is difficult for the multidimensional systems. This is because of the collisions integral which must be computed at each physical-space and velocity-space points. In many of the previous studies, the so-called Bhatnagar-Gross-Krook kinetic (BGK) model is the most used approximation to simplify the collisional operator J. In this case, the collisions are modelled by the relaxation of f towards the Maxwellian equilibrium distribution fM that is a function of the particle velocity vector c and fluid macroscopic quantities, that is, mass density ρ, hydrodynamic velocity v=(vx,vy), and temperature T. All macroscopic characteristics of the gas flow such asmass density,hydrodynamic velocity,hydrostatic pressure,pressure tensor, andheat flux vector can be related in the steady state to the distribution function f(r,c), respectively, as follows:(2)ρ=m∫fr,cdc,ρv=m∫cfr,cdc,p=m3∫C2fr,cdc,Pij=m∫CiCjfr,cdc,q=m2∫C2Cfr,cdc,where m,n, and C=c-v are, respectively, the molecular mass, number density, and peculiar velocity.


## 3. Continuum Description

The evolution equations for density, velocity, and internal energyε given by the conservation laws of mass, momentum, and energy can be written in the general form, if one neglects the external forces, as follows:(3a)∂tρ+div⁡ρv=0,(3b)∂tρv+div⁡ρvvT+pI+σ=0,(3c)∂tρε+12ρvi2+div⁡ρεv+12ρv2v+pv+σv+q=0,where p=ρRT=ρθ=2/3ρε is the hydrostatic pressure in which θ corresponds to the temperature in energy units. To close the above system, one needs to define the constitutive equations including the stress tensor σij=Pij-pδij and heat flux vector q with δij is the Kronecker delta.The Grad’s moment method is the best-known alternative of the CE expansion [20]. Multiplication Boltzmann equation by functions φA=m{1,ci,1/2C2,CiCj,1/2C2Ci} and subsequent integration over velocity space yields the 13-moment equations of Grad. To solve the closure problem of this set of equations, Grad constructed a phase density by an expansion of the Maxwellian into the Hermite polynomials. However, the discontinuous subshock structures problem is the most failure of this method. Recently, Struchtrup and his coauthors show [21–23] that the Grad’s approach accuracy can be further improved up O(Kn3) in the CE expansion by the regularization of the thirteen Grad’s moments. Note that (3a), (3b), and (3c) can be found easily by setting φA=m,mci and 1/2mC2, respectively, and using the invariance property of the collisions operator and right side of (1). The two other values ofφAlead to the nonlinear constitutive equations of stress tensor σij and heat flux vector qi given, respectively, in the so-called R13 set of equations by [19] (4)∂σij∂t+∂σijvk∂xk+45∂qi∂xj+2p∂vi∂xj_+2σki∂vj∂xk+∂mijk∂xk=-pμσij_,∂qi∂t+∂qivk∂xk+52p∂θ∂xi_+52σik∂θ∂xk+θ∂σik∂xk-σikθρ∂ρ∂xk-σijρ∂σjk∂xk+75qk∂vi∂xk+25qk∂vk∂xi+25qi∂vk∂xk+12∂Rik∂xk+16∂Δ∂xi+mijk∂vj∂xk=-23pμqi_.The indices in the angular brackets denote the symmetric trace-free parts of tensors. The condensed notation, with indices, used here can be found with more details in the appendix of [19]. The above equations contain the additional, nonlinear, quantities mijk, Rik, and Δ given by considering the appropriate terms of higher order moment equations as (see [24])(5a)Δ=-σklσklρ+6σklσklNSFρ+566qkqkNSFρ-12μpθ∂qk∂xk-θqk∂ln⁡p∂xk,(5b)Rij=-47σkiσjkρ+247σkiσjkNSFρ+19275qiqjNSFp-245μpθ∂qi∂xj-θqi∂ln⁡p∂xj,(5c)mijk=815σijqkNSFp+45qiσjkNSFp-2μpθ∂σij∂xk-θσij∂ln⁡p∂xk.If one neglects these quantities, the R13 set of equations is reduced to the well-known set of 13-moment equations of Grad. At the small values of Kn (Kn<1), the stress tensor and heat flux vector can be written in the CE expansion of Knudsen number powers as(6)σij=σij0+Knσij1+Kn2σij2⋯qi=qi0+Knqi1+Kn2qi2⋯.Note that the zeroth order of this expansion corresponds to Euler equations. An analysis of the order of Knudsen number in (4)shows that only the underlined terms are of first order in Knudsen number, while all others are of higher order [25]. Thus, the underlined terms are the NSF contributions which lead to [26](7a)σijNFS=-2μ∂vi∂xj,(7b)qiNFS=-κ∂θ∂xi.To solve the R13-moment equations one needs a set of boundary conditions for the normal components of the fluxes. Torrilhon and Struchtrup [27] show that vn, στn, qn, Rτn, mnnn, and mττn are the only normal components that can be prescribed. The subscripts τ and n denote the tangential and normal, pointing towards the gas, components of the tensors. The Maxwell accommodation model [28] for the boundary conditions in kinetic theory was used to derive boundary conditions for the R13 system, which link the moments of the gas in front of the wall to the tangential wall velocity vτw and wall temperature θw. The details of the construction procedure for the boundary conditions (BCs) for the R13 equations can be found in [27]; the result reads(8a)vn=0,(8b)στn=-χ2-χ2πθPVτ+15qτ+12mτnn,(8c)qn=-χ2-χ2πθ2PT-12PVτ2+12θσnn+115Δ+528Rnn,(8d)Rτn=χ2-χ2πθ6PTVτ+PθVτ-PVτ3-115θqτ-12θmτnn,(8e)mnnn=χ2-χ2πθ25PT-35PVτ2+75θσnn+175Δ-114Rnn,(8f)mττn=-χ2-χ2πθ15PT-45PVτ2+114Rττ+θσττ-15θσnn+1150Δ,where P=ρθ+1/2σττ-1/120Δ/θ-1/28Rττ/θ,Vτ=vτ-vτw and T=θ-θw are the velocity-slip and temperature-jump at the vicinity of the wall, respectively. The parameter χ is the Maxwell accommodation coefficient where the value is χ=1 for the full diffuse reflection and χ=0 for the specular one. The boundary conditions for NSF solution, first order of velocity-slip and temperature-jump, can be derived from the R13 ones by ignoring the higher order terms, thatis,Rnn, mτnn, and Δ:(9a)vn=0,(9b)στnNSF=-χ2-χPVτ-34μ∂θ∂xτ,(9c)qnNSF=-χ2-χ2PT-12PVτ2-μθ∂vn∂xn.In this case P reduces to P=ρθ+1/2σττ.


## 4. Numerical Scheme

The present scheme is based mainly on the algorithm developed recently by Rana et al. [24]. To solve the differential equations (3a)–(4) in the steady state, it is more convenient to rewrite the system in the matrix and dimensionless form as in [24] to facilitate their numerical resolution:(10)Aω∂ω∂x+Bω∂ω∂y+1KnPωω=0,where ω=[ρ,vx,vy,θ,qx,qy,σxx,σxy,σyy, Rxx,Rxy,Ryy, mxxx,mxyy,mxxy,myyy,Δ] is the vector of field variables and Aω and Bω are the coefficient matrices in x- and y-directions, respectively. The matrix Pω is the production matrix. These matrices are given in Appendix  A of [24]. The matrix equation (10) can be resolved using the finite difference approach in which the BCs (8a)–(8f), in the matrix form, are implemented. Thus, the cavity domain is discretised on a network,N×M, of nodes points. The spatial coordinates (xi,yj) of each node, including the boundary ones, are given by xi=iΔx and yj=jΔy with i=0,N+1 and j=0,M+1. Using the indices notation, (9a), (9b), and (9c) can be written for each node (i,j) as(11)Ai,j∂ωi,j∂x+Bi,j∂ωi,j∂y+1KnPi,jωi,j=0.i,j∈1,N×1,M,where ωi,j=ω(xi,yj) and the matrix are functions of ωi,j.The boundary conditions (8a)–(8f) are incorporated in matrix form, and linear extrapolation is used for boundary nodes at x-walls and y ones, that is, wall normal pointing in x- or y-direction:(12a)ω0,j=X1,j+2ω1,j-ω2,j+X1,jd+,(12b)ωN+1,j=XN,j-2ωN,j-ωN-1,j+XN,jd-j∈1,M,(12c)ωi,0=Yi,1+2ωi,1-ωi,2+Yi,1d+,(12d)ωi,M+1=Yi,M-2ωi,M-ωi,M-1+Yi,Md-i∈1,N,where the superscript signs denote the positive/negative directions of the wall normal. The matrices X and Y correspond to the coefficient matrix of field variable ω in (8a)–(8f)but the vectors Xd and Yd are related to surface inhomogeneity such as wall temperature T0 and velocity Uw. The expressions of these matrices are given in the appendix of [24]. By replacing the derivatives in (10) using a central difference scheme we obtain (13)Aijωi+1,j2Δx-Aijωi-1,j2Δx+Bijωi,j+12Δy-Bijωi,j-12Δy+1KnPi,jωi,j=0.i,j∈2,N-1×2,M-1.Thus, by combining this equation with the boundary conditions given by  (12a), (12b), (12c), and (12d)  and setting the indices i and j one gets the equations for different cavity walls. The solution convergence criterion is chosen to be ~O(Δx3,Δy3) which corresponds in our simulation to 10-6.For more convenience, the results are shown with the dimensionless variables using the following normalisations:X=x/L; Y=y/L; Vx=vx/v0; Vy=vy/v0; Qx,y=qx,y/(P0v0); T=T/T0; ρ-=ρ/ρ0, θ-=θ/θ0≡T; σ-xy=σxy/(ρ0θ0), where ρ0 is the reference mass density related the Knudsen number Kn and reference viscosity μ0 by the relation ρ0=μ0/KnLθ0 and P0=ρ0θ0 is the reference pressure.


## 5. Results and Discussion

In the present paper, we are interested in the study of the viscous and thermal behaviours of a dilute gas flow within the microcavity. The algorithm of this new approach is tested for this standard geometry for two crucial reasons. In the first hand, the cavity is the classical geometry encountered at the lot of devices in engineering area. The second reason is to compare with the previous study especially those of the one-dimension problem such as the classical Couette flow between two parallel plates. The rarefaction effects are discussed in the slip and early transition regimes. The streamlines of gas flow obtained using both NSF and R13 solutions forKn=0.01 and Kn=0.1 are sketched in Figures 2(a)–2(d). We see that the plots are almost similar and both solutions predict a main vortex provided by the walls motion. The secondary vortex observed at the low rarefied case Kn=0.01 in the cavity centre disappears at the rarefied one (Kn=0.1). To clarify the velocity variation more, especially, the longitudinal component, the profiles of Vx are plotted as a function of y-coordinate for different values of Kn, Figures 3(a) and 3(b). In the classical NSF solution, without velocity-slip and temperature-jump BCs, Vx is a linear function of y-coordinate between ±Uw/v0~±0.2. Using the first order of slip and jump (BCs) we prove the validity NSF solution for velocity profiles up Kn=0.1, slip regime, but above this value the solution is not sensitive to the rarefaction effects contrary to the new approach of R13. On the other hand, both approaches predict that the length of so-called Knudsen layer which corresponds to the nonequilibrium area near the walls is about 0.2L.Figure 2
(a, b) Streamlines obtained by NSF (a) and R13 (b) solutions forKn=0.01. (c, d) Streamlines obtained by NSF (c) and R13 (d) solutions for Kn=0.1.
(a) (b) (c) (d)Figure 3
(a, b)Vx profiles obtained by NSF and R13 solutions for Kn=0.01 (square), 0.05 (diamond), 0.1 (triangle), and 0.5 (circle).
(a) (b)The thermal behaviour of a gas is studied through the temperature and heat flux changes. Figures4(a) and 4(b) show the temperature profiles as a function of y-coordinate. As in the classical study of Couette flow problem [29, 30], the parabolic nature, thatis,Ty~Ay2, with positive coefficient A>0, of temperature profile is shown by both solutions for Kn=0.01. Note that the bearing observed in the temperature profiles curvature is due to the low speed or signal problem which is also the major difficulty encountered in the kinetic method like DSMC. With the increase of rarefaction effects, R13 solution predicts a strong inversion in the sharp of profiles A<0 especially for Kn=0.1; however, the NSF one is not well sensitive to the Knudsen value. This inversion can be explained by the prevalence of particle-wall collisions against to the particle-particle ones which take place at low rarefied cases.Figure 4
(a, b) Temperature profiles obtained by NSF (a) and R13 (b) solutions forKn=0.01 (square), 0.05 (diamond), and 0.1 (triangle).
(a) (b)The heat flux transfer without the temperature gradient is another important result provided by the present paper. Figures5(a) and 5(b) give x-component of heat flux vector as a function of y-coordinate. At the continuum limit, Kn=0.01, both NSF and R13 solutions give almost a constant value; Qx≈0. But, with the increasing of rarefaction degree the R13 approach, unlike the NSF one, leads to horizontal heat flux profiles compatible with the opposite motion of the top/bottom walls. This 1D plots of Qx can be interpreted as an extension of the parallel heat flux component expression given, for gas flow between two infinite plates in opposite motion, by [29](14)qx=B2dvxdydTydy,where B is a constant relative to the particles-interaction model. Note that this expression is obtained from the Burnett equations, second order approximation of CE expansion [31, 32]. At the continuum limit, x-velocity component and temperature profiles show that vxy~y and T(y)~y2. This agrees with the linear variation of qxy predicted by (14). Regarding Qx values given by both solutions (Figures 5(a) and 5(b)), we conclude also that NSF approach predicts almost constant parallel heat flux opposite to R13 one which gives a strong variation. For more details, a two-dimensional plot (Figures 6(a)–6(d)) is used to describe both heat flux transfer and viscous stress behaviours in the whole cavity domain. In Figures 6(a), 6(b), 6(c), and 6(d), the heatlines are sketched on the shear stress contours. These figures confirm the breakdown of NSF solution and the accuracy of R13 one for gas flow with no neglected rarefaction effects.Figure 5
(a, b)Qx profiles obtained by NSF (a) and R13 (b) solutions for Kn=0.01 (square), 0.05 (diamond), and 0.1 (triangle).
(a) (b)Figure 6
(a, b) Heatlines superimposed on the shear stressσ-x,y contours for Kn=0.05. (c, d) Heatlines superimposed on the shear stress σ-x,y contours for Kn=0.1.
(a) (b) (c) (d)


## 6. Conclusion

In this paper, we have analysed, numerically, the steady planar Couette flow of a dilute and monatomic gas beyond the scope of the classical Navier-Stokes description using velocity-slip and temperature-jump boundary conditions as a function of rarefaction degree. The behaviour of the gas in rectangular cavity is described by means of extended macroscopic approach R13. The results confirm that with slip and jump boundary conditions, the NS description can be improved. The difference between the results is due to the rarefaction effect which becomes important and the local equilibrium is broken down. The NSF equations cannot explain the variation of heat flux transfer andxy-component of shear stress tensor. Therefore, the R13 approach is more convenient to capture rarefaction effects heat transfer problem of gas flow in the slip and early transition regimes. Thus, this method can be considered as a good alternative to the kinetic approaches, like DSMC, which are expensive in the time cost.



---

*Source: 7619746-2016-12-14.xml*