1. Introduction
The dynamics of a dilute monatomic gas in terms of the single-particle distribution function is described by the Boltzmann equation (Chapman & Cowling Reference Chapman and Cowling1970; Cercignani Reference Cercignani1988). Unlike the continuum Navier–Stokes–Fourier hydrodynamics equation, the Boltzmann equation is a valid description even at highly non-equilibrium states (Mott-Smith Reference Mott-Smith1951; Liepmann, Narasimha & Chahine Reference Liepmann, Narasimha and Chahine1962; Cercignani Reference Cercignani1988), 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, Oran & Sinkovits Reference Oh, Oran and Sinkovits1997; Struchtrup Reference Struchtrup2004; Ansumali et al. Reference Ansumali, Karlin, Arcidiacono, Abbas and Prasianakis2007b). 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, Gross & Krook Reference Bhatnagar, Gross and Krook1954; Lebowitz, Frisch & Helfand Reference Lebowitz, Frisch and Helfand1960; Holway Reference Holway1966; Shakhov Reference Shakhov1968; Gorban & Karlin Reference Gorban and Karlin1994; Levermore Reference Levermore1996; Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000; Ansumali et al. Reference Ansumali, Arcidiacono, Chikatamarla, Prasianakis, Gorban and Karlin2007a; Singh & Ansumali Reference Singh and Ansumali2015; Agrawal, Singh & Ansumali Reference Agrawal, Singh and Ansumali2020). An important example is the Bhatnagar–Gross–Krook (BGK) model (Bhatnagar et al. Reference Bhatnagar, Gross and Krook1954), which states that the relaxation of the distribution function towards the Maxwell–Boltzmann (MB) form happens on 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 a monatomic gas is $2/3$. Thus, several other variations of the collision model such as the ellipsoidal statistical BGK (ES–BGK) model (Holway Reference Holway1966; Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000), the quasi-equilibrium models (Gorban & Karlin Reference Gorban and Karlin1994; Levermore Reference Levermore1996), the Shakhov model (Shakhov Reference Shakhov1968) and the Fokker–Planck model (Singh & Ansumali Reference Singh and Ansumali2015; Singh, Thantanapally & Ansumali Reference Singh, Thantanapally and Ansumali2016) are used as kinetic models with tunable Prandtl number. The ES–BGK model (Holway Reference Holway1966; Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000) 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, thus permitting 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 monatomic gases as they do not account for the internal molecular structure. However, many real gases such as nitrogen, oxygen and 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 _{b}$, which is crucial for a number of aerodynamic and turbomachinery engineering applications (von Backstrom Reference von Backstrom2008; Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015). The specific heat ratio predicted by the Boltzmann equation is that of a monatomic 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 Reference Wang-Chang and Uhlenbeck1951; Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015; Chapman & Cowling Reference Chapman and Cowling1970). 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 Reference Morse1964; Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000; Kataoka & Tsutahara Reference Kataoka and Tsutahara2004; Watari Reference Watari2007; Nie, Shan & Chen Reference Nie, Shan and Chen2008; Tsutahara et al. Reference Tsutahara, Kataoka, Shikata and Takada2008; Larina & Rykov Reference Larina and Rykov2010; Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015; Wang et al. Reference Wang, Yan, Li and Xu2017; Bernard, Iollo & Puppo Reference Bernard, Iollo and Puppo2019). 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 (Arima et al. Reference Arima, Taniguchi, Ruggeri and Sugiyama2012; Müller & Ruggeri Reference Müller and Ruggeri2013; Ruggeri & Sugiyama Reference Ruggeri and Sugiyama2015). A few BGK-like models have also been proposed for polyatomic gases which accept the Prandtl number as a tunable parameter (Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000; Brull & Schneider Reference Brull and Schneider2009).
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 degrees of freedom can be approximated by considering two single-particle distribution functions (one for the translational and another for the rotational degrees of freedom) whose dynamics are coupled to each other (Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000). However, recently it was pointed out that a simplified description in terms of a 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, Atif & Ansumali Reference Kolluru, Atif and Ansumali2020a). 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 time scales. The ratio of the two relaxation time scales 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 the modelling of the 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 monatomic gases 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 aims to be analytically and numerically tractable (Succi Reference Succi2001; Ansumali et al. Reference Ansumali, Karlin, Arcidiacono, Abbas and Prasianakis2007b; Atif et al. Reference Atif, Kolluru, Thantanapally and Ansumali2017; Atif, Kolluru & Ansumali Reference Atif, Kolluru and Ansumali2022). 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 monatomic gas ES–BGK model and can be implemented in the mesoscale framework such as the lattice Boltzmann method quite easily. This approach is distinctly different from and is more detailed than the existing approach in the lattice Boltzmann 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 Reference Kataoka and Tsutahara2004; Nie et al. Reference Nie, Shan and Chen2008; Chen et al. Reference Chen, Xu, Zhang, Li and Succi2010; Huang, Lan & Li Reference Huang, Lan and Li2020). In contrast, this model of polyatomic gases enlarges the set of microscopic degrees of freedom and models dynamics of rotational energy in an explicit manner.
The paper is organized as follows. Brief kinetic descriptions of monatomic and polyatomic gases are given in §§ 2 and 3 respectively. In § 4, we propose an extension to the ES–BGK model for polyatomic gases. The lattice Boltzmann formulation is described in § 5. The proposed model is numerically validated in § 6. Finally, in § 7, we discuss the outlook of the present work.
2. Kinetic description of a monatomic gas
The dynamics of dilute monatomic gases is well described by the Boltzmann equation in terms of the evolution of the single-particle distribution function $f$, where $f(\boldsymbol {x}, \boldsymbol {c}, t) \, {\rm d}\kern 0.06em \boldsymbol {x} \, {\rm d} \boldsymbol {c}$ is the probability of finding a particle within $(\boldsymbol {x}, \boldsymbol {x} + {\rm d}\kern 0.06em \boldsymbol {x})$, possessing a velocity in the range $(\boldsymbol {c}, \boldsymbol {c} + {\rm d} {\boldsymbol {c}})$ at a time $t$. The hydrodynamic variables are density $\rho (\boldsymbol {x}, t)$, velocity $\boldsymbol {u}(\boldsymbol {x}, t)$ and total energy $E(\boldsymbol {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:
where we define the averaging operator $\langle \phi _1(\boldsymbol {c}), \phi _2(\boldsymbol {c}) \rangle = \int _{-\infty }^{\infty } \phi _1(\boldsymbol {c}) \phi _2(\boldsymbol {c})\,{\rm d}\boldsymbol {c}$. In the co-moving reference frame with fluctuating velocity ${\boldsymbol {\xi }} = \boldsymbol {c} - \boldsymbol {u}$, the stress tensor $\sigma _{\alpha \beta }$ is the traceless part of the flux of the momentum tensor $\varTheta _{ij}\equiv \sigma _{\alpha \beta }+\rho \theta \delta _{\alpha \beta }$ and the heat flux $q_{\alpha }$ is the flux of the energy. Thus,
where the symmetrized traceless part $\overline {A_{\alpha \beta }}$ for any second-order tensor $A_{\alpha \beta }$ is $\overline {A_{\alpha \beta }} =( A_{\alpha \beta } + A_{ \beta \alpha } - 2 A_{\gamma \gamma } \delta _{\alpha \beta }/3 )/2$. The stress tensor and heat flux tensor are related to the pressure tensor $P_{\alpha \beta } = \left \langle c_{\alpha } c_{\beta }, f\right \rangle$ and energy flux $\langle c^2 c_\alpha /2,f \rangle$ respectively as
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 $\partial _t f + c_\alpha \partial _\alpha f = \varOmega (\,f,f)$, where the collision kernel $\varOmega (\,f,f)$ models the binary collisions between particles under the assumptions of molecular chaos (Chapman & Cowling Reference Chapman and Cowling1970; Cercignani Reference Cercignani1988; McQuarrie Reference McQuarrie2000). The nonlinear integro-differential Boltzmann collision kernel is often replaced by simpler models that should recover the following essential features such as the conservation laws, the $H$ theorem, and that the MB distribution is the zero of collision (Cercignani Reference Cercignani1988).
We briefly describe the two most widely used models, the BGK collision model and the ES–BGK model. The 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^{MB}$ (Bhatnagar et al. Reference Bhatnagar, Gross and Krook1954) as
This assumes that the process occurs at a single time scale $\tau$ corresponding to the mean free time. 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 monatomic ideal gas as $5/2$ (Chapman & Cowling Reference Chapman and Cowling1970):
Despite this defect of ${Pr}=1$, the BGK model is extremely successful as both a numerical and an analytical tool for analysis. The ES–BGK model (Holway Reference Holway1965) 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 the ES–BGK model is a quasi-equilibrium form of distribution derived by minimizing the $H$ function under the constraints of an additional condition of fixed stresses. The ES–BGK model uses the anisotropic Gaussian distribution $f^{ES}\equiv f^{Quasi}(\rho, \boldsymbol {u}, \lambda _{ij} )$:
where instead of pressure tensor a positive-definite matrix $\lambda _{ij} = (1-b) \theta \delta _{ij} + b \varTheta _{ij} \equiv \theta \delta _{ij} + b \sigma _{ij}$ is used with $-1/2 \leqslant b \leqslant 1$ as a free parameter, the range of $b$ being dictated by the positive definiteness of $\lambda ^{-1}_{ij}$ (Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000). The Chapman–Enskog analysis of this model yields
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.
3. Kinetic description of a polyatomic gas
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_{T} = \rho u^2/2 + 3\rho \theta _T/2$ and rotational component $E_{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(\boldsymbol {x}, \boldsymbol {c}, t, I)$ which defines the probability of finding a molecule with a position in the range $(\boldsymbol {x}, \boldsymbol {x}+ {\rm d}\kern 0.06em \boldsymbol {x})$ possessing a velocity in the range $(\boldsymbol {c}, \boldsymbol {c}+ {\rm d}\boldsymbol {c})$ in an internal energy parameter range $(I, I+ {\rm d}I)$ due to the additional degrees of freedom (Morse Reference Morse1964; Rykov Reference Rykov1975; Kuščer Reference Kuščer1989). The internal energy due to these additional degrees of freedom is defined by assuming a continuous variable in the internal energy space as $e_{int} = I^{2/\delta }.$
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 monatomic gas):
like its monatomic counterpart and the operator $\langle \langle \dots \rangle \rangle$ is defined as $\langle \langle \phi _1(\boldsymbol {c}, I), \phi _2(\boldsymbol {c}, I) \rangle \rangle = \int \int \phi _1(\boldsymbol {c}, I) \phi _2( \boldsymbol {c}, I) \,{\rm d}\boldsymbol {c} \,{\rm d}I$. For the reduced-order modelling, the distribution function $F(\boldsymbol {x}, \boldsymbol {c}, t, I)$ is often split into two coupled distribution functions $f_1(\boldsymbol {x}, \boldsymbol {c}, t)$ and $f_2(\boldsymbol {x}, \boldsymbol {c}, t)$ defined as $f_1(\boldsymbol {x}, \boldsymbol {c}, t) = \int F(\boldsymbol {x}, \boldsymbol {c}, I, t) \,{\rm d}I$ and $f_2(\boldsymbol {x}, \boldsymbol {c}, t) = \int F(\boldsymbol {x}, \boldsymbol {c}, I, t) I^{2/\delta } \,{\rm d}I$, where $f_1$ is related to the translational energy and $f_2$ to the rotational energy dynamics (Rykov Reference Rykov1975; Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000). The moments of reduced distribution $f_1(\boldsymbol {x}, \boldsymbol {c}, t)$ are then the same as the moments of the single-particle distribution function (2.1).
By construction, we have the zeroth moment of $f_2(\boldsymbol {x}, \boldsymbol {c}, t)$ as the rotational energy:
In other words, the temperature $\theta$ consists of contributions from the translational and rotational temperatures, and they follow the relation
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 to rotational energy.
As for a monatomic gas, the evolution equation for this distribution function $F(\boldsymbol {x}, \boldsymbol {c}, t, I)$ with collisional kernel $\varOmega (F,F)$ in the Boltzmann form is
which is consistent with the equipartition of energy at equilibrium (Pullin Reference Pullin1978; Kuščer Reference Kuščer1989). Similar to a monatomic gas, one defines the BGK collision kernel in terms of the two-particle distribution function for a polyatomic gas as (Brull & Schneider Reference Brull and Schneider2009)
with normalization factor $\varLambda _{\delta } = \int \exp (-I^{2/\delta }) \,\textrm {d}I$. Equation (3.4) with $\varOmega _{BGK}$ is written as two kinetic equations for the reduced distributions $f_1(\boldsymbol {x}, \boldsymbol {c}, t)$ and $f_2(\boldsymbol {x}, \boldsymbol {c}, t)$ by multiplying by $1$ and $I^{2/\delta }$ and then integrating over the internal energy variable as
This approach, where two reduced distributions are weakly coupled via temperature, recovers all the features of (3.4) and is widely adopted for polyatomic gases. Andries et al. (Reference Andries, Le Tallec, Perlat and Perthame2000) extended this approach via an extended ES–BGK collision kernel as
where $\lambda _{ij} = (1-\alpha ) [ (1 -b) \theta _{T} \delta _{ij}+ b \varTheta _{ij} ] + \alpha \theta \delta _{ij}$ with a generalized Gaussian $F^{ES}$:
with $\theta _{rel} = \alpha \theta + (1 - \alpha ) \theta _{R}$ and $Z_{ES} = 1/(1 - b + b \alpha )$. Similar to the monatomic 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 $f_1$ and $f_2$ is
A Chapman–Enskog analysis shows that the momentum equation yields the familiar compressible Navier–Stokes equation form as
with the shear and bulk viscosities
Similarly, an analysis of the translational and rotational heat flux dynamics at $O({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 ${Pr} = Z_q$, i.e. ${Pr} = 1$ for the BGK model and ${Pr} = 1/(1-b+b\alpha )$ for the ES–BGK model.
4. Reduced ES–BGK model for polyatomic gases
In the kinetic theory of gases, one often builds an extended moment system in terms of physically relevant lower-order moments (Grad Reference Grad1958). In this spirit of Grad's moment method, one may ask whether a reduced description for rotational degrees of freedom is feasible. It should be noted that the evolution equation of $f_1$ is only weakly coupled with the evolution of $f_2$ via $\theta _R$. An appropriate choice in the current context is a reduced description in terms of lower-order moments of second distribution $f_2$. For example, the rotational component can be modelled by the evolution of two scalars – rotational energy and its flux (which are the zeroth and the first moment of $f_2$). Such a class of reduced-order kinetic models might be better suited for large-scale hydrodynamic simulations. An extended BGK model for diatomic gases was formulated by Kolluru et al. (Reference Kolluru, Atif and Ansumali2020a) wherein the BGK collision model was coupled with the rotational part of energy (zeroth moment of $f_2$) which in itself follows an advection–relaxation equation.
We extend this approach by a generalized ES–BGK model for polyatomic gases with tunable Prandtl numbers where the collision term is a linear combination of ES–BGK and the BGK collision kernels. In this model, the ES–BGK term describes relaxation to a temperature $\theta _T$ over a time $\tau$ whereas the BGK collision kernel describes relaxation to a temperature $\theta$ over a time $\tau _1$. The kinetic equation of the unified model along with the evolution equation for the rotational energy is
with the form of heat flux due to internal degrees of freedom as
This model is a minimal extension of the monatomic ES–BGK model needed for modelling polyatomic gases which also recovers all important features such as the positivity, macroscopic limit and the $H$ theorem. Here, the Prandtl number is a tunable parameter due to the presence of two relaxation time scales, whereas the rotational part of the internal energy alters the specific heat ratio to that of a polyatomic gas. The model satisfies the $H$ theorem, thus ensuring convergence to a unique equilibrium state. A few important characteristic of the present model are as follows:
(i) Conservation laws. The mass and momentum conservation equations for the proposed model are obtained by taking the zeroth and first moments of $f_1$ evolution (4.1). The second moment signifying the translational energy evolution equation is
(4.3)\begin{equation} \partial_t E_T + \partial_{\beta} \left[ \left( E_T + \rho \theta_T \right)u_\beta + \sigma_{\beta\gamma}u_\gamma + q^{T}_{\beta} \right] = \frac{\rho}{\tau_1} \left( \frac{3}{2}\theta - \frac{3}{2}\theta_T \right), \end{equation}which when combined with the rotational energy equation shows that the total energy is conserved. This implies that the evolution equation for slow moments(4.4)\begin{equation} M_{slow} = \left\{\rho, \rho u_\alpha, \frac{1}{2} \rho u^2 + \frac{3+\delta}{2} \rho \theta \right\}, \end{equation}mass density, momentum density and total energy density are(4.5)\begin{equation} \left.\begin{gathered} \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,\\ \partial_t \left( E_T+E_R \right) + \partial_{\beta} \left[ \left( E_T+E_R+\rho\theta \right) u_\beta + \hat{\sigma}_{\beta\gamma}u_\gamma + q_{\beta} \right] = 0. \end{gathered}\right\} \end{equation}Thus, the conservation laws have the correct macroscopic form.(ii) $H$ theorem. For polyatomic gases, a part of entropy production should be due to rotational degrees of freedom. In the current model, as internal degrees of freedom are accounted for in a mean-field manner, similar to the Enskog equation one would expect that entropy contribution should only depend on rotational energy (Resibois Reference Resibois1978). Thus, we write a generalized $H$ function for polyatomic gas $H_1$ in Sackur–Tetrode form as a sum of Boltzmann part for monatomic contribution and rotational part $k \rho \ln \theta _R$ (Huang Reference Huang2009):
(4.6)\begin{equation} H_1 = H + k \rho \ln \theta_R, \end{equation}with $k$ being an unknown scale factor to be fixed later. On multiplying (4.1) with $\ln f$, and integrating over the velocity space, we obtain the evolution of $H$ as(4.7)\begin{equation} \partial_t H + \partial_\alpha J_\alpha^H = \varSigma^{ESBGK} + \frac{\tau}{\tau_1}\varSigma^{BGK} (\,f^{MB} (\rho, \boldsymbol{u}, \theta)) - \frac{3\rho}{2\tau_1} \frac{\theta-\theta_T}{\theta}, \end{equation}where $J_\alpha ^H$ is related to the entropy flux with $\varSigma ^{ESBGK}$ and $\varSigma ^{ BGK}$ being the entropy production due to the ES–BGK and the BGK terms respectively. The evolution of the rotational energy (second equation in (4.1)) can be rewritten as(4.8)\begin{equation} \partial_t \ln \theta_R + u_\alpha \partial_\alpha \ln \theta_R + \frac{2}{\delta \rho \theta_R} \partial_\alpha q^R_\alpha = \frac{ \theta - \theta_R }{\tau_1 \theta_R}. \end{equation}Multiplying the above equation with $\rho$ and exploiting continuity we obtain(4.9)\begin{equation} \partial_t \left( \rho \ln \theta_R \right) + \partial_\alpha \left( \rho u_\alpha \ln \theta_R + \frac{2}{\delta} \frac{q^R_\alpha}{\theta_R} \right) = \frac{\rho}{\tau_1} \frac{ \theta - \theta_R}{\theta_R} - \frac{2}{\delta} \frac{q^R_\alpha}{\theta_R^2} \partial_\alpha \theta_R. \end{equation}Thereby, adding (4.7) and (4.9) and using the form of $q^R_\alpha$ from (4.2), the evolution of $H_1$ is obtained as(4.10)\begin{equation} \partial_t H_1 + \partial_\alpha \left(J_\alpha^H + k \rho u_\alpha \ln \theta_R + k \frac{2}{\delta}\frac{q^R_\alpha}{\theta_R} \right) = \varSigma^{ESBGK} + \varSigma^{BGK} + \hat{\varSigma}, \end{equation}where the right-hand side is the net entropy production with contributions from the ES–BGK collision, the BGK collision and the rotational component of the model. Here,(4.11)\begin{equation} \hat{\varSigma} ={-}\frac{3\rho}{2\tau_1} \frac{\theta-\theta_T}{\theta} + \frac{k \rho}{\tau_1 } \frac{\theta - \theta_R }{\theta_R} + k \frac{2}{\delta}\kappa_R \left( \frac{\partial_\alpha \theta_R}{\theta_R} \right)^2. \end{equation}Similar to the standard BGK or ES–BGK case (Andries et al. Reference Andries, Le Tallec, Perlat and Perthame2000), the entropy production $\hat {\varSigma }$ in this model is also non-positive. This is achieved by choosing $k=-\delta /2$ and exploiting relation (3.3) to rewrite $\hat {\varSigma }$ as(4.12)\begin{equation} \hat{\varSigma} ={-} \frac{\delta}{2} \frac{\rho}{\tau_1} \frac{(\theta-\theta_R)^2}{\theta \theta_R} - \kappa_R \left( \frac{\partial_\alpha \theta_R}{\theta_R} \right)^2 \leq0. \end{equation}Hence, the proposed model satisfies the $H$ theorem.(iii) Hydrodynamics. In order to derive the hydrodynamic limit and the transport coefficients, the moments are typically categorized into fast $M_{fast}$ and slow $M_{slow}$ moments. The stress tensor and the heat flux constitute the relevant set of fast moments along with the translational and rotational temperatures as they are not conserved:
(4.13)\begin{equation} M_{fast} = \left\{\theta_T, \theta_R, \sigma_{\alpha \beta}, q_{\alpha} \right\}. \end{equation}The base state is obtained from zero of collision from (4.1) as(4.14a,b)\begin{equation} f=f^{MB}\implies \theta = \theta_T \quad {\rm and} \quad \theta = \theta_R. \end{equation}Thus, the fast moment can be expanded around the equilibrium values in a series as(4.15)\begin{equation} M_{fast} = M_{fast}(\,f^{MB}) + \tau M_{fast}^{(1)} + \cdots. \end{equation}In the Chapman–Enskog expansion, the time derivative of any quantity $\phi$ is expanded as(4.16)\begin{equation} \partial_t \phi = \partial_t^{(0)} \phi + \tau \partial_t^{(1)} \phi + {{O}} (\tau^2). \end{equation}The set of conservation laws (4.5) upon substituting the expansion of time derivative provide the definition of time derivative at ${{O}}(1)$ of slow variables as Euler equations:(4.17)\begin{equation} \left.\begin{gathered} \partial_t^{(0)} \rho + \partial_\alpha \left( \rho u_\alpha \right) = 0,\\ \partial_t^{(0)} \left( \rho u_\alpha \right) + \partial_\beta \left( \rho u_\alpha u_\beta + \rho \theta \delta_{\alpha \beta} \right) = 0,\\ \partial_t^{(0)} E + \partial_\beta \left( E u_\beta + \rho \theta u_\beta \right) = 0. \end{gathered} \right\} \end{equation}Thus, pressure evolution at $O(1)$ satisfies the adiabatic condition for a polyatomic gas:(4.18)\begin{equation} \left( \partial_t^{(0)} + u_{\beta} \partial_{\beta} \right) \left( \frac{p} {\rho^\gamma} \right) = 0,\quad \text{where } \gamma = \frac{5+\delta}{3+\delta}. \end{equation}Similarly, at order ${{O}}(\tau )$ we have(4.19)\begin{equation} \left.\begin{gathered} \partial_t^{(1)} \rho = 0, \\ \partial_t^{(1)} \left( \rho u_\alpha \right) + \partial_\alpha (\rho \theta^{(1)}_T) + \partial_\beta \sigma_{\alpha \beta}^{(1)} = 0,\\ \partial_t^{(1)} E + \partial_\beta \left( \sigma_{\alpha \beta}^{(1)} u_\alpha + \rho \theta^{(1)}_T u_\beta + q_\beta^{(1)}\right) = 0. \end{gathered}\right\} \end{equation}The expressions for $\rho \theta _T^{(1)}$ and $\sigma ^{(1)}_{\alpha \beta }$ can be obtained from the evolution equations of translational temperature (A13) and stress tensor (A15) respectively as (details of derivations in Appendix A)(4.20a,b)\begin{equation} \rho \theta_T^{(1)} ={-}\frac{2 \delta}{3(3+\delta)} \frac{\tau_1}{\tau} p \partial_\gamma u_\gamma, \quad \sigma^{(1)}_{\alpha \beta} ={-}\frac{2 p \overline{\partial_\beta u_\alpha}} {B-b}, \end{equation}where $B = 1 + \tau /\tau _1$. Substituting the above expressions in the momentum conservation equation, $O(\tau )$ hydrodynamics with shear viscosity $\eta$ and bulk viscosity $\eta _b$ is(4.21a,b)\begin{equation} \eta = \frac{p \tau}{B-b} {\quad} {\rm and} \quad \eta_b = \frac{2 \delta}{3(3+\delta)} p \tau_1. \end{equation}Similarly, translational thermal conductivity for the model is obtained from the translational heat flux evolution (A16). A Chapman–Enskog expansion of (A16) by substituting $q^T_\alpha = (q^T_\alpha )^{MB} + \tau {(q^T_\alpha )}^{(1)} + {{O}}(\tau ^2)$ at ${{O}}(1)$ yields(4.22)\begin{equation} \left(1 + \frac{\tau}{\tau_1} \right) {(q^T_\alpha)}^{(1)} ={-}\frac{5}{2} \rho \theta \partial_\alpha \theta. \end{equation}Thus, the translational thermal conductivity is $\kappa _T = 5 p \tau /(2B)$ which means the total thermal conductivity $\kappa = \kappa _T + \kappa _R$ with $k_r = \kappa _R/\kappa _T$ is $\kappa = {5 p \tau }/(2B) (1 + k_r )$. Hence,(4.23)\begin{equation} {Pr} = \frac{\eta C_p}{\kappa} = \left(\frac{B}{B-b} \right) \left(1+\frac{\delta}{5}\right) \left(\frac{1}{1+k_r} \right). \end{equation}To summarize, the parameter $\delta$ is determined directly from specific heat ratio $\gamma$. While the rotational relaxation time $\tau _1$ is fixed based on $\delta$ and bulk viscosity $\eta _b$. The translational relaxation time $\tau$ is then determined based on $\delta$, shear viscosity $\eta$ and $\tau _1$. The parameter $b$ is then fixed using $\tau$, $\tau _1$, ${Pr}$ and $k_r$ thereby imposing the target Prandtl number.
5. Discretizing via lattice Boltzmann method
In this section, we formulate a lattice Boltzmann scheme for solving the proposed model. Firstly, the velocity space is discretized into a discrete velocity set $\boldsymbol {c} = \{ \boldsymbol {c}_i, i=1 \dots N \}$ consisting of $N$ vectors. These vectors form the links of a lattice that should satisfy appropriate isotropy conditions (Succi Reference Succi2001; Atif, Namburi & Ansumali Reference Atif, Namburi and Ansumali2018). In particular, we validate the proposed kinetic model using a 41-velocity crystallographic lattice from Kolluru et al. (Reference Kolluru, Atif, Namburi and Ansumali2020b), which uses a body-centred cubic arrangement of grid points. The body-centred cubic lattice contains two simple cubic lattices offset by half the length of grid spacing in all directions for better spatial discretization (Namburi, Krithivasan & Ansumali Reference Namburi, Krithivasan and Ansumali2016). The weights for this lattice Boltzmann model are derived by imposing the following constraints on the discrete velocity set:
where $\theta _0$ is a reference lattice temperature. The discrete velocities and their corresponding weights for this RD3Q41 model are given in table 1.
The kinetic equation (4.1) in discrete in velocity space for populations $f_{1i}$ is
where
with moments of the discrete populations $f_{1i}$ defined as
To have a numerically efficient scheme for the $N$ coupled partial differential equations of (5.2), it is desirable to have large time steps, i.e. $\Delta t \gg \tau$. Upon integrating (5.2) along the characteristics and approximating the integral related to the collision term via the trapezoid rule, we obtain the implicit relation
which is made explicit by a transformation to an auxiliary population $g_{1i}(\boldsymbol {x},\boldsymbol {c},t ) = f_{1i} (\boldsymbol {x},\boldsymbol {c},t ) - ({\Delta t}/{2}) \varOmega _{i} (\boldsymbol {x},\boldsymbol {c},t )$. This implies the evolution equation for $g_{1i}$ is
with ${1}/{\tau ^*} = {1}/{\tau } + {1}/{\tau _1}$ and $\beta ^* = {\Delta t}/{(2 \tau ^* + \Delta t)}$. The moments of the auxiliary distribution $g_{1}$ are related to the moments of discrete populations $f_{1}$ as
To solve the second part of (4.1) that represents the internal energy, we write
by exploiting equations (3.3) and (4.2) and the continuity equation. The above equation is an advection–relaxation–diffusion equation that is solved by the steps detailed below:
(i) The relaxation equation
(5.9)\begin{equation} \frac{\partial \theta_{R}}{\partial t } = \frac{1}{\tau_1} \left(\theta - \theta_{R} \right) \end{equation}is first solved using the backward Euler method for a half-time step along with the relation in (3.3) to obtain(5.10)\begin{equation} \theta^{t+\Delta t/2}_{R} = \frac{1}{1+X} \theta^{t}_{R} + \frac{X}{1+X} \theta_T(\,f_{1}), \end{equation}where $X = 3 \Delta t / (2 \tau _1 (3+\delta ))$.(ii) The MacCormack scheme (MacCormack Reference MacCormack2003), which uses forward and backward differences for spatial derivatives in the predictor and corrector steps respectively, is used to solve the advection equation
(5.11)\begin{equation} \frac{\partial \theta_{R}}{\partial t } + u_\alpha \frac{\partial \theta_{R} }{\partial x_\alpha} = 0. \end{equation}(iii) The diffusion equation
(5.12)\begin{equation} \frac{\partial \theta_{R}}{\partial t} = \frac{2}{\rho \delta} \kappa_R \nabla^2 \theta_{ R} \end{equation}is then solved using the standard forward time centred space scheme to get $\theta ^{dif}_{R}$.(iv) Finally, the second part of relaxation is completed by an advance of $\theta ^{dif}_{R}$ by another half-time step $\Delta t/2$ leading to the final solution at $t+\Delta t$ as
(5.13)\begin{equation} \theta^{t+ \Delta t}_{R} = \frac{1}{1+X}\theta^{dif}_{R} + \frac{X}{1+X} \theta_T(\,f_{1}). \end{equation}
Note that the choice of the solver for the evolution equation of rotational energy is independent of the lattice Boltzmann solver used for solving $f_{1}$. The use of a scalar equation for rotational energy offers a significant computational advantage over using an additional distribution function, with memory usage and computational time reduced by a factor of 3–4. The difference between the two approaches lies in the treatment of internal degrees of freedom, but the accuracy remains the same for both formulations when the same flow solver is used (for low-Knudsen problems). The computational cost savings are entirely due to using a scalar variable instead of a full distribution function.
In contrast, solving kinetic equations with a $D3Q15$ velocity model requires at least five times more storage than advection–diffusion–reaction solvers, which typically work with 2–3 variables per grid point. We recall that both kinetic equation and advection–diffusion solvers are memory-bound; hence, the cost of floating point variables does not have a significant impact on the computational cost of either method. Instead, the number of memory read and write operations determines the computational cost.
A 15-velocity lattice Boltzmann model, for instance, requires 30 memory read/write operations per grid point for advection, and another 30 memory read/write operations and approximately 100 floating-point operations per grid point for the collision step. Thus, a simulation of the second kinetic equation requires 60 memory read/write operations and around 100 floating-point operations per grid point. Since most stencil-based finite-difference-type codes, including lattice Boltzmann codes, are memory-bound due to slower memory access in modern computers, the proposed model with advection–diffusion–relaxation would require only about 20–24 reads and writes per grid point. This leads to a 3–4 times reduction in computational time when using a scalar equation compared with a D3Q15 stencil.
In the next section, we validate the proposed numerical model by simulating a few benchmark problems related to acoustics, hydrodynamics and heat transfer such as propagation of an acoustic pulse, startup of a simple shear flow, thermal conduction and viscous heat dissipation.
6. Validation
In this section, we validate the model by contrasting simulation results with various benchmark results. As a first example, we consider a periodic domain $[-{\rm \pi},{\rm \pi} ]$ with $128 \times 4 \times 4$ lattice points to verify numerical sound speed. We initialize the domain with a pressure fluctuation of the form $p(x, t = 0) = p_0 (1 + \epsilon \cos (x))$ with $p_0 = \theta _0$. The pressure pulse is expected to reach the same state as the initial condition after one acoustic time period $(t_a)$. The $L_2$-norm of the pressure fluctuation is computed using the current state and initial state which is expected to be minimum when the waves are in phase. The number of time steps taken to achieve the least $L_2$-norm is used to compute $t_a$ and the speed of sound as $c_s = L/t_a$, where $L$ is the domain length. Thus, $\gamma = c^2_s/\theta _0$ is computed from the speed of sound and the lattice temperature. Here, we demonstrate the versatility of the model by simulating several real fluids by imposing the effective rotational degrees of freedom $\delta$ as given in table 2. We show in figure 1 that our model accurately recovers the specific heat ratio for various polyatomic gases, even for fractional (effective) rotational degrees of freedom. The proposed model remains accurate even for fractional rotational degrees of freedom, thereby achieving any target specific heat ratio values. In figure 2, we perform a grid convergence study for air and observe a second-order convergence. We perform additional validation studies by restricting our attention to diatomic gases with variable Prandtl number.
Next, we study the absorption of sound in a dissipative compressible medium. The presence of both viscosity and thermal conductivity leads to the dissipation of energy in the sound waves. For an emitted wave, the pressure perturbations $p'$ far away from the source decay during a finite time as (Landau & Lifshitz Reference Landau and Lifshitz1987)
where the dimensionless Landau number ${La}$ is (Ansumali, Karlin & Öttinger Reference Ansumali, Karlin and Öttinger2005)
Here, $\lambda$ is the ratio of bulk to shear viscosities and ${Kn}$ is the Knudsen number. The form of pressure perturbation shows that the wave profile is Gaussian-like at large distances and the width of the wave is proportional to $\sqrt {La}$ for a fixed domain length $L$.
To demonstrate the effectiveness of the lattice Boltzmann scheme, we perform a simulation at a fixed ${Kn}$ value of $10^{-3}$ on a domain of size $400 \times 400$ at Prandtl numbers $1.4$, $2, 5$ and $10$. We initialize the domain with a normal density perturbation of amplitude $0.001$ at the centre of the fluid of uniform density $1$ at rest. From (4.21a,b) and setting $\tau = (3/5) \tau _1$ one obtains $\lambda = {224}/(225 \,{Pr})$.
Using this relation and $\gamma = 7/5$, the Landau number ${La}$ is calculated as
The Landau numbers for the chosen set of parameters are listed in table 3. It is evident that $\rm {La}$ is a linear function of $1/{Pr}$ which suggests that the width of the Gaussian increases with a decrease in ${Pr}$. The pressure fluctuations far from the source of perturbation after $t=0.2 t_a$ for various Prandtl numbers are plotted in figure 3, where $t_a = L/c_s$ is the acoustic time scale. It is evident that as expected the width of the wave decreases with an increase in the Prandtl number.
Next, we investigate the propagation of an acoustic pulse in a diatomic gas with $\gamma =7/5$, where the isentropic speed of sound is $c_s = \sqrt {\gamma \theta }$. An axisymmetric pressure pulse is initialized at the centre of a domain of size $[-1,1]$ with $256\times 256\times 4$ grid points. The acoustic pulse is of the form
with $p_0 = \theta _0$, $\epsilon = 0.001$, $b = 0.1$, $\alpha = \ln (2)/b^2$ and $r=\sqrt {x^2 + y^2}$. For low amplitudes of pressure fluctuations and low viscosity, the exact form of the pressure fluctuation is known as the solution of the linearized Euler equations as (Tam & Webb Reference Tam and Webb1993)
where $\textrm {J}_0$ is the Bessel function of the first kind of zero order (Abramowitz & Stegun Reference Abramowitz and Stegun1965). Figure 4 shows that the pressure fluctuations from the simulation and the exact solution along the centreline of the $y$ axis are in agreement.
Next, we simulate the transient hydrodynamics in the startup of a simple shear flow between two flat plates separated by a distance $L$ on a grid of size $128 \times 64 \times 8$ with diffusive wall boundary condition (Ansumali & Karlin Reference Ansumali and Karlin2002) and periodicity in the other two directions. Here, the top plate is suddenly started with a velocity $u_w$ while the bottom plate remains stationary. The viscous effects play an important role in the development of the flow which is driven by momentum diffusion. Figure 5 contrasts the solutions at various diffusion times $t^* = t/(L^2/\nu )$ obtained from our simulations with the known analytical solution for the velocity (Pozrikidis & Jankowski Reference Pozrikidis and Jankowski1997):
As expected, the simulation recovers the analytical solution with good accuracy.
Next, we investigate the effects of thermal conduction by considering a set-up consisting of fluid confined in a square cavity of size $[L,L]$ with $128\times 128$ points and stationary walls. The top wall is maintained at a higher temperature $T_1$, while the other three walls are maintained at a temperature $T_0$ (${<}T_1$). Diffusive wall boundary conditions are applied in both directions. The analytical solution for the temperature profile at steady state is (Leal Reference Leal2007)
Figure 6 shows that the simulated temperature profiles along lines $x=0.1L$, $0.2L$ and $0.5L$ and along $y=0.25L$, $0.5L$ and $0.75L$ for a temperature difference of $0.1 \theta _0$ match well with the analytical solution.
Next, we validate our model for a thermal Couette flow problem to evaluate its capability in simulating viscous heat dissipation at various Prandtl numbers. We study the steady state of a flow induced by a wall at $y=H$ moving with a constant horizontal velocity $U_0$ and maintained at a constant elevated temperature $T_1$. The lower wall at $y=0$ is kept stationary at a constant temperature $T_0$ ($T_1 > T_0$). The analytical solution for the temperature profile for this set-up is (Bird et al. Reference Bird, Stewart, Lightfoot and Klingenberg2015)
where $\Delta T = T_1 - T_0$ is the temperature difference between the two walls and ${Ec} = U_0^2/(c_p \Delta T)$ is the Eckert number that represents the ratio of viscous dissipation to heat conduction with $c_p=7/2$ as the specific heat at constant pressure for a diatomic gas.
Simulations were performed for ${Pr} = 0.75$, 2.5, 5.0, $7.5$ and $10$ at Eckert number fixed at unity on a domain with $128$ grid points. The walls were maintained at temperatures $\theta _0+0.5\Delta \theta$ and $\theta _0-0.5\Delta \theta$ and plate velocity $U_0$ is chosen corresponding to a Mach number of $0.1$. Figure 7 compares the temperature profiles obtained analytically and via simulations and they are found to be in good agreement.
6.1. Transonic flows
A higher-order model (RD3Q167) capable of simulating transonic flows is employed (Hanumantharayappa et al. Reference Hanumantharayappa, Thantanapally, Namburi, Kumaran and Ansumali2021) for the following validation cases. The energy shells and their corresponding velocities with weights for this RD3Q167 model are given in table 4.
We present a Sod's test problem with a initial density ratio of $8$ and pressure ratio of $10$ on the left and right halves of the domain
on a grid of size $4096\times 4\times 4$ for this test case. In figure 8, we present the density, Mach number and pressure profiles at $\textrm {time}=0.1$ for specific heat ratios of helium ($\delta =0.03$), air ($\delta =1.96$) and methylal ($\delta =30.33$) as given in table 2 and contrast the solution with the exact Reimann solution.
To demonstrate the capability of the proposed model to simulate complex flows, we show a qualitative simulation of flow past a circular cylinder at a Reynolds number of 5000 defined based on diameter and free-stream velocity. The specific heat ratio for this simulation was chosen to be $1.22$ with additional rotational degrees of freedom $\delta = 6$. The Mach number for the present simulation is fixed at $0.45$ with 40 lattice points per diameter. We show the vorticity profile around the cylinder after 10 convection times in figure 9. It should be noted that the proposed formulation is quite stable even for such a severely under-resolved simulation.
7. Outlook
We have proposed a kinetic model for polyatomic gases with a tunable Prandtl number, by augmenting the ES–BGK model, an extension of the BGK model, at the level of the single-particle distribution function with an advection–diffusion–relaxation equation for the rotational energy. We show that close to the hydrodynamic limit, the internal degrees of freedom tend to be well represented just by rotational kinetic energy density and the proposed model recovers the compressible hydrodynamic equations of polyatomic gases as its macroscopic limit. It was shown that the transport coefficients of the model can be tuned for simulation of flows at different Prandtl numbers and specific heat ratios. A set of free parameters such as number of rotational degrees of freedom, translation relaxation time, rotational relaxation time, rotational thermal conductivity and a free parameter $b$ in the ES–BGK model can be used to tune the specific heat ratio, shear viscosity and bulk viscosity and set a target Prandtl number. This framework is general enough to deal with a more complex model of internal structures. We also demonstrated that the model respects the $H$ theorem. The simplicity of the model makes it suitable for lattice Boltzmann and other numerical implementations.
Acknowledgements
The support and resources provided by the PARAM Yukti Facility under the National Supercomputing Mission, Government of India, at the Jawaharlal Nehru Centre for Advanced Scientific Research are gratefully acknowledged.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Evolution equations
A.1. The BGK model
By taking appropriate moments of the Boltzmann BGK equation, we can see that the evolution of the stress and the heat flux are (Liboff Reference Liboff2003)
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\textrm {th}$-order moment involves $(M+1)\textrm {th}$ moment and thus forms an infinite-order-moment chain.
A.2. Two-population polyatomic model
For polyatomic gases, the energy equation gets an additional contribution from the rotational energy. In both the BGK and ES–BGK models, the evolution equations for the translational part of the energy and the rotational part of the energy $E_{R}=\delta \rho \theta _{R}/2$ are of the form
with $Z_E = 1$ for the BGK model and $Z_E = \alpha Z_{ES}$ for the ES–BGK model. The evolution equation for the total energy is written as the sum of (A2) as
where the relationship between translational and rotational temperatures (3.3) is used to show that energy is collisional invariant.
From (3.6) and (3.9), the stress evolution and the translational heat flux evolution equations in explicit form are
with $Z_q = 1$ for the BGK model and $Z_q = Z_{ES}$ for the ES–BGK model. Multiplying the rotational energy equation from (A2) with $u_\alpha$ and using velocity evolution we have
The rotational heat flux evolution can be obtained as the first moment of $f_2$ dynamics as
A.3. Proposed polyatomic model
We derive the evolution equations for kinetic energy, internal energy, pressure, stress, heat flux and translational and rotational temperatures for the proposed model. Using the momentum evolution equation (4.5), we obtain the evolution equation for $\rho u_\alpha u_\beta$ as
Evolution of kinetic energy obtained by taking the trace of the above equation is
Subtracting the above equation from the evolution equation of total energy from (4.5) gives the evolution equation for internal energy $e = (3+\delta )\rho \theta /2$ as
and the evolution equation of pressure $p = \rho \theta$ as
Using the pressure and continuity equation, the evolution equation for temperature $\theta$ is
From the evolution of kinetic energy (A8) and the translational temperature (4.3), the evolution for translational energy can be evaluated as
Using the continuity equation, the evolution for translational temperature $\theta _T$ is
For the evolution of the stress tensor $\sigma _{\alpha \beta }$, we multiply the kinetic equation (4.1) with $\xi _\alpha \xi _\beta$ and integrate over the velocity space to obtain
Thereafter, multiplying (A12) with $\delta _{\alpha \beta }$ and subtracting from the above equation, one obtains evolution of the stress tensor as
Similarly, the evolution equation for translational heat flux is obtained by multiplying the kinetic equation (4.1) with $\xi ^2 \xi _\alpha$ to obtain