Reduced kinetic model of polyatomic gases

The dynamics of a dilute monoatomic gas in terms of the single-particle distribution function is described by the Boltzmann equation (Chapman & Cowling, 1970; Cercignani, 1988). Unlike the continuum Navier–Stokes–Fourier hydrodynamics equation, the Boltzmann equation is a valid description even at highly non-equilibrium states (Mott-Smith, 1951; Liepmann et al., 1962; Cercignani, 1988), encountered in the presence of strong shock waves at high Mach number (ratio of flow speed to sound speed) and in a highly rarefied flow characterized by a large Knudsen number (ratio of the mean free path to characteristic length scale) (Oh et al., 1997; Struchtrup, 2004; Ansumali et al., 2007b). However, any analysis of the integro-differential Boltzmann equation is a formidable task even for the simplest problems. Thus, one often models the Boltzmann dynamics via a simplified collision term that converts the evolution equation to a partial differential equation (Bhatnagar et al., 1954; Lebowitz et al., 1960; Holway Jr, 1966; Shakhov, 1968; Gorban & Karlin, 1994; Levermore, 1996; Andries et al., 2000; Ansumali et al., 2007a; Lebowitz et al., 1960; Singh & Ansumali, 2015; Agrawal et al., 2020). An important example is the BGK model (Bhatnagar et al. 1954), which states that the relaxation of the distribution function towards the Maxwell–Boltzmann (MB) form happens in a time scale corresponding to the mean free time $\tau$ with the assumption that every moment of the distribution function relaxes at the same rate. The BGK model is quite successful in replicating qualitative features of the Boltzmann dynamics (collisional invariants, the zero of the collision, H-theorem, conservation laws, etc). However, the BGK model predicts the Prandtl number of the fluid to be unity, while the value predicted by the Boltzmann equation for monoatomic gas is $2/3$. Thus, several other variations of the collision model such as ES–BGK model (Holway Jr, 1966; Andries et al., 2000), the quasi-equilibrium models (Gorban & Karlin, 1994; Levermore, 1996), the Shakhov model (Shakhov, 1968), and the Fokker-Plank model (Singh & Ansumali, 2015; Singh et al., 2016) are used as kinetic models with tunable Prandtl number. The ES–BGK model (Holway Jr, 1966; Andries et al., 2000) is an elegant but simple improvement over the BGK model. This model assumes that the distribution function relaxes to an anisotropic Gaussian distribution within mean free time $\tau$. The anisotropic Gaussian in itself evolves towards the MB distribution with a second time scale. The presence of a second time scale as free parameter ensures that the time scales related to momentum and thermal diffusivity are independent and thus permits one to vary the Prandtl number in the range of $2/3$ to $\infty$.

Despite their success, the Boltzmann collision kernel and its aforementioned simplifications are limited to monoatomic gases as they do not account for the internal molecular structure. However, many real gases such as nitrogen, oxygen, or methane are polyatomic. At the macroscopic level, the internal molecular structure predominantly manifests in terms of modified specific heat ratio $\gamma$ and bulk viscosity $\eta_{\rm B}$, which is crucial for a number of aerodynamic and turbomachinery engineering applications (von Backstrom, 2008; Wu et al., 2015). The specific heat ratio predicted by the Boltzmann equation is that of a monoatomic gas ($\gamma=5/3$), whereas that of a diatomic gas is $7/5$.

Two-particle kinetic theory as an extension of the Boltzmann equation as expected correctly predicts the specific heat ratio for polyatomic gases along with heat conductivity and the bulk viscosity (Wang-Chang & Uhlenbeck, 1951; Wu et al., 2015; Chapman & Cowling, 1970). However, it is often not feasible to do any analysis on the Boltzmann-type equation for polyatomic gases. Therefore, several simplifications to model polyatomic gases have also been proposed. They essentially incorporate the rotational kinetic energy by decomposing the two-particle distribution function into two independent single-particle distribution functions (Morse, 1964; Andries et al., 2000; Tsutahara et al., 2008; Kataoka & Tsutahara, 2004; Watari, 2007; Nie et al., 2008; Larina & Rykov, 2010; Wu et al., 2015; Wang et al., 2017; Bernard et al., 2019). Furthermore, a thermodynamic framework and extensions thereof were developed for modelling highly nonequilibrium phenomena in dense and rarefied polyatomic gases where the Navier–Stokes–Fourier theory is no longer valid (Müller & Ruggeri, 2013; Ruggeri & Sugiyama, 2015; Arima et al., 2012). A few BGK like models have also been proposed for polyatomic gases which accept the Prandtl number as a tunable parameter (Andries et al., 2000; Brull & Schneider, 2009).

Hydrodynamic simulations for a realistic system require the development of reduced-order models to account for rotational degrees of freedom ideally without increasing the phase-space dimensionality. Indeed, the standard approach is to demonstrate that the two-particle distribution function describing the translational and rotational degree of freedom can be approximated by considering two single-particle distribution functions (one each for the translational and rotational degree of freedom) whose dynamics are coupled to each other (Andries et al., 2000). However, recently it was pointed out that a simplified description in terms of single-particle distribution function for the translational degree of freedom and a scalar field variable for rotational kinetic energy is sufficient for modelling the change in specific heat ratio for a dilute diatomic gas in the hydrodynamic limit (Kolluru et al., 2020a). This model supplemented the standard Boltzmann BGK equation with an advection-relaxation equation for the evolution of rotational energy. It preserved the correct conservation laws for diatomic gases in the hydrodynamic limit and satisfied the $H$ theorem. However, the model was restricted to diatomic gases and a Prandtl number $7/5$, limiting its application for heat transfer problems.

We propose a kinetic model of polyatomic gases to tune Prandtl number, specific heat ratio, and bulk viscosity in a physically transparent fashion. To do so, we write a new collision kernel which is a linear combination of the ES–BGK and BGK kernels that are locally relaxing to different temperatures at different timescales. The ratio of the two relaxation timescales is used to tune the Prandtl number. We couple the evolution of the single-particle distribution function (with this modified collision kernel) via an advection-diffusion-relaxation equation for the rotational energy. The rotational contribution to the internal energy alters the specific heat ratio to that of a polyatomic gas and allows to model bulk viscosity contribution arising out of the rotational degree of freedom. Such an extension of the ES–BGK model indeed reproduces the hydrodynamic behaviour of a polyatomic gas and also has a valid $H$ theorem. These minimal extensions of the ES–BGK model of monoatomic gas are constructed at a single-particle level for polyatomic gases and are phenomenological by construction. It is commensurate with the top-down modelling approach as developed in the context of the lattice Boltzmann models and aim to be analytically and numerically tractable (Succi, 2001; Ansumali et al., 2007b; Atif et al., 2017). The present model which requires only the solution of an advection-diffusion-relaxation equation along with the Boltzmann ES–BGK equation adds only a minor complexity over analogous monoatomic gas ES–BGK model and can be implemented in the mesoscale framework such as lattice Boltzmann (LB) method quite easily. This approach is distinctly different and is more detailed than the existing approach in the LB models where the effect of rotational degree of freedom is further coarse-grained and the correction needed to model specific heat ratio is directly added as a force term in the BGK collision model (Kataoka & Tsutahara, 2004; Nie et al., 2008; Chen et al., 2010; Huang et al., 2020). In contrast, this model of polyatomic gas enlarges the set of microscopic degrees of freedom and models dynamics of rotational energy in an explicit manner.

The manuscript is organized as follows: A brief kinetic description of monoatomic and polyatomic gases is given in Sections 2 and 3 respectively. In Section 4 we propose an extension to the ES–BGK model for polyatomic gases. The lattice Boltzmann formulation is described in Section 5. The proposed model is numerically validated in Section 6.

The dynamics of dilute monoatomic gases is well-described by the Boltzmann equation in terms of the evolution of the single-particle distribution function $f$, where $f({\bm{x}},{\bm{c}},t)\,d{\bm{x}}\,d{\bm{c}}$ is the probability of finding a particle within $({\bm{x}},{\bm{x}}+d{\bm{x}})$, possessing a velocity in the range $({\bm{c}},{\bm{c}}+d{\bf c})$ at a time $t$. The hydrodynamic variables are density $\rho({\bm{x}},t)$, velocity ${\bm{u}}({\bm{x}},t)$ and total energy $E({\bm{x}},t)=(\rho u^{2}+3\rho k_{B}T/m)/2$. Here onwards, we use a scaled temperature $\theta$ defined in terms of Boltzmann constant $k_{B}$ and mass of the particle $m$ as $\theta=k_{B}T/m$. The thermodynamic pressure $p$ and the scaled temperature $\theta$ are related via the ideal gas equation of state as $p=\rho\theta$. These hydrodynamic variables are computed as the moments of the single-particle distribution function

$\displaystyle\begin{split}\{\rho,\rho{\bm{u}},E\}=\left\langle\left\{1,{\bm{c}},\frac{c^{2}}{2}\right\},f\right\rangle,\end{split}$

(1)

where we define the averaging operator $\langle\phi_{1}({\bm{c}}),\phi_{2}({\bm{c}})\rangle=\int_{-\infty}^{\infty}\phi_{1}({\bm{c}})\phi_{2}({\bm{c}})d{\bm{c}}$. In the comoving reference frame with fluctuating velocity ${\bm{\xi}}={\bm{c}}$ - ${\bm{u}}$, the stress tensor $\sigma_{\alpha\beta}$ is traceless part of the flux of the momentum tensor $\Theta_{ij}\equiv=\sigma_{\alpha\beta}+\rho\theta\delta_{\alpha\beta}$ and the heat flux $q_{\alpha}$ is the flux of the energy. Thus,

$$\Theta_{ij}\equiv\left<\xi_{i}\xi_{j},f\right>,\;\;\sigma_{\alpha\beta}=\left<\overline{\xi_{\alpha}\xi_{\beta}},f\right>,\quad q_{\alpha}=\left<\xi_{\alpha}\frac{\xi^{2}}{2},f\right>,\quad$$

(2)

where the symmetrized traceless part $\overline{A_{\alpha\beta}}$ for any second-order tensor $A_{\alpha\beta}$ is $\overline{A_{\alpha\beta}}=\left(A_{\alpha\beta}+A_{\beta\alpha}-2A_{\gamma\gamma}\delta_{\alpha\beta}/3\right)/2$. The stress tensor and heat flux tensor are related to the pressure tensor $P_{\alpha\beta}=\left<c_{\alpha}c_{\beta},f\right>$ and energy flux $\left\langle c^{2}c_{\alpha}/2,f\right\rangle$ respectively as

$\displaystyle\begin{split}\sigma_{\alpha\beta}&=P_{\alpha\beta}-\rho u_{\alpha}u_{\beta}-\rho\theta\delta_{\alpha\beta},\qquad q_{\alpha}=\left\langle\frac{c^{2}c_{\alpha}}{2},f\right\rangle-u_{\alpha}\left(E+\rho\theta\right)-u_{\beta}\sigma_{\alpha\beta}.\end{split}$

(3)

We also define third-order moment $Q_{\alpha\beta\gamma}$, with the traceless part $\overline{Q_{\alpha\beta\gamma}}$ and the fourth-order moments divided into contracted part $R_{\alpha\beta}$ and the trace $R$

$$Q_{\alpha\beta\gamma}=\left<\xi_{\alpha}\xi_{\beta}\xi_{\gamma},f\right>,\quad\overline{Q_{\alpha\beta\gamma}}=Q_{\alpha\beta\gamma}-\frac{2}{5}\left(q_{\alpha}\delta_{\beta\gamma}+q_{\beta}\delta_{\alpha\gamma}+q_{\gamma}\delta_{\alpha\beta}\right),$$

(4)

$$R_{\alpha\beta}=\left<\xi^{2}\xi_{\alpha}\xi_{\beta},f\right>,\quad R=\left<\xi^{4},f\right>.$$

(5)

In the dilute gas limit, the time evolution of the distribution function is a sequence of free-flight and binary collisions well described by the Boltzmann equation

$\displaystyle\begin{split}\partial_{t}f+c_{\alpha}\partial_{\alpha}f=\Omega(f,f),\end{split}$

(6)

where the collision kernel $\Omega(f,f)$ models the binary collisions between particles under the assumptions of molecular chaos (Chapman & Cowling, 1970; Cercignani, 1988; McQuarrie, 2000). The nonlinear integro-differential Boltzmann collision kernel is often replaced by simpler models that should recover the following essential features (Cercignani, 1988):

1. (a)

   Conservation Laws: As the mass, momentum, and energy density of the particles are conserved during the elastic collision, any collision model must satisfy

   $$\langle\Omega(f,f),\{1,{\bm{c}},c^{2}\}\rangle=\{0,\mathbf{0},0\}.$$

   (7)

   Thus, the macroscopic conservation laws obtained by taking appropriate moments (with respect to $\{1,{\bm{c}},c^{2}\}$) of the Boltzmann equation

   $\displaystyle\begin{split}&\partial_{t}\rho+\partial_{\alpha}\left(\rho u_{\alpha}\right)=0,\\ &\partial_{t}\left(\rho u_{\alpha}\right)+\partial_{\beta}\left(\rho u_{\alpha}u_{\beta}+\rho\theta\delta_{\alpha\beta}+\sigma_{\alpha\beta}\right)=0,\\ &\partial_{t}E+\partial_{\alpha}\left[\left(E+p\right)u_{\alpha}+q_{\alpha}+\sigma_{\alpha\beta}u_{\beta}\right]=0,\end{split}$

   (8)

   are the same as those of the compressible hydrodynamics.

2. (b)

   Zero of collision: The collision term is zero if and only if the distribution function attains Maxwell–Boltzmann form, i.e.,

   $$\Omega(f,f)=0\iff f=f^{\rm{MB}},$$

   (9)

   where the Maxwell–Boltzmann distribution $f^{\rm MB}$ is

   $$f^{\rm MB}(\rho(f),{\bm{u}}(f),\theta(f))=\rho\,\left(\frac{1}{2\pi\theta}\right)^{3/2}\exp{\left(-\frac{(c-u)^{2}}{2\,\theta}\right)},$$

   (10)

   which ensures that the dynamics is consistent with the equilibrium thermodynamics.

3. (c)

   The H Theorem: The Boltzmann equation generalizes the second law of thermodynamics to nonequilibrium situations. Boltzmann defined a nonequilibrium generalization of the entropy known as the $H$ function (Cercignani, 1988)

   $$H[f]=\int(f\ln f-f)dc.$$

   (11)

   At equilibrium, the thermodynamic entropy is calculated as (Cercignani, 1988)

   $$S^{\rm eq}=-k_{\rm B}H[f^{\rm MB}]=-\rho k_{B}\left[\ln{\frac{\rho}{\left(2\pi\theta\right)^{3/2}}-\frac{5}{2}}\right].$$

   (12)

   The evolution of this $H$ function is

   $$\partial_{t}H+\partial_{\alpha}J_{\alpha}^{H}=\underbrace{\left\langle\Omega(f,f),\ln f\right\rangle}_{\Sigma}\leq 0,$$

   (13)

   where $J_{\alpha}^{H}$ is the entropy flux and $\Sigma$ is the negative of the entropy generation. Boltzmann demonstrated that the entropy generation is positive, hence, the $H$ function is nonincreasing in time (Cercignani, 1988). Thus, any consistent approximation for the Boltzmann collision kernel should also satisfy the same condition.

We briefly describe the two most widely used models, the BGK collision model and the ES–BGK model. The Bhatnagar–Gross–Krook (BGK) model, perhaps the simplest and most widely-used model of the Boltzmann collision kernel models the collision as a relaxation of the distribution function towards the equilibrium $f^{\rm{MB}}$ (Bhatnagar et al., 1954) as

$$\Omega_{\rm{BGK}}=\frac{1}{\tau}\left(f^{\rm{MB}}(\rho,{\bm{u}},\theta)-f\right).$$

(14)

This assumes that the process occurs at a single time scale $\tau$ corresponding to the mean free time. It is trivial to see that this model has the same collisional invariants as the Boltzmann kernel, hence, recovers the same conservation laws (Bhatnagar et al., 1954). The entropy production $\Sigma^{\rm BGK}$ due to the BGK model is written as

$\displaystyle\begin{split}\Sigma^{\rm{BGK}}=\left\langle\Omega_{\rm{BGK}},\ln f\right\rangle=\frac{1}{\tau}\int\left(f^{\rm{MB}}(\rho,{\bm{u}},\theta)-f\right)\ln\left(\frac{f}{f^{\rm MB}(\rho,{\bm{u}},\theta)}\right)d{\bm{c}}\leq 0.\end{split}$

(15)

Thus, like the Boltzmann equation, the BGK model also satisfies the $H$ theorem. By taking appropriate moments of the Boltzmann BGK equation, we can also see that the evolution of the stress and the heat flux are (Liboff, 2003)

$\displaystyle\begin{split}\partial_{t}\sigma_{\alpha\beta}+\partial_{\gamma}\left(\sigma_{\alpha\beta}u_{\gamma}\right)+\partial_{\gamma}\overline{Q_{\alpha\beta\gamma}}+2\overline{\sigma_{\gamma\beta}\partial_{\gamma}u_{\alpha}}+2p\overline{\partial_{\beta}u_{\alpha}}+\frac{4}{5}\overline{\partial_{\beta}q_{\alpha}}=-\frac{1}{\tau}\sigma_{\alpha\beta},\\ \partial_{t}{q}_{\alpha}+\frac{1}{2}\partial_{\beta}\left(\overline{R_{\alpha\beta}}+\frac{1}{3}R\,\delta_{\alpha\beta}\right)+\overline{Q_{\alpha\beta\gamma}}\partial_{\gamma}u_{\beta}+\partial_{\beta}\left(q_{\alpha}u_{\beta}\right)+\frac{7}{5}q_{\beta}\partial_{\beta}u_{\alpha}\\ +\frac{2}{5}q_{\alpha}\partial_{\beta}u_{\beta}+\frac{2}{5}q_{\beta}\partial_{\alpha}u_{\beta}-\frac{5}{2}\frac{p}{\rho}\partial_{\alpha}p-\frac{5}{2}\frac{p}{\rho}\partial_{\beta}\sigma_{\alpha\beta}-\frac{\sigma_{\alpha\beta}}{\rho}\partial_{\beta}p-\frac{\sigma_{\alpha\beta}}{\rho}\partial_{\eta}\sigma_{\beta\eta}=-\frac{{q}_{\alpha}}{\tau}.\end{split}$

(16)

The form of relaxation dynamics shows that the time-scales for the momentum diffusivity and thermal diffusivity are equal for the BGK model. These equations show that, like any equation of Boltzmann type, the dynamics of $M^{th}$ order moment involves $(M+1)^{th}$ moment and thus form an infinite order moment chain. The hydrodynamic limit is typically analysed via Chapman–Enskog expansion which allows evaluating the dynamic viscosity $\mu$ and thermal conductivity $\kappa$ for the model with the specific heat $C_{p}$ for a monoatomic ideal gas as $5/2$ is (Chapman & Cowling, 1970)

$$\mu=p\tau,\quad\kappa=\frac{5}{2}p\tau\quad\implies{\rm Pr}=\frac{\mu C_{p}}{k}=1.$$

(17)

Despite this defect of ${\rm Pr}=1$, the BGK model is extremely successful both as a numerical and an analytical tool for analysis. The ES–BGK model (Holway Jr, 1965) also describes the collision as simple relaxation process but unlike the BGK model, it overcomes the restriction on the Prandtl number. The extra ingredient for ES–BGK model is quasi-equilibrium form of distribution derived by minimizing the H–function (Eq.(11)) under the constraints of an additional condition of fixed stresses, which implies absolute minimization of

$$\Xi[f]=\int d{\bf c}\left[(f\ln f-f)+\alpha f+\beta_{j}c_{j}f+\gamma_{ij}c_{i}c_{j}f\right].$$

(18)

The solution to this minimization problem is an anisotropic Gaussian

$$f^{\rm Quasi}\left(\rho,\bm{u},\Theta_{ij}\right)=\frac{\rho}{\sqrt{\rm{det}[2\pi\Theta_{ij}]}}\exp\left({-\frac{1}{2}\xi_{i}\Theta_{ij}^{-1}\xi_{j}}\right).$$

(19)

Like the MB distribution, this has the same conserved moments as that of the single-particle distribution function but this also treats $\Theta_{ij}(f)$ as a conserved variable (Kogan, 1969) i.e.,

$$\{\rho(f^{\rm Quasi}),u_{\alpha}(f^{\rm Quasi}),\theta(f^{\rm Quasi}),\Theta_{ij}(f^{\rm Quasi})\}=\{\rho(f),u_{\alpha}(f),\theta(f),\Theta_{ij}(f)\}.$$

(20)

On this quasi-equilibrium manifold with stress as variable when the $H$ is minimum we have

$$H[f^{\rm Quasi}\left(\rho,\bm{u},\Theta_{ij}\right)]=\rho\,\ln{\left(\frac{\rho}{\sqrt{\rm{det}[2\pi\Theta_{ij}]}}\right)}-\frac{5}{2}\rho.$$

(21)

The ES–BGK model uses the anisotropic Gaussian distribution $f^{\rm ES}\equiv f^{\rm Quasi}\left(\rho,\bm{u},\lambda_{ij}\right)$

$$f^{\rm ES}\left(\rho,\bm{u},\lambda_{ij}\right)=\frac{\rho}{\sqrt{\rm{det}[2\pi\lambda_{ij}]}}\exp\left({-\frac{1}{2}\xi_{i}\lambda^{-1}_{ij}\xi_{j}}\right),$$

(22)

where instead of pressure tensor a positive definite matrix $\lambda_{ij}$

$$\lambda_{ij}=(1-b)\theta\delta_{ij}+b\Theta_{ij}\equiv\theta\delta_{ij}+b\sigma_{ij},$$

(23)

is used with $-1/2\leq b\leq 1$ as a free parameter, the range of $b$ is dictated by the positive definiteness of $\lambda^{-1}_{ij}$. For $0\leq b\leq 1$, $\lambda_{ij}$ is trivially positive as it is a convex combination of two positive definite matrices. The non-trivial range $-1/2\leq b<0$ is better understood from the eigenvalue analysis of $\lambda_{ij}$ (Andries et al., 2000). Let the eigenvalues of $\lambda_{ij}$ be $\Lambda_{i}$ and that of positive definite matrix $\Theta_{ij}$ be $\phi_{i}$ and thus $\phi_{1}+\phi_{2}+\phi_{3}=\Lambda_{1}+\Lambda_{2}+\Lambda_{3}=3\theta$. In terms of these eigenvalues, the matrix ${\bm{\lambda}}$ after suitable rotation can be rewritten as

$${\bm{\lambda}}=\begin{pmatrix}(1-b)\theta+b\phi_{1}&0&0\\ 0&(1-b)\theta+b\phi_{2}&0\\ 0&0&(1-b)\theta+b\phi_{3}\end{pmatrix},$$

(24)

which is a convex combination of diagonal matrices $\Psi_{1},\Psi_{2}$ and $\Psi_{3}$ as

$${\bm{\lambda}}=\left(\frac{1+2b}{3}\right)\underbrace{\begin{pmatrix}\phi_{1}&0&0\\ 0&\phi_{2}&0\\ 0&0&\phi_{3}\end{pmatrix}}_{\Psi_{1}}+\left(\frac{1-b}{3}\right)\underbrace{\begin{pmatrix}\phi_{2}&0&0\\ 0&\phi_{3}&0\\ 0&0&\phi_{1}\end{pmatrix}}_{\Psi_{2}}+\left(\frac{1-b}{3}\right)\underbrace{\begin{pmatrix}\phi_{3}&0&0\\ 0&\phi_{1}&0\\ 0&0&\phi_{2}\end{pmatrix}}_{\Psi_{3}}.$$

(25)

As in the diagonal representation, the non-zero components are the eigenvalues and one obtains the relationship between $\phi_{i}$ and $\Lambda_{i}$ as

$\displaystyle\begin{split}\Lambda_{i}=\left(\frac{1+2b}{3}\right)\phi_{i}+\left(\frac{1-b}{3}\right)\sum_{i\neq j}\phi_{j}.\end{split}$

(26)

Thus, the eigenvalues of $\lambda_{ij}$ are non-negative if the range of $b$ is restricted to $-1/2\leq b\leq 1$ as we also know that $\phi_{i}\geq 0$. In the ES–BGK model, the collisional relaxation is towards this anisotropic Gaussian distribution $f^{\rm ES}$ which itself attains the form of the Maxwellian at the equilibrium. The collision kernel in explicit form is

$$\Omega_{\rm ESBGK}=\frac{1}{\tau}\left(f^{\rm ES}\left(\rho,\bm{u},\lambda_{ij}\right)-f\right),$$

(27)

where it is evident that $b=0$ corresponds to the limit of the BGK equation and $b=1$ would imply that stress is conserved. For $b\neq 1$, the model has the same set of conservation laws as the BGK equation. As the stress and heat flux tensors follow the relation

$\displaystyle\begin{split}\sigma_{\alpha\beta}(f^{\rm ES})=b\sigma_{\alpha\beta}(f),\quad\left\langle\xi_{\alpha}\xi^{2},f^{\rm ES}\right\rangle={0},\end{split}$

(28)

thus, the evolution equations for stress tensor and heat flux are

$\displaystyle\begin{split}\partial_{t}\sigma_{\alpha\beta}+\partial_{\gamma}\left(u_{\gamma}\sigma_{\alpha\beta}\right)+\partial_{\gamma}\overline{Q_{\alpha\beta\gamma}}+2\overline{\sigma_{\gamma\beta}\partial_{\gamma}u_{\alpha}}+2p\overline{\partial_{\beta}u_{\alpha}}+\frac{4}{5}\overline{\partial_{\beta}q_{\alpha}}=-\left(\frac{1-b}{\tau}\right)\sigma_{\alpha\beta},\\ \partial_{t}{q}_{\alpha}+\frac{1}{2}\partial_{\beta}\left(\overline{R_{\alpha\beta}}+\frac{1}{3}R\,\delta_{\alpha\beta}\right)+\overline{Q_{\alpha\beta\gamma}}\partial_{\gamma}u_{\beta}+\partial_{\beta}\left(q_{\alpha}u_{\beta}\right)+\frac{7}{5}q_{\beta}\partial_{\beta}u_{\alpha}\\ +\frac{2}{5}q_{\alpha}\partial_{\beta}u_{\beta}+\frac{2}{5}q_{\beta}\partial_{\alpha}u_{\beta}-\frac{5}{2}\frac{p}{\rho}\partial_{\alpha}p-\frac{5}{2}\frac{p}{\rho}\partial_{\beta}\sigma_{\alpha\beta}-\frac{\sigma_{\alpha\beta}}{\rho}\partial_{\beta}p-\frac{\sigma_{\alpha\beta}}{\rho}\partial_{\eta}\sigma_{\beta\eta}=-\frac{{q}_{\alpha}}{\tau},\end{split}$

(29)

which shows that the momentum and thermal diffusivities are different in an ES–BGK model and the Prandtl number ${\rm Pr}=\mu c_{p}/\kappa$ is a free parameter. In particular, the Chapman-Enskog analysis of this model yields

$$\mu=\frac{p\tau}{1-b},\quad\kappa=\frac{5}{2}p\tau\quad\implies{\rm Pr}=\frac{1}{1-b}.$$

(30)

Thus, the free parameter $b$ in the anisotropic Gaussian allows one to vary the Prandtl number from $2/3$ to infinity. At $b=-1/2$, the Prandtl number predicted by the ES–BGK model is $2/3$, which matches with the value obtained from the Boltzmann equation, and when $b=0$, the model is equivalent to the BGK model. The thermal conductivity is fixed only by $\tau$ while the viscosity can be tuned via $b$ to obtain the required Prandtl number.

The $H$ theorem for this model was first proved by Andries & Perthame (2001) by showing that the entropy production $\Sigma^{\rm ESBGK}$ is non-positive. We briefly review the proof of the $H$ theorem for the ES-BGK model. The expression for the entropy production is

$\displaystyle\Sigma^{\rm ESBGK}=\left\langle\Omega_{\rm ESBGK},\ln f\right\rangle=\frac{1}{\tau}\left\langle\left(f^{\rm ES}\left(\rho,\bm{u},\lambda_{ij}\right)-f\right),\frac{\partial H}{\partial f}\right\rangle.$

(31)

The proof is built on the property of an arbitrary convex function $G(x)$ that

$$\frac{\partial G}{\partial x}\left(y-x\right)\leq G(y)-G(x)$$

(32)

using which we can write

$\displaystyle\Sigma^{\rm ESBGK}\leq\frac{1}{\tau}\left(H[f^{\rm ES}\left(\rho,\bm{u},\lambda_{ij}\right)]-H\left[f^{\rm Quasi}\left(\rho,\bm{u},\Theta_{ij}\right)\right]\right)+\frac{1}{\tau}\Delta H^{\rm Quasi},$

(33)

with the last term as $\Delta H^{\rm Quasi}=H[f^{\rm ES}\left(\rho,\bm{u},\Theta_{ij}\right)]-H[f]$ in the above equation is non-positive as $f^{\rm Quasi}\left(\rho,\bm{u},\Theta_{ij}\right)$ is by construction the minima of $H$. To prove that $H[f^{\rm ES}\left(\rho,\bm{u},\lambda_{ij}\right)]-H\left[f^{\rm Quasi}\left(\rho,\bm{u},\Theta_{ij}\right)\right]\leq 0$, using Eq.(21) we have

$\displaystyle H[f^{\rm ES}\left(\rho,\bm{u},\lambda_{ij}\right)]-H\left[f^{\rm Quasi}\left(\rho,\bm{u},\Theta_{ij}\right)\right]=\frac{1}{2}\rho\ln\left(\frac{{\rm det}[\Theta_{ij}]}{{\rm det}[\lambda_{ij}]}\right).$

(34)

Starting from the Brunn-Minkowsky inequality

$${\rm det}[aA+(1-a)B]\geq\left({\rm det}[A]\right)^{a}\left({\rm det}[B]\right)^{1-a},$$

(35)

relating the determinants of two positive matrices $A$ and $B$ and their convex combinations, we can show that ${\rm det}[\lambda_{ij}]\geq{\rm det}[\Theta_{ij}]$ (Horn & Johnson, 2012). This inequality along with Eq.(25) allows us to write

$${\rm det}[\lambda_{ij}]\geq\left({\rm det}[\Psi_{1}]\right)^{\frac{1+2b}{3}}\left({\rm det}[\Psi_{2}]\right)^{\frac{1-b}{3}}\left({\rm det}[\Psi_{3}]\right)^{\frac{1-b}{3}}.$$

(36)

However, from the definitions of $\Psi_{1},\Psi_{2}$, and $\Psi_{3}$ one can see that

$${\rm det}[\Psi_{1}]={\rm det}[\Psi_{2}]={\rm det}[\Psi_{3}]=\phi_{1}\phi_{2}\phi_{3}={\rm det}[\Theta_{ij}].$$

(37)

Hence, the total entropy production is non-positive, i.e., $\Sigma^{\rm ESBGK}\leq 0$.

The rotational degrees of freedom of a polyatomic gas manifest themselves at the continuum level in terms of change in specific heat ratio $\gamma$ and a non-zero bulk viscosity due to interaction among the translational component $E_{\rm T}=\rho u^{2}/2+3\rho\theta_{T}/2$ and rotational component $E_{\rm R}$ of energy. Thus, the rotational degrees of freedom need to be explicitly accounted for in any microscopic or kinetic description. Indeed, typically the kinetic descriptions are in terms of a two-particle distribution function $F({\bm{x}},{\bm{c}},t,I)$ which defines the probability of finding a molecule with a position in the range $({\bm{x}},{\bm{x}}+d{\bm{x}})$ possessing a velocity in the range $(\bm{c},\bm{c}+\bm{dc})$ with an internal energy $(I,I+dI)$ due to the additional degrees of freedom (Morse, 1964; Rykov, 1975; Kuščer, 1989). For a polyatomic gas with $\delta$ additional rotational degrees of freedom, the moments of this distribution function give the density, momentum, and total energy (with $\delta=0$ corresponding to a monoatomic gas)

$\displaystyle\begin{split}\{\rho,\rho{\bm{u}},E_{T}+E_{R}\}=\left\langle\left\langle\left\{1,{\bm{c}},\frac{c^{2}}{2}+I^{2/\delta}\right\},F\right\rangle\right\rangle,\end{split}$

(38)

like its monoatomic counterpart and the operator $\left<\left<,\right>\right>$ is defined as

$$\left\langle\left\langle\phi_{1}({\bm{c}},I),\phi_{2}({\bm{c}},I)\right\rangle\right\rangle=\int\int\phi_{1}({\bm{c}},I)\phi_{2}({\bm{c}},I)d{\bm{c}}dI.$$

(39)

For the reduced-order modeling, the distribution function $F({\bm{x}},{\bm{c}},t,I)$ is often split into two coupled distribution functions $f_{1}({\bm{x}},{\bm{c}},t)$ and $f_{2}({\bm{x}},{\bm{c}},t)$ defined as

$\displaystyle\begin{split}f_{1}({\bm{x}},{\bm{c}},t)=\int F({\bm{x}},{\bm{c}},I,t)dI,\quad\quad f_{2}({\bm{x}},{\bm{c}},t)=\int F({\bm{x}},{\bm{c}},I,t)I^{2/\delta}dI,\end{split}$

(40)

where $f_{1}$ is related to the translational energy and $f_{2}$ with the rotational energy dynamics (Rykov, 1975; Andries et al., 2000). The moments of reduced distribution $f_{1}({\bm{x}},{\bm{c}},t)$ are then same as the moments of single-particle distribution function

$\displaystyle\{\rho,\rho{\bm{u}},E_{T}\}=\left\langle\left\{1,{\bm{c}},\frac{c^{2}}{2}\right\},f_{1}\right\rangle.$

(41)

By construction, we have the zeroth moment of $f_{2}({\bm{x}},{\bm{c}},t)$ as the rotational energy

$\displaystyle E_{R}=\left\langle\left\langle\frac{c^{2}}{2}+I^{2/\delta},F\right\rangle\right\rangle-\left\langle f_{1},\frac{c^{2}}{2}\right\rangle=\left\langle\left\langle F,I^{2/\delta}\right\rangle\right\rangle\equiv\left\langle f_{2},1\right\rangle=\frac{\delta}{2}\rho\theta_{R}.$

(42)

In other words, the temperature $\theta$ consists of contributions from the translational and rotational temperatures, and they follow the relation

$$\theta=\left(\frac{3}{3+\delta}\right)\theta_{T}+\left(\frac{\delta}{3+\delta}\right)\theta_{R},$$

(43)

and in thermodynamic equilibrium the equipartition of energy requires $\theta_{R}=\theta_{T}$. The heat flux for a polyatomic gas is $q_{\alpha}=q^{T}_{\alpha}+q^{R}_{\alpha}$ where $q^{T}_{\alpha}$ is the translational heat flux and $q^{R}_{\alpha}$ is an additional heat flux due rotational energy. The rotational heat flux and a stress like quantity (second moment) are defined as

$$q^{\rm R}_{\alpha}=\int d\xi f_{2}\xi_{\alpha},\quad\sigma^{R}_{\alpha\beta}=\int d\xi f_{2}\xi_{\alpha}\xi_{\beta}-\rho\theta^{2}\delta_{\alpha\beta}.$$

(44)

Like the monoatomic gas, the evolution equation for this distribution function $F({\bm{x}},{\bm{c}},t,I)$ with collisional kernel $\Omega(F,F)$ in the Boltzmann form is

$\displaystyle\begin{split}\partial_{t}F+c_{\alpha}\partial_{\alpha}F=\Omega(F,F),\end{split}$

(45)

which is consistent with the equipartition of energy at equilibrium (Pullin, 1978; Kuščer, 1989). Similar to the monoatomic gas, one defines the BGK collision kernel in terms of the two-particle distribution function for a polyatomic gas as (Brull & Schneider, 2009)

$\displaystyle\begin{split}\Omega_{\rm{BGK}}&=\frac{1}{\tau}\left(F^{\rm MB}(\rho,\bm{u},\theta,I)-F\right),\\ F^{\rm MB}(\rho,\bm{u},\theta,I)&=\frac{\rho\Lambda_{\delta}}{\left(2\pi\theta\right)^{3/2}\theta^{\delta/2}}\exp{\left(-\left(\frac{({\bm{c}}-{\bm{u}})^{2}}{2\theta}+\frac{I^{2/\delta}}{\theta}\right)\right)},\end{split}$

(46)

with normalisation factor $\Lambda_{\delta}=\int\exp{(-I^{2/\delta})dI}$. Equation (45) with $\Omega_{\rm{BGK}}$ is written as two kinetic equations for the reduced distributions $f_{1}({\bm{x}},{\bm{c}},t)$ and $f_{2}({\bm{x}},{\bm{c}},t)$ by multiplying with $1$ and $I^{2/\delta}$ and then integrating over the internal energy variable as

$\displaystyle\begin{split}\partial_{t}f_{1}+c_{\alpha}\partial_{\alpha}f_{1}&=\frac{1}{\tau}\left(f_{1}^{\rm MB}(\rho,\bm{u},\theta)-f_{1}\right),\\ \partial_{t}f_{2}+c_{\alpha}\partial_{\alpha}f_{2}&=\frac{1}{\tau}\left(\frac{\delta}{2}\theta f_{1}^{\rm MB}(\rho,\bm{u},\theta)-f_{2}\right).\end{split}$

(47)

This approach, where two reduced distributions are weakly coupled via temperature, recovers all the features of Eq.(46) and is widely adopted for polyatomic gases. Andries et al. (2000) extended this approach via an extended ES–BGK collision kernel as

$$\Omega_{\rm ESBGK}(F)=\frac{Z_{\rm ES}}{\tau}\left(F^{\rm ES}(\rho,u,\lambda_{ij},\theta_{\rm rel})-F\right),$$

(48)

where $\lambda_{ij}=(1-\alpha)\left[(1-b)\theta_{T}\delta_{ij}+b\Theta_{ij}\right]+\alpha\theta\delta_{ij}$ with a generalized Gaussian $F^{\rm ES}$

$$F^{\rm ES}(\rho,u,\lambda_{ij},\theta_{\rm rel})=\frac{\rho\Lambda_{\delta}}{\theta_{\rm rel}^{\delta/2}\sqrt{\rm{det}[2\pi\lambda_{ij}]}}\exp\left({-\frac{1}{2}\xi_{i}\lambda^{-1}_{ij}\xi_{j}}-\frac{I^{2/\delta}}{\theta_{\rm rel}}\right),$$

(49)

with $\theta_{\rm rel}=\alpha\theta+(1-\alpha)\theta_{R}$ and $Z_{\rm ES}=1/(1-b+b\alpha)$. Similar to the monoatomic ES–BGK model, the parameter $b$ is used to tune the Prandtl number, while parameter $\alpha$ is used to tune the bulk viscosity coefficient independently. The reduced description which generalizes the ES–BGK model in terms of the $f_{1}$ and $f_{2}$ is

$\displaystyle\begin{split}\partial_{t}f_{1}+c_{\alpha}\partial_{\alpha}f_{1}&=\frac{Z_{\rm ES}}{\tau}\left(f^{\rm ES}(\rho,u,\lambda_{ij})-f_{1}\right),\\ \partial_{t}f_{2}+c_{\alpha}\partial_{\alpha}f_{2}&=\frac{Z_{\rm ES}}{\tau}\left(\frac{\delta}{2}\theta_{\rm rel}f^{\rm ES}(\rho,u,\lambda_{ij})-f_{2}\right).\end{split}$

(50)

A common feature between Eqs. (47) and (50) is that the kinetic equation for translation distribution function $f_{1}$ is coupled with the kinetic equation for the rotational distribution function $f_{2}$ via rotational temperature $\theta_{R}$ only. At this point, it might be instructive to analyse the moment chains of the BGK and ES–BGK systems. As both the collision kernels conserve mass and momentum, the evolution equations for density and momentum are of the same form as that of monoatomic gas. In particular

$\displaystyle\begin{split}\partial_{t}\rho+\partial_{\alpha}(\rho u_{\alpha})&=0,\\ \partial_{t}(\rho u_{\alpha})+\partial_{\beta}\left(\rho u_{\alpha}u_{\beta}+\rho\theta\delta_{\alpha\beta}+\hat{\sigma}_{\alpha\beta}\right)&=0,\end{split}$

(51)

where $\hat{\sigma}_{\alpha\gamma}=\sigma_{\alpha\gamma}+\rho\left(\theta_{T}-\theta\right)\delta_{\alpha\gamma}$ is the modified stress tensor and the velocity evolution is

$$\partial_{t}u_{\alpha}+u_{\beta}\partial_{\beta}u_{\alpha}+\frac{1}{\rho}\partial_{\beta}\left(\rho\theta\delta_{\alpha\beta}+\hat{\sigma}_{\alpha\beta}\right)=0.$$

(52)

For polyatomic gases, the energy equation gets an additional contribution from the rotational energy. In both the BGK and ES–BGK models, the evolution equation for translational part of the energy and the rotational part of the energy $E_{R}=\delta\rho\theta_{R}/2$ are of the form

$\displaystyle\begin{split}\partial_{t}E_{T}+\partial_{\alpha}\left[\left(E_{T}+\rho\theta\right)u_{\alpha}+q_{\alpha}^{T}+u_{\gamma}\hat{\sigma}_{\alpha\gamma}\right]&=\frac{Z_{E}}{\tau}\frac{3\rho}{2}\left(\theta-\theta_{T}\right),\\ \partial_{t}E_{R}+\partial_{\alpha}\left(E_{R}u_{\alpha}+q^{R}_{\alpha}\right)&=\frac{Z_{E}}{\tau}\frac{\delta\rho}{2}\left(\theta-\theta_{R}\right),\end{split}$

(53)

with $Z_{E}=1$ for the BGK model and $Z_{E}=\alpha Z_{\rm ES}$ for the ES–BGK model. The evolution equation for the total energy is written as the sum of Eqs. (53) as

$$\partial_{t}\left(E_{T}+E_{R}\right)+\partial_{\alpha}\left[\left(E_{T}+E_{R}+\rho\theta\right)u_{\alpha}+q_{\alpha}+u_{\gamma}\hat{\sigma}_{\alpha\gamma}\right]=0,$$

(54)

where the relationship between translational, rotational temperatures (Eq.(43)) is used to show that energy is collisional invariant. From Eqs. (47) and (50), the stress evolution and the translational heat flux evolution equations in explicit form are

$\displaystyle\begin{split}\partial_{t}\sigma_{\alpha\beta}&+\partial_{\gamma}\left(u_{\gamma}\sigma_{\alpha\beta}\right)+\partial_{\gamma}\overline{Q_{\alpha\beta\gamma}}+\frac{4}{5}\overline{\partial_{\beta}q^{T}_{\alpha}}+2\rho\theta_{T}\overline{\partial_{\beta}u_{\alpha}}+2\overline{\partial_{\gamma}u_{\alpha}\sigma_{\gamma\beta}}=-\frac{1}{\tau}\sigma_{\alpha\beta},\\ &\partial_{t}q^{T}_{\alpha}+\partial_{\beta}\left(u_{\beta}q^{T}_{\alpha}\right)+\overline{Q_{\alpha\beta\gamma}}\partial_{\beta}u_{\gamma}+\frac{1}{2}\partial_{\beta}R_{\alpha\beta}+\frac{7}{5}q^{T}_{\beta}\partial_{\beta}u_{\alpha}+\frac{2}{5}q^{T}_{\alpha}\partial_{\eta}u_{\eta}+\frac{2}{5}q^{T}_{\beta}\partial_{\alpha}u_{\beta}\\ &-\frac{5}{2}\theta_{T}\partial_{\alpha}(\rho\theta_{T})-\frac{\sigma_{\alpha\beta}}{\rho}\partial_{\beta}(\rho\theta_{T})-\frac{5}{2}\theta_{T}\partial_{\kappa}\sigma_{\kappa\alpha}-\frac{\sigma_{\alpha\beta}}{\rho}\partial_{\kappa}\sigma_{\kappa\beta}=-\frac{Z_{q}}{\tau}q^{T}_{\alpha},\end{split}$

(55)

with $Z_{q}=1$ for BGK and $Z_{q}=Z_{\rm ES}$ for the ES–BGK model. Multiplying rotational energy equation from Eq.(53) with $u_{\alpha}$ and using Eq.(52) we have

$\displaystyle\begin{split}\partial_{t}\left(E_{R}u_{\alpha}\right)+\partial_{\beta}\left(E_{R}u_{\alpha}u_{\beta}\right)+u_{\alpha}\partial_{\beta}q^{R}_{\beta}+\frac{E_{R}}{\rho}\partial_{\beta}\left(\rho\theta\delta_{\alpha\beta}+\hat{\sigma}_{\alpha\beta}\right)=\frac{Z_{E}}{\tau}\frac{\delta}{2}\rho u_{\alpha}\left(\theta-\theta_{R}\right),\end{split}$

(56)

using which the rotational heat flux evolution obtained as first moment of $f_{2}$ dynamics is

$$\partial_{t}q^{R}_{\alpha}+\partial_{\beta}\left(u_{\beta}q^{R}_{\alpha}\right)+q^{R}_{\beta}\partial_{\beta}u_{\alpha}+\partial_{\beta}\sigma^{R}_{\alpha\beta}+\partial_{\alpha}\left(\frac{\delta}{2}\rho\theta^{2}\right)-\frac{\delta}{2}\theta_{R}\partial_{\beta}\left(\rho\theta\delta_{\alpha\beta}+\hat{\sigma}_{\alpha\beta}\right)=-\frac{Z_{q}}{\tau}q^{R}_{\alpha},$$

(57)

where $\sigma^{R}_{\alpha\beta}=\int f_{2}\xi_{\alpha}\xi_{\beta}$. The Eqs. (51) and (54) form the compressible Navier–Stokes–Fourier equations for a polyatomic gas. A Chapman-Enskog analysis shows that the Eq.(51) can be written in familiar compressible Navier–Stokes equation form as

$$\partial_{t}(\rho u_{\alpha})+\partial_{\beta}(\rho u_{\alpha}u_{\beta})+\partial_{\alpha}p-\partial_{\beta}\left(2\eta\overline{\partial_{\beta}u_{\alpha}}+\eta_{b}\partial_{\kappa}u_{\kappa}\delta_{\alpha\beta}\right)=0$$

(58)

with the shear and bulk viscosities

$$\eta=p\tau,\quad\quad\frac{\eta_{b}}{\eta}=\frac{2\delta}{3(3+\delta)Z_{E}}.$$

(59)

Similarly, an analysis of the translational and rotational heat flux dynamics at $O({\rm Kn})$ leads to $q^{T}_{\alpha}=-\kappa_{T}\partial_{\alpha}\theta,\,q^{R}_{\alpha}=-\kappa_{R}\partial_{\alpha}\theta,$ with the translational and rotational thermal conductivities $\kappa_{T}=5p\tau/(2Z_{q}),\,\kappa_{R}={\delta p\tau}/(2Z_{q})$ respectively. The effective thermal conductivity $\kappa=\kappa_{T}+\kappa_{R}={(5+\delta)p\tau/(2Z_{q})}$ Thus, the Prandtl number is ${\rm Pr}=Z_{q}$ i.e., ${\rm Pr}=1$ for BGK model and ${\rm Pr}=1/(1-b+b\alpha)$ for the ES–BGK model.