The NR Method for Polyatomic Gases
Abstract
In this paper, we propose a numerical regularized moment method to solve the Boltzmann equation with ESBGK collision term to simulate polyatomic gas flows. This method is an extension to the polyatomic case of the method proposed in [9], which is abbreviated as the NR method in [8]. Based on the form of the Maxwellian, the Laguerre polynomials of the internal energy parameter are used in the series expansion of the distribution function. We develop for polyatomic gases all the essential techniques needed in the NR method, including the efficient projection algorithm used in the numerical flux calculation, the regularization based on the Maxwellian iteration and the order of magnitude method, and the linearization of the regularization term for convenient numerical implementation. Meanwhile, the particular integrator in time for the ESBGK collision term is put forward. The shock tube simulations with Knudsen numbers from up to are presented to demonstrate the validity of our method. Moreover, the nitrogen shock structure problem is included in our numerical experiments for Mach numbers from to .
Keywords: polyatomic ESBGK model; moment method; NR method
1 Introduction
The kinetic theory has long been playing an important role in the rarefied gas dynamics. As a mesocopic theory standing between the fluid dynamics and the molecular dynamics, the kinetic theory is built on the basis of the Boltzmann equation, which uses a distribution function to give a statistical description of the distribution of microscopic particles’ velocities. In 1940s, Grad [11] proposed the idea using the Hermite expansion to approximate the distribution function, and a 13moment theory was given in detail in [11]. Recently, based on the idea of Grad, systems with large numbers of moments together with their numerical schemes are considered in [7, 9, 8], where some regularizations inspired by [21, 19] are also taken into account. In [8], the numerical regularized moment method is abbreviated as the NR method. However, all these works concentrate only on the monatomic gases, and in this paper, we will develop the NR method for the polyatomic case.
The study to apply the moment method to polyatomic gases can be traced back to McCormack [15], where a 17moment model was proposed. As far as we know, the most recent polyatomic extension of Grad’s 13moment equations is the work of Mallinger [14], whose system contains only 14 moments. In both models, a great amount of work is devoted to the deduction of the collision terms. In order to generalize the moment theory to large number of moments, we prefer a BGKlike simplified collision operator. As in the monatomic case, the simplest BGK model fails to give correct heat conduction, and for polyatomic gases, it also gives incorrect relaxation collision number, resulting in qualitative errors in temperatures compared with the Boltzmann equation [3]. Possible alternatives include the Rykov model [18] and the ESBGK model [4], which incorporate physical Prandtl number and relaxation collision number into the collision term. In this work, our investigation is restricted to the ESBGK model.
For polyatomic gases, besides the velocities of microscopic particles, an additional nonnegative ordinate representing the energy of internal degrees of freedom appears in the distribution function. Thus, in order to expand the distribution function into series, the basis functions are chosen as a combination of Hermite polynomials and Laguerre polynomials with proper translation and scaling based on the macroscopic velocity and translational and rotational temperatures of the gas. By considering the coefficients of the basis functions as moments, a system with infinite number of moment equations is derived. A moment closure is then followed to truncate the system with infinite equations and get a system with only finite equations. The framework for the moment closure is the same as [9]:

the Maxwellian iteration is applied to determine the order of magnitude for each moment;

by dropping higher order terms, the truncated moments are expressed by moments with lower orders;

for easier numerical implementation, the expression is linearized around a Maxwellian.
However, the details of the Maxwellian iteration are significantly different. In the polyatomic case, the iteration is much more complicated than the monatomic case because of the existence of both translational and rotational temperatures in the basis functions, and the process should be conducted carefully. Moreover, for the ESBGK model, analysis on the moments of the Gauss distribution also increases the complexity. Fortunately, the final result remains a similar form as simple as in [9].
As to the numerical scheme, the general framework in [8] is applicable. A split scheme is applied to divide the transportation part and the collision part, and the transportation part is processed by a finite volume method. Recalling that a special “projection” introduced in [7, 8] is required in the calculation of numerical fluxes, we further develop this technique to the polyatomic case in this paper. Meanwhile, the polyatomic ESBGK collision term can no longer be solved analytically as the BGK operator [7], and the CrankNicolson scheme is applied to ensure the unconditional numerical stability. Our numerical experiments show that our scheme correctly converges to the solution of the Boltzmann equation as the number of moments increases. The distinction between BGK and ESBGK models, together with the relation between monatomic and polyatomic cases, is illustrated by the numerical results of shock tube problems. Also, we apply the NR method to the nitrogen shock structure problem, and the results are comparable to the experimental data.
The rest of this paper is arranged as follows: in Section 2, a brief review of the polyatomic ESBGK model is given. In Section 3, the polyatomic NR method is introduced comprehensively, and in Section 4, a number of numerical experiments are carried out to validate our algorithm. As a summation, some concluding remarks are given in Section 5. Finally, some involved calculations are collected in the appendix for better readability to the body matter.
2 The ESBGK Boltzmann equation for polyatomic gases
The ESBGK model for polyatomic gases, which gives correct NavierStokes heat conduction compared with the BGK model, has been deduced in [4, 6]. The polyatomic ESBGK Boltzmann equation reads
(2.1) 
where denotes the molecule distribution, which is a positive function with and . The parameters , and stand for the time, spatial position and microscopic molecule velocity respectively, and is an internal energy parameter. In the right hand side of (2.1), is the Prandtl number, is the pressure, and denotes the viscosity coefficient. is a generalized Gaussian defined as
(2.2) 
Here is the total number of molecular internal degrees of freedom, and is the gas constant. The density and the macroscopic velocity are related to the distribution function through
(2.3) 
and is a relaxation temperature
(2.4) 
where is the relaxation collision number. For polyatomic gases, three temperatures are used frequently, including the translational temperature , the internal temperature , and the equilibrium temperature . They are defined by
(2.5)  
(2.6)  
(2.7) 
And the pressure is obtained from the ideal gas law:
(2.8) 
Now it only remains to define and :
(2.9) 
where
(2.10) 
and stands for the identity matrix.
3 The Nr method for polyatomic ESBGK equation
In this section, we are going to extend the NR method proposed in [7, 9] to the polyatomic case, which includes the following steps:

The distribution function is expanded into a series with specially selected basis functions.

A system with infinite number of moment equations is deduced.

The moment system is truncated at a certain place and made closed by regularization.

The regularization term is linearized in order to simplify the numerical implementation.

The numerical method is carried out following [8].
The details are introduced in the following five subsections.
3.1 Spectral representation of the velocity space
In the NR method for the monatomic gases, the Hermite polynomials have been employed to construct the basis functions of the velocity space, since Hermite polynomials are orthogonal over the region . For the polyatomic distribution function, since , we use the Laguerre polynomials, which are orthogonal over the region , as the basis functions in the ordinate . Thus the basis function has the following form:
(3.1) 
where is a multiindex, and
(3.2)  
(3.3)  
(3.4) 
Some properties of the Laguerre polynomials and the Hermite polynomials can be found in Appendix A. With equation (3.1), the distribution function is expanded as
(3.5) 
Let us consider the general case when , and have no relation with the distribution function . Using the orthogonality of the Laguerre and Hermite polynomials, we can deduce that
(3.6a)  
(3.6b)  
(3.6c)  
(3.6d) 
If is the macroscopic velocity and , are the translational and internal temperatures for the distribution , using (2.3), (2.5), (2.6) and (3.6), we conclude
(3.7) 
If (3.7) is satisfied, then (3.5) is called as a normal representation of . If (3.5) is not a normal representation, then the density, momentum and translational and internal energies can be easily calculated through (3.6). For a normal representation, we have
(3.8) 
where is defined in (2.10).
3.2 The moment equations for the ESBGK model
In this section, we are going to derive equations for the moment set . The general strategy is to substitute (3.5) into (2.1), and then match the coefficients of the same basis functions. For the left hand side of (3.5), the process is similar as that in [9], and the detailed derivation can be found in Appendix B. Suppose has the following expansion:
(3.9) 
Then the analytical expressions of the moment equations are obtained as
(3.10) 
where is taken as zero when or any of the components of is negative, and is defined in (3.2).
Now we focus on the relation between and . The expression of (2.2)—(2.10) and the equalities under normal representation (3.7) and (3.8) show that are functions of , , and with . Thus the system (3.10) is closed for and , which means only the expressions of and are needed. The following results are trivial:
(3.11a)  
(3.11b)  
(3.11c) 
where (3.11c) comes from the physical meaning of the relaxation collision number. Moreover, since has the form
(3.12) 
where
(3.13) 
and the basis functions have the same structure (see (B.1))
(3.14) 
we can conclude
(3.15) 
where is independent of . From (3.11a) and (3.11c), we find
(3.16) 
Until now, what remains is to work out the expressions of for . This needs some involved calculation with details in Appendix C. The final result is in a recursive form as
(3.17) 
where such that , and is taken as zero when . With (3.11b), one can easily observe that when is odd.
As the end of this subsection, we give the equations for the velocity and the translational and internal temperatures , . The equation for can be obtained by substituting with in (3.10), and the result is
(3.18) 
The equation for can be obtained by substituting with , , , and then summing up all the three equations. The result is
(3.19) 
where
(3.20) 
Similarly, if we set and , then we have
(3.21) 
3.3 Truncation and closure with regularization
Since the moment system (3.10) contains an infinite number of equations and cannot be used for computation, we need to choose a finite set from them as the governing equations of our method. However, due to the existence of the last term in the second line of (3.10), the resulting moment system will be unclosed, which leads to the “closure problem” for the moment method.
We first consider the truncation of the spectral expansion. In general, we can choose two nonnegative integers and , and use the moment set as the finite subset. Such choice well retains the Galilean invariance since and only couple with each other in the collision term. We will postpone the discussion of the relation between and . Below we use to denote the index set such that the set contains all the moments appearing in the final moment equations.
Now we are going to make the system closed. The simplest way is to set in (3.10) if , which leads to the Gradtype moment equations for polyatomic gases. Mallinger’s work [14] has generalize the Grad moment equations to the polyatomic case. Here, we follow [17, 19, 9] and use the Maxwellian iteration together with the order of magnitude method to give a more reasonable closure. The procedure of Maxwellian iteration is constructed by rearranging (3.10) and adding superscripts representing the iteration steps:
(3.22) 
where is considered as a small parameter and satisfies
(3.23) 
In the first line of (3.22), and are defined as
(3.24)  
(3.25) 
The reason why the iteration scheme for is special is that are functions of , . Note that (), , () and () remain invariant during the iteration, and according to (3.19) and (3.21), and evolves as follows:
(3.26)  
(3.27) 
where also remains invariant during the iteration. The iteration starts with
(3.28) 
In order to simplify the notation, we define the following vectors:
(3.29) 
and rewrite (3.22) and (3.26) as
(3.30)  
(3.31)  
(3.32) 
where is a linear operator, and is a vector of linear operators. Each of ’s components has the following form:
(3.33) 
where , which can be considered as a “constant” during the iteration. is a vector, whose components can be expressed as
(3.34) 
where are also constants. Now we are ready to carry out the Maxwellian iteration.
The first step of iteration In the first step, the formulae for the translational and internal temperatures can be written as
(3.35) 
where and . It is easy to find
(3.36) 
Thus
(3.37) 
Now, it can be easily deduced from (3.30) that is nonzero if and only if , or , . Precisely, we have
(3.38) 
and
(3.39) 
where and has an order of magnitude . The meaning of the subscripts of is the same as .
The general progress In general case, we have
(3.43) 
and
(3.44)  
(3.45)  
(3.46)  