A Closed-Form Solution of a Kinetic Integral Equation for Rarefied Gas Flow in a Cylindrical Duct

A spectral method based on Hermite cubic splines expansions combined with a collocation scheme is used to develop a solution for the vector form integral S-model kinetic equation describing rarefied gas flows in cylindrical geometry. Some manipulations are made to facilitate the computational treatment of the singularities inherent to the kernel. Numerical results for the simulation of flows generated by pressure and thermal gradients, Poiseuille and thermal-creep problems, are presented.


Introduction
Research in the field of rarefied gas dynamics (RGD) involving the analysis and description of micro and nanoflows has been fundamental for the development of micro devices and new technologies, specially in the field of micro-electro-mechanical systems (MEMS) [1]- [5].
The Boltzmann equation is the fundamental mathematical model for dealing with problems in RGD [6].Investigations of its solution and associated kinetic models have found much attention in the last few years [7]- [13].Numerical and analytical approaches have been developed, mainly in Cartesian coordinates [14] [15].
The solution of Boltzmann-like problems is even more challenging in multidimensional and complex geometries, particularly, if one seeks for analytical approaches.In cylindrical coordinates, to overcome such difficulties, the integral form of the kinetic equations has been employed.In this context, several approaches have been proposed in the past years [16]- [19], particularly related to the solution of the BGK model [20], as, for example, the use of classical Legendre polynomials expansions [21].A special approach used in an earlier work of Ferziger [22], is based on a transformation proposed by Mitsis [23], which allows us to recast the problem in cylindrical coordinates into a much simpler problem in plane Cartesian coordinates.The transformation, however, is restricted to some specific (simpler) physical cases.
Mitsis' transformation [23] was used by Siewert and Valougeorgis [24] to develop a solution for the integral form of the S-model [25] kinetic equation, which is a vector equation, while the BGK kinetic equation is a scalar problem.It has been addressed in the literature, however, that the S-model is more appropriate than the BGK model to deal, for instance, with nonisothermal flows [24] [26].
In this work, we use an spectral approach based on Hermite cubic splines expansions along with a collocation scheme to solve the vector integral form of the S-model kinetic equation for flow in cylindrical tube.Numerical results are obtained for the Poiseuille and thermal-creep problems.As in other analytical approaches [24] the solution is explicit, in terms of the spatial variable.However, the technique proposed in this work, that has been also effective to deal with the BGK integral equation [27], may be extended to deal with the problems which include reflexive boundary conditions.

The Formulation
The initial formulation of the problems we treat in this work, is the integro-differential formulation associated with the kinetic S-model [24], for describing the behavior of gas flows in a cylindrical tube of radius R , with no variation in the axial direction, such that the distribution of particles and the boundary conditions depends on the spatial variable ( ) , written in dimensionless units, and the particle velocity vector c , also given in dimensionless units, expressed in the cylindrical coordinates.Still, in the initial formulation, the flow is subject to reflexive boundary conditions at r R = .In Ref. [24], the authors follow a previous derivation regarding to the BGK model [28] to derive in details the integral form of the equation for the S-model, which describes the flow of a rarefied gas in a straight cylindrical duct of radius R , which will be our starting point of this work.Since the derivations were done in details in [24] [28], we do not repeat it here.We start, already, from considering the integral vector equation [24] ( ) ( ) ( ) ( ) . In Equation (1) the inhomogeneous term Γ is given by ( ) where the constants 1 k and 2 k are defined so that 1 1 k = and 2 0 k = for the so-called problem of Poiseuille flow while 1 0 k = and 2 1 k = for the thermal-creep problem.Still, the kernel of integral equation is where and Here ( ) 0 I x and ( ) 0 K x are the modified Bessel functions of zero order, first and second kind [29], respectively.Furthermore, in Equation (3), ( ) with ( ) ( ) Now, if we define [24], we can evaluate the quantities of physical interest for the gas flow, such as the velocity profile [24] ( ) ( )  G (11) and the heat-flow profile ( ) ( ) ( ) Consequently, we write the particle-flow rate as ( ) and the heat-flow rate ( ) The Equations (11) to (14) are used for both problems, Poiseuille and thermal-creep, depending on the specified values of the constants 1 k and 2 k which define the expression used for the source term given in Equation (2).
Finally, we note that, in Equation (1), the term brings the information from the boundary condition, and it is not really an inhomogeneous term when reflexive boundary conditions are applied.For this reason, the only case considered in Ref. [24] was the diffuse reflection, where, at the end ( ) r = B 0 .At this point, in this work, we also disregard this term in our formulation and will comment on that later on in this text.

The Solution
Firstly, we rewrite Equation (1) as, ( ) ( ) ( ) and for the case with no specular reflection, ( ) . We then propose a solution to Equation (15) in terms of a truncated expansion of Hermite cubic splines functions [30], as where whose components (constants) a α and b α are to be determined.In regard to the splines functions ( ) first of all, we consider there to be 1 (see Appendix for details).
If we now substitute Equation ( 16) into Equation ( 15) we obtain which still can be expressed in the form and At this point, we introduce a collocation scheme, such that, if we evaluate Equation ( 19) at the collocation points 2 , 0,1, 2, , , we end up with + linear system to the constants a α and b α used in expansion given by Equation ( 16).

Computational Aspects
The procedure described in Section 3 may be considered quite straightforward.However the solution of the linear system given by Equations ( 23), depends on the explicit evaluation of the integrals , which are written in terms of the kernel, Equations (3) to (9), which clearly includes singularities that introduce numerical difficulties.In this way, a series of computational steps have to be taken into account.
As a first step, as used in previous works [24] [31], for computational reasons, the modified Bessel functions are rewritten as ˆê and e , such that Equations ( 4) and ( 5) are now expressed as for < y x and for > y x.To be consistent, from now on, we consider the previously given integral equation kernel, Equation (3), expressed in the form Continuing, if we look back to Equations ( 6) to ( 9), we can see that with and where, again, the (constant components) matrices 0 Δ , 2 Δ , 4 Δ are defined in Equations ( 7) to (9).More importantly, in writing the kernel as above, it is to point out that in these definitions, the only expression which effectively involves the term ( ) which induces a singular behavior, is the first one in the right hand-side.We note, however, that this same expression was part already of the BGK model kernel we had to deal in a previous work [21].
Finally, for future use, later on in this text, we can also express the kernel of the integral equation, in a matrix form, , , , , , following, of course, the previous definitions given in Equations ( 7) to (9) and Equations ( 29) to (31).Here, ( and ( ) ( ) Going back to the definition of the linear system we need to solve, Equations (23), we now focus our attention in the definition of the functions given in Equations ( 20) and (21).
is the support of the spline function ( ) , in other words, In an analogous form, In fact, to be more specific, we can express the vector , in terms of its components ( ) ( ) ( ) and , where a α and b α are the constants to be defined in the proposed expansion, Equation ( 16), obtained as the solution of the linear system Equations ( 23).The terms ( ) , defined by Equations ( 39) and (40), is rewritten in a matrix form , and .
Here ( ) , 0, Considering the explicit derivations above, at this point, we can rewrite the linear system defined in Equations ( 23), in the following form ( ) , p q i h x α are given by the expressions shown in Equations ( 41)-( 46) and ( 50)-(51).Still, the terms 1 B and 2 B are defined by the components of the source term Γ , so that depending on the problem, as mentioned earlier: 1 1 k = and 2 0 k = for the Poiseuille flow; 1 0 k = and 2 1 k = for the thermal-creep problem.
The linear system given by Equations ( 52) and ( 53), is of order 2 2 L + (or 4 4) M + , which solution, the constants a α and b α , allow us to evaluate the quantities of physical interest such as velocity profile ( ) ( ) and the heat flow profile ( ) ( )

Numerical Quadrature Schemes
Searching for a general procedure, as much as possible, the basic approach for evaluating the integrals involved in this derivation was based on the use of the Gauss-Legendre quadrature scheme [29] with nodes n µ and weights n ω defined in the interval [ ] 1,1 − .In this form, convenient mappings, according to the integration interval, were used.
In this way, we define and we consider a quadrature scheme with i t N nodes and weights for evaluating the integrals defined in Equations ( 27) to (29) as At this point, however, a special remark has to be done in regard to the evaluation of Equation (60).While the general procedure described above works well as long as Ry Rx − is not small-in fact, using one obtains results which agree in seven digits with results obtained by the software Maple-the case Ry Rx → has to be treated more carefully, if ones search for more accurate results as a general case.To deal with that, we subtract (and add) from Equation (29) the expression ( ) ( ) In the resulting integral expression, we introduce the change of variables and, subsequently, once more, we add (and subtract) the integral term ( ) To be clear, when Ry Rx → (we considered typically 0.05 Ry Rx − ≤ ) we evaluate ( ) ( ) The reasons for doing that is, firstly, the use of these additional functions results in a smoother function when approaching 0 τ = .In addition to that, ( ) 1 , S y x can be expressed as

E y x Rx y x S y x E x y Ry y x
where ( ) Ê x is the elliptic integral of the first order, that can be evaluated numerically by applying the algorithm given by [29].Still, 2 S , can be evaluated by using Maple, such that 2 1.4020222280980703 S = .The integral ( ) , S y x , given by Equation (67), can be then evaluated without much difficulty by applying the scheme of Gauss-Legendre quadrature, by mapping [ ] .This procedure is, in fact, analogous to the one we used to deal with the BGK model [21] [27] and it was implemented once we seek for increasing the accuracy of the final results, as we will comment later on in this text.
To all other integrations processes needed in this formulation, more specifically, the expressions given by Equations ( 46) and (51) as well as Equations ( 13) and ( 14), linear mappings from the integration interval to the interval [ ] 1,1 − were performed, to the use of the Gauss-Legendre quadrature scheme.Clearly, in the case of Equations ( 46) and (51) the dependence of the integration limits on the interval definition of the splines functions have to be taken into account to improve the computational results.In this way, we evaluate these equations as With this, we conclude the evaluation of the elements of the system given by Equations ( 52) and (53).

Numerical Results
As input data in our program we need to specify the radius of the duct, R , the order of the quadrature schemes ( ) N N N and the number M of knots (associated to the splines functions) which will define L in the linear system given by Equations ( 52) and (53).We implemented a FORTRAN program and subroutines from LINPACK [32] were used to solve the linear system.The solution was then used to obtain the physical quantities of interest, the velocity and heat-flow profiles, Equations ( 57) and (58), and subsequently Equations ( 13) and (14).
As a first test, to check and to have confidence in our approach, we reproduced all the results listed in Ref. [24], for the velocity and heat-flow profiles (case 2 R = ), as well as velocity slips and particle flow rates for values of R ranging from 0.001 = R to 100 R = .We do not repeat the results here since they were already reported in [24].
As we mentioned earlier, the approach used in that work can be very accurate, since it works with a Cartesian coordinates problem associated with a very accurate discrete ordinates solution, after using Mitsis' transformation [23].In order to get agreement in six to seven digits, for all cases, with Siewert and Valougeorgis' results (listed with seven digits) [24], we used 300 M = , r N = results can also be reproduced, in general, with two to three digits of agreement with [24], particularly for < 1 R .We then investigated additional cases.To provide the results we list in Table 1 and Table 2, we used the first set of parameters listed above.In regard to Table 2, we remark that, the Onsager reciprocity relation [33] [34], T P U Q = is verified.Still, the results listed in Table 2 have some agreement with Ref. [26], within the accuracy   declared in that paper (0.5%).In addition to provide accurate results, it is important to remark that, the approach proposed in this work may be used to simulate the case which includes specular reflection at the surface ( ) ( ) , as we show in a subsequent work, and, a better description of the boundary effects is a relevant issue when it comes to the modeling of microflows.

Conclusion
In this work, we have used a classical spectral approach, based on the use of Hermite cubic splines combined with a collocation scheme, to obtain a solution, in a closed-form, to the S-model kinetic integral vector equation relevant to describe rarefied gas flow in cylindrical coordinates.The main purpose of this work was to verify the performance of this approach as an extension of an earlier work related to a simpler integral model equation, the BGK model, trying to derive, as much as possible, a general approach which can be now, in a subsequent work, extented to include the treatment of reflexive boundary conditions.The numerical results obtained showed good agreement with accurate results available in the literature.

10 R
= .If we use, for instance, lower order approximations as 10

Table 1 .
Velocity and heat flow profiles in a duct with R = 1.