Pressure-driven Nitrogen Flow in Divergent Microchannels with Isothermal Walls

Gas flow and heat transfer in confined geometries at micro and nano scales differ considerably from those at macro-scales, mainly due to nonequilibrium effects such as velocity slip and temperature jump. The nonequilibrium effects enhance with a decrease in the characteristic length-scale of the fluid flow or the gas density, leading to the failure of the standard Navier-Stokes-Fourier (NSF) equations in predicting thermal and fluid flow fields. The direct simulation Monte-Carlo (DSMC) method is employed in the present work to investigate pressure-driven nitrogen flow in divergent microchannels with various divergence angles and isothermal walls. The thermal fields obtained from numerical simulations are analysed for different inlet-to-outlet pressure ratios (1.5 $\leq \Pi \leq$ 2.5), tangential momentum accommodation coefficients and Knudsen numbers (0.05 $\leq \mathrm{Kn} \leq$ 12.5), covering slip to free-molecular rarefaction regimes. The thermal field in the microchannel is predicted, heat-lines are visualised, and the physics of heat transfer in the microchannel is discussed. Due to the rarefaction effects, the direction of heat flow is largely opposite to that of the mass flow. However, the interplay between thermal and pressure gradients, which are affected by geometrical configurations of the microchannel and applied boundary conditions, determines the net heat flow direction. Additionally, the occurrence of thermal separation and cold-to-hot heat transfer (also known as anti-Fourier heat transfer) in divergent microchannels is explained.


Introduction
Gas flow in micro-and Nano-channels with non-uniform cross-sections offers opportunities to develop small devices with novel applications (see for example, [1,2,3,4,5,6]). Understanding heat and fluid flow in microchannels is essential to engineer novel micro-electromechanical systems (MEMS) [7,8,9]. However, this is a challenging task since the rate of collisions between the gas molecules and solid walls reduces with decreasing the characteristic length scale of the fluid flow or the gas density, affecting the random movement of molecules, which is commonly known as nonequilibrium effects [10,11]. Nonequilibrium molecular transport processes cause velocity slip and temperature jump at walls and influence thermal and fluid flow fields [12,13]. The Knudsen number (Kn), which is the ratio of the molecular mean free path λ to a characteristic length scale L, is often employed as an indicator of deviation from the equilibrium condition. The higher the Knudsen number, the higher the deviation from the equilibrium condition. Four different rarefaction regimes have been defined based on the Knudsen number [14], and are commonly named as free molecular (Kn > 10), transition (10 −1 ≤ Kn ≤ 10), slip (10 −3 ≤ Kn ≤ 10 −1 ) and continuum (Kn ≤ 10 −3 ) regimes.
Micro-devices usually operate in the slip and transition rarefaction regimes at the standard pressure and temperature [12]. It is widely acknowledged that the standard Navier-Stokes-Fourier (NSF) equations fail to predict thermal and fluid flow fields accurately when the nonequilibrium effects are significant [15,16,17]. The Boltzmann equation governs fluid flow at micro-scales in the entire range of Knudsen number [18]. However, direct numerical solution of the Boltzmann equation is generally expensive with available computational capabilities and are limited to simple applications because of the high dimensionality of the Boltzmann equation and the complexity of the collision integral. Moreover, the application of the Boltzmann equation is limited to dilute gasses due to the assumption of binary intermolecular collisions [19]. The direct simulation Monte-Carlo (DSMC) method is a numerical technique that has been developed based on the kinetic theory to approximate the solution of the Boltzmann equation [18].
Zheng et al. [20] compared the numerical predictions obtained from DSMC simulations with the solutions of the NSF equations for a Poiseuille gas flow between two parallel plates, and reported quantitative and qualitative differences between the predictions. Deviations between numerical predictions obtained using the DSMC method and the NSF equations increases with increasing the Knudsen number [20,21,22,23]. Using the regularised 13-moment equations, Torrilhon and Struchtrup [24] achieved a better agreement with the DSMC results up to Kn ≈ 1 in comparison with the NSF-based solutions. Varoutis et al. [25] studied Poiseuille gas flow in channels with finite lengths using the DSMC method, and stated that the nonequilibrium effects are the most significant factor affecting rarefied gas flow in channels.
The majority of published literature on subsonic internal gas flows at micro-scales are limited primarily to flow passages with uniform cross-sections [26,27,28]. Moreover, previous studies have mostly focused on characterising pressure loss, mass flow rate, and velocity field [29]. There is a crucial lack of detailed investigations concerning the thermal field analysis of gas flow in micro-and Nano-channels with non-uniform cross-sections in which nonequilibrium effects play an important role [30]. Using the NSF equations with the first-order slip boundary condition, Varade et al. calculated temperature variations in divergent [31] and convergent [32] microchannels for both continuum and slip flow regimes. Using the same approach, Hemadri et al. [33,34] calculated temperature distribution in convergent microchannels in the slip and early-transition rarefaction regimes. Milićev and Stevanović [35] proposed an analytical solution to the one-dimensional NSF equations to describe steady-state pressure-driven isothermal gas flow in microchannels with variable cross-sections. Ohwada et al. [36] stated that heat and mass flow are in opposite directions for Poiseuille flow between two parallel plates for Kn > 10 −1 while they are in the same direction adjacent to the walls at low Knudsen numbers. It should be noted that only the tangential heat flow profiles were considered in their study while both normal and tangential heat fluxes contribute to the net heat flow. John et al. [37] numerically studied the influence of pressure ratio and surface accommodation coefficient on the thermal field of a Poiseuille gas flow between two parallel plates. The DSMC method has been recently employed to investigate temperature variations in divergent microchannels [38,39,40]. Characterising the thermal field in microchannels with non-uniform cross-sections is emerging for micro-and Nano-system engineering and further investigations are essential.
Despite the extensive interest in studying fluid flow at the micro-and Nano-scales, numerical investigations of gas flow under nonequilibrium conditions are deficient in thermal analysis. Previous studies on thermal field analysis in microfluidics systems are mostly limited to the channels with uniform cross-sections (see for instance, [41,42,37,43,44,45,46,47,48]). Therefore, there is an indispensable need for detailed investigations on heat and fluid flow at the micro-and Nano-scales under nonequilibrium conditions. The goal of the present study is to enhance our understanding of heat and fluid flow in divergent microchannels, which are applicable in micro-devices such as micro-pumps, micro-actuators and micro-thrusters. Numerical simulations based on the DSMC method are performed to explore the effects of divergence angle, inlet-to-outlet pressure ratio, rarefaction, and tangential momentum accommodation coefficient on thermal and fluid flow fields. Heat and fluid flow under nonequilibrium conditions are considerably different from its equilibrium counterpart, providing the motivation to conduct the present study.

Problem description
A probabilistic numerical approach based on the direct simulation Monte-Carlo method is employed to study steady-state, pressure-driven, Poiseuille gaseous nitrogen flow in divergent microchannels shown schematically in Figure 1. The divergence angle of the channel (φ) is defined as follows [39]: where, H i and H o are the height of the channel at the inlet and the outlet, respectively. The length of the channel (L) is twenty times larger than its inlet height (i.e. L = 20 × H i = 8 × 10 −6 m). The problem is assumed to be two-dimensional, which is a valid assumption for channels with considerably large depths [43,49]. The molecular properties of gaseous nitrogen are summarised in Table 1, where d p is molecular diameter, m p molecular mass, ω the viscosity-temperature index, T ref the reference temperature in the viscosity-temperature relation, and µ 0 the viscosity at the reference temperature. Figure 1: Schematic of the divergent microchannel considered in the present study. The hatched region was considered for calculations.

Fluid flow direction
Due to the symmetric configuration of the channel and the flow field, one-half of the physical domain (the hatched region in Figure 1) is considered for calculations. Heat and fluid flow in the channel are described in a Cartesian coordinate system. The temperature of the solid walls (T w ) and the gas at the inlet (T i ) is constant equal to 300 K. The Knudsen number is defined based on the inlet height, i.e. Kn = λ/H i . The inlet-to-outlet pressure ratio (Π) is defined as p i /p o and the gas pressure at the channel inlet (p i ) is calculated based on the Knudsen number and the ideal gas law. The effects of tangential momentum accommodation coefficient (α) on thermal and fluid flow fields are also studied in the present work. The tangential momentum accommodation coefficient ranges between zero (specular reflection) and one (fully-diffused reflections) depending on internal gas molecule structure, the molecular mass of the gas, and surface material [50,51].

Numerical procedure
An extended version of the dsmcFoam solver [52], which has been developed within the framework of an open-source solver (OpenFOAM [53]), was employed in the present study. This solver has been meticulously validated for a variety of benchmark cases [54,55,52] including subsonic Poiseuille micro-flows. The variable hard sphere (VHS) model was selected to treat the binary inter-molecular collisions. The energy exchange between both translational and rotational modes was allowed. The no-time-counter (NTC) [18] collision sampling model was chosen. Additionally, the transient adaptive sub-celling (TAS) scheme [56] was utilised to divide collision cells every time-step according to the number of particles inside each collision cell to have at least two simulator particles in each subcell. This results in the selection of collision partners that are located within a relative distance smaller than the local molecular mean free path.
Oran et al. [57] suggested a cell size of λ/3 to capture the gradients in flow fields. In the case of low-speed micro-flows, larger grids might be employed in the stream-wise direction due to insignificant flow gradients [58,59,60,61]. After a grid independence test, a grid with cell sizes of λ/4 and λ/2 in the normal and the stream-wise directions, respectively, was selected. To capture the flow gradients for the cases with Kn > 2.5 × 10 −1 , the grid size was the same as that for the cases with Kn = 2.5 × 10 −1 . All simulations were initialised with 50 particles per cell to minimise statistical relations between particles [62] and to have at least 15 particles in the computational cells adjacent to the channel's outlet. The time-step was selected sufficiently smaller than the local mean collision time so particles stay inside a cell for multiple time-steps. The results of the grid, time-step size and particle-per-cell independence tests are presented in Figure 3. To return the relevant macroscopic data, the results were sampled over 3 × 10 6 time-steps after achieving the steady-state condition. The sample size was sufficiently large to let the solution procedure continues even after equating the inlet and outlet mass flow rates; leading to a converged solution with negligible (less than 1%) statistical noises [63,64,65]. Each simulation was carried out in parallel on four cores of an Intel Core i7-3520M processor, and took roughly 200 h to complete.

Model validation and verification
The Poiseuille argon flow in a microchannel with a uniform cross-section was considered to verify the reliability of the present model. The microchannel length is 15 µm and its height is 1 µm. Gaseous argon enters the channel with a constant temperature of 300 K, and the inlet-to-outlet pressure ratio Π was set to 3. The channel walls were assumed to be fully diffuse (i.e. α = 1), and their temperature was set to 300 K. Four different cases with different Knudsen numbers are studied, which covers the slip and transition rarefaction regimes. For each case about 1.5 × 10 6 DSMC simulator particles were utilised. The predicted pressure and Mach number (Ma) distributions along the microchannel centreline are shown in Figure 2. The Mach number is defined as follows: where, V is the fluid velocity vector, γ the ratio of specific heat of nitrogen at a constant pressure to its specific heat at a constant volume, R the specific gas constant, and T the temperature. The results are in reasonable agreement with the data reported by White et al. [65]. Further comparison of the results obtained from the present model with experimental, analytical and numerical data was performed by the authors for gaseous nitrogen flow in divergent microchannels, and can be found in [38,39,31]. Figure 3 shows the influence of grid size, number of simulator particles per cell (PPC), and time-step size on the results obtained from the present DSMC simulations. For this study, nitrogen flow in a microchannel with uniform cross-section (φ = 1), Kn i = 0.1 (transition rarefaction regime), and the inlet-to-outlet pressure ratio Π of 2.5 was considered. To determine an appropriate timestep size ∆t within which the molecular movement and collision are distinguishable, the following correlations were employed: where, ∆t c is the mean collision time, ∆t t the mean transit time and κ b the Boltzmann constant. The time-step size ∆t was defined as a fraction of the minimum value of ∆t c and ∆t t (i.e. ∆t = ξ · min (∆t c , ∆t t ) and ξ ≤ 1). For the range of parameters studied in the present work, the results are practically insensitive to the grid size, number of simulator particles per cell, and time-step size. x

Results and Discussion
The influence of the microchannel divergence angle on variations of temperature and Mach number along the microchannel centreline is shown in Figure 4. For φ < 2.5, the gas temperature decreases along the channel because of the augmentation of the gas velocity towards the channel outlet. The increase in kinetic energy of the gas results in a reduction of its internal energy and thus its temperature. The gas velocity remains almost constant along the channel with φ = 2.5 and the gas temperature hardly changes along the channel, except in regions close to the inlet and outlet boundaries. Further increase in the divergence angle φ results in a decrease of gas velocity along the channel, which is attributed to the increased cross-sectional area of the channel as well as the increased shear stress at the channel walls [39,31]. The kinetic energy of the gas converts into internal energy as the fluid velocity decreases towards the channel outlet for the cases with φ > 2.5, leading to an increase in gas temperature. Additionally, the bulk gas temperature reduces with an increase in divergence angle of the channel because of the increased bulk gas velocity and the reduced rate of molecular collision. Because of the development of a viscous boundary layer close to the channel inlet [39,65], the gas flow accelerates in the entrance region; however, the flow decelerates after passing this region moving towards the channel outlet to conform to the pressure prescribed at the outlet. Close to the channel outlet (x/L > 0.9), the gas temperature drops that is attributed to the rapid gas expansion in that region, which is known as 'the expansion cooling phenomenon' [66,20]. A similar behaviour close to the channel outlet has been reported by others [28,37,45,43,44] for Poiseuille micro-flows in microchannels with uniform cross-sections. The overall heat flow direction, for the problem considered in the present work, is determined by both the thermal and pressure gradients in the microchannel [67]. Heat flow direction and contours of dimensionless temperature (T /T i ) obtained from the DSMC simulations are shown in Figure 5 for microchannels with different divergence angles at Kn i = 0.05 and Π = 2.5. Heat flows generally towards the channel inlet except in regions adjacent to the walls in the Knudsen layer, where the heat flow is dominated by viscous heating. For the cases with φ > 2.5, the amount of heat induced by viscous dissipation decreases along the microchannel due to the reduction of wall shear stress and slip velocity [39], decreasing the influence of viscous dissipation on net heat flow direction in regions close the walls. In contrast, the amount of heat generated by viscous dissipation increases as the gas moves along the channels with φ ≤ 2.5, increasing its influence on net heat flow direction in regions adjacent to the walls. Accordingly, it can be argued that the heat flow is dominated by the pressure gradient (known as a rarefaction effect [67]) in microchannels with a divergence angle φ ≤ 2.5. Due to the rarefaction effects, heat flow from cold to hot regions is observed, which is not predictable using the NSF equations because of the neglect of high-order rarefaction terms [68]. Heat flow in the central region of the channel is dominated by the mass flow driven by the pressure gradient. However, the magnitude of thermal gradient increases with an increase in divergence angle φ, affecting the net heat flow direction. For the cases with φ > 2.5, a thermal separation occurs at certain location in the channel beyond which heat flows towards the channel outlet. Thermal separation for pressure-driven gas flow in microchannels of variable cross-sections is attributed to the enhanced contribution of thermal gradient to the total heat transfer [69]. The position of thermal separation point moves towards the inlet with an increase in the divergence angle of the channel. It should be noted that no hydrodynamic separation was observed for the cases considered in the present study [39,38]. These observations can also be explained through the asymptotic analysis of the Boltzmann equation [67] for small Knudsen numbers (i.e. slip and early-transition regimes). However, this type of analysis, which is based on the continuum theory, fails to predict heat transfer accurately at high Knudsen numbers [38].  Numerical predictions of the gas temperature and Mach number along the channel centreline are presented in Figure 6 for different inlet Knudsen numbers Kn i at φ = 7 and Π = 2.5. The bulk gas velocity along the channel and variations in the gas velocity decreases with an increase in the Knudsen number. However, this reduction in the velocity magnitude and its variations becomes insignificant for Kn i ≥ 1. Variation of gas velocity along the channel and the molecular collision frequency decreases with an increase in Knudsen number, limiting the variation of gas temperature along the channel. The magnitude of temperature drop close to the channel outlet also decreases with an increase in the inlet Knudsen number. Moreover, the predicted gas temperatures close to the channel inlet are slightly higher than the inlet gas temperature T i , which is because of the augmentation of temperature jump with increasing the Knudsen number [37].  Figure 6: The effect of rarefaction on the profiles of (a) gas temperature and (b) Mach number along the microchannel centreline (φ = 7.00 and Π = 2.5). The inlet temperature (Ti) is used to non-dimensionalise the predicted temperatures.
The influence of the Knudsen number on the thermal field and the heat flow direction in a divergent channel with φ = 7 and Π = 2.5 is shown in Figure 7. The temperature gradient in the channel decreases with an increase in Knudsen number, decreasing the contribution of the Fourier heat transfer in net heat flow. According to Figure 7 (c) and (d), at sufficiently large Kn i the heat flow in the channel is merely towards the channel inlet that demonstrates the domination of the pressure gradient in the net heat flow direction, intensifying the counter-gradient cold-to-hot heat transfer (also known as anti-Fourier heat transfer). Moreover, the total heat flow rate decreases with increasing the Knudsen number. The distance between the thermal separation point and the channel outlet increases with a decrease in Kn i . The bulk gas velocity and temperature in the channel decrease with a decrease in the inlet-tooutlet pressure ratio Π, as shown in Figure 8. Increasing the inlet-to-outlet pressure ratio increases the gas kinetic energy and decreases its internal energy. The thermal field as well as the heat lines are shown in Figure 9 for channels with a divergence angle of φ = 7 and Knudsen number of 0.05. The magnitude of thermal gradient close to the channel outlet increases with decreasing the pressure ratio Π, which enhances the contribution of the thermal gradient to the net heat flow; this affects the net heat flow direction and leads to the occurrence of thermal separation farther away from the channel outlet. The heat generated in the Knudsen layer due to viscous dissipation diminishes with reducing the inlet-to-outlet pressure ratio Π, decreasing the contribution of thermal gradients to the net heat flow close to the channel walls.  The influence of the tangential accommodation coefficient (α) on the profiles of temperature and velocity along divergent channels is shown in Figure 10 for Π = 2.5 and Kn i = 0.05. For all the cases considered in the present work, decreasing the tangential accommodation coefficient weakens the interactions between the gas molecules and the channel walls, leading to an increase in the gas velocity. A similar behaviour has also been observed both experimentally and numerically for rarefied gas flows in microchannels with uniform cross-sections [37,51,50,70]. As expected, the augmentation of gas velocity in the channel results in an increase in the kinetic energy and a decrease in internal energy and thus the gas temperature. The results presented in Figure 10 show that the gas reaches supersonic velocities when α = 0 (i.e. specular reflection); this happens for the case with φ = 7.0 only when α = 0.5 and does not happen when α = 1 (i.e. fully-diffused reflection).  The results presented in Figure 10 show that heat and fluid flow patterns in the channel are sensitive to the value of the tangential accommodation coefficient α, particularly for low values of the Knudsen numbers. The number of gas molecules contributing to the adiabatic exchange of energy increases with decreasing the value of the tangential accommodation coefficient α, increasing the contribution of thermal gradients to the net heat flow. Reducing the value of the tangential accommodation coefficient to α = 0 as an extreme condition, it is observed that the heat flow converges to a sink point, where pressure gradients balance thermal gradients. Heat lines are directed towards the channel outlet in regions where thermal gradients dominate the net heat flow.

Conclusions
Two-dimensional numerical simulations based on the direct simulation Monte-Carlo (DSMC) method were performed to predict thermal and flow fields of pressure-driven nitrogen flow in divergent microchannels. A wide range of Knudsen numbers was studied covering the slip to the free molecular rarefaction regimes. Moreover, the influence of microchannel divergence angle, inlet-to-outlet pressure ratio and incomplete surface accommodation were investigated. The following conclusions are drawn based on the present study.
the influence of thermal gradient on net heat flow direction enhances with an increase in the divergence angle of the microchannel and can lead to the occurrence of thermal separation even when hydrodynamic separation is absent. The contribution of the pressure gradient in determining the heat flow direction enhances with an increase in the Knudsen number. At sufficiently large Knudsen numbers, the heat flow is merely towards the channel inlet. Cold-to-hot heat transfer (also known as anti-Fourier heat transfer) was observed in the microchannel as a result of nonequilibrium effects. Viscous dissipation in the Knudsen layer diminishes with a decrease in applied pressure ratio, which results in the domination of the pressure gradient influence on heat flow close to the channel walls.
Further studies may focus on realising the effects of the molecular weight and structure, or different gas mixtures on thermal and fluid flow fields inside micro-and Nano-channels with with complex geometries.