Boltzmann collision operator for polyatomic gases in agreement with experimental data and DSMC method

This paper is concerned with the Boltzmann equation based on a continuous internal energy variable to model polyatomic gases with constant specific heats. We propose a family of models for the collision kernel and evaluate the nonlinear Boltzmann collision operator to get explicit expressions for transport coefficients like shear and bulk viscosities, thermal conductivity, depending on the collision kernel parameters. This model is shown to contain as a special case the collision kernel used in the direct simulation Monte Carlo method with the variable hard sphere cross section. Then, we show that it is possible to choose parameters in such a way that we recover various physical phenomena, in particular, experimental data for the shear viscosity, Prandtl number and the ratio of bulk and shear viscosities at the same time.


Introduction
In the present paper, we study kinetic modeling of polyatomic gases. The focus is on the Boltzmann equation with continuous internal energy variable introduced in [1][2][3]. The approach is based on a specific parametrization of the molecular collision process which makes the model accessible both from the rigorous analytical and computational points of view. We examine in detail the idea of this parametrization by studying the Larsen-Borgnakke DSMC procedure to calculate internal energy exchanges during collisions [4]. Our attention is restricted to calorically perfect gases that have the constant specific heat with respect to temperature (sometimes called polytropic). In applications, this includes rarefied gas flows of linear molecules at room temperature. The polyatomic assumption influences the form of the Boltzmann collision operator which contains a weight factor depending on a parameter interpreted as number of molecular internal degrees of freedom. For such a model, the rigorous existence and uniqueness theory were established in [5] in the case of space homogeneity.
In the case of a single monoatomic gas, the Boltzmann equation can be used to derive continuum models able to capture strong non-equilibrium processes in rarefied gas flows, improving classical models such as the Navier-Stokes-Fourier system. This goal can be achieved via the moment method [6,7] or within extended thermodynamics [8]. Although polyatomic gas models are highly relevant in applications, there are only partial results of continuum models derived from the Boltzmann equation that address moment equations of some particular order, see [9] and references therein. The system of 14-moment equations with two hierarchies, momentum-and energy-like, was firstly presented for dense gases in [10]. The same structure was derived in [11] for rarefied gases starting from the Boltzmann model. The hierarchy of moment equations for polyatomic gases was analyzed by the order-of-magnitude method in [12] where also kinetic boundary conditions are investigated. When restricted to non-equilibrium processes in which shear stresses and heat conduction are neglected, the 6-moment model evokes and is one of the rare systems that admit an exact entropy principle [13]. The next goal was to obtain physically relevant transport coefficients by evaluating the nonlinear collision operator or computing production terms. The main question that arises is the modeling of cross section or collision kernel. First results were concerned with the constant cross section or Maxwell molecules [14], the relative speed alone [11,15] or the relative speed combined multiplicatively with the internal energy [16]. They are all of limited agreement with theoretical restrictions and experimental data. These results were improved in [17,18] using the collision kernel proposed in the rigorous analysis [5]. In particular, [17] shows an agreement with the theoretical value of the Prandtl number given by Eucken's formula [19] with an error of up to 12%. This error was eliminated in [18] by studying 17-moment equations and by introducing frozen collisions that do not change microscopic internal energies and consequently enriching the collision kernel from [5] with additional free parameters.
In this paper, we focus on the experimental data for transport coefficients and the Prandtl number. As a matter of fact, experimental data for the Prandtl number are even lower than data obtained from the Eucken formula and for some calorically perfect gases cannot be covered by the collision kernel model used in [18]. This is the motivation to propose in this paper the extended collision kernel model capable to recover, at the same time, experimental data for the shear viscosity exponent and the Prandtl number as much as the bulk viscosity in agreement with [20]. We also show that this extended model contains the cross-sectional model frequently used in the DSMC method [4,21].
The rest of the paper is organized as follows. Section 2 introduces calorically perfect gases for which the Boltzmann equation is described in Sect. 4. The microscopic dynamics of collisions of binary polyatomic molecules is studied in Sect. 3 bringing the comparison with Larsen-Borgnakke DSMC procedure as well. The macroscopic system of 17-moment equations and transport coefficients are derived in Sects. 5 and 6, respectively. We then propose the model for the collision kernel in Sect. 7 for which we firstly show agreement with the model frequently used in DSMC method in Sect. 8. The next goal of Sect. 9 is to provide the fit of collision kernel parameters such that for certain gases that behave as calorically perfect over a significant temperature range (around room temperature) the experimentally measured shear viscosity exponent is recovered as much as various formulas for the Prandtl number, including the experimentally observed value. Additionally, the agreement with the bulk viscosity is provided.

Physical background: calorically perfect gases
Thermally perfect gases (sometimes called non-polytropic) are ideal gases that obey the thermal equation of state and have temperature-dependent specific heats [22]. The hydrodynamic pressure p is connected with temperature T by where ρ is the gas mass density, m the molecule mass and k the Boltzmann constant k = 1.38 · 10 −23 J/K. The specific internal energy, denoted by e, depends on temperature T through its relation to the specific heat at the constant volume here expressed in the units of kg·J/K. The specific heat c v (T ) is an experimentally measured quantity, and its explicit functional dependence on T can be withdrawn from [23,24].
A special case is the calorically perfect gas, also called polytropic, for which c v (T ) = const and (2) reduces to where we took e 0 = 0 [25]. This paper is concerned with calorically perfect gases, and in our setting we introduce the number of internal degrees of freedom δ > 0, which is a constant, but not necessarily an integer, related to c v as follows Therefore, in our notation (3) reads Note that in some literature [9] a parameter α > −1 is used instead of δ, and they are connected through δ = 2(α + 1).

Study of collisions
At the microscopic level, it is known that the kinetic energy of a polyatomic molecules pair is not conserved during their collision. The idea of the continuous kinetic approach [1][2][3] is to associate a polyatomic molecule, apart from the velocity v ∈ R 3 , also the microscopic internal energy I ≥ 0 aimed to capture all peculiarities related to the collision. Consider two colliding polyatomic molecules of mass m with the velocity-internal energy (v, I ) and (v * , I * ) that change to (v , I ) and (v * , I * ) due to the collision. During this interaction, momentum and total energy are conserved mv + mv * = mv + mv * , Defining the center-of-mass frame of reference via the center of mass velocity and relative velocity, laws (6) can be rewritten These equations are parametrized with a spherical angle σ ∈ S 2 and two dimensionless collision weights r, R ∈ [0, 1] according to the procedure explained in detail in [1], for instance. Roughly speaking, the idea is to split the total energy E into translational energy of the relative motion of the colliding particles and total internal energy using the weight R ∈ [0, 1] and then to associate those energies to particles with the parameters σ ∈ S 2 and r ∈ [0, 1], respectively. The collision rules read Conservation laws (6) allow to define primed parameters enabling to see pre-post change of variables as an involution. For later convenience, we introduce translational and internal energies scaled with the total energy as follows Parametrization (9)-(10) yields In the next section, we explain the Larsen-Borgnakke DSMC algorithm for post-collision energy sampling.

Larsen-Borgnakke DSMC procedure of collision energy exchange
Within the direct simulation Monte Carlo (DSMC) method [21], perhaps the most widely used procedure to calculate internal energy exchanges is the Larsen-Borgnakke model [4].
A brief overview of this procedure is presented below in the current setting of a single internal energy variable [26]. Let δ > 0 be the number of particles' internal degrees of freedom, and let ξ tr ≥ 3 be the effective number of translation degrees of freedom, which is dependent on the specific choice of the total cross-sectional model [21]. The procedure is explained in terms of an algorithm for sampling post-collisional energies of the particles.
In the case, an internal energy exchange does occur; the total collisional energy E from (8) is computed and then re-distributed in accordance with the number of degrees of freedom. With the notation (12), firstly the fraction of the total energy E is assigned to the internal energy I of the first of the two particles by sampling i from the Beta distribution Beta δ 2 ; δ+ξ tr 2 . After the first particle's internal energy has been chosen, the remaining collisional energy is used to get the second particle's internal energy. Namely, sample φ from the Beta distribution φ ∼ Beta δ 2 ; ξ tr 2 , and set i * = (1 − i )φ. Finally, the remaining energy e tr = 1 − i − i * is the translation energy of the relative motion of the particles scaled with E. This procedure will be compared in detail to the Boltzmann collision operator in Sect. 8.

Kinetic model: the Boltzmann equation
In the continuous kinetic setting, the behavior of polyatomic gases is described by the distribution function f := f (t, x, v, I ) ≥ 0 which depends on macroscopic variables, time t ≥ 0 and space position x ∈ R 3 and on microscopic variables, velocity v ∈ R 3 and internal energy I ∈ [0, ∞). This internal energy phase space can be replaced by internal state phase space, as explained in [18]. If the distribution functions change due to particles' collisions, an efficient way to describe its evolution is to write the Boltzmann equation where Q( f, f )(v, I ) is the collision operator. Its structure is influenced by the Jacobian of microscopic prepost collision transformation (15) and the physical requirement on the model to capture macroscopic behavior (5). The Jacobian of the transformation (v, v * , I, I * , r, σ, R) → (v , v * , I , I * , r , σ , R ) defined by relations (10)- (11) is computed in [3,5] where the second equality is obtained by taking the ratio of |u | = √ R E/m and its counterpart |u| = R E/m. We take the form of the collision operator introduced in [1] and used in the rigorous mathematical theory [5], as well as for moment equations of various order and computational purposes [17,18]. Here, standard abbreviations are used f := f (v , I ), , and the primed quantities are defined in (10). The collision kernel B = B(v, v * , I, I * , r, σ, R) ≥ 0 satisfies microreversibility assumptions corresponding to the pre-post change and an interchange of colliding molecules The form of H δ is influenced by the weighted factor in the collision operator (16) and Jacobian (15) which ensures the microreversibility of the model, as explained below.

Microreversibility
In the continuous kinetic approach [1], it is shown that the microreversibility of the collision kernel (17) together with the Jacobian (15) yields the following invariance written in terms of parameters r and R which can be seen as a translation of the (second) quantic Fermi's golden rule, in the spirit of [27]. This notion can be related to the one used in DSMC method that does not deal with r, R-parametrization. Instead on the collision kernel, the focus is on the cross section E, I, I * → I , I * representing the cross section of a collision of particles with the total collision energy E = m 4 |u| 2 + I + I * as in (8) and precollisional internal energies I , I * that results in a post-collisional internal energies I , I * . Namely, the following microreversibility condition on the cross section is imposed [28] |u| (I I * ) δ/2−1 E, I, I * → I , I * dσ dI * dI * dI dI dv * dv = The two concepts can be related. Indeed, for the sake of simplicity, instead of (I , I * ) we can consider scaled energies (i , i * ) defined in (13) and then observe that the change of variables (i , i * ) → (r, R) has the Jacobian Consequently, the collision kernel B and the cross section are related via This relation is used in Sect. 8 in the concrete example of the collision kernel and cross section.

Collision operator weak form and equilibrium
The invariance property (19) ensures a well defined collision operator weak form. Namely, for any suitable test function ψ(v, I ), This weak form provides conservation laws at the macroscopic level from the Boltzmann equation and the H-theorem with the equilibrium distribution function [1] where c = v − U is peculiar velocity with respect to the macroscopic velocity U and is the gamma function. This equilibrium distribution gives for the total energy of the gas consistent to a calorically perfect gas as introduced in (5).

Moment equations
In this paper, we consider 17-moment equations derived from the Boltzmann equation, namely the system of equations for equilibrium variables: mass density ρ, gas velocity U , energy density ρe that satisfy conservation laws, and non-equilibrium variables: the trace of the pressure tensor P and its deviatoric part σ , also called stress, and translational and internal heat fluxes q and s, ⎛ with c = v − U . As usual in kinetic gas theory, we assume (25) to hold also for non-equilibirum distributions and define the temperature of the gas through it. The internal energy density can be split into translational and internal parts, ρe = ρe tr + ρe int which gives rise to define a translational and internal temperature by such that δ+3 The non-equilibrium pressure P can be related either to the sum of equilibrium pressure p in (1) and dynamical pressure , or equivalently to the translational part of the internal energy density ρe tr , that is, translational temperature T tr [29]. More precisely, Using (25) in non-equilibirum, the last equation of (28) implies In this paper, we choose to express P in terms of pressures p and . The equation for p is derived from the equation for the conserved energy ρe, while the equation for follows from the one for P.
Thus, integrating the Boltzmann equation (14) against test functions from (26) yields moment equations for the corresponding densities where the non-convective fluxes are defined as ⎛ while the production terms read ⎛ and One possibility to close this system is to provide an approximative distribution function and to specify the collision kernel, which is the standard procedure used in moment method [7] or extended thermodynamics [9].
In this paper, we choose the distribution function as in [18], where f M is the local Maxwellian of the shape (24), yielding the following expressions for the unknown fluxes (32), Moreover, for the choice of the distribution function (34), the linearized production terms (33) are of the form [18] Explicit expressions for coefficients P (0) , P s are provided by evaluating the collision operator for a certain choice of the collision kernel. Computations are implemented in computer algebra software [30] and are used in this paper as well.

Transport coefficients and Prandtl number
The standard procedure of Maxwellian iteration [6,9] for moment Eq. (29) provides expressions for transport coefficients-shear viscosity μ, bulk viscosity ν and thermal conductivity κ [18], in terms of coefficients of the production terms (35), In the polytropic regime (4), Prandtl number is defined by while the Eucken formula reduces to Thus, the evaluation of the collision operator yielding explicitly computed production terms (35) enables to formulate the model function for the Prandtl number by using its definition (37) and models for the transport coefficients (36) This model function depends on the choice of the collision kernel. The idea of our approach is to propose a model for the collision kernel enriched with parameters such that it is possible to simultaneously recover experimental data for the shear viscosity dependence upon temperature, experimentally measured value of the Prandtl number and expected behavior of the ratio between bulk and shear viscosities. Experimental data are withdrawn from the database [21,23,31,32] for gases that behave as calorically perfect on certain temperature range around room temperature.

Collision kernel models
In Sect. 7.1, we firstly analyze the collision kernel proposed in the rigorous mathematical analysis [5], already used in a modified version in [18]. As this model provides a limited agreement with experimental data, we propose its improvement in Sect. 7.2.

Initial model
This model was originally proposed in [5] as an example that satisfies the imposed assumption on the collision kernel allowing to solve the Cauchy problem for the space homogeneous Boltzmann Eqs. (14)- (16). It was used firstly in [17] to compute transport coefficients from 14-moment system and production terms for the pressure tensor and the (total) heat flux, which led to the Prandtl number showing the error of about 10% with respect to the Eucken formula for some choice of gases. The agreement with the Eucken formula is reached in [18] with the model of 17 moments implying essentially different formula for the heat conductivity and by introducing frozen collisions. Frozen collisions are the ones that do not change particles' internal energies, i.e., I = I and I * = I * , also known as elastic in the DSMC community. In [18], it is proposed to convexly split the collision operator (16) into a pure polyatomic (non-frozen) and a frozen operator using the parameter ω ∈ [0, 1], or equivalently in (16) to consider the following model for the collision kernel, where u = v − v * is the relative velocity,û is its normalized vector and E is the total energy (8). The angular part of the collision kernel is taken to be constant, i.e., b(û · σ ) = K , following the variable-hard-sphere (VHS) model [21]. The constant serves to achieve consistency of the pure polyatomic collision operator (16) with the monoatomic limit. The parameter ζ ≥ 0 is related to the shear viscosity exponent, as it will be shown in the upcoming sections, while the parameters η ≥ 0 and η f ≥ 0 control the influence of internal energies during collision.
In [18], this model is shown to provide agreement of the Prandtl number (37) with the Eucken formula (38) for the following choice of the collision kernel parameters Since the value of ζ = 0 is nonphysical, in [18] a generalized Eucken formula for the Prandtl number is proposed, referred to as the frozen formula, Obviously, it reduces to the Eucken formula (38) for ζ = 0. However, this model is of limited agreement with experimentally observed values of the Prandtl number, which motivates us to improve it further.

Extended model
The frozen formula for the Prandtl number (42) is the lower limit of possible Pr-numbers in the initial model (40). However, it is known that experimental data are below this curve, which makes the model (40) inadequate and motivates us to propose its improvement as follows. With dimensionless energies introduced in (12) and for the parameters ω, ζ and K δ defined as in the initial model (40), and additionally forζ ,ζ f ≥ 0 andη,η f ≥ −1/2, we propose the following model The crucial improvement relies in the new set of parameters and the construction which allows negative influence of the internal energy, while keeping the collision kernel B still non-negative. Indeed, as 0 ≤ i, i * ≤ 1, ζ ≥ 0 and r, R ∈ [0, 1], and ifη is negative, i.e., −1/2 ≤η ≤ 0, the fist term in (43) is still non-negative, The similar computation can be done for the second term implying the non-negativeness of the collision kernel B.
Note that the extended model (43) reduces to the initial model (40) for the notation (12) and the following choice of parametersζ and by restricting the range ofη,η f ≥ 0.

Comparison with collision kernel models used in DSMC
The DSMC method also uses the concept of frozen collisions. More precisely, according to the Larsen-Borgnakke model, when two particles are selected for a collision, an internal energy relaxation event occurs with a probability of p int ; otherwise, the collision is frozen (or elastic) and the internal energy states of the particles remain unchanged. This probability p int of an internal energy exchange occurring during a collision can be related to the macroscopic collision number Z int , defined as the ratio of the times τ int (time it takes for the system to reach equilibrium between internal and translational modes) and τ tr (mean collision time). Namely, it was shown in [33] that in order to reproduce a macroscopic collision number Z int , p int should be defined as where the number of internal degrees of freedom δ > 0 and the effective number of translation degrees of freedom ξ tr ≥ 3 are introduced in Sect. 3.2. Thus, approximately every 1/ p int collisions, an internal energy exchange occurs. Note that the correct choice of p int ensures that the Larsen-Borgnakke model can adequately reproduce the macroscopic relaxation rates of the internal degrees of freedom [28] , although it does not capture state-specific transition probabilities. The next goal is to rewrite the cross section in the DSMC formulation and to translate it to the continuous model of the present paper. The total cross section is the sum of the frozen and non-frozen one, weighted by the corresponding probabilities defined by (44) Here n-f E, I, I * → I , I * is the cross section of the transition of internal energies I and I * to I and I * given a total collisional energy E. Thus, the first part corresponds to frozen collisions with no internal energy exchange (which occur with a probability of 1 − p int ) and the second one corresponds to non-frozen collisions (which occur with a probability of p int ).
We translate the Larsen-Borgnakke DSMC procedure to our formulation. On one side, the internal energy cross section n-f E, I, I * → I , I * is usually assumed [21] to be related to the probability of obtaining post-collisional energies I , I * given the total collision energy E. Note that the reference to energies I, I * is irrelevant, since the sampling steps presented in Sect. 3.2 indicate that the Larsen-Borgnakke model does not explicitly depend on I and I * -only the total collision energy E affects the post-collisional values I , I * . The exact relation can be obtained by taking into account scaled energies i , i * and e tr = 1 − i − i * introduced in (12) and sampling procedure presented in Sect. 3.2. Indeed, if the scaled energies are understood as random variables, then they follow multi-variable Beta distribution with the probability density function where is the usual Gamma function. The random variables describing internal energy i of the first particle, internal energy i * of the second particle and translational energy of their relative motion and Beta[ ξ tr 2 , δ], respectively. The internal energy cross section is defined as n-f E, I, I * → I , where the shape ofˆ n-f (|u|) will be assumed later.
In the special case of δ = 2, (46) can be further simplified tõ with e tr being the scaled post-collision translational energy (12). It can be seen that in the case of δ = 2, the probability does not depend explicitly on the internal energies and is a function of only the total pre-collision energy and post-collision translational energy.
On the other side, in the DSMC method, when working with the analytical expressions for the cross section, the usual assumption is to model both f andˆ n-f with the standard variable hard sphere (VHS) model [21], with possibly two different constants C f VHS and C VHS , respectively, where s visc is the viscosity exponent. We thus have that the number of effective translational degrees of freedom [21] is From expressions (45), (47) and (49), we can construct the cross section based on the DSMC model satisfying microreversibility assumption (20) and relate it to the models presented in the previous section. We find In order to be consistent with the notation used in the previous section, we note that the parameters of the present DSMC model are related to the parameters of the models introduced in Sect. 7 as follows The collision kernel in DSMC is related to the cross section via which in our case reads, specifying (51), The aim is to translate this collision kernel used in the DSMC method into the collision kernel used in the continuous approach. For that, the non-frozen part is rewritten using parametrization (13), the change of measure (21) and taking Value of δ connected with the experimentally measuredĉ v at 300K as given in (4). Value of ζ according to the fit with experimental data for the shear viscosity μ according to (57) and the reference value μ 0 of the shear viscosity in the units of µPa·s at temperature 300K and pressure p = 0.092bar which is equal to K in the monoatomic case corresponding to the limit δ = 0. Here K is the constant as defined in (40). Indeed, Now we can translate the non-frozen part of the DSMC model (53) into the model used in the continuous approach introduced in Sect. 7, which corresponds to both models (40) and (43) for η = 0,η = 0, respectively, and is in consistency to the relation (22). The transformation of the frozen part of the DSMC model is done by analogy. Therefore, we can conclude that the Larsen-Borgnakke model is more restrictive, as it does not account neither for the influence of the pre-collisional internal energies I , I * on the post-collisional internal energies in the case of non-frozen collisions, nor for the potential influence of the pre-collisional internal energies I , I * on the frozen scattering cross section, which they might have via influence on molecular diameters [34].

Collision kernel parameters' fit with experimental data
The goal of this section is to fit the parameters of the proposed collision kernel model (43) in order to recover experimentally observed shear viscosity exponent and the value of the Prandtl number, as much as the ratio of bulk and shear viscosities in agreement with [20], for specific gases that behave as calorically perfect over a significant temperature range, also listed in [24]. Partial results are shown for the collision kernel model (40) also studied in [18].

Polytropic gases
Looking at experimental data in the section of gas phase thermochemistry data in [23] for c p at constant pressure of 1 bar in the units of J/(mol K), and transferring it to the dimensionless formĉ p = c p /R, where R = 8.314 J/(mol K) is the universal gas constant, we get experimental data forĉ v =ĉ p − 1. We consider a gas to be calorically perfect if its relative change is below 5% over a significant temperature range that starts with the room temperature of 300K. Gases are listed in Table 1.

Shear viscosity exponent
It is known that the shear viscosity μ dependence on temperature T can be modeled as a power law function [35],   Table 2 that correspond to the viscosity exponent s visc from Table 1 (color figure online) for some reference temperature T 0 and viscosity μ 0 . In our analysis, we take T 0 = 300K, and for μ 0 we take the value of μ in the units of µPa·s at temperature 300K and pressure p = 0.092bar. The model of μ computed via production terms as described in (36) has the same functional form as (56), making possible the relation between collision kernel parameter ζ and the viscosity exponent s visc , namely Taking experimental data for μ [23,24], and performing the fit with the model function (56) give the value of parameter ζ using (57). Results are shown in Table 1. It is worthwhile to mention that results are in a good agreement with the ones from [35],

Prandtl number
For a specific choice of the gas, parameters δ and ζ are fixed. The wealth of our collision kernel models relies on the fact that they provide additional parameters, namely ω, η, η f in the initial model (40) and ω,η,η f ,ζ ,ζ f in the extended model (43) that can be fit in order to capture some additional physical phenomena. Let us first consider the Prandtl number. The model function for the Prandtl number is explained in (39) and depends on the parameters of the collision kernel. The goal of this section is to verify validity of the collision kernel by checking the agreement of the Prandtl number model function with the given data for the Prandtl number. We consider three possibilities: the experimentally measured value, Eucken, and frozen formulas.
The experimentally measured value of the Prandtl number can be obtained by using its definition (37) and experimental data for the heat conductivity κ at T = 300K, together with δ and μ 0 from Table 1.  The Eucken formula for the Prandtl number in our setting is given by relation (38) and is fully determined only by specifying δ.
Additionally, we take into account the generalization of the Eucken formula proposed in [18] for the purely frozen case obtained by taking ω = η f = 0 in the initial model (40) and ω =η f = 0 in the extended model (43). This frozen formula for the Prandtl number depends on δ and ζ and provides the same value for both models. The crucial difference between the two models is that the frozen formula is the limiting curve (or the lower bound) for the initial model (40), while the extended model covers the area even below this curve thanks to the specific range of parameters related to the internal energy. The difference is clearly observed in Figs. 1 and 2, where the frozen curve is presented as blue dashed line.
These three values of the Prandtl number-experimentally observed, according to Eucken and frozen formulas-are presented in Table 2 for calorically perfect gases from Table 1. The next goal is to fit parameters of the collision kernel models (40) and (43) in order to recover these given data for the Prandtl number, taking its model function (39) that depends on the parameter of the collision kernel model.
The initial model (40) has been already studied in [18]. It is shown that the agreement with the Eucken formula (38) can be attained and that the frozen formula (42) is the lowest possible set of values for the Prandtl number reachable by this model. Here, we illustrate the model in Fig. 1. We plot the Prandtl number (37) after selecting various values of the parameter s visc from Table 1 representing different gases. The possible range, shown as yellow background, is obtained by varying the free model parameters 0 ≤ ω ≤ 1 and 0 ≤ η, η f ≤ 2, after, for the simplicity of the presentation, we fix the degrees of freedom δ = 2.
Additionally, we consider experimental data for the Prandtl number. Already from Table 2 one can notice that the experimental value of the Prandtl number is sometimes lower than the Eucken formula, and even below the frozen curve for H 2 or outside the (yellow) region of reachability of the initial model. This is also displayed in Fig. 1.
Agreement with all considered measured value of the Prandtl number is attained for the extended model (43). Namely, for the given δ and viscosity exponent s visc as in Table 1, it is possible to adjust free parameters ω,η,η f ,ζ andζ f such that the experimental data from Table 2 are recovered.
The key point lies in additional parameters,η andη f ,ζ andζ f , allowing to extend the range of attainable values of the Prandtl number even bellow the frozen curve. This is possible becauseη andη f can take the negative vales, i.e.,η,η f ≥ −1/2, and moreover the model function for the Prandtl number has non-monotone behavior with respect toζ andζ f for other parameters fixed. For the fixed δ and s visc , the lowest possible value of the Prandtl number is reached for the frozen case ω = 0, the lowest influence of the internal energy expressed asη f = −1/2, and then by minimizing the Prandtl number model function (39) with respect tô ζ f . The illustration is provided in Fig. 2. The range of Prandtl number, shown as the yellow background, is   Table 4. Blue points represent experimental values for the pair (ratio of bulk and shear viscosities ν μ , Prandtl number Pr) (color figure online) obtained by varying the free parameters of the model 0 ≤ ω ≤ 1, − 1 2 ≤η,η f ≤ 2 and 0 ≤ζ ,ζ f ≤ 2, after fixing δ = 2 for the sake of simplicity. In particular, the limiting points of the model are shown as orange dots.

Bulk viscosity
Apart from the agreement with the experimental value of Prandtl number in Sect. 9.3 provided by the extended model for the collision kernel (43), our next aim is to study the ratio of the bulk and shear viscosities. As explained in (36), its model can be obtained from the evaluation of the Boltzmann collision operator ν μ We consider the extended model (43), and therefore, the ratio (58) depends on the collision kernel parameters ω,η,η f ,ζ ,ζ f , apart from δ and ζ (or s visc by (57)). The goal is to fit these parameters in order to simultaneously recover experimental data for the Prandtl number and the ratio of the bulk and shear viscosities, presented in Tables 2 and 3, respectively. The fit can be obtained firstly by fixing δ and ζ as in Table 1, and then by solving the system of equations for remaining parameters ω,η,η f ,ζ andζ f requiring (i) the model function for the Prandtl number (39) to be equal to its experimental value and (ii) the model function for the ratio (58) to coincide with its experimental value.
This system does not have a unique solution and one possible choice is given in Table 4. First, values for η,ζ andζ f are arbitrary chosen. Then the system is solved numerically for the parameters ω andη f using Mathematica software. Let us highlight that the freedom in the choice of parametersη,ζ , andζ f could be used to match additional physical transport properties.
It is worthwhile to mention that beyond recovering experimental data for the ratio of bulk and shear viscosities, extended model also provides a wide range for its values. In Fig. 3, the curves represent parametric plots of the Prandtl number that corresponds to the viscosity ratio after varying 0 ≤ ω ≤ 1 and setting other Viscosity exponent s visc and number of internal degrees of freedom δ are given in Table 1 parameters to be as given in Table 4. Additionally, blue points represent experimental values for the pair bulk and shear viscosities ratio-Prandtl number, recovered for certain ω. Note that for ω → 0 the viscosity ratio has the tendency to take larger values. The reason for this behavior lies in the fact that P (0) is proportional to ω, since the frozen contribution to this production term is zero [18].

Conclusions
A new collision kernel for polyatomic gases with continuous internal energy has been proposed. It is shown that it contains as a limiting case the VHS/Larsen-Borgnakke model, widely used in the DSMC community to model internal energy transfer. The system of 17-moment equations is used to derive expressions for transport properties of interest, namely shear viscosity, Prandtl number, and ratio of bulk to shear viscosities. The proposed collision kernel is used to evaluate the transport properties in the polytropic regime for several polyatomic gases, and it shown to be able to reproduce all relevant parameters simultaneously.