Computational Study of Rarefied Gas Flow and Heat Transfer in Lid-driven Cylindrical Cavities

The gas flow characteristics in lid-driven cavities are influenced by several factors, such as cavity geometry, gas properties, and boundary conditions. In this study, the physics of heat and gas flow in cylindrical lid-driven cavities with various cross-sections, including fully or partially rounded edges, is investigated through numerical simulations using the direct simulation Monte Carlo (DSMC) and the discrete unified gas kinetic scheme (DUGKS) methods. The thermal and fluid flow fields are systematically studied for both constant and oscillatory lid velocities, for various degrees of gas rarefaction ranging from the slip to the free-molecular regimes. The impact of expansion cooling and viscous dissipation on the thermal and flow fields, as well as the occurrence of counter-gradient heat transfer (also known as anti-Fourier heat transfer) under non-equilibrium conditions, are explained based on the results obtained from numerical simulations. Furthermore, the influence of the incomplete tangential accommodation coefficient on the thermal and fluid flow fields is discussed. A comparison is made between the thermal and fluid flow fields predicted in cylindrical cavities and those in square-shaped cavities. The present work contributes to the advancement of micro/nano-electromechanical systems (MEMS/NEMS) by providing valuable insights into rarefied gas flow and heat transfer in lid-driven cavities.

of gas behavior in lid-driven cylindrical cavities to design and optimize functional devices and systems that are both efficient and sustainable.
The present work focuses on describing the physics of heat and gas flow in cylindrical lid-driven cavities of various cross-sections, including those with fully or partially rounded edges, utilizing the direct simulation Monte Carlo (DSMC) and the discrete unified gas kinetic scheme (DUGKS) methods.
The effects of both constant and oscillatory lid velocities on the thermal and fluid flow fields are described. The manuscript is organized as follows: Section 2 presents the cylindrical cavity geometries and boundary conditions, Section 3 describes the numerical techniques used, Section 4 presents the results of the thermal and gas flow fields, and finally, the conclusion is provided in Section 5.

Problem description
The present work investigates the thermal and fluid flow behavior of argon, a monoatomic gas, in two different cavities, as depicted in Figure 1. The first cavity, designated as P-Cavity, is a cylindrical cavity with a flat top lid, with the horizontal lid located at y = R·sin(π/4), where R represents the radius of the cavity, as illustrated in Figure 1(a). The second cavity, referred to as C-Cavity (see Figure 1(b)), is a cylindrical cavity with a 90° arc lid. Both cavities have a diameter of 1 µm and are centered at the coordinates O(x, y) = (0, 0 µm). The solid boundaries of the cavity are maintained at a constant temperature of Tw = 300 K. The gas flow within the cavities is induced by the harmonic movement of the lid.
Specifically, the lid oscillates at a constant frequency , and the magnitude of its velocity is defined as follows: where, UM is the maximum amplitude of the oscillating lid velocity magnitude, and t is the time. 5 Thermal and gas flow behavior in cavities with constant lid velocity are also considered in the present work as a special case, where the oscillation frequency is equal to zero. The oscillatory gas flow behavior in cavities can be characterized by the cavity geometry, the Knudsen number and the Strouhal number. The Knudsen number (Kn) and the Strouhal number (St) are defined as follows [40]: where, λ is the molecular mean free path, D is the cavity diameter, KB is the Boltzmann constant, Tw is the wall temperature, m is the molecular mass, and vm is the most probable molecular velocity that is determined as follows [40]: The computational grid configuration used in the simulations is also shown in Figure 1 for different cavity geometries.

Direct simulation Monte Carlo (DSMC)
The DSMC method is a useful computational approach for simulating rarefied gas flows as it can accurately capture the behavior of individual gas molecules, even at low densities. The present work utilizes the direct simulation Monte Carlo (DSMC) method, a probabilistic, particle-based approach for simulating the kinetic behavior of gas particles. This method is used to approximate the solution of the Boltzmann kinetic equation that is defined as follows [41][42][43]: where, f and f1 are the velocity distribution functions, ξ is the molecular velocity vector, r is the position vector, t is the time, is the number density, * , 1 * are the corresponding velocity distribution after collision, cr is the relative velocity vector of pre-and post-collision, c is the collision cross-section, and Ω is the unit solid angle. In the DSMC method, the following splitting scheme is applied to the distribution function at tk to obtain the solution at tk+1 [18]; [ + , ( + ), ( where, is the time step and operators ,ℎ and are the DSMC numerical algorithms approximating the collision and free molecular motion terms in the Boltzmann equation respectively. In this context, the Maxwell-type diffuse-specular condition is applied to model gas-solid surface interactions as collisions between gas molecules and surface material atoms [44]. The tangential accommodation coefficient (αt) was employed to represent the degree of specular-diffuse reflection of the surface, with values ranging from 0 to 1. A fully diffuse surface (αt = 1) is assumed to result in reflected 7 molecules that are in a Maxwell-Boltzmann equilibrium distribution. Conversely, a fully specular boundary (αt = 0) is equivalent to a binary collision of hard spheres, where the particles retain their velocity in the tangent direction but change direction in the normal direction. The probability density function and distribution function are given by [45]: where, v is the velocity vector, nw is the unit vector normal to the surface, α is the accommodation coefficient, ℛ is the gas constant, 3 is the generalized Dirac delta function, nr is the number density, and Tr is the surface temperature. The subscripts i and r indicate the incident and reflected parameters respectively. The effects of tangential accommodation coefficient αt on the thermal and fluid flow fields in lid-driven cylindrical cavities are studied in the present work for the cases with Kn = 1. In the present study, the direct simulation Monte Carlo (DSMC) simulations were conducted utilizing the open-source solver, dsmcFoam+ [46]. The dsmcFoam+ solver has been widely employed to model rarefied gas flow and heat transfer at micro and nano scales [15,[47][48][49][50][51]. The working fluid selected for this investigation is argon, with a molecular mass of m = 6.63×10 −26 kg and a molecular diameter of

Discrete unified gas kinetic scheme (DUGKS)
The discrete unified gas kinetic scheme (DUGKS) is a deterministic method that is efficient and accurate for simulating rarefied gas flows with a wide range of Knudsen numbers, including those with complex geometries, while DSMC is a stochastic method that is particularly well-suited for simulating rarefied gas flows with high Knudsen numbers but can become computationally expensive with large number of simulated molecules. Thus, the choice between DUGKS and DSMC depends on the specific requirements of the problem at hand. DUGKS shares the advantages of the unified gas kinetic scheme (UGKS) and the lattice Boltzmann method [52].
DUGKS is a numerical method for solving the kinetic equation of rarefied gas flows. The method is based on the unified gas kinetic scheme (UGKS) and the Bhatnagar-Gross-Krook (BGK)-Shakhov model, which includes both the Navier-Stokes-Fourier and the Burnett equations as special cases. The 9 DUGKS method discretizes the UGKS model using a finite-volume method and a two-step time-marching scheme. The resulting algorithm is highly efficient and accurate and is able to handle a wide range of rarefied gas flow problems, including those with complex geometry and boundary conditions [23].
Additionally, the DUGKS method has been shown to be more robust and less sensitive to numerical errors than other methods for solving the kinetic equation of rarefied gas flows [53]. Details of the DUGKS method are reported comprehensively elsewhere [19,20,23,53], and thus are not repeated here.
An open-source DUGKS solver, dugksFoam [53], was employed to construct DUGKS simulations in the present work. A grid size of 0.2λ with 28 Gauss-Hermite discrete velocity points was employed in the simulations, and the time-step size Δt was adapted automatically based on the Courant-Friedrich-Lewy number (CFL = ξ·Δt / Δx) to achieve a CFL less than 0.5, where ξ is the molecular velocity magnitude, and Δx is the computational grid size.

Model verification
A comprehensive study was performed to assess the impact of the computational grid size, time , where U is the gas velocity magnitude, γ is the ratio of the specific heats, ℛ is the gas constant, and T is the temperature). As seen in Figure 2(a-c), among the various grid sizes tested (0.5λ, 0.1λ, and 0.02λ), a grid size of 0.1λ was determined to be sufficient to accurately capture the thermal and gas flow characteristics. 10 The time-step size in DSMC simulations should be small enough to decouple the molecular movement from molecular collision. To determine the optimal time-step size (Δt), the mean collision time (Δtc) was calculated as follows [54]: where, λ is the molecular mean free path, KB is the Boltzmann constant, T is the temperature, and m is the molecular mass. The time-step size in the simulations was set to 0.01Δtc, based on the results presented in Figure 2(d-f). The influence of PPC on the numerical predictions is shown in Figure 2(g-i), leading to the selection of a PPC value of 30 for the DSMC simulations. Figure 3 shows the effect of computational grid size on the numerical predictions obtained from the DUGKS simulations. The results showed that a grid size of 0.2λ is adequate for accurately simulating the thermal and gas flow characteristics using the DUGKS model.  14 The selection of the model was based on its performance in specific conditions, with the DUGKS model being employed for low Knudsen number values of 10 −2 and 10 −1 due to its computational efficiency and accuracy in modelling oscillatory flows, as demonstrated in the literature [55][56][57]. The DSMC model was employed for higher Knudsen number values of 10 −1 , 1, and 10 due to its superior performance in accurately modelling gas-surface interactions.

Results and discussion
In  The dynamic viscosity of a gas is generally related to the mean free path of the molecules and their collision frequency. As the mean free path increases with increasing the Knudsen number, the collision frequency decreases, which leads to a decrease in the dynamic viscosity. However, the precise relationship between the Knudsen number and the dynamic viscosity can be complex and depends on factors such as gas composition and temperature [58]. At very low Knudsen numbers (i.e. in the continuum regime), the dynamic viscosity is independent of the Knudsen number and depends only on the gas properties and temperature [17]. As the Knudsen number increases and the gas becomes more rarefied, the dynamic viscosity decreases due to the reduced frequency of molecular collisions. Accordingly, an increase in the Knudsen number results in a decrease in the rate of momentum transfer through the fluid and a subsequent decrease in the shear stress. Figure 6 shows the predicted location of the vortex center as a function of Knudsen number for cavities with different geometries. It appears that the center of the vortex in the P-Cavity moves away from the moving boundary and shifts slightly towards the cavity center with increasing the Knudsen number. Similar behavior has been reported for rarefied gas flow in lid-driven square cavities [5,28,29,56]. As the Knudsen number increases from 10 −2 to 10 −1 , the vortex center in the C-Cavity moves towards the moving boundary and shifts towards the cavity center. However, the vortex center moves   The results indicate that an increase in the Knudsen number results in an increase in the velocity slip at the solid boundaries. However, for both cavity shapes, the slip velocity remains nearly constant when the Knudsen number increases from 1 to 10. When the Knudsen number exceeds 10 −1 , the velocity field exhibits symmetry about the y-axis crossing the cavity center. In lid-driven cavities, the velocity difference between the moving lid and the adjacent gas generates a velocity gradient near the lid, resulting in shear stress. The magnitude of shear stress is directly proportional to the velocity difference between the moving lid and the adjacent gas. Therefore, the magnitude of shear stress is greatest in regions close to the lid and decreases moving away from the lid. However, the magnitude of shear stress is significantly reduced as the Knudsen number increases, which is consistent with the results reported by John et al. [29] for rarefied gas flow in lid-driven square cavities. This can be attributed to the reduced frequency of molecular collisions and the decrease in the rate of momentum transfer through the gas. John et al. [29] for Kn = 1.

19
The behavior of gas molecules in proximity to solid surfaces has a significant impact on the macroscopic behavior of the flow [60,61]. When the Knudsen number is low (Kn ≪ 1), molecular collisions occur frequently, and the gas behaves as a continuum. As a result, viscous forces dominate the molecular interactions with solid surfaces. However, when the Knudsen number increases, molecular collisions occur less frequently, and the molecular interactions with solid surfaces are dominated by the ballistic collisions of gas molecules with the surface. The behavior of gas molecules near solid surfaces is significantly influenced by the tangential accommodation coefficient αt, which is a dimensionless quantity that ranges between 0 (fully specular reflection) and 1 (fully diffusive reflection) [38,[62][63][64]. The value of the tangential accommodation coefficient is generally dependent on the type of gas and the conditions of the wall surface [64][65][66][67]. The tangential accommodation coefficient defines the ratio of the momentum accommodation coefficient for tangential momentum to that for normal momentum and plays a critical role in describing the ability of gas molecules to transfer momentum with a solid surface and determining the slip velocity in rarefied flows. Figure 8 shows the normalized velocity profiles for different cavity shapes and tangential accommodation coefficients at Kn = 1. When the tangential accommodation coefficient αt is close to unity, which corresponds to diffusive reflection, the slip velocity at the moving lid is small. A decrease in the tangential accommodation coefficient results in an increase in the slip velocity at the moving lid, which in turn lead to a decrease in the shear stress near the lid. This decrease in shear stress can significantly affect the flow behavior, such as the development of vortices and the mixing of fluids. Figure 9 shows the change in the locus of the vortex center in cylindrical cavities as a function of the tangential accommodation coefficient. A decrease in the tangential accommodation coefficient results in an increase in the slip velocity at the moving lid. This increase in slip velocity leads to a decrease in the shear stress near 20 the lid, affects the pressure field in the cavity and results in a less efficient momentum transfer [68].
Hence, the vortex center moves downwards and towards the geometric center of the cavity with decreasing the tangential accommodation coefficient. The results suggest that changes in the tangential accommodation coefficient lead to larger variations in the vortex center's location in cylindrical cavities than in square cavities.   [30] for Kn = 1.   Figure 11 shows the time evolution of the velocity field in the P-Cavity for various lid oscillation frequencies. At St = 1, the gas flow in the cavity is in sync with the lid oscillation, and the streamlines are closed. However, as the frequency of the lid oscillation increases, the phase lag relative to the lid also increases. This causes the gas flow far away from the lid in the cavity to become out of phase with the lid, resulting in open streamlines. Figure 12 shows the evolution of the velocity field in the C-Cavity for different Strouhal numbers, and it exhibits a similar behavior to the P-Cavity. The profiles of normalized velocity in cylindrical cavities for different Strouhal numbers are shown in Figure 13. As the lid oscillation frequency increases, the slip velocity at the lid also increases, resulting in a larger magnitude of shear stress in regions close to the lid. The slip velocity at the moving lid decreases by 26.8% in the P-Cavity and by 39.6% in the C-Cavity with an increase in the Strouhal number from 0 to 10. However, the slip velocities at the lid of a square cavity are generally greater than those observed in cylindrical cavities, as shown in Figure 13. The results from Figure 13(a) reveal that the gas horizontal velocity decreases as moving away from the lid. This decrease in velocity is more pronounced at higher Strouhal numbers, leading to an increased velocity gradient and hence an increase 24 in shear stress near the lid at higher lid oscillation frequencies. For the vertical velocity, its magnitude decreases as the Strouhal number increases until it reverses direction at St = 5, corresponding to the downward bending of streamlines. Moreover, the maximum and minimum velocity values move towards both sides of the cavities as the Strouhal number increases.  Figure 14 shows the predicted thermal field in cylindrical cavities with different shapes and Knudsen numbers, driven by a constant lid velocity. A non-uniform temperature distribution is observed in the cavities, with the left side experiencing a temperature decrease and the right side experiencing a temperature increase. When the Knudsen number is low, significant temperature variations occur only in the regions close to the corners of the moving lid, while the temperature remains relatively constant in the substantial part of the cavity. However, as the Knudsen number increases, temperature variations extend throughout the cavity. The results indicate that an increase in the Knudsen number leads to an 25 increase in the temperature difference induced in the cavity. Moreover, the temperature difference in the C-Cavity is found to be larger than that in the P-Cavity.

Constant lid velocity
As the moving lid drags the gas along, it generates a gas flow within the cavity that creates regions of high and low pressure in the cavity. In regions of lower pressure, the gas molecules experience a reduction in kinetic energy and subsequently slow down. This reduction in kinetic energy results in a decrease in the gas temperature, a phenomenon referred to as expansion cooling. The magnitude of the expansion cooling effect in a lid-driven cavity flow depends on the Knudsen number. At low Knudsen numbers (Kn ≪ 1), the gas behaves like a continuum fluid, and the expansion cooling effect is negligible.
However, as the Knudsen number increases, the gas becomes more rarefied, and the expansion cooling effect becomes more significant. Additionally, the C-Cavity experiences a more pronounced expansion cooling effect than the P-Cavity.
Interactions between gas molecules produce frictional forces, resulting in the dissipation of mechanical energy into thermal energy, known as viscous dissipation. This process plays a crucial role in altering the thermal field in rarefied lid-driven cavity flow with isothermal walls. The effect of viscous dissipation increases with increasing Knudsen number because the gas molecules are less likely to equilibrate with each other and transfer energy through collisions. Therefore, the mechanical energy input to the system from the moving lid is more likely to be dissipated into thermal energy through viscous dissipation, rather than being transferred to the gas molecules through collisions. Maximum viscous dissipation occurs where the largest magnitude of velocity gradient exists, which is predicted to be the top right corner, where the maximum shear stress occurs, making it the hottest region in the cavity. The magnitude of the velocity gradient generated in the C-Cavity is larger than that in the P-Cavity, resulting in a greater temperature rise in the former than the latter (see Figure 7(a)). A more notable temperature 26 difference exists in the C-Cavity as compared to the P-Cavity, due to the greater expansion cooling effect and higher temperature increase caused by viscous dissipation in the former relative to the latter. The results presented in Figure 14 demonstrate that anti-Fourier heat transfer, also known as coldto-hot heat transfer, occurs in the rarefied lid-driven cylindrical cavities even in cases where the Knudsen number is low (i.e., Kn = 10 −2 ). This phenomenon is attributed to the non-equilibrium nature of the rarefied gas, which cannot be described accurately by the Fourier heat conduction constitutive law based on the continuum theory [31,48,70]. The effects of rarefication on the heat transfer behavior are accounted for in the gas kinetic theory [45]. In a rarefied fluid, the heat transfer from a hot region to a cold region is governed by the temperature gradient, while the heat transfer from a cold region to a hot region is governed by the gradient of shear stress, which is proportional to the second derivative of velocity with respect to space [70]. The results indicate that the contribution of the shear stress gradient to total heat transfer dominates that of the temperature gradient. Furthermore, an increase in the Knudsen 27 number leads to an increase in the contribution of shear stress gradient to total heat transfer, resulting in predominantly cold-to-hot heat transfer, as depicted by the heat flux stream-traces in Figure 14. The presence of sharp corners in the P-Cavity causes a sudden change in the velocity field close to the corners, leading to an increase in the second derivative of velocity in this region. As a result, anti-Fourier heat transfer is more significant in the P-Cavity compared to the C-Cavity.
The predicted thermal field in cylindrical cavities with different shapes, employing various values for the tangential accommodation coefficients, is shown in Figure 15. The gas flow in the cavities is driven by a constant lid velocity and the Knudsen number is set to 1. The results reveal that changes in the tangential accommodation coefficient significantly affect the temperature distribution and the direction of heat flux in the cavity. A decrease in the tangential accommodation coefficient results in a reduction in the maximum temperature in the cavity and an increase in the minimum temperature, leading to a more uniform temperature distribution. Moreover, the heat flux stream-traces indicate that the heat flows predominantly from the moving lid towards the stationary bottom wall, with a decrease in the value of the tangential accommodation coefficient. This observation is attributed to the reduced effect of expansion cooling in thermal energy transfer in the cavity. A decrease in the tangential accommodation coefficient results in an increase in the number of molecules contributing to elastic energy exchange with the walls, leading to less wall shear stress and lower viscous heat dissipation [30].
28 Figure 15. Contours of the temperature overlayed by heat flux stream-traces for different cavity shapes and tangential accommodation coefficients (αt). The velocity magnitude of the cavity lid Ulid is constant, and the Knudsen number is equal to 1.

Oscillatory lid velocity
The contours of the temperature overlayed with heat flux stream-traces in cylindrical cavities are shown in Figure 16 for different Strouhal numbers at t/ts = 0. For these cases, the Knudsen number is equal to 10 −1 . The results indicate that an increase in the lid oscillation frequency leads to an increase in the temperature difference induced in the cavity. In both cavity geometries, cold-to-hot heat transfer predominates when the Strouhal number equals 1. Interestingly, the heat flux stream-traces in the C-Cavity with St = 1 reveal the existence of a vortex heat flow structure near the cavity center. When the Strouhal number is increased to 2, it appears that the hot-to-cold heat transfer dominates in significant parts of the cavity, except in the vicinity of the corners of the lid. As the Strouhal number increases, the heat flow structure within the cavity becomes more intricate. Nevertheless, in substantial regions of the cavity, the hot-to-cold heat transfer continues to be dominant, except for the regions located proximate to the lid. The dominance of hot-to-cold heat transfer, particularly in the central region of the cavity, can be 29 attributed to the increased viscous heat dissipation and, consequently, the elevated temperature gradient in the cavity. In regions proximate to the lid, the contribution of shear stress gradient to total heat transfer supersedes that of the temperature gradient, leading to the domination of cold-to-hot heat transfer near the cavity lid.

Conclusions and future directions
The rarefied gas flow and heat transfer in cylindrical lid-driven cavities with various cross-sections were studied numerically using both the direct simulation Monte Carlo (DSMC) and the discrete unified gas kinetic scheme (DUGKS) methods for a wide range of Knudsen numbers. The present study focused on understanding the effects of cavity geometry, degree of gas rarefaction, and boundary conditions on the thermal and fluid flow fields for constant and oscillatory lid velocities. The results were also compared with those obtained in square-shaped cavities.
The results obtained from the simulations suggest that the geometry of the cavity plays a significant 30 role in influencing the characteristics of the thermal and gas flow fields. Specifically, the average gas velocity and induced temperature difference in cylindrical cavities with fully rounded edges (C-Cavity) are generally greater than that in cylindrical cavities with partially rounded edges (P-Cavity). Furthermore, it was observed that the expansion cooling and viscous dissipation are more pronounced in the C-Cavity compared to that in the P-Cavity. The results also revealed that anti-Fourier heat transfer, also known as cold-to-hot heat transfer, occurs in rarefied lid-driven cylindrical cavities with constant lid velocity, even when the Knudsen number is low (i.e., Kn = 10 −2 ). Moreover, anti-Fourier heat transfer is more significant in the P-Cavity compared to the C-Cavity.
The results indicate that changes in the tangential accommodation coefficient significantly affect the temperature distribution and the direction of heat flux in the cavity. A decrease in the tangential accommodation coefficient leads to a reduction in the maximum temperature in the cavity and an increase in the minimum temperature, resulting in a more uniform temperature distribution.
It was found that lid oscillations can reverse the heat flow direction from cold-to-hot to hot-to-cold, particularly in the central regions of the cavity. The outcomes of this study provide valuable insights into the behavior of rarefied gas flow and heat transfer in lid-driven cylindrical cavities with various cross-sections and can inform the design and optimization of relevant systems and devices.
To improve the computational efficiency of simulating rarefied gas flows in complex geometries, future research could investigate the development of hybrid numerical methods that combine the strengths of different approaches. Since the present work is limited to two-dimensional cylindrical cavities, further studies could explore three-dimensional rarefied gas flow instabilities in lid-driven cavities with complex geometries. Experimental investigations could also be conducted to verify the accuracy of the numerical simulations and shed light on the practical implications of the results for relevant 31 engineering applications.