Modelling Thermally Induced Non-Equilibrium Gas Flows by Coupling Kinetic and Extended Thermodynamic Methods

Thermally induced non-equilibrium gas flows have been simulated in the present study by coupling kinetic and extended thermodynamic methods. Three different types of thermally induced gas flows, including temperature-discontinuity- and temperature-gradient-induced flows and radiometric flow, have been explored in the transition regime. The temperature-discontinuity-induced flow case has shown that as the Knudsen number increases, the regularised 26 (R26) moment equation system will gradually loss its accuracy and validation. A coupling macro- and microscopic approach is employed to overcome these problems. The R26 moment equations are used at the macroscopic level for the bulk flow region, while the kinetic equation associated with the discrete velocity method (DVM) is applied to describe the gas close to the wall at the microscopic level, which yields a hybrid DVM/R26 approach. The numerical results have shown that the hybrid DVM/R26 method can be faithfully used for the thermally induced non-equilibrium flows. The proposed scheme not only improves the accuracy of the results in comparison with the R26 equations, but also extends their capability with a wider range of Knudsen numbers. In addition, the hybrid scheme is able to reduce the computational memory and time cost compared to the DVM.


Introduction
The advent of micro-electro-mechanical systems (MEMS) and the associated fabrication technologies has inspired a renewed impetus in understanding thermally driven flows [1][2][3][4]. Typical examples include thermal transpiration [5,6], radiometric flow [7,8] and Knudsen pumps [5,9]. Since radiometers could also work as Knudsen pumps, transporting gas without any moving parts, they can not only be developed as radiometric actuators for spacecraft attitude control system [10], but also be used for micro-scale gas chromatography and gas separation applications [11,12]. Due to their important applications in industry, many experimental and numerical investigations of thermally induced flow have been carried out [13][14][15].
Numerically, the behaviour of thermally induced non-equilibrium gas flows can be described and modelled from either a microscopic or macroscopic point of view. The Boltzmann equation is the fundamental model for non-equilibrium flows at the microscopic level, which uses a molecular velocity distribution function (VDF) to describe the system state. From a historical perspective, two major categories of approaches have been developed to solve the Boltzmann equation. One is a stochastic approach, such as the direct simulation Monte Carlo (DSMC) method developed by Bird [16], and the in which, t and x i = (x, y, z) are the temporal and spatial coordinates, respectively, and any subscript i, j, k represents the usual summation convention. The pressure p is related to the temperature T by the ideal gas law p = ρRT, where R is the gas constant. However, the stress tensor, σ ij , and the heat flux vector, q i , in the set of Equation (1) are unknown. The classic way to close this set of equations is through a CE expansion of the molecular distribution function in terms of Knudsen, Kn, around the Maxwellian [21] to obtain the constitutive relationships for σ ij and q i . In the method of moments, Grad [23] derived their governing equations from the Boltzmann equation as: in which, µ is the dynamic viscosity of the gas. The collision constants, A σ and A q , are determined by the molecular collision model. The high moments, m ijk , R ij and ∆, represent the difference between the true value of the higher moments and their corresponding approximation with the truncated distribution function, f G , at the third order in Hermite polynomials. In Grad's original method [23], such deviations were omitted, so that m ijk = R ij = ∆ = 0 and the set of Equations (1) and (2) are well known as Grad's 13 moment equations (G13). Struchtrup and Torrilhon [24] and Struchtrup [31] regularised the G13 moment equations by applying a CE-like expansion and an order-of-magnitude approach and obtained the algebraic constitutive expressions of m ijk , R ij and ∆ in terms of the derivatives of lower order moments. The regularised G13 moment equations are denoted as the R13 moment equations. Although the R13 moment equations improve the performance of the G13 moment equations significantly, they cannot provide sufficient accurate description of the Knudsen layer [30,32]. To remedy the deficiency of the R13 equations, the governing equations of the high-order moment quantities m ijk , R ij , ∆ that can be derived from the Boltzmann equation are employed in the present study. They are [28,29]: where the non-linear source terms M ijk , ij , ℵ are listed in Appendix A. Similarly, the higher-order moments, φ ijkl , ψ ijk and Ω i in the set of Equation (3) represent the difference between the true value of the higher moments and their corresponding approximation with the truncated distribution function, f G , at the fourth order in Hermite polynomials. A CE-like expansion was employed to obtain the following constitutive relationships [28]: The closed set of Equations (4) is known as the R26 moment equations. The values of the collision constants, , A Ω2 and A Ω3 depend on the molecular collision model adopted and represent the relaxation time-scale for each moment. They are given in Table 1 for the case of the Shakhov model [33] as employed in the present study. The Prandtl number, Pr, is set to be 2/3 for monoatomic gases.
To apply the extended thermodynamic equations to flows in confined geometries, appropriate wall boundary conditions are required to determine a unique solution. Macroscopic wall boundary conditions for confined flows were obtained based on Maxwell's kinetic wall boundary condition and a fifth-order approximation of the VDF in Hermite polynomials and they are listed in [28]. However, as a truncated VDF is used in the derivation of the set of boundary conditions at a wall where the state of non-equilibrium is strong, the accuracy near the wall is reduced. As a result, it hampers the capability of the moment method. To increase the accuracy of the solution while maintaining a low computational cost, a hybrid algorithm [34], which couples the thermodynamic equations with the kinetic equation in the wall boundary layer, is adopted in the present study.

Kinetic Equation and the Shakhov Model
From the microscopic view of point, the behaviour of a gas can be described by the kinetic equation: where f = f (t, x i , C i ) is the VDF of gas molecules and C i is the molecular velocity. The mean relaxation time, τ, is evaluated from The Shakhov model [33] is adopted so that the equilibrium VDF, f eq , is given by Maxwell's kinetic wall boundary condition [27] is used in association with the kinetic Equation (5). It states that a fraction α of gas molecules undergoes diffusive reflection with a Maxwellian distribution f w M at the temperature of the wall T w while the remaining fraction (1 − α) will be reflected specularly. In a frame in which the coordinates are attached to the wall, with n i the normal unit vector of the wall pointing towards the gas such that all molecules with C i n i < 0 are incident upon the wall and molecules with C i n i ≥ 0 are emitted by the wall, Maxwell wall boundary condition can be expressed by with in which, u w i is the wall velocity and ρ w is the density of the thermalised molecules determined to ensure that no molecules accumulate on the wall.
The modelled Boltzmann equation can be solved alone for the whole flow domain and the solution can be served as the benchmark data. In the hybrid algorithm employed in the present study,

Hybrid Algorithm of Coupling the Moment Equations and the Kinetic Equation (HYBR)
A hybrid algorithm [34] is used to balance the efficiency and the accuracy of a solution, which couples the moment equations and the kinetic equation. Taking a thermally induced cavity flow as an example, we can divide the computational domain into two sub-domains, as shown in Figure 1. In the near-wall subdomain where the non-equilibrium effects are strong, the kinetic equation in association with the Maxwell's diffusive wall boundary condition is applied, so that accurate description of the flow field can be achieved. The kinetic Equation (5) is solved by the discrete velocity method (DVM) [17,35]. The thickness of the kinetic layer is represented by l. The second-order upwind is used to discretise the spatial derivative in equation. An iteration method is employed as detailed in [35]. As the gas is away from the wall, the R26 moment equations are employed to improve the efficiency. The set of the R26 moment Equations (1)-(3) is solved by a pressure-based numerical algorithm [36] for weakly compressible and low-speed flows. It has been successfully applied in the study of pressure-driven Poiseuille flow [28,37], thermal transpiration flow [5] and gas flows in porous media [38,39]. Once the distribution function, f , is obtained from Equation (19), its moments with respect to the molecular velocity, , C can be determined. For example, the density, ρ , and the momentum, For convenience, the intrinsic or peculiar velocity is introduced as so that the moments with respect to i c can be conveniently calculated. A set of N moments are then used to describe the state of the gas through Any moment can be expressed by its trace and traceless part [31]. For example, the pressure tensor can be separated as follows: where ij δ is the Kronecker delta, is the deviatoric stress tensor. The angular brackets are used to denote the traceless part of a symmetric tensor. The temperature, T , is given by thermal energy density as The heat flux vector, i q , is defined as: Once the distribution function, f , is obtained from Equation (5), its moments with respect to the molecular velocity, C, can be determined. For example, the density, ρ, and the momentum, ρu i can be obtained from ρ = f dC and ρu i = C i f dC.
For convenience, the intrinsic or peculiar velocity is introduced as so that the moments with respect to c i can be conveniently calculated. A set of N moments are then used to describe the state of the gas through Any moment can be expressed by its trace and traceless part [31]. For example, the pressure tensor can be separated as follows: where δ ij is the Kronecker delta, p = p kk /3 is the pressure, and σ ij = p ij is the deviatoric stress tensor. The angular brackets are used to denote the traceless part of a symmetric tensor. The temperature, T, is given by thermal energy density as The heat flux vector, q i , is defined as: Furthermore, the Grad's moments, m ijk , R ij , ∆, φ ijkl , ψ ijk and Ω i can be evaluated from their definitions by [28].
These macroscopic quantities not only describe the boundary layer accurately but also provide the boundary information for the R26 moment equations solved away from the wall.
On the other hand, the VDF can be approximated by different order of Hermite polynomials using moments [40] as, where H

(n)
A is the Hermite function, and a (n) A is the corresponding coefficient [40]. The local Maxwellian distribution function is given by With the moments available in the R26 moment system, the fifth-order expansion of VDF in Hermite polynomials f (5) can be expressed by [28]: The reconstructed VDF at the coupling interface serves as the boundary condition for the kinetic equation from the bulk flow. In this way, the kinetic equation and the moment equations supply boundary information at the interface for each other and iterate between them until a converged solution is reached.

Entropy and H-Theorem
Two important features of the solutions to the Boltzmann equation are that the distribution function, f , is non-negative and that the solution must satisfy the H-theorem [41]: In the method of moment, the VDF is truncated into the fourth-or fifth-order accuracy; hence, it may not satisfy the H-theorem. Struchtrup and Torrilhon have proved that the linearised R13 equations naturally fulfil the H-theorem [26]. In this section, we will explore the validity of the hybrid DVM/R26 method based on its entropy. For equilibrium flow, the value of H can be calculated with the equilibrium VDF so that, where η is the thermodynamic or equilibrium entropy given by [23] η = 3 2 Rln p ρ 5/3 (22) and e o is an entropy constant equal to (3/2)(ln2π + 1). Therefore, H can be regarded as the entropy for non-equilibrium flows [23,41]. For a homogenous system, the generalised entropy H never decreases with time. In the R26 moment equations, with the VDF truncated at the fifth-order accuracy in Hermite polynomials, Gu and Emerson [28] derived an approximated entropy equation for R26 equations where the flow is not far from equilibrium. An alternative way to evaluate the entropy is by directly integrating the reconstructed VDF. In the dimensionless form, H, it becomes, ln f (5) dC. (23) in which, f (5) represents the dimensionless form of the reconstructed VDF f (5) as where ρ re f and T re f are the reference density and temperature, respectively.

Numerical Test Cases
In this section, we apply the foregoing hybrid DVM/R26 (HYBR) method to simulate several kinds of thermally induced non-equilibrium flows. In all of the cases, both the DVM and the hybrid DVM/R26 method share the same spatial meshes, and the gas medium is modelled as an argon gas. All of the wall boundaries are treated as diffusive walls. The viscosity is obtained from Sutherland's law [42]: where, the Sutherland's constant S for argon is S = 144 K and the reference viscosity µ 0 = 21.25 × 10 −6 Pa · s at the temperature T 0 = 273 K.

2D Thermal Cavity Flow Induced by Temperature Discontinuity
The first case is the flow in a square cavity induced by temperature discontinuities at the cavity boundaries. The geometric configuration is sketched in Figure 1 and the origin of the x, y coordinates sits at the centre of the cavity. The characteristic flow length is defined as the side length of the square cavity, which is L 0 = 10 −5 m. The temperature on the top wall is maintained at T h = 400 K, while on other walls it is maintained at a lower temperature T c = 200 K. The reference temperature is set to be T re f = 300 K and the reference pressure p re f is determined by the Knudsen number, Kn, defined by in which, the reference mean free path λ can be calculated from the initial reference pressure p re f and viscosity µ re f at T re f by We consider three Knudsen numbers in this case, i.e., Kn = 0.1, 0.5 and 1. For all of the cases, the discrete velocity space is discretised in the range of −6 2RT re f , 6 2RT re f 3 with 64 × 64 × 24 non-uniform points. The physical space is meshed with 101 × 101 non-uniform points. For the hybrid DVM/R26 method, we apply 5 grid points near the wall boundary where the ratio of the thickness of the kinetic layer l to the characteristic length L 0 is l/L 0 = 3.64%. The temperature contours and the streamlines, as well as the system entropies at different Knudsen numbers are shown in Figure 2.
On the left side of each plot are the results from the DVM solution of the kinetic equation. The hybrid DVM/R26 results are presented on the right side of each plot.
As indicated in Figure 2, the overall agreement between the DVM and the hybrid DVM/R26 results, especially in terms of the temperature field, T = T/T re f , is good. When Kn = 0.1, the gas molecules travel from the upper hot wall directly to the bottom wall, and there are two more vortices near the left and right walls. As the Knudsen number increases, two vortices closer to the upper wall shrink. The vortices near the side walls begin to dissolve the vortices on the bottom wall; as a result, the side wall vortices grow larger, and the bottom wall vortices become smaller from Kn = 0.5 to 1. Both the DVM and the hybrid DVM/R26 method can capture these vortices accurately. When Kn = 1, the hybrid DVM/R26 method slightly overpredicts the size of the vortices near the bottom wall and slightly underpredicts the size of the vortices near the side walls. Figure 2b shows that the entropy contours calculated from the DVM and hybrid DVM/R26 methods agree well with each other. The entropy at the bottom wall is higher than that in the upper wall. That is because the gas molecules travel away from the upper wall, and they accumulate near the bottom wall region. When Kn = 0.5, the hybrid DVM/R26 method slightly overpredicts the entropies near the side wall. As the Knudsen number further increases, i.e., Kn = 1, the vortices near the side walls and bottom wall merge together, which drive more gas molecules from the upper and lower walls to the centre of the cavity, causing the entropy there to be larger than that near the bottom wall. Both the DVM and the hybrid DVM/R26 methods can reproduce the entropies accurately.
The normalised velocity profiles along the centre lines of the cavity are shown in Figure 3a-f. From these figures, we can see that the hybrid DVM/R26 results agree well with the DVM results when the Knudsen number is below 0.5. In contrast, the R26 moment equations in association with their wall boundary conditions underpredict the absolute value of maximum velocity by about 9.5% and 28.5% at Kn = 0.1 and 0.5, respectively. The hybrid scheme improves the accuracy of the results in comparison with the original R26 moment equations. When Kn = 1, it is very difficult for the R26 moment equations to find a converged solution, and the hybrid DVM/R26 method overpredicts the velocity by about 9.5% at the centre of the cavity.  when the Knudsen number is below 0.5. In contrast, the R26 moment equations in association with their wall boundary conditions underpredict the absolute value of maximum velocity by about 9.5% and 28.5% at Kn = 0.1 and 0.5, respectively. The hybrid scheme improves the accuracy of the results in comparison with the original R26 moment equations. When Kn = 1, it is very difficult for the R26 moment equations to find a converged solution, and the hybrid DVM/R26 method overpredicts the velocity by about 9.5% at the centre of the cavity. Since we only applied 5 grid points in the computational kinetic layer in the hybrid DVM/R26 method, the computational costs can be significantly reduced. For this case, all tests are done on a single processor. The computational cost, in terms of the computational memory and time cost, has been given in Table 2. The convergence criterion for the steady-state is defined by ( ) where n and n + 1 stand for the n-th and (n + 1)-th iterations. Since we only applied 5 grid points in the computational kinetic layer in the hybrid DVM/R26 method, the computational costs can be significantly reduced. For this case, all tests are done on a single processor. The computational cost, in terms of the computational memory and time cost, has been given in Table 2. The convergence criterion for the steady-state is defined by where n and n + 1 stand for the n-th and (n + 1)-th iterations. It is clear to see that the hybrid DVM/R26 method has the ability to save tremendous memory usage by about 70.1% in comparison with the DVM, and thus save the computational time cost, especially at the low Knudsen numbers. As the Knudsen number increases, the convergence rate of implicit DVM also increases rapidly. Therefore, our newly developed hybrid DVM/R26 method is suitable for Kn ≤ 1, especially when the computational domain is much larger than the near-wall region. When Kn > 1, the implicit DVM is fast enough to get the steady-state solutions.

Radiometric Flow
In this case, we investigate another thermally induced flow, i.e., radiometric flow, which is generated by a small plate with differentially heated sides placed in a chamber. The force acting on the small plate is called radiometric force known to be the driven mechanism of the radiometer [1]. The flow configuration as well as the hybrid arrangement is sketched in Figure 4. The hot small plate with a dimension of 3.81 × 0.95 cm 2 sits at the geometric centre of the chamber (x = 0, y = 0). The temperatures of the left and right surface of the plate are kept at T h = 419 K and T c = 384 K, respectively. The upper and lower sides of the plate are maintained at 400K. The size of the outer chamber is 45 × 45 cm 2 , and the temperature is kept at T w = 300 K. The reference length and the reference temperature are defined as the height of the plate L 0 = 3.81 cm and the temperature of the outer chamber T re f = 300 K, respectively. A non-uniform mesh with 63,800 cells is used and it is refined near the surface of the plate. The discrete velocity space is discretised in the range of  The temperature fields T = T/T re f and streamlines, as well as the entropy fields predicted by both the DVM and hybrid DVM/R26 methods are presented in Figure 5. For the case of Kn = 0.1 and 0.5, the overall agreements between the DVM and the hybrid DVM/R26 results are very good in terms of both temperature and velocity fields. The hybrid method slightly underpredicts the size of the vortex near the right side of the hot plate at Kn = 0.1. In terms of the system entropy, the hybrid DVM/R26 method is able to reproduce the accurate entropies with the truncated VDF f (5) at Kn = 0.1 and 0.5. It slightly underpredicts the entropy by about 0.3%. When Kn = 1, the gas is well into the transition regime and the VDF is far away from the equilibrium, the hybrid scheme slightly overpredicts the entropy by about 0.3%. Shown in Figure 5a,c,e are four vortices generated near four corners of the hot plate. The two vortices near the right side of the plate become larger and the other two vortices become smaller as the Knudsen number increases. In each sub-figure, upper and lower half are the results using the DVM and Hybrid DVM/R26 method, respectively. Figure 6 presents the dimensionless normal pressure (normal stress) difference between the left and right side of the plate along the vertical direction, defined in Equation (39), which is the main contribution to the radiometric force as having been analysed [7], where the subscript nn stands for the normal component of the stress tensor relative to the wall. At the lower Knudsen number, i.e., Kn = 0.1, the left/right pressure difference takes larger value near the top and bottom of the plate and lower value near the centre of the plate. A good agreement can be found between the DVM and the hybrid DVM/R26 method at Knudsen number below a value of 0.5. When Kn = 1, the distribution of the normal pressure difference is nearly uniform along the plate surface in the vertical direction, and the hybrid DVM/R26 method overpredicts the pressure difference by about 10%. For this case, all tests are done on a single processor, and the computational memory and time cost are listed in Table 3. The convergence criterion for the steady-state is defined by Equation (41). We can see that the hybrid DVM/R26 method has the ability to save the computational memory and time cost. It is because we only use 5 grid points in each kinetic layer in the hybrid DVM/R26 scheme, so that the Boltzmann model equation is solved only in a small region.  Figure 6 presents the dimensionless normal pressure (normal stress) difference between the left and right side of the plate along the vertical direction, defined in Equation (29), which is the main contribution to the radiometric force as having been analysed [7], where the subscript nn stands for the normal component of the stress tensor relative to the wall. At the lower Knudsen number, i.e., Kn = 0.1, the left/right pressure difference takes larger value near the top and bottom of the plate and lower value near the centre of the plate. A good agreement can be found between the DVM and the hybrid DVM/R26 method at Knudsen number below a value of 0.5. When Kn = 1, the distribution of the normal pressure difference is nearly uniform along the plate surface in the vertical direction, and the hybrid DVM/R26 method overpredicts the pressure difference by about 10%. For this case, all tests are done on a single processor, and the computational memory and time cost are listed in Table 3. The convergence criterion for the steady-state is defined by Equation (28). We can see that the hybrid DVM/R26 method has the ability to save the computational memory and time cost. It is because we only use 5 grid points in each kinetic layer in the hybrid DVM/R26 scheme, so that the Boltzmann model equation is solved only in a small region.

2D Thermal Cavity Flow Induced by the Temperature Gradients
The last case is the thermal cavity flow induced by temperature gradients at wall, which is also a benchmark case to evaluate the accuracy of the numerical scheme. The computational domain is

2D Thermal Cavity Flow Induced by the Temperature Gradients
The last case is the thermal cavity flow induced by temperature gradients at wall, which is also a benchmark case to evaluate the accuracy of the numerical scheme. The computational domain is 10 −5 × 10 −5 m 2 square partitioned by structured rectangular mesh as shown in Figure 7. The left and right walls are maintained at constant temperature T C = 263 K. At the top and bottom walls, we introduce a linearly increasing temperature (from T C = 263 K to T H = 283 K) in left half of domain, and a linearly decreasing temperature (from T H to T C ) in the right half. All the walls are treated as diffusive boundaries. The reference mean free path λ is calculated from the initial uniform density. The reference temperature T re f and the characteristic length L 0 are set to be 273K and 10 −5 m, respectively.

2D Thermal Cavity Flow Induced by the Temperature Gradients
The last case is the thermal cavity flow induced by temperature gradients at wall, which is also a benchmark case to evaluate the accuracy of the numerical scheme. The computational domain is   Three cases corresponding to Kn = 0.1, Kn = 0.5 and Kn = 1.0 are computed. For all of these cases, the spatial space is discretised with 101 × 101 uniform points, and the discrete velocity space is discretised in the range of −6 2RT re f , 6 2RT re f 3 with 64 × 64 × 24 non-uniform points. For the hybrid DVM/R26 method, 5 grid points are used in each kinetic layer near the wall boundary with l/L 0 = 3.64%. The entropy in the hybrid DVM/R26 method is calculated from Equation (23). The comparison of the DVM and the hybrid DVM/R26 method on temperature-gradient-induced thermal cavity flow are shown in Figure 8. In addition, the non-dimensional temperature profiles T = T/T re f along the vertical and horizontal centre lines are presented in Figure 9.
comparison of the DVM and the hybrid DVM/R26 method on temperature-gradient-induced thermal cavity flow are shown in Figure 8. In addition, the non-dimensional temperature profiles ref T T T = along the vertical and horizontal centre lines are presented in Figure 9.
For the cases of Kn = 0.1, 0.5 and 1.0, the hybrid DVM/R26 results are compared with DVM solutions. The overall agreement between the two approaches, especially in terms of the temperature field, is very good. As indicated in Figure 8a,d,g, four vortices are generated with two of them rotating counter-clockwise at the lower left and upper right of the cavity, and another two vortices rotating clockwise at the upper left and lower right of the cavity. As a consequence, the maximum and minimum stresses appear at the centre of clockwise and counter-clockwise vortices, respectively. The absolute values of velocities near the edge of the vortices are higher than that in the centre of the vortices. Both the DVM and the hybrid DVM/R26 methods have the ability to capture these four vortices and flow parameters accurately. It is found in Figure 9 that from the regions near solid walls to the cavity centre, the gas temperature increases along horizontal lines, while it decreases along vertical lines. The maximum temperature value decreases as the degree of rarefaction increases. It is because both the collisions among gas molecules and the interactions between hot wall and gas molecules become weak when the gas is far away from the equilibrium state. In terms of the system entropy, it is very interesting to note that, unlike the other flow properties, the distributions of entropies are totally different between the lower and higher Knudsen numbers. We have tracked the iteration history of the system entropy to investigate its fundamental characteristics in detail, as shown in Figure 10. As indicated in Figures 8c,f,i and 10, for all of the cases, the maximum value of entropy first appears near the centre of the upper and lower walls. Meanwhile, the minimum value of entropy first occurs at the four corners. When the Knudsen number is above 0.5, the variations of the overall entropy distribution contours with respect to the iterations are small. Therefore, the basic shapes of the steady-state entropy contours are similar to those of their initial contours. However, when Kn = 0.1, the maximum values of entropies first appear near the centre of upper and lower walls, and then they move to the centre of the cavity. In this period, the energy transfers from the hot upper and lower walls to the centre, and produces the system entropy, as indicated in Figure 10b. After the left and right cold walls absorb the energy, the values of entropy near the side walls increase subsequently. With strong gas-wall interactions and intensive molecular collisions, more energy can be transferred directly from the centre of the upper and lower hot walls to the side walls, which lead to higher values of entropies near the side wall. For the cases of Kn = 0.1, 0.5 and 1.0, the hybrid DVM/R26 results are compared with DVM solutions. The overall agreement between the two approaches, especially in terms of the temperature field, is very good. As indicated in Figure 8a,d,g, four vortices are generated with two of them rotating counter-clockwise at the lower left and upper right of the cavity, and another two vortices rotating clockwise at the upper left and lower right of the cavity. As a consequence, the maximum and minimum stresses appear at the centre of clockwise and counter-clockwise vortices, respectively. The absolute values of velocities near the edge of the vortices are higher than that in the centre of the vortices. Both the DVM and the hybrid DVM/R26 methods have the ability to capture these four vortices and flow parameters accurately. It is found in Figure 9 that from the regions near solid walls to the cavity centre, the gas temperature increases along horizontal lines, while it decreases along vertical lines. The maximum temperature value decreases as the degree of rarefaction increases. It is because both the collisions among gas molecules and the interactions between hot wall and gas molecules become weak when the gas is far away from the equilibrium state.
In terms of the system entropy, it is very interesting to note that, unlike the other flow properties, the distributions of entropies are totally different between the lower and higher Knudsen numbers. We have tracked the iteration history of the system entropy to investigate its fundamental characteristics in detail, as shown in Figure 10. As indicated in Figure 8c,f,i and Figure 10, for all of the cases, the maximum value of entropy first appears near the centre of the upper and lower walls. Meanwhile, the minimum value of entropy first occurs at the four corners. When the Knudsen number is above 0.5, the variations of the overall entropy distribution contours with respect to the iterations are small. Therefore, the basic shapes of the steady-state entropy contours are similar to those of their initial contours. However, when Kn = 0.1, the maximum values of entropies first appear near the centre of upper and lower walls, and then they move to the centre of the cavity. In this period, the energy transfers from the hot upper and lower walls to the centre, and produces the system entropy, as indicated in Figure 10b. After the left and right cold walls absorb the energy, the values of entropy near the side walls increase subsequently. With strong gas-wall interactions and intensive molecular collisions, more energy can be transferred directly from the centre of the upper and lower hot walls to the side walls, which lead to higher values of entropies near the side wall. For this case, both the DVM and the hybrid DVM/R26 method are simulated on a single processor. The computational costs for this case are listed in Table 4. The convergence criterion for the steady-state is defined by Equation (43). As expected, the hybrid DVM/R26 method is able to save the memory cost by about 68.5%, and thus reduce the computational time subsequently.

Conclusions
In the present study, the application of coupling macro-and microscopic approaches for simulating thermally induced non-equilibrium flows has been explored. The R26 moment equation system is employed at the macroscopic level, meanwhile the Boltzmann model equation associated with the DVM are used to describe the gas dynamics at the microscopic level. Three types of thermally induced flows have been investigated with different Knudsen numbers and the results have been validated using DVM results. The simulation results show that the hybrid DVM/R26 approach can be faithfully used for thermally induced non-equilibrium low-speed flows. Since we only solve the Boltzmann model equation in the near-wall regions, tremendous computational memory and time can be saved in comparison with the DVM. The entropy fields also show that the reconstructed VDF ( ) 5 f is able to yield accurate results when the Knudsen number is less than unity. It is also interesting to find that, unlike the other flow parameters, the distributions of system entropy present totally different characteristics between the lower and higher Knudsen numbers in the temperature gradient cases.  For this case, both the DVM and the hybrid DVM/R26 method are simulated on a single processor. The computational costs for this case are listed in Table 4. The convergence criterion for the steady-state is defined by Equation (28). As expected, the hybrid DVM/R26 method is able to save the memory cost by about 68.5%, and thus reduce the computational time subsequently.

Conclusions
In the present study, the application of coupling macro-and microscopic approaches for simulating thermally induced non-equilibrium flows has been explored. The R26 moment equation system is employed at the macroscopic level, meanwhile the Boltzmann model equation associated with the DVM are used to describe the gas dynamics at the microscopic level. Three types of thermally induced flows have been investigated with different Knudsen numbers and the results have been validated using DVM results. The simulation results show that the hybrid DVM/R26 approach can be faithfully used for thermally induced non-equilibrium low-speed flows. Since we only solve the Boltzmann model equation in the near-wall regions, tremendous computational memory and time can be saved in comparison with the DVM. The entropy fields also show that the reconstructed VDF f (5) is able to yield accurate results when the Knudsen number is less than unity. It is also interesting to find that, unlike the other flow parameters, the distributions of system entropy present totally different characteristics between the lower and higher Knudsen numbers in the temperature gradient cases. Funding: This research was funded by EPSRC Programme Grant EP/N016602/1 "Nano-Engineered Flow Technologies: Simulation for Design across Scale and Phase" and James Weir Foundation.