PARAMETRIC STUDY ON PROPULSION PERFORMANCE OF MICROTUBES

The pressure-driven rare¦ed gas §ow of polyatomic gases through short tubes in a wide range of the Knudsen number is numerically investigated. The downstream over the upstream pressure ratio is taken very close to zero. Such §ows are characterized by low Reynolds numbers and high viscous losses and, therefore, short circular microtubes may be used instead of typical micronozzles. The main computed quantities include the §ow rate, the discharge coe©cient, the thrust, and the impulse factor which are provided in terms of the gas rarefaction and the tube dimensionless length. Based on the above, a parametric study on the propulsion characteristics of microtubes is provided. Furthermore, a comparison between corresponding polyatomic and monoatomic results is performed and the e ̈ect of the internal degrees of freedom on the results is investigated.


INTRODUCTION
In several practical applications including micro propulsion nozzles in high altitude, the §ow may be in the whole range of the Knudsen number and modeling should be based on deterministic techniques solving the Boltzmann equation or, alternatively, on the stochastic direct simulation Monte-Carlo (DSMC) method.Micronozzles are often used as low-thrust propulsion systems in order to produce accurate orbital maneuvers in microsatellites.Therefore, a systematic study of the gas §ow in such devices is needed to determine the optimal geometry and design.It is well known that at low Reynolds numbers, the viscous losses in micronozzles become large enough making the concept of a nozzle expansion useless and micronozzles can be replaced by short circular tubes.Rare¦ed monoatomic gas §ows through capillaries due to pressure gradients have been extensively studied in the past [13].The corresponding work in polyatomic gases is commonly handled based on the DSMC method, while similar simulations based on kinetic modeling are limited.
In the present work, the rare¦ed gas §ow of polyatomic gases through short circular tubes due to pressure gradients is modeled via the polyatomic Holway model of the Boltzmann equation [4].This model holds the entropy inequality and ensures good agreement with experimental data concerning heat transfer [5,6] and steady condensation [7] con¦gurations.The purpose of the present work is to numerically investigate polyatomic gas expansion into very low pressures in a wide range of the Knudsen number and to compute the deduced §ow rate, discharge coe©cient, thrust, and impulse factor in terms of §ow and geometric parameters as well as to examine the e¨ect of the internal degrees of freedom by comparison with corresponding monatomic §ows.

FLOW CONFIGURATION AND BASIC PARAMETERS
Consider two large reservoirs A and B which are connected by a tube of radius R and ¦nite length L. The polyatomic gas in the containers far from the tube is in equilibrium at pressures P A and P B , with P B /P A = 0.01.The walls and the gas in the container areas far from the tube are maintained at the same temperature T 0 .The computational domain consists of the two large computational areas, which correspond to the upstream and downstream reservoir including an intermediate area which contains the tube.The §ow con¦guration and the computational domain are shown in Fig. 1.The reference quantities are R, T 0 , and P A , while the most probable speed, de¦ned as υ 0 = 2k B T 0 /m with m and k B denoting the molecular mass and the Boltzmann constant, respectively, is taken as the reference velocity.The §ow is characterized by the ratio L/R and the rarefaction parameter de¦ned as where µ 0 is the reference viscosity at temperature T 0 .It is noted that the rarefaction parameter is proportional to the inverse Knudsen number, with the limiting values at δ 0 = 0 and δ 0 → ∞ corresponding to the free molecular and hydrodynamic limits, respectively.In the temperature range where the e¨ects of vibrational degrees of freedom can be neglected, the problem may be modeled by the Boltzmann equation for a gas of rigid rotators.The investigation is based on the description of the state of a polyatomic gas using the distribution function f ( r, z, ξ, e), which is a function of the spatial coordinates r and z, the molecular velocity ξ, and the rotational motion energy e.Then, the macroscopic quantities of practical interest are obtained by the moments of the distribution function.Furthermore, it is convenient to introduce the dimensionless independent variables z = z/R, r = r/R, and c = ξ/υ 0 as well as the dimensionless macroscopic quantities: where N is the number density; P is the pressure; u is the velocity vector; and T is the total (thermodynamic) temperature.The subscripts tr and rot denote the translational and rotational parts, while the parameter j is the number of rotational degrees of freedom (j = 0 refers to monoatomic molecules, j = 2 to diatomic and linear polyatomic molecules, and j = 4 to nonlinear polyatomic molecules).The dimensionless §ow rate, thrust, impulse factor, and discharge coe©cient, which are the overall quantities of practical interest, are de¦ned, respectively, as Here, ' M , F t , and I sp are the dimensional mass §ow rate, thrust, and impulse factor, respectively; g r = 9.81 m/s 2 is the gravity accelaration; and γ = (5 + j)/(3 + j) is the ratio of the speci¦c heats of the gas.The subscript ¢ex£ denotes the values at the exit of the tube (z = L/R).Also, the local Mach number is given by where c s = γk B T 0 /m is the speed of sound and |u| = u 2 r + u 2 z is the magnitude of u.

KINETIC MODELING
The e¨ort of numerically solving the Boltzmann equation is signi¦cantly reduced by substituting its collision term with reliable kinetic models.The well-known model introduced by Holway [4] is implemented.The H-theorem can be proved in a straightforward manner for this model following the arguments leading to analogous proof of the BhatnagarGrossKrook (BGK) model.It is obvious that the dependency of the distribution function f on the energy e of the rotational motion signi¦cantly increases the computational e¨ort.However, for the speci¦c problem under consideration, the computational e¨ort is reduced by eliminating, based on a projection procedure, the e component of energy by introducing the following reduced distributions [8]: Then, for the present §ow problem, the Holway model may be written in the dimensionless form as where Here, the hard-sphere intermolecular potential has been applied.The parameter 0 < Z −1 ≤ 1 indicates the rotational collisions as a fraction of the total number of collisions.It is noted that as the parameter Z → ∞, the ¦rst equation in (1) reduces to the kinetic BGK-model for a monoatomic gas [9].The dimensionless macroscopic quantities are expressed in terms of the functions g and h as Next, to close the problem, the formulation of the boundary conditions is provided.Purely di¨use boundary re §ections are considered at the walls and symmetry is imposed on r = 0.At the open boundaries, for the incoming particles, a Maxwellian distribution is imposed based on the local values of the pressure and temperature assuming zero bulk velocity.The part of distribution function corresponding to molecules exiting the domain is a part of the solution.Based on the above, the boundary conditions in dimensionless form can be written as upstream reservoir: (3) downstream reservoir: solid walls: The superscripts ¢+£ denote the outgoing distribution from a surface.The parameter n w is calculated by the no penetration condition at the walls.The set of integrodi¨erential Eqs.(1) with the boundary conditions (3)(5) are solved numerically discretizing in the physical space by the control volume approach and in the molecular velocity space by the discrete velocity method.The macroscopic quantities are computed by GaussLegendre quadrature in the velocity magnitudes and trapezoidal rule in the polar angles.The number of polar angles θ ∈ [0, 2π] is set equal to 160, while the numbers of c p and c z are set equal to 16 each, resulting in a total of 16 × 16 × 160 = 40,960 discrete velocities.Parallelization in velocity space is introduced.Also, to further reduce the computational e¨ort, grid re¦nement is applied in the physical space (r, z).Initially, the physical mesh is uniformly distributed with only 10 intervals per unit length in each direction.The simulations are performed on this grid, and after convergence, they are repeated on re¦ned meshes, where the number of intervals at each physical direction has been doubled.Four levels of re¦nement are introduced resulting in 80 × 80 nodes in the physical space.The discretized equations are solved in an iterative manner.The macroscopic quantities are assumed and then the kinetic Eqs. ( 1) are solved to deduce the reduced distribution functions, which are introduced into the moment Eqs.(2) to yield the updated values of the macroscopic distributions.The iteration process is terminated when the convergence criteria with t denoting the iteration index and K the number of nodes in the physical space, is ful¦lled, while the termination parameter is set to ε = 10 −9 .This discretization ensures grid independent results up to several signi¦cant ¦gures.The implemented algorithm has been extensively applied in previous works to solve with considerable success heat transfer con¦gurations [10] and nonlinear §ows through short tubes due to pressure and temperature gradients [3,11].

RESULTS AND DISCUSSION
The calculations have been carried out in the range of the rarefaction parameter δ 0 from 0 to 10, i. e., in the free molecular and transition regimes and for L/R = 1 and 5.It is noted that typically, the parameter Z varies from 1 to 5 and the choice of Z = 3 for the problem under question is reasonable.The presented results have been obtained for purely di¨use boundary conditions and the Hard Table 1 Dimensionless §ow rate W in terms of the rarefaction parameter and tube length to radius ratio  Sphere model with the upstream and downstream domains being 15 × 15 unit lengths.Tabulated results are presented for the dimensionless §ow rate, thrust, impulse factor, and discharge coe©cient as well as plotted results for the distribution of various macroscopic quantities.
In Table 1, the dimensionless §ow rate W for j = 0, 2, 3 (j = 0 refers to monoatomic gases) is given.Clearly, the e¨ect of the internal degrees of freedom on the gas §ow rate is very small for all values of the rarefaction parameter and for both L/R = 1 and 5.However, it is noted that for δ 0 = 1 and 10, W is decreased as j is increased.Also, W is increased as the length of the channel is decreased and the rarefaction parameter is increased.More speci¦cally, the §ow rate for δ 0 = [0.1,1] increases very slowly and then more rapidly for δ 0 = [1,10].Additional simulations have been performed with Z = 6 for δ 0 = 1 and L/R = 1 and 5 showing very small e¨ect on the §ow rates.
In Table 2, the variation of the dimensionless thrust F t in terms of the rarefaction parameter δ 0 and the ratio L/R is presented.The thrust is increased as the rarefaction parameter δ 0 is increased and the ratio L/R is decreased.It is clear that the propulsion e©ciently is increased as the tube length is decreased.Similarly to the §ow rates, the rotational degrees of freedom j and the parameter Z have a very small e¨ect on the values of F t as the rarefaction parameter δ 0 is increased.As expected, no e¨ect is observed in the free molecular regime (δ 0 = 0) since the transferred momentum in not a¨ected due to the absence Table 3 Dimensionless impulse factor Isp in terms of the rarefaction parameter and tube length to radius ratio   of the intermolecular collisions and the energy exchange between rotational and translational degrees of freedom.
In Table 3, the impulse factor I sp is shown.As the §ow becomes more rare¦ed, I sp is decreased.As the rarefaction parameter is increased, the increment of the rotational degrees of freedom leads to a very small increment of the impulse factor.This is well expected since the impulse factor is de¦ned as the ratio of the thrust over the §ow rate with the former one increasing and the latter one decreasing as j is increased.
In Table 4, the discharge coe©cient C d is presented.The discharge coe©cient C d decreases by increasing the tube ratio L/R, while for ¦xed L/R, C d is increased as δ 0 is increased.In addition, as the rotational degrees of freedom are increased from 0 to 2 and then to 3, the coe©cient C d is increased.This is due to the fact that the ratio of the speci¦c heats of the gas is decreased as j is increased, taking also into account that the §ow rates of the two types of gases are about the same.Overall, it may be concluded that the propulsion characteristics of polyatomic gas expansion through microtubes are slightly improved compared with the corresponding ones in the case of monoatomic gases.
In Fig. 2, the distributions of the Mach number along the symmetry axis r = 0 for L/R = 1 and r = 5 at δ 0 = 10 are shown.The Mach number far upstream is almost zero and is increased in the region just before the tube, while after the tube, it is rapidly decreased.It is seen that as the number of the internal degrees of freedom is increased, the Mach number is increased due to the decrease of the ratio of the speci¦c heat while the magnitude of the velocity vector in the two types of gas is almost the same.
In Fig. 3, the distributions of the dimensionless axial velocity, pressure, and temperature along the symmetry axis r = 0 for δ 0 = 0.1 and 10 with L/R = 1 are shown.In Fig. 4, the corresponding results for L/R = 5 are presented.Starting with the pressure variation, it is seen that far upstream, it is equal to one, then it rapidly decreases through the tube, and ¦nally, after the tube, it gradually approaches the far downstream conditions.As expected, the axial velocity has the same behavior with the Mach number.The maximum value of the velocity is increased as δ 0 is increased.The axial velocity and the pressure pro¦les in polyatomic gases (j = 2 and 3) are quantitatively very close to the corresponding pro¦les for monoatomic gases.The temperature equals unity in most of the domain, while inside the tube, it is decreased.The minimum value of the temperature distribution is decreased as the rarefaction of the gas is decreased and the ratio L/R is increased.In the case of polyatomic gases, the translational τ tr and total τ tot temperatures have the same qualitative behavior with the temperature of the monoatomic gas.The rotational temperature τ rot is maintained almost constant in the whole domain for small δ 0 , but as the rarefaction level of the gas is decreased, it is also decreased in the same way as the translational and total temperatures.It is noted that the increase of the dimensionless temperatures above one at the right open boundary for δ 0 = 10, both in Figs. 3 and 4, is contributed to the size of the downstream computational domain which is not adequately large.Increasing further the size of the downstream container will eliminate this behavior but it will greatly increase the computational e¨ort, while as shown in [12], it is expected to have Distributions of the dimensionless axial velocity and temperatures in the radial direction at the middle (z = L/(2R)) of the tube are shown in Fig. 5 for 2 ¡ 2; and 3 ¡ j = 3), and temperatures (c, 1 ¡ τ , j = 0; 2 ¡ τtr; 3 ¡ τrot; 4 ¡ τtot; black signs refer to j = 2 and blank signs to j = 3) for δ0 = 0.1 (left column) and δ0 = 10 (right column) with L/R = 5 along the symmetry axis δ 0 = 1 and 10 with L/R = 1.As expected, the velocities have a parabolic type shape with minimum and maximum values at the wall and at the center of the tube, respectively.The velocity pro¦les of diatomic gases (j = 2) are almost identical with the corresponding monoatomic pro¦les.The corre- sponding temperature pro¦les are also shown.In all cases, a temperature drop across the tube (radial direction) is observed.For δ 0 = 1, the translational temperature of a diatomic gas is close to the corresponding temperature of a monoatomic gas, while the rotational temperature is kept almost constant.For δ 0 = 10, the translational temperature of a diatomic gas is higher than the temperature of a monoatomic gas, while the rotational temperature is not constant anymore and it is reduced moving from the wall towards the center of the tube.

CONCLUDING REMARKS
The characteristic parameters of short tubes operating as propulsion systems in the case of polyatomic gases have been computed implementing kinetic modeling.Solving the Holway kinetic model subject to di¨use boundary conditions, the §ow rate, the thrust, the impulse factor, and the discharge coe©cient as well as the distributions of the macroscopic quantities with practical interest have been obtained.It has been found that the e¨ect of the rotational degrees of freedom on the macroscopic quantities is small except for the case of temperature distributions.It may be concluded that the overall propulsion ef-¦ciency in the case of polyatomic gases compared to the one in monoatomic gases is slightly improved.In addition, it has been demonstrated that polyatomic kinetic modeling may be applied as an alternative approach to DSMC.The implemented parametric study may be useful in the design of micronozzles operating under rare¦ed conditions with applications in cold gas propulsion systems [13,14].

Figure 1
Figure 1 Flow con¦guration and computational domain

Table 2
Dimensionless thrust Ft in terms of the rarefaction parameter and tube length to radius ratio

Table 4
Discharge coe©cient C d in terms of the rarefaction parameter and tube length to radius ratio