Kinetic modeling of polyatomic heat and mass transfer in rectangular microchannels

The present study aims at estimating the heat and the mass transfer coefficients in the case of the polyatomic gas flows through long rectangular microchannels driven by small and large pressure (Poiseuille flow) and temperature (Thermal creep flow) drops. The heat and mass transfer coefficients are presented for all gas flow regimes, from free molecular up to hydrodynamic ones, and for channels with different aspect ratios as well as for various values of translational and rotational Eucken factors. The applied values of the Eucken factors were extracted based on the Rayleigh-Brillouin experiments and the kinetic theory of gases. The numerical study has been performed on the basis of a kinetic model for linear and non-linear gas molecules considering the translational and rotational degrees of freedom. The solution of the obtained system of the kinetic equations is implemented on the Graphics Processing Units (GPUs), allowing the reduction of the computational time by two orders of magnitude. The results show that the Poiseuille mass transfer coefficient is not affected by the internal degrees of freedom and the non-dependence of the previous observed deviations with the experimental data on the molecular nature of the gas molecules is confirmed. However, the study shows that the deviation between monatomic and polyatomic values of the mass transfer coefficient in the thermal creep flow is increased as the gas rarefaction is decreased, and for several polyatomic gases met in practical applications in the temperature range from 300 to 900 K might reach 15%. In addition, the effect of the internal degrees of freedom on the heat transfer coefficient is found to be rather significant. The polyatomic heat transfer coefficients are obtained essentially higher than the monatomic ones, with the maximum difference reaching about 44% and 67% for linear and non-linear gas molecules. In view of the large differences between monatomic and polyatomic gases, the present results may be useful in the design of technological devices in which the thermal creep phenomenon plays a dominant role.


Introduction
One of the very critical factors for the design of new technological devices, as well as for the improvement of existing ones, is the study of the internal gas flows in the whole range of the gas rarefaction. Some indicative examples include the vacuum industry (vacuum pumps and gas separators [1][2][3][4]), high altitude aerodynamics (mono-and bi-propellant thrusters and resistojets [5,6]) and the Micro Electronic Mechanical Systems (MEMS) industry (microresonators, comb-drive sensors and gauges [7][8][9]). More information on rarefied gas dynamics applications may be also found in [10,11]. The flow in the aforementioned devices can be categorized according to the gas rarefaction degree introducing the dimensionless Knudsen number (Kn), which is defined as the ratio of the mean free path of the gas molecules to the characteristic length of the examined system. When the Knudsen number is significantly low (typically less than 0.1) the flow behavior can be captured properly by the wellknown Navier-Stokes-Fourier approach supplemented by the velocity slip and temperature jump boundary conditions. However, when the Knudsen number reaches higher values and the gas molecules mean free path is comparable to the characteristic dimensions of the system, the flow behavior in these systems cannot be properly captured by the typical Navier-Stokes-Fourier approach and must be described based on kinetic theory of gases, as described by the integrodifferential Boltzmann equation [12,13] or reliable kinetic model equations.
Over the past decades, a significant amount of research has been conducted studying internal single gas flows in long channels over the whole range of the gas rarefaction on the basis of the kinetic theory of gases where the flow is considered on a molecular level. From the 60's until nowadays, the temperature and pressure driven monatomic gas flows through channels of various cross sections have been investigated by many researches applying the linear form of the original Boltzmann equation [13] or the Bhatnagar-Gross-Krook-Welander (BGKW) [14], and Shakhov (S) [15] kinetic models. Since the literature survey is very extensive, here only an indicative list of these research works is provided [16][17][18][19][20][21][22]. In addition, the temperature and pressure driven flows in micro-channels have also been studied based on the regularized 13-moment equations flows in the transition regime and slip regime [23,24]. Very recently, in [25], the monatomic temperature driven flow in a plane channel has been simulated by applying physics-informed neural network (PINNs)-based approach showing a very good agreement with the available numerical data in the literature. Pressure and temperature driven gaseous flows through microchannels have also been studied experimentally [26][27][28][29][30][31][32]. These experimental studies focus mainly on the correct measurement of mass flow rates and velocity profiles for different gases. The provided results have been used extensively to compare and validate numerical methods.
The available research work on polyatomic gas flows is rather limited and is focused on the flows through parallel plates and tubes (i.e. channels of circular cross-section). The polyatomic heat and mass transfer phenomena appearing in the polyatomic gas flows through parallel plates have been investigated in [33,34]. Corresponding results for channels with circular cross section were presented for both complete (diffuse reflection) [35] and arbitrary accommodation [36]. The aforementioned research on polyatomic gas flows has been conducted by use of the Hansen and Morse model [37] of the linearized Wang Chang and Uhlenbeck equation [38]. Moreover, the Rykov model [39] as well as its modified version proposed by Wu [40], were used in the study of polyatomic gas flows in flat channels and tubes due to small pressure and/or temperature gradients [41,42]. The study of polyatomic gas flows in long tubes was further extended in the case of large temperature difference for a wide variety of polyatomic gases in [43], while in [44] the effect of the vibrational degrees of freedom on the heat flow rates is investigated based on the Holway kinetic model [45]. Overall, in the aforementioned studies, it was deduced that the internal degrees of freedom play a significant role on both heat and mass transfer coefficients in the case of the temperature driven flows. The effect of the internal structure of the gas molecules can reach 30-40% in the case of mass transfer coefficients and 80% in the case of heat transfer coefficients, with this percentage depending on the type of the gas, the working temperature, the flow geometry and the gas rarefaction.
Based on the aforementioned large deviations between monatomic and polyatomic modeling and the high applicability of the polyatomic flows in long rectangular crosssection channels in practice, the aim of the present work is twofold: i) to present data for heat and mass transfer coefficients in the case of polyatomic flows in rectangular channels considering a wide range of flow parameters, and ii) to systematically study the effect of the internal degrees of freedom on the mass and heat transfer coefficients performing comparisons with their corresponding monatomic values. To the best of the authors' knowledge the polyatomic values of the mass and heat transfer coefficients in the case of rectangular channels of various aspect ratios are not available in the literature. Once these coefficients become available the study of the polyatomic gas flows in a wide range of the flow parameters can be performed in a very computationally efficient way. For that reason in this work the obtained results are provided in tabulated form as supplementary material in order to allow for an easy calculation of the pressure and temperature driven flows of various polyatomic gases in the whole range of the gas rarefaction under small and large pressure and temperature gradients. Furthermore, based on the numerical data of the mass and heat transfer coefficients presented in this work for polyatomic flows driven by small temperature/ pressure drops, the study is extended in the case of the flows caused by large pressure-and temperaturedifferences using a methodology previously proposed in the literature for monatomic gases. Several polyatomic kinetic models and numerical methods have been proposed in the literature for solving the Boltzmann equation. A detailed description on the various kinetic models as well as on the numerical methods can be found in Refs. [13,46]. In [47][48][49][50], a comparative study among kinetic models showed that the Rykov model is a solid alternative to describe accurately the mass and heat transfer phenomena simultaneously in the whole range of the gas rarefaction providing results very close to the corresponding experimental measurements. The Rykov model satisfies the collision invariants of mass, momentum and energy and recovers simultaneously all the transport coefficients. A shortcoming of the Rykov model is that the H-theorem has not been proved for it yet, but this does not affect the robustness of the model in practical applications [47][48][49][50]. Thus, in the present study, the kinetic solution is obtained based on the Rykov model subject to Maxwell diffuse boundary conditions by the discrete velocity method. The solution of the kinetic equations has been accelerated using Graphics Processing Units (GPUs).
The rest of this paper is organized as follows: In Sect. 2, the problem description and the kinetic formulation are provided. In Sect. 3, the applied numerical approach is analyzed in detail as well as a description of the adopted parallelization strategies in Compute Unified Device Architecture (CUDA) for solving the kinetic equations on GPUs is provided. Section 4, the obtained numerical results are presented in detail. Finally, Sect. 5 summarizes the main conclusions.

Problem description & kinetic formulation
The flow configuration considered in the present study is illustrated in Fig. 1. Consider a microchannel of length L having a rectangular cross section of height H and width W. The upstream and downstream ends of the microchannel are maintained at constant pressure and temperature P 1 , T 1 and P 2 , T 2 respectively. Without loss of generality, it is assumed that P 1 < P 2 and T 1 < T 2 . As a result of the imposed pressure and temperature gradient a flow is formed in the longitudinal direction z . The length of the microchannel L is considered large enough with respect to its height H, i.e. L > > H, so that the end-effects can be ignored. Based on this last assumption, the local pressure and temperature drop at any cross section of the microchannel may be considered small, even in the case of externally imposed large pressure and temperature drops. Thus, based on the local pressure and temperature conditions, at any cross section of the microchannel the flow can be modeled by applying the linear kinetic theory and the combined flow problem can be decomposed into two independent subproblems: in the first problem only flow caused by the pressure difference is studied (Poiseuille flow), while the second problem is focused on the flow due to the temperature difference (thermal creep flow). The local dimensionless pressure and temperature gradients in the longitudinal direction are defined as where P G and T G are the pressure and temperature of the gas in the flow direction z . The rarefaction level of the gas at a certain position z =z 0 in the flow direction is determined by the rarefaction parameter defined as [51] where P 0 and T 0 are the pressure and temperature of the gas at z =z 0 , R G is the individual gas constant, defined as the ratio of the Boltzmann constant k B to the molecular mass M G , and μ 0 = μ (T 0 ) is the dynamic viscosity at T 0 . The limiting cases of δ 0 = 0 and δ 0 → ∞ correspond to the collisionless and hydrodynamic regimes respectively. Alternatively, the gas rarefaction can be expressed via the Knudsen number (Kn 0 ) which in the case of Hard-Sphere molecules is linked to the gas rarefaction parameter as follows: Kn 0 = √ /(2δ 0 ). When the study involves polyatomic gases and for moderate gas temperatures, besides the translational degrees of freedom, the rotational degrees of freedom have to be considered. A kinetic model which can describe properly the behavior of the diatomic gases in the whole range of the gas collisionality by taking into account the translational and rotational degrees of freedom of gas molecules has been proposed by Rykov in 1975 [39]. Recently, the Rykov model has been generalized to polyatomic gases including linear and non-linear molecules by Wu et al. [40]. Comparisons between numerical and experimental data showed the ability and robustness of the Rykov model to describe properly the heat and mass transfer mechanism in the whole range of the gas rarefaction [47][48][49][50]. In the absence of external forces of any kind, and assuming steady state conditions the linearized form of the Rykov model can be read as [41,43] (1) where h 0 (x,ỹ, ξ p , θ, z ) and h 1 (x,ỹ, ξ p , θ, z ) are the unknown perturbation functions, ξ = (ξ p , θ, z ) is the molecular velocity vector, with ξ p , θ and z being its three cylindrical coordinates, ũ (x,ỹ) is the bulk velocity, while q tr (x,ỹ) and q rot (x,ỹ) are the heat fluxes due to the translational and rotational degrees of freedom of the gas molecules. The single variable j describes the number of the rotational degrees of freedom taking only two values, i.e. 2 and 3 for linear and non-linear molecules, respectively. The macroscopic quantities ũ , q tr , and q rot are expressed as moments of the unknown distribution functions h 0 and h 1 as where the local Maxwellian f G is defined as The determination of the parameters C 1 and C 2 in Eqs. (4) and (5) presupposes the proper description of the thermal conductivities due to translational and rotational motion of the gas molecules by the kinetic model. More specifically, applying the Chapman-Enskog expansion [52] to the kinetic model equations, it can be shown that the rotational and translational thermal conductivity can be obtained in the following form [40,53] : Alternatively, the model parameters C 1 and C 2 are related to the translational and rotational Eucken factors, denoted as f tr and f rot respectively, as follows: It is worth mentioning that, for C 1 = 1 (f tr = 5/2) the Rykov model is transformed into the Shakhov kinetic model, while for C 1 = 0 (f tr = 5/3) into the BGK kinetic model for (10) f tr = 5 3 − C 1 , the monatomic gases. At this point, the following dimensionless quantities are introduced: The kinetic formulation includes the system of two 5D kinetic Eqs. (3)-(7), i.e. 3D in the velocity space and 2D in the physical space. Since the computational cost becomes high, further mathematical processing of the kinetic equations in the direction of the reduction of the computational cost is worth the effort. Following the projection procedure proposed by Chu [54], the dependence on the molecular velocity component in the longitudinal flow direction z can be eliminated by introducing the following reduced distribution functions: After introducing the dimensionless quantities given in Eqs. (12) and the aforementioned reduced functions φ 0 x, y, c p , , φ 1 x, y, c p , and φ 2 x, y, c p , into the kinetic Eqs. (3)-(5), the following reduced kinetic system is obtained: The applied projection process converts the initial 5D problem into a 4D problem (2D in the molecular velocity space and 2D in the physical space).
Introducing the same projection functions φ 0 , φ 1 and φ 2 and the dimensionless quantities c, u, q tr , and q rot into the macroscopic quantities given in Eq. (6), the corresponding macroscopic quantities u, q tr , and q rot are expressed in terms of φ 0 , φ 1 and φ 2 as follows: The kinetic formulation is completed when appropriate boundary conditions along the walls of the channel are formulated. The particles are reflected from the walls following diffuse scattering. Thus, the boundary conditions in each of the four channel walls can be read as h 0 = h 1 = 0 , which in terms of the reduced distributions are written as φ 0 = φ 1 = φ 2 = 0.
Once the macroscopic velocity, the translational heat flux and the rotational heat flux fields are known, the mass transfer coefficients G p , and G T can be calculated as follows: The heat transfer coefficients Q p and Q T are calculated as the summation of their individual parts due to the translational and rotational degrees of freedom: The total mass G and heat flow Q coefficients can be obtained as It is noted that the Poiseuille coefficients (i = P) are obtained for X P , X T = (1, 0) , while the thermal creep coefficients (i = T) for X P , X T = (0, 1) . From Eqs. (14)- (17) it is readily deduced that in the Poiseuille flow the rotational part of the total heat coefficient is zero, i.e. Q rot P = 0 . In the next section, all the details with respect to the applied numerical are presented.
It should be mentioned that depending on the characteristic vibrational temperature of the gas molecules and the desired level of the computational accuracy, the contribution of the vibrational degrees of freedom on the total heat transfer coefficient should also be taken into account [44,46]. The present study was carried out in a wide range of the working gas temperature, gas rarefaction, and translational and rotational Eucken factors. The obtained data for the heat and mass transfer coefficients are presented in tabulated (21) G = −G P X P + G T X T , Q = Q P X P − Q T X T .
form and are representative for several polyatomic gases met very often in several practical applications. The effect of the internal degrees of freedom on the transfer coefficients has been studied with special attention given to three non-polar gases namely N 2 , CO 2 , and CH 4 covering a range of the gas temperature from 300 to 600 K. This temperature range is well below the characteristic vibrational temperatures of the gases, so that the contribution of the vibrational degrees of freedom can be considered effectively small. Based on the gas characteristic data presented in the book of [55] for several polyatomic gases, for the diatomic molecules there is only one characteristic vibrational temperature and in the case of N 2 is 3400 K, while for molecules consisting of more than one atom, several characteristic vibrational temperatures exist, one for each vibrational mode. The lowest characteristic vibrational temperature in the case of CO 2 (four in total) and CH 4 (nine in total) is 960 K and 1880 K respectively. It is noted that N 2 , CO 2 , and CH 4 gases are of great practical importance mainly due to their applicability in several thermal and chemical processes. For instance, N 2 is used as a test gas for microstructure devices for thermal and chemical processes applications [56], CO 2 is used in absorption processes [57,58], and CH 4 is used as a test gas for gas sensor with microchannels [59], as well as a fuel for combustion in microchannels [60] or even in the heat transfer process with application to cryogenic micro channel heat exchangers [61]. Furthermore, for these three gases explicit expressions for the translational and rotational Eucken factors extracted from experimental measurements are available in the literature [62][63][64].

Numerical method & GPU implementation
The numerical solution of the system of the 4D kinetic equations presented in the Sect. 2 is obtained numerically by applying the well-known deterministic Discrete Velocity Method (DVM) in the 2D molecular velocity space and a second order finite scheme in the physical 2D space. More specifically, the continuum molecular speed space c p ∈ [0, ∞) is replaced by a discrete set of molecular speeds c p = c p,1 , c p,2 , … , c p,N C , which are taken to be roots of the Legendre polynomial of order N c accordingly mapped from [−1, 1] to [0, ∞) , while the continuum polar angle space is replaced by a set of N t discrete angles, i.e. Θ = 1 , 2 , … , N t , uniformly distributed in [0, 2π] . In the physical space, the rectangular domain is divided into N x × N y elements, resulting in N x + 1 and N y + 1 grid points in the x and y directions. Then, for each one of these N x × N y elements the kinetic Eqs. (14) are replaced by their corresponding discretized form, which is obtained by integrating both sides of the Eqs. (14) with respect to x direction first and the resulting equations with respect to y direction, using as integral limits the physical limits of the current element, i.e. in the x direction: x min = x i and x max = x i + Δx i with i = 1, … , N x and in the y direction:y min = y j and y max = y j + Δy j with j = 1, … , N y . It is noted that, the spatial derivatives in the x and y directions are eliminated by performing the corresponding integrations analytically, while all the remaining integrations on the left-and right-hand sides of the Eqs. (14) are approximated by applying the trapezoidal rule of second order. For instance, the final discretised expression for the distribution function φ 0 is read as follows: where the indices m = 1, … , N c and n = 1, … , N t represent the discrete velocity magnitudes and angles respectively. In Eq. (22), the superscripts with the signs + and − represent the local position of the corner points of the current numerical element, e.g. i+j+ 0,m,n = 0 x i + Δx i , y j + Δy j , c p,m , n . Once the system of the discretized kinetic equations is solved at each computational node, the integrals of the macroscopic moments (17) are estimated using the current values of the distribution functions φ 0 , φ 1 and φ 2 and applying the Gauss-Legendre quadrature for the velocity magnitude c p and the trapezoidal rule for the angle θ. Then, the updated values of the macroscopic quantities are used as new estimations for the next iteration. The iteration process between the kinetic equations and the corresponding moments of the distribution functions is terminated, when the following convergence criterion: with the pr superscript referring to the computed quantities in the previous iteration step, is satisfied. The double integrals (18)- (20) in the definitions of the mass G i and heat Q i transfer coefficients have been calculated using the trapezoidal rule in both x and y directions. The accuracy of the applied numerical method has been verified in the previous studies [43,49,65], and its reliability to describe heat and mass transfer phenomena in the whole range of the gas collisionality has been confirmed. The results presented in the next section have been obtained with N x × N y = 128 × 128 × (H∕W) , N c = 32 , and N t = 128 .
Numerical tests showed that by doubling independently the numerical parameters N x , N y , N c , and N t , the values of the mass and heat coefficients do not change more than 0.1%. In the recent years the use of Graphics Processing Units (GPUs) has gained significant attention to speedup scientific computations. Since the literature survey on the analysis of the CUDA model is very extensive only a brief description is provided herein. For more details about the CUDA model the reader is referred to Refs. [66,67]. CUDA is a new computing architecture that offloads data in parallel and computes intensive tasks on GPU. A kernel, i.e. a subroutine running on a GPU, is launched with a grid, which consists of multiple threads and blocks. Many threads compose a block, and all blocks together consist the grid. Threads are executed in groups of 32 threads called warps, following the single-instruction, multiple thread execution model [66]. Considering the hardware architecture, instructions are processed on streaming processors (SPs), also known as CUDA cores. Many SPs are organized in streaming multiprocessors (SMs) that share control logic and an instruction cache. Threads of a warp execute the same instruction concurrently, unless a conditional statement forces some of them to follow different path. Thus, with the GPU consisting of thousands of cores for handling multiple tasks simultaneously the reduction of computational time can be achieved.
As the rarefaction degree of the gas is decreased the convergence of the numerical code becomes noticeably slower and the computational cost increases significantly. For instance, for δ 0 = 1 the computational time required is of the order of hours, while for δ 0 = 50 it amounts to several days. In [68,69], it has been shown that solving the kinetic equations on GPUs the computational time can be reduced by orders of magnitude. In the present work, the recently proposed GPU algorithm for the DVM method proposed by Zhu et al. [69], utilizing memory reduction techniques in the velocity and physical spaces, has been adopted in order to speed up the construction of an extensive database for the heat and mass transfer coefficients. In every iteration the kinetic equations are solved independently for each velocity component (c p , θ). The three parts of the code, namely the computation of the distribution functions, the macroscopic quantities and the relative error between iterations, are transformed into separate kernels. More specifically, for each discrete velocity the spatial numerical grid sweeping for the distribution functions is performed from the grid index 1 to the grid index N x in the x-direction and from N y to 1 and 1 to N y in the y-direction for the discrete angles that belong in the fourth and first quarter respectively, and subsequently from N x to 1 in the x-direction and from N y to 1 and 1 to N y in the y-direction for the third and second quarter respectively. Therefore, four kernels were created for the calculation of the distribution functions per quadrant, with each one being called with N t /4 blocks and N c threads. As a result, the dimensions of the arrays in the velocity space needed for the calculation of the distribution functions φ 0 , φ 1 and φ 2 are reduced from (N c , N t ) to N c , N t ∕4 . This reduces significantly the global memory capacity requirements. However, at large values of H/W in order to further reduce the demand for available global memory in GPU, the spatial domain in y-direction is divided into smaller subdomains, which are solved in sequential order with the boundary values of each subdomain being used as boundary conditions for the adjacent subdomains. Examining the memory usage, for N x × N y = 128 × 128 × (H∕W) , N c = 32 , and N t = 128 , the amount of the GPU global memory allocated by the numerical algorithm after applying the memory reduction technique in the velocity space (only N c × N t ∕4 points are stored for each of φ 0 , φ 1 and φ 2 ) for the aspect ratios H/W = 1, 5, 10 and 20 was about 0.4 GB, 2 GB, 4 GB, and 8 GB respectively, instead of about 1.6 GB, 8 GB, 16 GB and 32 GB needed, respectively, in the case that all the N c × N t points are stored in the global memory. Thus, for small and moderate aspect ratios H/W, i.e. H/W < 5, the memory reduction technique in the velocity space allows for running the numerical code in the most of the modern GPUs (having at least 12 GB), while for H/W > 5 the greater need for available GBs in the global memory can be further reduced by applying the memory reduction technique in the physical space too. The macroscopic moments are updated after each distribution function kernel. The kernel of the macroscopic quantities is called with a grid of N x × N y thread blocks each one with N c threads. Each thread block is in charge of calculating the moments of the distribution function at a single point of the numerical grid and the threads of the block work combinedly using the fast shared memory of the CUDA model to perform the binary reductions. Further details with respect to the use of the shared memory and advanced reduction techniques can be found in the books [66,67]. Finally, a kernel for the estimation of the relative error has been developed based on the well-known GPU parallel reduction techniques [66].
Although in [69] a systematic study of the GPU kinetic modeling implementation has been performed, it is still interesting to examine the achieved speedup for the present problem. The results of the present study were obtained using three GPUs, namely the TESLA V100, A100, and K80, released by NVIDIA ® in 2014, 2017, 2020 respectively. The V100 is equipped with 80 SMs (2560 double precision cores) having a maximum clock rate of 1.53 GHz and a peak theoretical double precision floating-point performance of 7.8 TFLOPs. The GPU has 32 Gb of global memory, with a memory bandwidth up to 900 Gb/s, and 48 kB of shared memory. Similarly, the K80 and A100 have 1.37 and 9.7 TFLOPs double precision floating-point performance, respectively. Additionally, the K80 consists of two GK210 GPUs each one with 13 SMs (832 double precision cores) and a global memory of 12 GB with a 240.6 GB/s memory bandwidth, while the A100 is equipped with one GA100 GPU consisting of 108 SMs (3456 double precision cores) and 40 GB of global memory with 1555 GB/s maximum bandwidth. The parallel code was compiled using the CUDA ™ version 11. Sequential code computations were performed on the central processing unit (CPU) Intel ® Xeon ® Gold 6130 (22 MB Cache, 2.10 GHz). In Fig. 2 the achieved speedup in the case of H/W = 1 for one numerical iteration, defined as the ratio of the execution time of the sequential code on one CPU core to the corresponding time of the parallel code on GPU, for various combinations of the numerical parameters N x , N y , N c , and N t is presented. It is noted that, both parallel and serial codes have been compiled without using optimization flags. It is obtained that by doubling the number of N c × N t points in the velocity space the speedup is increased 1.5-3 times. Moreover, in most cases for a given number of N c × N t points, by doubling the number of nodes N x × N y in the physical space, the speedup is increased about 2 times. In addition, the A100 card provides speedups of about 2.5-4 and 1.3-1.8 times higher comparing to the ones achieved by the K80 and V100 cards, respectively. For H/W = 1 with N x × N y × N c × N t = 128 × 128 × 32 × 128 , the obtained speedup of K80, V100 and A100 is about 41, 105 and 147, respectively. For the more demanding case H∕W = 10 with N x × N y × N c × N t = 128 × 1280 × 32 × 128 the achieved speedup when running the code on V100 GPU increases even more, to about 300.

Results and discussion
In this section, the numerical results for the heat and mass transfer coefficients are presented covering a wide range of the involved physical parameters, namely δ 0 , H/W, C 1 (or f tr ) and C 2 (or f rot ). In the present work three aspect ratios H/W have been considered, namely, H/W = 1 which corresponds to the limiting case of channels of square cross section, H/W = 10 which corresponds to moderate aspect ratio and H/W = 20 which is the case of a narrow channel. For the later aspect ratio H/W = 20 the numerical data nearly approach the limiting case of the parallel plates (H/W → ∞) in a wide range of the gas rarefaction [33,39,41]. Explicit expressions for the translational and rotational Eucken factors have been proposed by Mason and Monchick as well as by Uribe et al. in [62] and [63], respectively. An analysis, shows that in the temperature range from 300 to 600 K the data on the translational and rotational Eucken factors extracted from [62] and [63] are in good agreement (within 5%) provided that the same value of the rotationaltranslational collision number and the rotational-energy diffusion coefficient is adopted [63]. It is noted that the translational Eucken factor is obtained essentially lower than its well-known monatomic value of 5∕2 . Very recently, a novel methodology has been proposed by Wu et al. [64] for the determination of the translational Eucken factor based on the Rayleigh-Brillouin experiments. The extracted data by Wu et al. [64] for N 2 are found to be in good comparison with that obtained previously by the theory of Mason and Monchick. Analyzing the data for the translational Eucken factor given in [62][63][64], in a temperature range of 300 K to 900 K and for certain polyatomic gases, namely N 2 , O 2 , CO 2 , and CH 4 , it is found that C 1 ∈ [0.65, 0.9] and C 2 ∈ [0.18, 0.30] . It is noted that the numerical data on the mass and heat transfer coefficients are presented in a wide range of the translational and rotational Eucken factors in order to build a database which represents a relatively large number of polyatomic gases in a wide range of the working temperatures.
In the subsection 4.1, the effect of the internal degrees of freedom on the mass and heat transfer coefficients in conjunction to the aspect ratio (height to weight ratio) is discussed in detail. Also, comparisons to the corresponding previously published numerical data for some limiting cases are performed. In the subsection 4.2, the work is extended in the case of large pressure and temperature differences, while comparisons with the corresponding experimental data are also included.

Polyatomic heat and mass transfer coefficients
In order to validate the numerical code, in Fig. 3, the present Rykov values of the mass transfer coefficients G p and G T for j, C 1 = (0, 1) are compared to that previously reported by Sharipov [51] using the Shakhov model. In addition, in the same figure the corresponding analytical values of these coefficients in the slip and free-molecular regimes are also illustrated for comparison purposes. The analytical expressions for the monatomic values of G p and G T can be found in the chapter 13 of [51]. As it is seen, the present estimations of G p and G T are in very good agreement with those calculated by Sharipov for all considered aspect ratios H∕W , namely H∕W = 1 , 10 and 20, covering a wide range of the gas collisionality, i.e. δ 0 = 10 −3 , 40 . The numerical results are also in very good agreement with the corresponding analytical ones for δ 0 = 0 and δ 0 = 100, with the maximum deviation remaining always less than 0.4% demonstrating the validity of the numerical code. Furthermore, it was found that the reciprocal Onsager relation G T = Q P [70] is satisfied for all examined cases with accuracy of at least four significant figures.
The presentation of the numerical results continues focusing on the effect of the rotational degrees of freedom on the mass flow rate. Figure 4 shows the relative deviation between monatomic and polyatomic mass transfer coefficients, for both Poiseuille flow (i = P) and thermal creep (i = T) in terms of the gas rarefaction degree. Since the mass transfer coefficients depend only on the model parameter C 1 , the results are given for various values of the model parameter C 1 that are representative for several polyatomic gases in the temperature range from 300 to 600 K. As it is seen, the Poiseuille mass transfer coefficient G p depends slightly on the model parameter C 1 and therefore on the internal structure of the gas molecules. It seems that the deviation between monatomic and polyatomic is zero in the free molecular and hydrodynamic regimes, while for intermediate values of δ 0 ∈ [0.1, 10] the deviation curve is Gaussian with its maximum being about 0.5, with this value being decreased as both C 1 and H/W are increased. Therefore, in the free molecular (δ 0 = 0) and hydrodynamic (δ 0 → ∞) regimes, the coefficient G p can be calculated by the corresponding analytical expression for monatomic gases given in the Sects. 13.1.4 and 13.1.2, respectively, of [51]. In the transition regime the values of G p can be calculated based either on the Shakhov model [51] or the present polyatomic Rykov modeling. On the contrary, the thermal creep mass transfer coefficient G T is obtained to be significantly dependent on the rotational degrees of freedom in a wide range of the gas rarefaction. More specifically, for all examined values of C 1 and for small values of δ 0 the relative deviation ΔG T (%) is almost zero, while as δ 0 is increased the relative deviation between monatomic and polyatomic values is also increased up to about δ 0 = 20, while for δ 0 > 20 the relative deviation curve becomes almost independent of δ 0 . Thus, in the free molecular regime the analytical monatomic values of G T , previously reported in Sect. 13.1.4 of [51], hold valid in the case of polyatomic gases, while in the transition and hydrodynamic regimes the present data should be used especially for small values of the translation Eucken factor (or for small values of the model parameter C 1 ). In the hydrodynamic regime applying the thermal slip boundary conditions in [16] and introducing the thermal slip coefficient expression proposed for polyatomic gases given in [71], it can be shown that the coefficient G T in the case of polyatomic gases and for δ 0 → ∞, assuming diffuse scattering of the molecules on the walls, can be estimated as: The present numerical data for G T are in good agreement with the corresponding ones provided by Eq. (24) as δ 0 → ∞. Regarding the effect of the aspect ratio H/W in conjunction with the internal degrees of freedom, it seems that, for any given value of δ 0 in the range from 0.1 to 20, as H/W is increased the deviation between monatomic and polyatomic values of G T is also increased, with this deviation becoming more pronounced at small values of the translation Eucken factor. In the hydrodynamic regime, the relative deviation curve becomes independent of the aspect ratio H/W coefficient, since, as it is easily shown from Eq. (24), both monatomic and polyatomic values of G T are independent of H/W.
The effect of the internal degrees of freedom on the heat transfer coefficients Q tr T and Q rot T due to the thermal creep flow is depicted in Figs. 5 and 6, respectively. In Fig. 5  , between monatomic and polyatomic values are plotted in terms of δ 0 . Similarly to G T , the translational heat transfer coefficient Q tr T , due to the translational degrees of freedom, is independent of C 2 (or f rot ) and depends solely on the model parameter C 1 (or f tr ). From Fourier's law it can be shown that, for large values of δ 0 the translational heat transfer coefficient Q tr T has the following form Equation (25), shows that for δ 0 → ∞ the coefficient Q tr T is independent of the aspect ratio H/W, while it is directly proportional to the translational Eucken factor and (24) G T = 0.45f tr δ 0 .
(25) Q tr T = 3f tr 2δ 0 . inversely proportional to δ 0 . From the present kinetic modeling, it is observed that the translational Q tr T is increased as δ 0 decreases reaching its maximum value in the free molecular regime, i.e. 9G P /4. In the free molecular regime the solution for Q tr T becomes independent of the parameter C 1 (or f tr ), and it is identical with that obtained for monatomic gases. It is noted that, the relative deviations between polyatomic and monatomic values of Q tr T present the same quantitative and qualitative behavior with respect to δ 0 , H/W and C 1 , with that previously mentioned for G T . In Fig. 6, the heat transfer coefficient due to the rotational motion of the gas, namely Q rot T , is shown. From kinetic Eqs. (14)-(16), it can be easily shown that Q rot T depends only on the rotational Eucken factor f rot via the model parameter C 2 . Values for Q rot T are presented for three indicative values of C 2 , as well as for linear (j = 2) and non-linear (j = 3) gas molecules. The effect of H/W on the rotational heat transfer coefficient Q rot T becomes negligible as δ 0 increases, while for small and moderate values of δ 0 the aspect ratio H/W has a strong impact on Q rot T . More specifically, for δ 0 → ∞ the rotational heat transfer coefficient Q rot T can be read as follows: while for δ 0 = 0 as It is noted that G P strongly depends on H/W [51]. In free molecular regime where the effect of H/W on Q rot T is more dominant, it is found that an increase of H/W by 10 times leads to an increase of Q rot T by a factor of about 2.37. In the transition regime, where the solution depends also on C 2 (or f rot ), the quantification of the effect of H/W should be performed considering specific gases. With respect to the model parameter C 2 (or f rot ), it was found that for both linear and non-linear molecules the heat transfer coefficient Q rot T is an increasing function of C 2 under moderate and small gas rarefaction. Comparing the values of Q rot T obtained for the two types of molecules (linear and non-linear), it is seen that the non-linear molecules are characterized by a higher Q rot T comparing to that obtained for linear molecules. In the free molecular regime the ratio of Q rot T of the non-linear molecules to that of the linear molecules is exactly equal to 1.5, while the determination of this ratio in the transition and hydrodynamic regimes presupposes the selection of specific gas molecules in order to consider the influence of the difference in their rotational Eucken factor f rot on the data.  As it is mentioned above, in order to quantify the deviations between monatomic and polyatomic values of the mass and the heat transfer coefficients of the thermal creep someone should consider specific gases. In Tables 1 and 2, the mass and the total heat transfer coefficients of N 2 , CO 2 and CH 4 are shown, respectively, for square (H/W = 1) and wide rectangular (H/W = 10) channels in a wide range of δ 0 . In addition, the corresponding monatomic values of the coefficients for f tr = 5/2 (C 1 = 1) are also included for comparison purposes. The used values of the translational and rotational Eucken factors for N 2 and CO 2 are those extracted from Rayleigh-Brillouin scattering experiments in [64] at 336.6 K and 296.5 K, respectively, while the corresponding values for CH 4 at 300 K have been calculated using the explicit expressions for the thermal conductivity given in [63]. It is noted that, since the translational Eucken factor of polyatomic gases is an increasing function of temperature, with its limiting value at high temperature approaching the monatomic value of f tr = 5/2, the highest deviations between its monatomic and polyatomic vales are expected to be observed at low temperatures. As the rarefaction level of the gas is decreased, the polyatomic values of G T are obtained to be lower than the corresponding monatomic ones, with the maximum relative deviation in the case of N 2 , CO 2 and CH 4 reaching about 3%, 10% and 7%, respectively at δ 0 = 40. Regarding the total heat transfer coefficient Q T = Q rot T + Q tr T , the differences between monatomic and polyatomic gas modeling are more pronounced. In the gases consisting of linear molecules, i.e. N 2 , and CO 2 , the coefficient Q T is increased about 44% at δ 0 = 0.01 and 26%-32% at δ 0 = 40 with respect to its monatomic values. The corresponding deviations in the case of a gas with non-linear molecules (CH 4 ) become 67% and 49% at δ 0 = 0.01 and 40, respectively. In the free molecular regime, the rotational part of the total coefficient Q T is about 44% and 66% of the translational part for linear (N 2 and CO 2 ) and non-linear (CH 4 ) gas molecules respectively, while in the hydrodynamic regime they are reduced to about 36-39% and 60%, respectively. The effect of H/W on the deviations between monatomic and polyatomic kinetic values of Q T becomes more remarkable at δ 0 = 1, with the noticeably high relative deviations observed at δ 0 = 1 for (H/W = 10) (about 40%, 36% and 59% for N 2 , CO 2 , and CH 4 , respectively) being further increased about 1.7%-3% for (H/W = 1).

The case of large pressure and temperature differences
Under the assumption of the long rectangular channel, i.e. L > > H, even in the case of large pressure and temperature drops on the channel ends, the local pressure and temperature gradients are small. In the case of the flows driven by large temperature and pressure differences, the rarefaction parameter varies in the flow direction, i.e. δ = δ (z), and the mass flow rate can be estimated based on the methodology proposed in [72]. According to this methodology, the mass flow rate of a flow through a channel under large temperature and pressure differences can be estimated based on the mass conservation principle. Let us consider a flow through a rectangular duct with its temperature and pressure being T 1 and P 1 at z = 0 , and T 2 and P 2 at z = L . Without loss of generality, it is assumed that P 1 < P 2 and T 1 < T 2 . At this stage it is useful to introduce the following reduced mass flow rate: where Ṁ is the dimensional mass flow rate. Following the proposed methodology in [72], after substituting Eqs. (18) and (21) into Eq. (28), the following ordinary differential equation can be obtained: where z =z∕L , Π(z) = P(z)∕P 1 and (z) = T(z)∕T 1 . Since the reduced mass flow rate G′ is not known a priori the solution of the ordinary differential Eq. (29) can be obtained in an iterative manner as follows: an initial value of the parameter G′ is assumed, and then the Eq. (29) is solved applying the Euler method with the boundary conditions Π(0) = 1 and (0) = 1 at z = 0. After the solution of the initial-value problem, the value of the pressure is estimated at z = 1, i.e. Π(1) = P(1)∕P 1 . In the next step, the calculated value of Π(1) is compared with the actual known value of P 2 ∕P 1 , and in the case that they are not consistent the Eq. (29) is re-solved assuming an updated value of G′. The iteration process is repeated until the absolute difference between the actual and the estimated value of the pressure is less than 10 -7 . At the end, the numerical value of G′ is obtained and the dimensional mass flow can be estimated easily by Eq. (28). With respect to the aforedescribed numerical procedure for calculating the mass flow rate for a given value of the ratios T 2 ∕T 1 and P 2 ∕P 1 the following points should be noted: (i) The rarefaction parameter δ (z) varies in the flow direction taking the values δ 1 and δ 2 at z = 0 and z = 1 respectively. The values of δ 1 and δ 2 are known and can be calculated by Eq. (2). Assuming a viscosity dependence on temperature according to the Inverse Power Law (IPL) model [73] the rarefaction parameter δ (z) can be linked to the local pressure and temperature as follows: where ω is the IPL coefficient (or viscosity index) with its two limiting values being 0.5 and 1 for the Hard-Sphere (HS) and Maxwell interactions. (ii) The mass transfer coefficients G P and G T are functions of both the model parameter C 1 and the local rarefaction parameter δ (z). Hence, based on the database of the coefficients presented in the previous Sect. 4.1 the local values of G P and G T can be estimated at the local value of δ (z), provided that the temperature distribution in the flow direction is known (iii) The temperature distribution (z) in the flow direction is obtained by solving the one-dimensional heat conduction differential equation in the following form [74] : where ̂ = T 2 ∕T 1 and B = T 1 1 − T 1 2 ∕ 1 T 1 − 1 T 2 ,with 1 and 2 being the thermal conductivity of the channel wall at the temperatures T 1 and T 2 , respectively. Values of the thermal conductivities 1 and 2 can be found in [75] for several alloys for a wide range of temperature. An analysis on the thermal conductivity data in [75] shows that in the temperature range from 300 to 900 K and for several alloys the coefficient B lies between 0.3 and 0.7. It is noted that for the limiting value of B → 0 a linear temperature distribution is assumed, i.e.
In Fig. 7, a comparison with the experimental data available for N 2 in [76] is performed in terms of the channel conductance in liter per second (l/s), defined as [27] (31) where R u = 8.31446 J K −1 mol −1 is the universal gas constant, M G = 28.0134 amu is the atomic mass of N 2 and T 0 = 273.15 K is the reference temperature. In the experimental measurements a channel with square cross section of length L = 1.28 m and height H = 0.01589 m has been tested. Hence, the hypothesis of the long channel, i.e. L > > H, adopted in the present study is valid. The experimental measurements have been performed in the experimental TRANSFLOW (Transitional Flow range experiments) facility which was set up at the Karlsruhe Institute of Technology (KIT) in Germany. A detailed description of the TRANSFLOW facility can be found in [77]. Since the flow is driven solely by the pressure difference between the channel ends, the corresponding numerical solution for the conductance is obtained by the solution of the Eq. (29) assuming d ∕dz = 0 and (z) = 1 supplemented by the Eqs. (28), (30) and (32). In Fig. 7, the experimental conductance data are compared to the corresponding numerical monatomic and polyatomic values of conductance in a range of the average rarefaction parameter δ av , i.e. δ av = (δ 1 + δ 2 )/2. The pressure ratio P 2 ∕P 1 varies from 26 to 173. As it seen, the experimental data and the present numerical calculations are in very good agreement. It is noted that, the polyatomic values of the conductance are obtained to be almost identical with the corresponding monatomic ones in the whole range of the gas rarefaction. This is justified by the fact that the mass flow rate in the pure pressure driven flows is almost unaffected by the internal degrees of freedom as it has been discussed in the Sect. 4.1. The observed deviations (always within 10%) may be explained by the back scattering phenomenon discussed in [76] and as it is obtained by the present analysis they are not related to the diatomic nature of N 2 .
It is worth mentioning that, due to the deviations between polyatomic and monatomic values of G T mentioned in the Sect. 4.1, corresponding large deviations are expected in the temperature driven flows under large temperature differences. Thus, it is interesting to quantify the corresponding deviations in the case of P 2 ∕P 1 = 1 and large values of T 2 ∕T 1 . Table 3 presents the values of the reduced flow rate G′ of N 2 and CO 2 for T 1 = 300 K, T 2 ∕T 1 = [2,3] , δ 1 = [0.01, 50] , H/W = [1,10], assuming HS viscosity values ( = 0.5) and linear variation of the temperature (B → 0) . In addition, the corresponding monatomic reduced flow rates G′ are also shown for comparison purposes. In the modeling of N 2 and CO 2 flows the temperature dependence of the translational Eucken factor has been taken into account and the values have been calculated using the explicit expressions proposed by Uribe et al. in [63]. As it is seen, for small values of the rarefaction parameter δ 1 , i.e. δ 1 = 0.01, the values of G′ obtained for N 2 and CO 2 are almost identical with that obtained by the Shakhov model for monatomic gases. However, for δ 1 > 0.01, the reduced flow rates G′ of N 2 and CO 2 are obtained lower than the corresponding monatomic ones, with this difference becoming more pronounced in the hydrodynamic regime and reaching about 5% and 10% at δ 1 = 50.
To extend the analysis beyond the limiting Hard-Sphere viscosity model, in Table 4, the corresponding values of G′ of N 2 and CO 2 based on the IPL viscosity model are shown. The values of the IPL coefficient for N 2 and CO 2 are taken equal to 0.71 and 0.87 respectively, in order to yield a very good fit (within 3%) for the viscosity experimental data given in [78] in the whole examined range of temperatures from 300 to 900 K. In addition, in order to investigate the effect of the temperature distribution on the reduced flow rate G′ the following two cases are considered: B → 0 and B = 0.5 . Comparing the results shown in Table 4 for N 2 and CO 2 with B → 0 to the corresponding ones given in Table 3, the influence of the type of the molecular interaction is studied. As it is seen, for δ 1 > 1 the deviation between the HS and IPL values is always higher than 2%, with the IPL values being always higher than the corresponding HS ones, reaching at δ 1 = 50 about 9% and 17% at T 2 ∕T 1 = 2 and 15% and 30% at T 2 ∕T 1 = 3 for N 2 and CO 2 , respectively. It is noted that similar conclusions have also been drawn in the case of temperature driven flows through long tubes in [43]. Next, focusing on the data for B → 0 and B = 0.5 , presented in  Table 4, it is deduced that the values of G′ assuming linear temperature distribution (B → 0) are always higher than those based on a non-linear temperature variation (B = 0.5) according to the thermal properties of the channel walls. From the data in Table 4, we conclude that the effect of the temperature distribution on the reduced flow rate G′ is significant for moderate and large values of δ 1 with the maximum deviation between the cases B → 0 and B = 0.5 reaching about 3-4% and 7-7.5% at T 2 ∕T 1 = 2 and 3, respectively. It is noted that, in the pure thermal creep flow despite the fact that the pressure difference between the two ends of the channel is zero, i.e. P 2 ∕P 1 = 1 , the pressure distribution varies in the flow direction. In Fig. 8 the dimensionless pressure distributions of N 2 and CO 2 along a channel with square cross section (H/W = 1) are compared to the corresponding monatomic ones for T 2 ∕T 1 = 2 and 3 assuming typical values of δ 1 . As it is seen, the pressure distribution does not remain constant, with its variation becoming more pronounced as T 2 ∕T 1 is increased. Initially the dimensionless pressure increases from its inlet value (Π(0) = 1) up to a certain value between z = 0.4 and 5 and then it is decreased to reach the value at the exit, i.e. Π(1) = 1 . The variation of the pressure along the channel is increased as δ 1 is decreased and the gas becomes more rarefied. It is also observed that, near the free molecular regime, i.e. δ 1 = 0.01, all the pressure distributions coincide, while for intermediate values of δ 1 , i.e. δ 0 = 1, they depart from each other with the pressure of the monatomic gas being higher than the corresponding polyatomic ones. Similar observations can be made in the case of a rectangular channel with high aspect ratio, i.e. H/W = 10.
Simulations have also been performed in the case of a flow driven by both pressure and temperature drops. It was found that the deviations between monatomic and polyatomic values of the reduced flow rate G′ are negligible. It seems that, for moderate and large values of the rarefaction parameter δ 1 the deviations between monatomic and polyatomic values of G′ observed in the thermal creep are eliminated by the Poiseuille flow which seems to be the dominant phenomenon in the case of the combined flow. In Table 5, the monatomic values of G′ for the combined flow for HS model at P 2 ∕P 1 = 100 with B → 0 are presented for completeness purposes. The negative values of G′ indicate a flow formed form the high pressure reservoir to low pressure reservoir (see Fig. 1). It is also observed that as the temperature is increased and the contribution of the thermal creep in the overall flow is also increased the reduced flow rate is decreased. Overall, the values of G′ in the combined flow are essentially higher than the corresponding ones in the thermal creep, and increase as the rarefaction of the gas is decreased.

Concluding remarks
This work presents a numerical investigation of the pressureand thermally-induced non-polar polyatomic gas flows in long rectangular cross-section channels, paying special attention on the differences between polyatomic and monatomic gases.
The investigation was carried out based on the linear Rykov kinetic model, which takes into account the translational and the rotational degrees of freedom of the gas molecules, and assuming diffuse scattering of the gas at the channel walls. The kinetic solution of the problem has been accelerated on Graphics Processing Units (GPUs), allowing the reduction of the computational time by two orders of magnitude. According to the assumption of the long channel the flow problem has been decomposed into two separate sub-problems: one in which the flow driven by the pressure drop, the so-called Poiseuille flow, and another with the flow being caused by the temperature drop, the so-called thermal creep flow. As for the mass transfer coefficients, in the Poiseuille flow, it was found that the internal degrees of freedom do not affect the mass transfer coefficient and the previously published data for monatomic gases can be applied in any polyatomic gas. In the thermal creep flow, the effect of the internal structure of the gas molecules on both mass and heat transfer coefficients may be significant. The results show that the deviation between monatomic and polyatomic values of the thermal creep mass flow rate is negligibly small in the free molecular regime. However, as the gas rarefaction is decreased and the flow enters the transition regime the deviation becomes a function of the translational Eucken factor and the aspect ratio of the channel, while as the rarefaction is further decreased, i.e. the so-called hydrodynamic regime, the observed deviations in the mass and heat transfer coefficients represent the corresponding deviations in the translational Eucken factor. In the case of the thermal creep under small or large temperature differences between the channel ends, it was found that in the hydrodynamic regime the mass flow rates of N 2 , CO 2 and CH 4 differ from the corresponding monatomic ones about 3%, 10% and 7%, respectively. These differences may be further increased for polyatomic gases with lower values of the translational Eucken factor.
As for the heat transfer coefficients, it is shown that the rotational degrees of freedom have a dominant effect on them. In the Poiseuille flow, the rotational heat transfer coefficient vanishes and the total heat transfer coefficient becomes independent of the rotational degree of freedom and equal to the mass transfer coefficient in the thermal creep. In the thermal creep flow the total heat transfer is divided into its translational and rotational parts. The polyatomic values of the translational heat transfer coefficient are obtained to be lower than the corresponding monatomic ones, with this difference being increased as both the gas rarefaction and the translational Eucken factor are decreased. However, in polyatomic gases the additional contribution to the total heat transfer caused by the rotational motion of the gas molecules leads to significantly higher values of the heat transfer coefficient compared to that obtained for monatomic gases. In the case of N 2 and CO 2 flows, the total heat transfer coefficient may differ, depending on the level of the gas rarefaction, about 26%-44% from the corresponding monatomic values, with this difference being 44%-67% for CH 4 flow.
Regarding the lateral wall effect on the transfer coefficients, it is concluded that for a specific translational Eucken factor and gas rarefaction, the increase of the aspect ratio causes an increase in the differences between monatomic and polyatomic values of the mass and translational heat transfer coefficients in the thermal creep flow. The influence of the aspect ratio on the rotational heat flow rates is increased as the gas rarefaction is increased. An increase of the aspect ratio by 10 times leads to an increase of rotational heat transfer coefficient by a factor of about 2.37.
In addition, the Poiseuille and thermal creep flows under large pressure and temperature ratios have been investigated. The results show that as the gas rarefaction is decreased the reduced flow rates based on the IPL viscosity model are always higher than the corresponding HS ones. The effect of the temperature distribution on the reduced flow rate is large for moderate and small gas rarefaction and can reach about 7%.
The present calculated values of the transport coefficients, which were not available previously in the literature for rectangular ducts and here are provided in a tabulated form in the accompanying supplementary material, may support the design of technological devices in which the thermal creep phenomenon plays a dominant role in the formation of the flow.

Conflict of interest
The authors and the corresponding author state that there is no confict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.