Nonlinear oscillatory rarefied gas flow inside a rectangular cavity

The nonlinear oscillation of rarefied gas flow inside a two-dimensional rectangular cavity is investigated on the basis of the Shakhov kinetic equation. The gas dynamics, heat transfer, and damping force are studied numerically via the discrete unified gas-kinetic scheme for a wide range of parameters, including gas rarefaction, cavity aspect ratio, and oscillation frequency. Contrary to the linear oscillation where the velocity, temperature, and heat flux are symmetrical and oscillate with the same frequency as the oscillating lid, flow properties in nonlinear oscillatory cases turn out to be asymmetrical, and second-harmonic oscillation of the temperature field is observed. As a consequence, the amplitude of the shear stress near the top-right corner of the cavity could be several times larger than that at the top-left corner, while the temperature at the top-right corner could be significantly higher than the wall temperature in nearly the whole oscillation period. For the linear oscillation with the frequency over a critical value, and for the nonlinear oscillation, the heat transfer from the hot to cold region dominates inside the cavity, which is contrary to the anti-Fourier heat transfer in a low-speed rarefied lid-driven cavity flow. The damping force exerted on the oscillating lid is studied in detail, and the scaling laws are developed to describe the dependency of the resonance and antiresonance frequencies (corresponding to the damping force at a local maximum and minimum, respectively) on the reciprocal aspect ratio from the near hydrodynamic to highly rarefied regimes. These findings could be useful in the design of the micro-electro-mechanical devices operating in the nonlinear-flow regime.


I. INTRODUCTION
Oscillatory gas motion at the micro-and nanoscale has attracted significant attention due to the development of microelectro-mechanical systems (MEMS) including the microaccelerometer, inertial sensor, and actuators [1].Since the surface-area-to-volume ratio in MEMS devices is large, surface forces such as the damping exerted on the oscillating parts by the gas are important for the design of moving microdevices [2].
To understand and quantify the damping, the gas rarefaction effect, caused by the small characteristic length scale of the MEMS, should be taken into account.The degree of gas rarefaction is normally categorized by the Knudsen number (Kn), which is defined as the ratio of the mean free path of gas molecules to the characteristic flow length.Most MEMS devices operate in the slip (10 −3 Kn 0.1) and early transition regimes (0.1 Kn 1) [3,4], where the Navier-Stokes equations fail, and the kinetic theory approach should be adopted in study of rarefied gas dynamics [5,6].
In general, the DSMC and DVM are robust approaches for the simulation of rarefied gas flows in respectively high-speed and low-speed flow regimes.However, most of the aforementioned studies were focused on the oscillatory gas flows in the transition regime (0.1 Kn 10) [5,7,9,10,15,18,20,21].This is mainly due to the well-known intrinsic limitation that the computational time step and spatial mesh size in these two methods should be smaller than the mean collision time and the mean free path of gas molecules, respectively, resulting in an extremely expensive computations in the near hydrodynamic regime, i.e., Kn 0.1.Special slip boundary treatments were introduced to improve the predictive capability in studying the oscillatory rarefied flow in the slip and early transition regimes [3,13,22].However, these methods can only capture some of the bulk properties [12].
The 1D oscillatory Couette flow has been studied extensively, while the 2D oscillatory flow received less attention.Wu et al. has studied the oscillatory rarefied gas flow [17] and sound propagation in a 2D rectangular cavity [19] by using the fast spectral method (FSM) for the linearized Boltzmann equation [23,24].It was reported that the damping force on the oscillating lid in a 2D cavity could be even smaller than that in 1D oscillatory Couette flow, due to the antiresonance in gas [17].Also, the damping force and sound speed near the oscillating lid are strongly affected by the aspect ratio [19].In these works, the velocity amplitude of the oscillating lid is assumed to be sufficiently small, so that the oscillatory flow can be regarded as a linear system.However, the nonlinear oscillation may be possible; for example, the oscillation frequency of MEMS device could reach 110 GHz, and the oscillation amplitude could reach several nanometers, so that the oscillation speed could be 330 m/s which is comparable or even larger than the sonic speed [16,[25][26][27][28].The nonlinear oscillation may dramatically change the flow behavior and damping.Besides, high oscillatory frequency could further aggravate the degree of rarefaction [25].Therefore, to accurately predict the gas damping in nonlinear oscillatory microsystems, the nonlinearity and compressibility of the flow should be taken into account.
In fact, Tsuji [16] et al. and Aoki et al. [27] recently applied the Boltzmann model equation to study the nonlinear motion of a rarefied gas caused by a plate oscillating in its normal direction (1D problem).As far as the authors are aware of, however, direct investigation of the nonlinear oscillation in lid-driven cavity flows from the hydrodynamic to free molecular regimes is absent.Little effort has been devoted to understanding the heat transfer in nonlinear rarefied oscillatory flow, which plays a significant role in MEMS devices [16,27,29,30].
This work is devoted to studying the nonlinear oscillatory rarefied gas flow in a 2D rectangular cavity on the basis of the gas kinetic equation.We introduce the Shakhov model in Sec.II and the numerical scheme in Sec.III.Numerical results of the gas dynamics, heat transfer, and damping force in both linear and nonlinear oscillations over a wide range of gas rarefaction, cavity aspect ratio, and oscillation frequency are presented, compared, and discussed in Sec.IV, which is followed by the conclusions in Sec.V.

II. PROBLEM FORMULATION
We consider a rarefied gas flow in a 2D rectangular cavity driven by a lid at y = H .The lid oscillates harmonically in the x direction with a frequency ω, see Fig. 1.The velocity of the oscillating lid is given by where U 0 is the amplitude of the oscillating velocity, and t is the time.The other three walls at x = 0, x = L, and y = 0 are fixed, and all four walls are isothermal with a fixed temperature T w .The problem considered is characterized by the aspect ratio of the cavity A, the Mach number Ma, the Strouhal number St, and the Knudsen number Kn, which are, respectively, defined by where γ is the specific-heat ratio, v m = √ 2RT w is the most probable molecular velocity with R being the specific gas constant, and λ is the mean free path of gas molecules, which is related to the shear viscosity μ of the gas as where p = ρRT w is the pressure and ρ is the density.In this work, the hard-sphere molecular model is used, where the shear viscosity of the gas is determined by the gas temperature in the following way: The gas kinetic theory is used to investigate how the flow responds to variations of the Knudsen number, the Strouhal number, the Mach number, and the aspect ratio of the cavity.To this end, the Shakhov equation is adopted, which has been widely used to describe the dynamics of monotonic gases [31].In the absence of external force, it takes the form of where f (x,ξ ,t) is the velocity distribution function of gas molecules at the position x = (x,y,z) and time t, with ξ = (ξ x ,ξ y ,ξ z ) being the molecular velocity.f S is the reference distribution function expressed by the Maxwellian distribution function f eq and a heat flux correction term: where c = ξ − U is the peculiar velocity with U being the macroscopic flow velocity, q = 1 2 cc 2 f dξ is the heat flux, and T is the temperature of the gas.The collision time τ in Eq. ( 5) is related to the viscosity μ and the pressure p = ρRT by τ = μ/p.The Maxwellian distribution function f eq is given by The conservative variables W ≡ (ρ,ρU,ρE) T are calculated from the velocity moments of the velocity distribution function: where ψ = (1,ξ , 1 2 ξ 2 ) T .Note that the gas temperature is related to the total gas energy as ρE = 1  2 ρU 2 + 3 2 ρRT .
Since only the 2D problem is considered in this work, two reduced velocity distribution functions are introduced to cast the three-dimensional molecular velocity space into 2D space [32]: For convenience, in what follows we denote ξ = (ξ x ,ξ y ) and x = (x,y).The governing equations for the two reduced distribution functions can be deduced from Eq. ( 5) as where the reduced reference distribution functions are g S = g eq + g Pr and h S = h eq + h Pr , with h eq = RT g eq , (11b) ) It is clear that the governing equations for g and h in Eq. ( 10) can be expressed in the following compact form: where the symbol φ is used to represent g or h.

III. NUMERICAL METHOD
The discrete unified gas-kinetic scheme (DUGKS) is used to solve the Shakhov equation [33,34], which is essentially an explicit finite-volume method.With the intrinsic coupling of molecular collision and transport processes in determining of the flux across the cell interface, the computational time step and mesh size are not limited by the mean collision time and mean free path of gas molecules, respectively, so that the multiscale-flow physics can be efficiently and self-adaptively captured from the hydrodynamic to the free-molecular-flow regimes.In addition, DUGKS has second-order accuracy in both temporal and spatial discretizations.The capacity of DUGKS has been evaluated for gas flows by comparing with the lattice Boltzmann method [35][36][37][38], the spectral method [38,39], the DVM [40], and the DSMC [41,42].It has also been successfully applied to study the phonon transport [43], turbulent flows [38,39,44], and thermally induced nonequilibrium flows [42].
To solve Eq. ( 12), the computational domain is first divided into control cells (the trapezoidal cells are applied in this study).By integrating Eq. ( 12) in a cell V j (centered at x j ) from time t n to t n+1 , and by using the midpoint rule for the temporal integration of the convection term and the trapezoidal rule for the collision term, one obtains the following evolution equation [33]: where t = t n+1 − t n , and is the flux across the cell interface, |V j | and ∂V j are the volume and surface of cell V j , n is the outward unit vector normal to cell interface, and φ n j = V j φ(x,ξ ,t n )d x/|V j | is the cell average of the distribution function.We have also introduced two auxiliary functions: φ = φ − ( t/2) and φ+ = φ + ( t/2) .
Since the mass, momentum, and energy are conserved during the collision, macroscopic variables can be calculated from ) and the heat flux from In practical computations, the evolution of φ is tracked according to Eq. ( 13), instead of the original distribution functions φ, to avoid the implicit computations.
The key procedure in updating φ is to evaluate the flux F in Eq. ( 14), which is solely determined by the gas distribution function φ n+1/2 (x f ,ξ ).To do so, in DUGKS, Eq. ( 12) is integrated along the characteristic line within a half time step s = t/2 [33], where time integration of the collision term is approximated by the trapezoidal rule.Again, to remove the implicitness in Eq. ( 17), we introduce two other distribution functions: Then Eq. ( 17) is expressed explicitly as Therefore, once φ+,n (x f − ξ s,ξ ) is obtained, the distribution function φn+1/2 (x f ,ξ ) can be determined from Eq. ( 19).Note that the conservative variables and heat flux can also be evaluated from φ as and which means that W (x f ,t n+1/2 ) and φ S,n+1/2 (x f ,ξ ) can be obtained from φn+1/2 (x f ,ξ ) directly.Therefore, the original distribution function across the cell interface can be calculated from Eq. ( 18): In this work, φ+,n ( where σ j is the slope of φ+ in the cell j , which is computed by the central difference method.Note that σ j can also be approximated by some gradient limiters for problems with discontinuities [41].The flux can then be computed from Eq. ( 22).Finally, φ at the cell center can be updated according to Eq. ( 13).
In numerical simulations, the continuous molecular velocity space is discretized into a finite discrete velocity set {ξ i } [32].Distribution functions such as g and h are defined at a number of discrete velocity points as gi and hi .Then, the appropriate quadrature rule is adopted to calculate the moments: , where i is the weight for the corresponding quadrature rule.Note that the time step of explicit time marching in the DUGKS is solely determined by the Courant-Friedrichs-Lewy (CFL) condition, i.e., t = η x min /ξ max , where η is the CFL number, x min is the minimum grid spacing, and ξ max is the maximum discrete velocity.

IV. RESULTS AND DISCUSSION
The numerical simulations are performed covering a wide range of rarefaction: Kn = 10, 1, 0.1, 0.01, and 0.001.Two Mach numbers, i.e., Ma = 0.01 and 1.2, are chosen to represent the linear and nonlinear oscillations, respectively.Various aspect ratios and Strouhal numbers are considered.
In the simulation, all the cavity walls are assumed to be fully diffuse; that is, the velocity distribution of the reflected gas molecules at the walls is Maxwellian with the wall temperature and velocity: with where x w is a point at the wall located at a cell interface, and n is the outward unit normal vector at the wall pointing to the cell.Note that the damping force, which is the average shear stress acting on the oscillating lid: is an important parameter in the design of MEMS devices.
According to the transformation introduced in Sec.III, the shear stress P xy is calculated as with U and V being macroscopic flow velocities along the horizontal and vertical directions, respectively.

A. Grid convergence and model accuracy
When Ma = 0.01, the 2D continuous molecular velocity space (ξ x ,ξ y ∈ [−4,4]) is discretized by the trapezoidal rule with 32 × 32 nonuniform grid points [24,45], except that 48 × 48 discrete points are used when Kn = 10.When Ma = 1.2, the molecular velocity space (ξ x ,ξ y ∈ [−6,6]) is represented by 48 × 48 nonuniform grid points.In terms of the spatial discretization, a set of nonuniform meshes are adopted with N x × N y grid points, and the mesh resolution is gradually refined from the cavity center to wall boundaries.The location of a control volume center (x i ,y j ) is generated by in which a is a constant that determines the mesh distribution.The larger a is, the smaller the mesh size near the walls.Here a in the x and y directions are set to be 2.5 and 5.0, respectively.Therefore, for the cases with cavity aspect ratio A of 0.5,1 and 2, N x × N y = 31 × 61, 61 × 61, and 121 × 61 nonuniform meshes are adopted, respectively.The independence of the results with respect to the discretizations of molecular velocity space and spatial space has been confirmed for the given conditions.
Note that DUGKS has a distinguished performance in robustness [35,37].For instance, for the case of Ma = 0.01, Kn = 0.1, St = 2, and A = 1, the change of the amplitude of the average shear stress [Eq.(26)] on the oscillating lid is less than 0.03% when the CFL number varies from 0.01 to 0.8.Therefore, a relatively large CFL number can be used to save the computational time.In all the simulations below, the CFL number η ≈ 0.5 is set to satisfy n t = π , n ∈ Z + .
The primary reason of adopting the DUGKS is that, as opposed to the traditional DVM, the grid size in the DUGKS does not have to be smaller than the mean free path in the near-hydrodynamic regime [40].For example, for the case of Kn = 0.01, A = 1, with St = 2 and 10, the changes in the amplitude of the average shear stress on the driven lid in the simulation adopting 61 (41 and 21) mesh points in one direction are less than 0.4% (1% and 6%), compared with the results of 121 points.With 61 grid points, the average mesh size is about twice the mean free path.
The capacity of DUGKS has been well demonstrated in the previous works.Particularly, the validation for high-Machnumber flow has been confirmed by comparing with the DSMC method [41].Here we only validate the code by comparing the amplitude of the average shear stress along the oscillating lid with the results of the FSM solutions of the linearized Boltzmann equation [17].As can be seen from Table I, the maximum relative difference in the amplitude of the average shear stress between these two methods is less than 1% for a wide range of Strouhal number from 2 to 30, which indicates that the DUGKS for the Shakhov model equation is appropriate to describe the oscillatory rarefied gas flow.

B. Flow characteristics and heat transfer
The velocity profiles on the oscillating lid at various Knudsen and Strouhal numbers are investigated when the cavity aspect ratio is A = 1 and the Strouhal number is St = 2.We are interested in the results that have already reached the periodic "steady-state;" that is, the solution at the next oscillation period will be exactly the same as the previous period.
Figure 2 presents the horizontal velocity on the oscillating lid when t = 0 for Ma = 0.01 and Ma = 1.2, with the Knudsen number from 10 to 0.001.When Ma = 0.01, the flow is in the linear oscillation region, and the flow velocity is symmetric along the central vertical line x = 0.5L; see Fig. 2(a).As the Knudsen number decreases, the maximum horizontal velocity on the lid increases towards the limit U = U 0 ; this is comprehensible because the smaller the Knudsen number, the smaller the slip velocity is.Also, as Kn decreases, the velocity profile becomes more and more flattened across the center of the lid; this is probably due to the fact that, when Kn is large, the Knudsen layer reaches to the center of the cavity, such that the velocity profile changes rapidly.When Ma = 1.2, the flow is in the nonlinear oscillation regime, and the velocity profiles are not symmetrical any more; see Fig. 2(b).Instead, when Kn is small (Kn = 0.001, 0.01, and 0.1), the maximum flow velocity occurs near the top-right corner of the 2D cavity, while when Kn is large (Kn = 1 and 10), the maximum flow velocity occurs close to the cavity center.
The time evolution of the temperature near the right-top corner of the cavity, e.g., at (x,y) = (0.95L, 0.95H ), is plotted in Fig. 3 when Kn = 0.1 and 0.01, Ma = 0.01 and 1.2, and St = 2.The results with other aspect ratios and Strouhal numbers are similar.When Ma = 0.01, it is seen from Figs. 3(a) and 3(c) that the average temperature variation is zero, while when Ma = 1.2, it is seen from Figs. 3(b) and 3(d) that, most of the time, the top-right corner is heated, which may lead to the burning of this corner.Also, the amplitude of temperature variation for Ma = 1.2 is about four orders of magnitude larger than that for Ma = 0.01; that is, the variation is roughly proportional to Ma 2 .With the assistance of the Fourier transform analysis, we note that the temperature evolves with a single frequency the same as the oscillation frequency of the lid when Ma = 0.01, but the second-harmonic frequency emerges when Ma = 1.2, which is a strong signature of the nonlinear interaction.The temperature variation is well-fit by a sinusoidal function as (29) where the coefficients are given in Table II.Third harmonic or even higher harmonic oscillations are possible when the Mach number increases further.
The influence of the amplitude and frequency of the oscillating lid on the heat transfer is also investigated.Figure 4 shows the temperature contours overlaid by heat flux streamlines during the first half period when Ma = 0.01 and 1.2, Kn = 0.1, St = 2, and A = 1.The temperature variation T − T w and flow velocities in the next half period has the same magnitude as in the first half period, but with reversed signs.As shown in Fig. 4(a), when Ma = 0.01 and t = 0, the moving lid has the maximum velocity, and the cold and hot flow fields appear near the top-left and -right corners of the cavity, respectively.The temperature variation satisfies the symmetry relation T (x,y) − T w = T w − T (L − x,y), so that the temperature along the center vertical line x = 0.5L is always equal to the wall temperature.As the oscillation of the lid continues, i.e., the lid velocity decreases, the hot and cold regions move toward each other, and finally swap their positions when the driving velocity reverses; see the top row of Fig. 4. In addition, the temperature variation decreases with the oscillating speed of the lid.When Ma = 1.2, the  29).Note that time t has been normalized by the oscillation period 2π/ω.temperature and heat flux asymmetrically spread inside the cavity, which is contrary to the case of Ma = 0.01 where their distributions are symmetric about the line x = 0.5L.In this nonlinear oscillation case, when t = 0, the temperature near the top-right corner of the cavity rise rapidly, while that in the top-left corner only decreases slightly; see Fig. 4(e).As time goes on, the hot region gradually expands toward the left vertical wall and eventually swap positions with the cold region; see Fig. 4(h).
The velocity and frequency of the oscillating lid change the behavior of heat transfer.It is recognized that, for the low-speed rarefied lid-driven cavity flow, the heat could be transferred from the cold to hot region [41,46]; see Fig. 5(a), where the Strouhal number is zero.When the lid oscillates, we find from Fig. 4(a) that the direction of heat flux for Ma = 0.01 is predominantly from the hot to the cold region at the beginning of a oscillation period, except in the two top corners.When the driving velocity deceases to zero, a heat-transfer circuit emerges near the upper part of the cavity; see Fig. 4(c).When Ma = 1.2, the hot-to-cold heat transfer occurs in almost the entire cavity.To assess the influence of the oscillation frequency on the heat-transfer property, we present the temperature contours and the heat flux streamlines at t = 0, Ma = 0.01 and 1.2, and St = 0, 0.5, and 1 in Fig. 5.As shown, when Ma = 0.01, the direction of heat flux for St = 0 is generally from the cold to hot flow regions, as is the case of St = 0.5.However, when St = 1, two heattransfer circuits appear respectively in the upper and lower parts of the cavity, suggesting a transition of the dominant mechanism in heat transfer.Finally, as observed from Fig. 4(a), the hot-to-cold transfer is already prevailing inside the cavity when St = 2; this situation does not change when St increases further.

C. Damping force on oscillating lid
We first try to investigate the behavior of the damping force at the limit of ω → ∞.In this case, binary collisions are negligible [13] because the high oscillation frequency ω is much larger than the mean collision frequency.As a consequence,  Eq. ( 5) is equivalent to the collisionless Boltzmann equation: Clearly, the absence of the nonlinear collision term means that the flow system is linear; that is, all the flow properties oscillate around their equilibrium values with the same frequency ω of the oscillation lid.Thus, the velocity distribution function can be expressed as where the equation for f is independent of time:  Integrating Eq. ( 32) with respect to the x, one can obtain where f = A 0 f dx/A is the average distribution function.Note that Eq. ( 33) is exactly the same as Eq.(3.3) given by Wu et al. [17] but derived without the assumption of small Mach number.Therefore, for the linear case at high frequency, the right-hand-side of Eq. ( 33) can be neglected.Hence, together with the boundary condition (24) the analytical solution of the amplitude of the average shear stress on the oscillating lid can be obtained, which is equal to 1/ √ π ≈ 0.564.However, for the nonlinear case, the right-hand-side of Eq. ( 33) cannot be neglected due to the asymmetrical distributions on the two vertical walls, which results in the asymmetrical distributions of the velocity distribution function; see Fig. 2(b) for the velocity profile.Unfortunately it is difficult to derive the analytical solution in this situation, but the numerical results below show that the asymptotic value of the damping force when St → ∞ is very close to that of the linear case.
We now numerically solve the Shakhov model at different Knudsen numbers, when St = 2, and Ma = 0.01 and 1.2, respectively.Results of the shear stress along the oscillating lid are plotted in Fig. 6.It is observed that the shear stress on the lid decreases with the Knudsen number for both linear and nonlinear oscillations, because the corresponding larger velocity at lower Kn should exert a smaller shear stress on the lid.For the linear oscillation, distributions of the shear stress are symmetrical along the line x = 0.5L, the same as that of the flow velocity in Fig. 2(a).However, for the nonlinear oscillation, distributions of the shear stress are not symmetrical along the center vertical line anymore, and the shear stress adjacent to the right wall is larger than that close to the left wall.For example, when Kn = 10, the shear stress near the right wall is about four times that near the left wall.
We continue to investigate the variation of the damping force on the oscillating lid as a function of Strouhal number.The results for the linear and nonlinear oscillations at various Knudsen numbers and cavity aspect ratios are plotted in Fig. 7.Note that the case of A = ∞ corresponds to the 1D oscillatory Couette flow, where the damping force always increases monotonically with Strouhal number.However, the damping force changes nonmonotonically in the 2D cavity.For example, when Kn = 0.1 and Ma = 0.01, as the Strouhal number increases, the damping force first increases to a local maximum, then decreases to a local minimum, and finally increases toward the limiting value of 1/ √ π ; see Fig. 7(e).Similar behavior is observed in the nonlinear oscillation case, except that the limiting value of the damping force when St → ∞ for Kn = 10 and 1 is slightly larger than the analytical value 1/ √ π , while for Kn = 0.1 and 0.01, the limiting value of the damping force is smaller than the analytical solution.Generally speaking, the maximum deviation of the limiting value of the damping force when St → ∞ between Ma = 1.2 and 0.01 is less than 7%.
Intuitively, due to the presence of the two vertical walls, the damping force in the 2D oscillatory cavity flow should be larger than that of the 1D oscillating Couette flow.This is indeed the case when the oscillation frequency is small.However, as St increases, the damping force in the 2D case could be even smaller than that of the 1D oscillatory Couette flow.For instance, see Fig. 7(a), when Kn = 10, Ma = 0.01 and A = 2, the damping force in the oscillatory cavity flow is larger than that of the 1D oscillatory Couette flow when St is less than 1, but above this value, the damping force first decreases to a minimum at St = 1.4,and then increases toward the same limits for both the 1D and 2D cases.The underlying physics leading to the minimum in shear stress will be discussed further in the Sec.IV D below.
When Kn = 0.01, the variation of amplitude of the average shear stress with respect to the Strouhal number is more complicated than those at larger Knudsen numbers.This situation has not been studied in Ref. [17] because the adopted iterative numerical method can hardly obtain the steady-state solution due to the strong binary collision.However, the explicit DUGKS can handle this problem easily.It is found from Figs. 7(g) and 7(h) that, when the cavity aspect ratio A is small, the damping force first increases, then decreases, and eventually increases almost linearly, when St increases.However, when A is large, the variation of the damping force with respect to the Strouhal number complicates: there are even two local peaks and valleys in the damping force when A = 2. Similar phenomena are observed for Kn = 0.001 (not shown here).
Finally, the variation of the damping force | Pxy | on the oscillating lid as a function of the Mach number is shown in Fig. 8.As expected, the normalized damping force increases monotonically with the Mach number, which is fitted by a three-order polynomial function, with the coefficients listed in Table III.As the Mach number changes from 0.01 to 1.8, for Kn = 0.1 and 1 at St = 2, the damping force is increased by more than 20%.However, the variation of the damping force for St = 10 is weaker than that of St = 2.This is comprehensible because, at large Knudsen numbers, the collision term in the Shakhov model is negligible so that a linear behavior is expected; see Eq. (30).

D. Scaling law for resonance frequency and aspect ratio
The drop (rise) in the amplitude of the average shear stress in Fig. 7 can be interpreted qualitatively by the gas antiresonance (resonance) inside the cavity [17].For the free molecular flow, the molecules leaving the top lid with the velocities v nearly parallel to the top lid, hitting the right vertical wall, then  being reflected and hitting the left vertical wall, and finally returning to the point from which they left, should have traveled a distance about 2L.Thus, 2L ≈ vδt, (34) where δt is the traveling time.Therefore, if then molecules leaving and hitting the top lid should have the same (opposite) phases.When the velocity distribution functions for molecules leaving and coming back to the oscillating lid are in-phase, the horizontal velocity of the gas near the oscillating lid will be maximum, however, the shear stress defined in Eq. ( 27) will be minimum, because the molecules leaving and coming back to the lid have opposite y-component velocity.The antiresonance and resonance refer to the states where the shear stress exerting on the oscillating lid are minimum and maximum, respectively.Equations ( 34) and ( 35), together with Eq. ( 2), give the resonance and antiresonance Strouhal numbers as and For the linear oscillation, perturbation of the temperature inside the cavity could be neglected, so that the molecular velocity v ≈ v m (= √ 2RT w ).Therefore, the dominant resonance Strouhal number St r ≈ π/2A and the antiresonance Strouhal number St a ≈ π/A can be obtained by setting n = 1 in Eqs.(36) and (37).On the other hand, for the nonlinear case, as shown in Fig. 3, the average temperature inside the cavity is larger than the reference value T w , so the most probable molecular velocity v = √ 2RT > v m .As a consequence, St a and St r should be larger than that of the linear case.
Note that these analyses are for the free molecular flow.As the Knudsen number decreases, binary collisions of molecules become more intense.Therefore, free transport of the gas molecules is more difficult, so the real traveling time is larger than 2L/v, and the obtained St a may be smaller than the theoretical values in Eqs.(36) and (37).
To conclude, the antiresonance (resonance) Strouhal number should be related to the degree of gas nonlinearity and rarefaction.On the other hand, the ability of transferring energy from the oscillating lid is boosted by frequent collisions at a smaller Knudsen number, such that the flow inside the cavity could have sufficient kinetic energy to oscillate against the lid.Therefore, the resonance is difficult to develop for the highly rarefied flow, and it is only visible for the small Knudsen number, see Figs. 7(g) and 7(h).
Figure 9 presents the linear scaling law to describe the dependency of the antiresonance Strouhal number on the reciprocal cavity aspect ratio for both linear and nonlinear cases at Kn = 0.1 and 1.The solutions predicted by the FSM of the linearized Boltzmann equation [17] are also included for comparison.The numerical results for Kn = 10 are not shown because they almost overlap those for Kn = 1.It is found that, in the linear case, our simulation results are in excellent agreement with the fitting functions obtained from the FSM solutions.In addition, as expected, the predicted antiresonance Strouhal number St a in the linear oscillation is overall smaller than that of the nonlinear case.
In addition, as we can see from Figs. 7(g) and 7(h) that for Kn = 0.01, both resonance and antiresonance appear, even though their strengths are quite weak relative to the background value.For the first time, we reveal that the linear scaling law for the resonance (antiresonance) Strouhal number and the inverse cavity aspect ratio at Ma = 0.01 and 1.2 with Kn = 0.01 in   36) and (37).Furthermore, the fitting lines of the resonant Strouhal number for both linear and nonlinear cases almost overlap each other, suggesting that the velocity amplitude of the lid has little effect on the resonance frequency in the early slip regime.

V. SUMMARY
We have investigated the oscillatory rarefied gas flow in a 2D rectangular cavity on the basis of the Shakhov equation.Both the linear and nonlinear oscillations were numerically studied by using the discrete unified gas-kinetic scheme.The effects of the gas rarefaction, oscillation frequency, and aspect ratio of the cavity were investigated.It has been found that the flow properties, including the flow velocity, temperature, shear stress, and heat flux, are asymmetrically distributed inside the cavity for the nonlinear oscillation, which are different from the linear case where these properties are symmetrical along the vertical centerline of the cavity.Interestingly, it was noticed that the hot-to-cold heat transfer could be dominant inside the linear oscillatory rarefied cavity, which is in contrast to the well-known anti-Fourier (cold-to-hot) heat transfer in a low-speed lid-driven rarefied cavity.For the nonlinear oscillatory flow, the hot-to-cold heat transfer occurs in almost the entire cavity, which is independent of the oscillation frequency.It is also noted that double-frequency oscillation of the temperature evolution is possible in nonlinear oscillation, where the temperature of the top-right corner is significantly higher than the wall temperature in almost the whole oscillation period.
We have also investigated the damping force (i.e., the average shear stress exerted on the oscillating lid).In the linear oscillation case, the amplitude of the shear stress is symmetrical along the center vertical line, while in the nonlinear oscillation case it is asymmetrical: when Ma = 1.2, the shear stress near the top-right corner of the cavity could be four times larger than that near the top-left corner.
One of the features of the oscillatory cavity flow is that the damping force exerted on the oscillating lid has local dips and peaks when the oscillation frequency changes.This is due to the antiresonance and resonance of rarefied gas flows, respectively.Linear scaling laws for the antiresonance frequency and the inverse aspect ratio of the cavity are established theoretically for the both linear and nonlinear oscillatory flows from the near hydrodynamic to highly rarefied regimes, which are then verified numerically.In particular, the linear scaling law for resonance frequency is also presented in the early slip regime, where the resonance is easy to develop.
The present study provides a useful guideline to avoid damping damage induced by the nonlinear oscillation for design of MEMS devices.

FIG. 1 .
FIG. 1. Schematic diagram of the oscillatory flow in a rectangular cavity.The coordinate origin is located at the bottom-left corner of the cavity.

FIG. 2 .
FIG. 2. The horizontal velocity on the oscillating lid when ωt/2π = 0, St = 2, A = 1, and (a) Ma = 0.01 and (b) Ma = 1.2.Along the direction of the arrow, the Knudsen number is respectively 10, 1, 0.1, 0.01, and 0.001 for each line.The velocity has been normalized by the velocity amplitude U 0 of the oscillating lid.

FIG. 4 .
FIG. 4. Temperature contours and streamlines of heat flux inside the cavity at different times, when Kn = 0.1, St = 2, A = 1, and the Mach numbers are Ma = 0.01 (top row) and 1.2 (bottom row).

FIG. 5 .
FIG.5.Temperature contours overlaid by the heat flux streamlines inside the cavity for St = 0, 0.5 and 1 (from let to right column), and Ma = 0.01 (top row) and 1.2 (bottom row), with Kn = 0.1 and A = 1, when the horizontal velocity of the lid is maximum along the positive direction (ωt/2π = 0).

FIG. 6 .
FIG.6.The magnitude of shear stress on the oscillating lid when t = 0, St = 2, and (a) Ma = 0.01 and (b) Ma = 1.2.Along the direction of the arrow, Kn is respectively 10, 1, 0.1, 0.01, and 0.001 for each line.Note that the shear stress has been normalized by ρ 0 RT w U 0 /v m .

FIG. 9 .
FIG.9.The antiresonance Strouhal number St a , at which the amplitude of the average shear stress at the lid is minimum, as a linear function of the inverse aspect ratio 1/A for (a) Kn = 0.1 and (b) Kn = 1.The results of the linearized Boltzmann equation solved by the FSM[17] are also included for comparison.
FIG.10.The resonance (antiresonance) Strouhal numbers, at which the amplitude of the average shear stress at the lid is a local maximum (minimum), as a linear function of the inverse aspect ratio 1/A when Kn = 0.01, in both the linear (Ma = 0.01) and nonlinear (Ma = 1.2) oscillations.

Fig. 10 ,
Fig. 10, which is in reasonable agreement with the theoretical St r (St a ) in Eqs.(36) and(37).Furthermore, the fitting lines of the resonant Strouhal number for both linear and nonlinear cases almost overlap each other, suggesting that the velocity amplitude of the lid has little effect on the resonance frequency in the early slip regime.

TABLE I .
[17]arison of the amplitude of the average shear stress on the oscillating lid at different values of the Strouhal number St, when Kn = 0.1, Ma = 0.01, and A = 1.Here the reference solutions are from the linearized Boltzmann equation solved by the FSM[17].

TABLE III .
Fitting the average shear stress on the oscillating lid by a 3 Ma 3 + a 2 Ma 2 + a 1 Ma + a 0 .