Hostname: page-component-586b7cd67f-gb8f7 Total loading time: 0 Render date: 2024-11-26T05:19:06.631Z Has data issue: false hasContentIssue false

Reduced kinetic model of polyatomic gases

Published online by Cambridge University Press:  12 May 2023

Praveen Kumar Kolluru
Affiliation:
Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore 560064, India
Mohammad Atif
Affiliation:
Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore 560064, India
Santosh Ansumali*
Affiliation:
Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore 560064, India
*
Email address for correspondence: [email protected]

Abstract

Kinetic models of polyatomic gas typically account for the internal degrees of freedom at the level of the two-particle distribution function. However, close to the hydrodynamic limit, the internal (rotational) degrees of freedom tend to be well represented just by rotational kinetic energy density. We account for the rotational energy by augmenting the ellipsoidal statistical Bhatnagar–Gross–Krook (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. This reduced model respects the $H$ theorem and recovers the compressible hydrodynamics for polyatomic gases as its macroscopic limit. As required for a polyatomic gas model, this extension of the ES–BGK model not only has the correct specific heat ratio but also allows for three independent tunable transport coefficients: thermal conductivity, shear viscosity and bulk viscosity. We illustrate the effectiveness of the model via a lattice Boltzmann method implementation.

Type
JFM Papers
Copyright
© The Author(s), 2023. Published by Cambridge University Press

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:

(2.1)\begin{equation} \{ \rho, \rho \boldsymbol{u}, E \} = \left\langle \left\{ 1, \boldsymbol{c}, \frac{c^2}{2} \right\}, f \right\rangle, \end{equation}

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,

(2.2ac)\begin{equation} \varTheta_{ij}\equiv \left\langle \xi_i \xi_j, f\right\rangle,\quad \sigma_{\alpha \beta} = \left\langle \overline{\xi_{\alpha} \xi_{\beta}}, f\right\rangle, \quad q_{\alpha} = \left\langle\xi_{\alpha} \frac{\xi^2}{2}, f\right\rangle, \end{equation}

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

(2.3a,b)\begin{equation} \sigma_{\alpha \beta} = P_{\alpha \beta} - \rho u_{\alpha} u_{\beta} - \rho \theta \delta_{\alpha \beta}, \quad 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{equation}

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

(2.4)\begin{equation} \varOmega_{BGK} = \frac{1}{\tau} (\,f^{MB}(\rho,\boldsymbol{u}, \theta) - f). \end{equation}

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):

(2.5a,b)\begin{equation} \mu = p \tau ,\quad \kappa = \frac{5}{2} p \tau \quad \implies {Pr} = \frac{\mu C_p}{k} = 1. \end{equation}

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} )$:

(2.6)\begin{equation} f^{ES}\left(\rho, \boldsymbol{u}, \lambda_{ij} \right) = \frac{\rho}{\sqrt{{\rm{det}}[2 {\rm \pi}\lambda_{ij}]}} \exp\left({-\frac{1}{2} \xi_i \lambda^{{-}1}_{ij} \xi_j} \right), \end{equation}

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

(2.7a,b)\begin{equation} \mu = \frac{p \tau}{1-b} ,\quad \kappa = \frac{5}{2} p \tau\quad \implies {Pr} = \frac{1}{1-b}. \end{equation}

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):

(3.1)\begin{equation} \{ \rho, \rho \boldsymbol{u}, E_{T} + E_{R} \} = \left\langle\left\langle \left\{ 1, \boldsymbol{c}, \frac{c^2}{2} + I^{2/\delta} \right\}, F \right\rangle\right\rangle, \end{equation}

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:

(3.2)\begin{equation} 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. \end{equation}

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

(3.3)\begin{equation} \theta = \frac{3}{3+\delta} \theta_T + \frac{\delta}{3+\delta} \theta_R, \end{equation}

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

(3.4)\begin{equation} \partial_t F + c_\alpha \partial_\alpha F = \varOmega(F,F), \end{equation}

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)

(3.5)\begin{equation} \left.\begin{gathered} \varOmega_{BGK} = \frac{1}{\tau} (F^{MB}(\rho,\boldsymbol{u}, \theta ,I) - F ),\\ F^{MB}(\rho,\boldsymbol{u}, \theta , I) = \frac{\rho \varLambda_{\delta} }{\left(2 {\rm \pi}\theta \right )^{3/2} \theta^{\delta/2}} \exp{\left(-\left(\frac{(\boldsymbol{c}-\boldsymbol{u})^2}{2 \theta} + \frac{I^{2/\delta}}{\theta} \right)\right)}, \end{gathered}\right\} \end{equation}

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

(3.6)\begin{equation} \left.\begin{gathered} \partial_t f_1 + c_\alpha \partial_\alpha f_1 = \frac{1}{\tau} (\,f_1^{MB}(\rho,\boldsymbol{u}, \theta) - f_1),\\ \partial_t f_2 + c_\alpha \partial_\alpha f_2 = \frac{1}{\tau} \left( \frac{\delta}{2}\theta f_1^{MB} (\rho,\boldsymbol{u}, \theta) - f_2 \right). \end{gathered}\right\} \end{equation}

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

(3.7)\begin{equation} \varOmega_{ESBGK}(F) = \frac{ Z_{ES}}{\tau} (F^{ES} (\rho, u, \lambda_{ij}, \theta_{rel} ) - F), \end{equation}

where $\lambda _{ij} = (1-\alpha ) [ (1 -b) \theta _{T} \delta _{ij}+ b \varTheta _{ij} ] + \alpha \theta \delta _{ij}$ with a generalized Gaussian $F^{ES}$:

(3.8)\begin{equation} F^{ES}(\rho, u, \lambda_{ij}, \theta_{rel} ) = \frac{\rho \varLambda_{\delta} }{\theta_{rel}^{\delta/2} \sqrt{{\rm{det}}[2 {\rm \pi}\lambda_{ij}]}} \exp\left({-\frac{1}{2} \xi_i \lambda^{{-}1}_{ij} \xi_j} - \frac{I^{2/\delta}}{\theta_{rel}} \right), \end{equation}

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

(3.9)\begin{equation} \left.\begin{gathered} \partial_t f_1 + c_\alpha \partial_\alpha f_1 = \frac{ Z_{ES}}{\tau} (\,f^{ES}(\rho, u, \lambda_{ij}) - f_1),\\ \partial_t f_2 + c_\alpha \partial_\alpha f_2 = \frac{Z_{ES}}{\tau} \left( \frac{\delta}{2} \theta_{rel} f^{ES}(\rho, u, \lambda_{ij}) - f_2 \right). \end{gathered}\right\} \end{equation}

A Chapman–Enskog analysis shows that the momentum equation yields the familiar compressible Navier–Stokes equation form as

(3.10)\begin{equation} \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, \end{equation}

with the shear and bulk viscosities

(3.11a,b)\begin{equation} \eta = p \tau, \quad \frac{\eta_b}{\eta} = \frac{2 \delta}{3(3+\delta)Z_E}. \end{equation}

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

(4.1) \begin{equation} \left.\begin{gathered} \partial_t f_1 + c_\alpha \partial_\alpha f_1 = \frac{1}{\tau} (\, f^{ES}(\rho ,\boldsymbol{u}, \theta_T \delta_{\alpha \beta} + b \sigma_{\alpha \beta} ) - f_1) + \frac{1}{\tau_1} (\, f^{MB}(\rho, \boldsymbol{u}, \theta ) - f_1 ),\\ \partial_t\left( E_R \right) + \partial_\alpha ( E_R u_\alpha + q^{R}_\alpha ) = \frac{1}{\tau_1} \left( \frac{\delta}{2} \rho\theta - E_R \right), \end{gathered}\right\} \end{equation}

with the form of heat flux due to internal degrees of freedom as

(4.2)\begin{equation} q^{R}_\alpha ={-} \kappa_R \partial_\alpha \theta_R. \end{equation}

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:

  1. (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.
  2. (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.
  3. (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:

(5.1)\begin{equation} \left.\begin{gathered} \sum w_i = 1, \quad \sum w_i c_{ix}^2 = \theta_0, \quad \sum w_i c_{ix}^4 = 3\theta_0^2, \\ \sum w_i c_{ix}^2 c_{iy}^2 = \theta_0^2, \quad \sum w_i c_{ix}^4 c_{i}^2 = 21\theta_0^3, \\ \sum w_i c_{ix}^2 c_{i}^4 = 35\theta_0^3, \quad \sum w_i c_{i}^8 = 945\theta_0^4, \end{gathered}\right\}\end{equation}

where $\theta _0$ is a reference lattice temperature. The discrete velocities and their corresponding weights for this RD3Q41 model are given in table 1.

Table 1. Velocities and their corresponding weights for the RD3Q41 model with $\theta _0 = 0.2948964908710633$.

The kinetic equation (4.1) in discrete in velocity space for populations $f_{1i}$ is

(5.2)\begin{equation} \partial_t f_{1i}\left(\boldsymbol{x},\boldsymbol{c},t \right) + c_\alpha \partial_\alpha f_{1i}\left(\boldsymbol{x},\boldsymbol{c},t \right) = \varOmega_{i} (\boldsymbol{x},t), \end{equation}

where

(5.3) \begin{equation} \varOmega_{i}\left(\boldsymbol{x},t \right) = \frac{1}{\tau} [ f^{ES}_{1i}(\rho ,\boldsymbol{u}, \theta_T \delta_{\alpha \beta} + b \sigma_{\alpha \beta} ) - f_{1i} ] + \frac{1}{\tau_1} [f_{1i}^{Eq}(\rho, \boldsymbol{u}, \theta) - f_{1i}], \end{equation}

with moments of the discrete populations $f_{1i}$ defined as

(5.4ac)\begin{align} \rho(\,f_{1}) = \sum_i f_{1i}, \quad \rho u_\alpha(\,f_{1}) = \sum_i f_{1i} c_{i \alpha}, \quad \theta_T (\,f_{1}) = \frac{1}{3\rho(\,f_{1})} \left( \sum_i f_{1i} c^2_{i} - \rho {\boldsymbol{u}}^2 \right). \end{align}

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

(5.5)\begin{equation} f_{1i}(\boldsymbol{x} + \boldsymbol{c} \Delta t, t + \Delta t) = f_{1i}( \boldsymbol{x}, t) + \frac{\Delta t}{2} \left[ \varOmega_i (\boldsymbol{x},t) + \varOmega_i (\boldsymbol{x} + \boldsymbol{c}t,t+\Delta t) \right], \end{equation}

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

(5.6)\begin{equation} g_{1i}\left(\boldsymbol{x}+ {\boldsymbol{c}}_i \Delta t, t+\Delta t \right) = g_{1i}\left(\boldsymbol{x},t \right) \left( 1 - 2\beta^* \right) + 2 \tau^* \beta^* \varOmega_i\left(\boldsymbol{x},\boldsymbol{c},t \right), \end{equation}

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

(5.7)\begin{gather} \left.\begin{gathered} \rho(g_{1}) = \rho(\,f_1), \quad \boldsymbol{u}(g_1) = \boldsymbol{u}(\,f_1), \quad \theta_T(g_1) = \theta_T(\,f_1) + \frac{\delta \Delta t}{2\tau_1 (3 +\delta )} \left( \theta_T(\,f_1) - \theta_R \right) , \\ \sigma_{\alpha \beta}(g_1) = \sigma_{\alpha \beta}(\,f_1) \left(1 + \frac{\Delta t}{2 \tau^*} - \frac{\Delta t}{2\tau} b \right). \end{gathered}\right\} \end{gather}

To solve the second part of (4.1) that represents the internal energy, we write

(5.8)\begin{equation} \frac{\partial \theta_{R}}{\partial t } + u_\alpha \frac{\partial \theta_{R} }{\partial x_\alpha} = \frac{1}{\tau_1} \left(\theta - \theta_{R} \right) + \frac{2}{\rho \delta} \kappa_R \nabla^2 \theta_{R}, \end{equation}

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:

  1. (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 ))$.
  2. (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}
  3. (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}$.
  4. (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.

Figure 1. Specific heat ratio in simulating sound propagation in different gases. The line represents the reference value with $\delta$ the number of effective rotational degrees of freedom for various gases as listed in table 2.

Figure 2. Grid convergence study showing a second-order convergence for air.

Table 2. Specific heat ratios of real fluids (Green & Southard Reference Green and Southard2019) and their effective rotational degrees of freedom.

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)

(6.1)\begin{equation} p'(r,t) \propto \left({La} \, rL \right)^{{-}1/2} \exp\left(- \frac{\left(r - c_s t \right)^2}{2 {La} \,rL } \right), \end{equation}

where the dimensionless Landau number ${La}$ is (Ansumali, Karlin & Öttinger Reference Ansumali, Karlin and Öttinger2005)

(6.2)\begin{equation} {La} = {Kn} \left(\frac{4}{3} + \lambda \right) + \frac{Kn}{Pr} \left(\gamma - 1 \right). \end{equation}

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

(6.3)\begin{equation} {La} = {Kn} \left(\frac{4}{3} +\frac{314}{225} \frac{1}{Pr} \right). \end{equation}

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.

Figure 3. (a) Pressure perturbations versus Pr with (b) a zoomed plot. The width of the Gaussian wave increases with a decrease in ${Pr}$.

Table 3. Variation of ${La}$ with ${Pr}$ at ${Kn} = 10^{-3}$.

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

(6.4)\begin{equation} p(x,y,t=0) = p_0 (1 + \epsilon {\rm e}^{-\alpha r^2}), \end{equation}

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)

(6.5)\begin{equation} p'(x,y,t) = p_0 \times \frac{\epsilon}{2 \alpha} \int_0^{\infty} \exp\left(\frac{-\xi^2}{4 \alpha}\right) \cos(c_s \xi t) {\rm J}_0(\xi r) \xi\,{\rm d}\xi, \end{equation}

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.

Figure 4. Comparison of the pressure fluctuation along the centreline at time $t^*$ from lattice Boltzmann simulation (points) and exact solution (line).

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):

(6.6)\begin{equation} u^* = \frac{u}{u_w} = \frac{y}{L} - \frac{2}{\rm \pi} \sum_{k=1}^\infty \left[ \frac{1}{k} \exp\left( - k^2 {\rm \pi}^2 \frac{\nu t}{L^2}\right) \sin\left(k {\rm \pi}\left(1 - \frac{y}{L} \right) \right) \right]. \end{equation}

As expected, the simulation recovers the analytical solution with good accuracy.

Figure 5. Transients in a planar Couette flow.

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)

(6.7)\begin{equation} \frac{T - T_0}{T_1 - T_0} = \frac{2}{\rm \pi} \sum_{n=1}^{\infty} \frac{({-}1)^{n+1}+1}{n} \sin\left({n{\rm \pi} x}\right) \frac{\sinh(n{\rm \pi} y)}{\sinh(n {\rm \pi})}. \end{equation}

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.

Figure 6. Temperature profiles at $y=0.25L$, $0.5L$ and $0.75L$ (a) and at $x=0.1L$, $0.2L$ and $0.5L$ (b) in a two-dimensional heated cavity.

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)

(6.8)\begin{equation} \frac{T - T_0}{\Delta T} = \frac{y}{H} + \frac{Pr \,Ec}{2} \frac{y}{H} \left(1 - \frac{y}{H} \right), \end{equation}

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.

Figure 7. Temperature profiles at steady state (symbols) compared with the analytical solution (lines) at varying Prandtl numbers.

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.

Table 4. Velocities and their corresponding weights for the RD3Q167 model with $\theta _0 = 0.7374708021487686$.

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

(6.9a,b)\begin{equation} \begin{pmatrix} \rho_L \\ u_L \\ p_L \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix},\quad \begin{pmatrix} \rho_R \\ u_R \\ p_R \end{pmatrix} = \begin{pmatrix} 0.125 \\ 0 \\ 0.1 \end{pmatrix} \end{equation}

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.

Figure 8. Density, Mach number and pressure after $\textrm {time}=0.1$ for Sod's test problem.

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.

Figure 9. Vorticity profile around a circular cylinder after 10 convection times.

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)

(A1)\begin{gather} \left.\begin{gathered} \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} + 2 p \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{gathered}\right\} \end{gather}

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

(A2)\begin{equation} \left.\begin{gathered} \partial_t E_T + \partial_{\alpha} [\left( E_T + \rho \theta \right) u_\alpha + q_\alpha^{T} + u_\gamma \hat{\sigma}_{\alpha \gamma}] = \frac{Z_E}{\tau} \frac{3 \rho}{2} \left( \theta - \theta_T \right) ,\\ \partial_t E_R + \partial_\alpha (E_R u_\alpha + q^R_\alpha) = \frac{Z_E}{\tau} \frac{\delta \rho }{2} \left( \theta - \theta_R \right), \end{gathered}\right\} \end{equation}

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

(A3)\begin{equation} \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, \end{equation}

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

(A4)\begin{gather} \left.\begin{gathered} \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{gathered}\right\} \end{gather}

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

(A5)\begin{gather} \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{gather}

The rotational heat flux evolution can be obtained as the first moment of $f_2$ dynamics as

(A6)\begin{align} \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(\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. \end{align}

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

(A7)\begin{equation} \partial_t \left( \rho u_\alpha u_\beta \right) + \partial_\gamma (\rho u_\alpha u_\beta u_\gamma) + u_\alpha \partial_\gamma ( \rho \theta_T \delta_{\beta \gamma} + \sigma_{\beta \gamma} ) + u_\beta \partial_\gamma (\rho \theta_T \delta_{\alpha \gamma} + \sigma_{\alpha \gamma} ) = 0 . \end{equation}

Evolution of kinetic energy obtained by taking the trace of the above equation is

(A8)\begin{equation} \partial_t \left(\tfrac{1}{2}\rho u^2\right) + \partial_\gamma \left( \tfrac{1}{2} \rho u^2 u_\gamma \right) + u_\beta \partial_\gamma (\rho \theta_T \delta_{\beta \gamma} + \sigma_{\beta \gamma} ) = 0. \end{equation}

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

(A9)\begin{equation} \partial_t e + \partial_\beta \left( e u_\beta\right) + \partial_\beta q_\beta + \sigma_{\beta \gamma} \partial_\gamma u_\beta + \rho \theta_T \partial_\gamma u_\gamma = 0 \end{equation}

and the evolution equation of pressure $p = \rho \theta$ as

(A10)\begin{equation} \partial_t p + \partial_\beta \left( p u_\beta\right) + \left(\frac{2}{3+\delta}\right) \left(\partial_\beta q_\beta + \sigma_{\beta \gamma} \partial_\gamma u_\beta + \rho \theta_T \partial_\gamma u_\gamma \right) = 0. \end{equation}

Using the pressure and continuity equation, the evolution equation for temperature $\theta$ is

(A11)\begin{equation} \partial_t \theta + u_\beta \partial_\beta \theta + \left(\frac{2}{(3+\delta)\rho}\right) \left( \partial_\beta q_\beta + \sigma_{\beta \gamma} \partial_\gamma u_\beta + \rho \theta_T \partial_\gamma u_\gamma \right) = 0. \end{equation}

From the evolution of kinetic energy (A8) and the translational temperature (4.3), the evolution for translational energy can be evaluated as

(A12)\begin{align} \partial_t \left( \frac{3 \rho \theta_T}{2} \right) + \partial_{\beta} \left( \frac{3 \rho \theta_T}{2} u_\beta \right) + \rho \theta_T \partial_\beta u_\beta + \sigma_{\beta\gamma} \partial_\beta u_\gamma + \partial_\beta q^{T}_{\beta} = \frac{\rho}{\tau_1} \left(\frac{3}{2}\theta - \frac{3}{2}\theta_T \right). \end{align}

Using the continuity equation, the evolution for translational temperature $\theta _T$ is

(A13)\begin{equation} \partial_t \theta_T + u_\beta \partial_{\beta} \theta_T + \frac{2}{3} \theta_T \partial_\beta u_\beta + \frac{2}{3 \rho} \sigma_{\beta\gamma} \partial_\beta u_\gamma + \frac{2}{3 \rho} \partial_\beta q^{T}_{\beta} = \frac{1}{\tau_1} \left(\theta - \theta_T \right). \end{equation}

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

(A14)\begin{align} & \partial_t (\rho \theta_T \delta_{\alpha \beta}+ \sigma_{\alpha \beta}) + \partial_\kappa Q^{T}_{\alpha \beta \kappa} + \partial_\kappa \left(u_\kappa (\rho \theta_T \delta_{\alpha \beta} + \sigma_{\alpha \beta} )\right) + (\rho \theta_T \delta_{\kappa \beta} + \sigma_{\kappa \beta} ) \partial_\kappa u_\alpha \nonumber\\ &\quad + (\rho \theta_T \delta_{\kappa \alpha} + \sigma_{\kappa \alpha}) \partial_\kappa u_\beta = \frac{1}{\tau} \left( b -1 \right) \sigma_{\alpha \beta} +\frac{1}{\tau_1} \left( \rho \theta \delta_{\alpha \beta} - \rho \theta_T \delta_{\alpha \beta} - \sigma_{\alpha \beta} \right). \end{align}

Thereafter, multiplying (A12) with $\delta _{\alpha \beta }$ and subtracting from the above equation, one obtains evolution of the stress tensor as

(A15)\begin{align} &\partial_t (\sigma_{\alpha \beta}) + \partial_\gamma \left(u_\gamma \sigma_{\alpha \beta} \right) + \partial_\gamma \overline{Q^{T}_{\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}}\nonumber\\ &\quad = \left(\frac{1}{\tau}\left(b-1\right) - \frac{1}{\tau_1} \right) \sigma_{\alpha \beta}. \end{align}

Similarly, the evolution equation for translational heat flux is obtained by multiplying the kinetic equation (4.1) with $\xi ^2 \xi _\alpha$ to obtain

(A16) \begin{align} & \partial_t q^{T}_\alpha + \partial_\beta (u_\beta q^{T}_\alpha ) + \overline{Q^{T}_{\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 \nonumber\\ &\quad - \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} ={-}\left(\frac{1}{\tau} + \frac{1}{\tau_1} \right) q^{T}_\alpha. \end{align}

References

Abramowitz, M. & Stegun, I.A. 1965 Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, vol. 55. Courier Corporation.Google Scholar
Agrawal, S., Singh, S.K. & Ansumali, S. 2020 Fokker–Planck model for binary mixtures. J. Fluid Mech. 899, A25.CrossRefGoogle Scholar
Andries, P., Le Tallec, P., Perlat, J.-P. & Perthame, B. 2000 The Gaussian–BGK model of Boltzmann equation with small Prandtl number. Eur. J. Mech. (B/Fluids) 19 (6), 813830.CrossRefGoogle Scholar
Ansumali, S., Arcidiacono, S., Chikatamarla, S.S., Prasianakis, N.I., Gorban, A.N. & Karlin, I.V. 2007 a Quasi-equilibrium lattice Boltzmann method. Eur. Phys. J. B 56 (2), 135139.CrossRefGoogle Scholar
Ansumali, S., Karlin, I.V., Arcidiacono, S., Abbas, A. & Prasianakis, N.I. 2007 b Hydrodynamics beyond Navier–Stokes: exact solution to the lattice Boltzmann hierarchy. Phys. Rev. Lett. 98 (12), 124502.CrossRefGoogle Scholar
Ansumali, S. & Karlin, I.V. 2002 Kinetic boundary conditions in the lattice Boltzmann method. Phys. Rev. E 66 (2), 026311.CrossRefGoogle ScholarPubMed
Ansumali, S., Karlin, I.V. & Öttinger, H.C. 2005 Thermodynamic theory of incompressible hydrodynamics. Phys. Rev. Lett. 94 (8), 080602.CrossRefGoogle ScholarPubMed
Arima, T., Taniguchi, S., Ruggeri, T. & Sugiyama, M. 2012 Extended thermodynamics of dense gases. Contin. Mech. Thermodyn. 24 (4–6), 271292.CrossRefGoogle Scholar
Atif, M., Kolluru, P.K. & Ansumali, S. 2022 Essentially entropic lattice Boltzmann model: theory and simulations. Phys. Rev. E 106, 055307.CrossRefGoogle ScholarPubMed
Atif, M., Kolluru, P.K., Thantanapally, C. & Ansumali, S. 2017 Essentially entropic lattice Boltzmann model. Phys. Rev. Lett. 119, 240602.CrossRefGoogle ScholarPubMed
Atif, M., Namburi, M. & Ansumali, S. 2018 Higher-order lattice Boltzmann model for thermohydrodynamics. Phys. Rev. E 98, 053311.CrossRefGoogle Scholar
von Backstrom, T.W. 2008 The effect of specific heat ratio on the performance of compressible flow turbo-machines. In ASME Turbo Expo 2008, pp. 2111–2117. ASME.CrossRefGoogle Scholar
Bernard, F., Iollo, A. & Puppo, G. 2019 BGK polyatomic model for rarefied flows. J. Sci. Comput. 78 (3), 18931916.CrossRefGoogle Scholar
Bhatnagar, P.L., Gross, E.P. & Krook, M. 1954 A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94 (3), 511.CrossRefGoogle Scholar
Bird, R.B., Stewart, W.E., Lightfoot, E.N. & Klingenberg, D.J. 2015 Introductory Transport Phenomena, vol. 1. Wiley.Google Scholar
Brull, S. & Schneider, J. 2009 On the ellipsoidal statistical model for polyatomic gases. Contin. Mech. Thermodyn. 20 (8), 489508.CrossRefGoogle Scholar
Cercignani, C. 1988 The Boltzmann equation. In The Boltzmann Equation and its Applications, pp. 40–103. Springer.CrossRefGoogle Scholar
Chapman, S. & Cowling, T.G. 1970 The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases. Cambridge University Press.Google Scholar
Chen, F., Xu, A., Zhang, G., Li, Y. & Succi, S. 2010 Multiple-relaxation-time lattice Boltzmann approach to compressible flows with flexible specific-heat ratio and Prandtl number. Europhys. Lett. 90 (5), 54003.CrossRefGoogle Scholar
Gorban, A.N. & Karlin, I.V. 1994 General approach to constructing models of the Boltzmann equation. Physica A 206 (3–4), 401420.CrossRefGoogle Scholar
Grad, H. 1958 Principles of the kinetic theory of gases. In Thermodynamik der Gase/Thermodynamics of Gases, pp. 205–294. Springer.CrossRefGoogle Scholar
Green, D.W. & Southard, M.Z. 2019 Perry's Chemical Engineers’ Handbook. McGraw-Hill Education.Google Scholar
Hanumantharayappa, M.N., Thantanapally, C., Namburi, M., Kumaran, V. & Ansumali, S. 2021 LES/DNS of flow past T106 LPT cascade using a higher-order LB model. AIAA Paper 2021-3485.CrossRefGoogle Scholar
Holway, L.H. Jr. 1965 Kinetic theory of shock structure using an ellipsoidal distribution function. Rarefied Gas Dyn. 1, 193.Google Scholar
Holway, L.H. Jr. 1966 New statistical models for kinetic theory: methods of construction. Phys. Fluids 9 (9), 16581673.CrossRefGoogle Scholar
Huang, K. 2009 Introduction to Statistical Physics. Chapman and Hall/CRC.CrossRefGoogle Scholar
Huang, R., Lan, L. & Li, Q. 2020 Lattice Boltzmann simulations of thermal flows beyond the Boussinesq and ideal-gas approximations. Phys. Rev. E 102 (4), 043304.CrossRefGoogle ScholarPubMed
Kataoka, T. & Tsutahara, M. 2004 Lattice Boltzmann model for the compressible Navier–Stokes equations with flexible specific-heat ratio. Phys. Rev. E 69 (3), 035701.CrossRefGoogle ScholarPubMed
Kolluru, P.K., Atif, M. & Ansumali, S. 2020 a Extended BGK model for diatomic gases. J. Comput. Sci. 45, 101179.CrossRefGoogle Scholar
Kolluru, P.K., Atif, M., Namburi, M. & Ansumali, S. 2020 b Lattice Boltzmann model for weakly compressible flows. Phys. Rev. E 101 (1), 013309.CrossRefGoogle ScholarPubMed
Kuščer, I. 1989 A model for rotational energy exchange in polyatomic gases. Physica A 158 (3), 784800.CrossRefGoogle Scholar
Landau, L.D. & Lifshitz, E.M. 1987 Fluid mechanics. Translated from the Russian by J.B. Sykes and W.H. Reid. Course Theor. Phys. 6, 300303.Google Scholar
Larina, I.N. & Rykov, V.A. 2010 Kinetic model of the Boltzmann equation for a diatomic gas with rotational degrees of freedom. Comput. Maths Math. Phys. 50 (12), 21182130.CrossRefGoogle Scholar
Leal, L.G. 2007 Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes, vol. 7. Cambridge University Press.CrossRefGoogle Scholar
Lebowitz, J.L., Frisch, H.L. & Helfand, E. 1960 Nonequilibrium distribution functions in a fluid. Phys. Fluids 3 (3), 325338.CrossRefGoogle Scholar
Levermore, C.D. 1996 Moment closure hierarchies for kinetic theories. J. Stat. Phys. 83 (5–6), 10211065.CrossRefGoogle Scholar
Liboff, R.L. 2003 Kinetic Theory: Classical, Quantum, and Relativistic Descriptions. Springer Science & Business Media.Google Scholar
Liepmann, H.W., Narasimha, R. & Chahine, M.T. 1962 Structure of a plane shock layer. Phys. Fluids 5 (11), 13131324.CrossRefGoogle Scholar
MacCormack, R. 2003 The effect of viscosity in hypervelocity impact cratering. J. Spacecr. Rockets 40 (5), 757763.CrossRefGoogle Scholar
McQuarrie, D.A. 2000 Statistical Mechanics. University Science Books.Google Scholar
Morse, T.F. 1964 Kinetic model for gases with internal degrees of freedom. Phys. Fluids 7 (2), 159169.CrossRefGoogle Scholar
Mott-Smith, H.M. 1951 The solution of the Boltzmann equation for a shock wave. Phys. Rev. 82 (6), 885.CrossRefGoogle Scholar
Müller, I. & Ruggeri, T. 2013 Rational Extended Thermodynamics. Springer Science & Business Media.Google Scholar
Namburi, M., Krithivasan, S. & Ansumali, S. 2016 Crystallographic lattice Boltzmann method. Sci. Rep. 6, 27172.CrossRefGoogle ScholarPubMed
Nie, X., Shan, X. & Chen, H. 2008 Thermal lattice Boltzmann model for gases with internal degrees of freedom. Phys. Rev. E 77 (3), 035701.CrossRefGoogle ScholarPubMed
Oh, C.K., Oran, E.S. & Sinkovits, R.S. 1997 Computations of high-speed, high Knudsen number microchannel flows. J. Thermophys. Heat Transfer 11 (4), 497505.CrossRefGoogle Scholar
Pozrikidis, C. & Jankowski, D. 1997 Introduction to Theoretical and Computational Fluid Dynamics, vol. 675. Oxford University Press.Google Scholar
Pullin, D.I. 1978 Kinetic models for polyatomic molecules with phenomenological energy exchange. Phys. Fluids 21 (2), 209216.CrossRefGoogle Scholar
Resibois, P. 1978 H-theorem for the (modified) nonlinear Enskog equation. J. Stat. Phys. 19 (6), 593609.CrossRefGoogle Scholar
Ruggeri, T. & Sugiyama, M. 2015 Rational Extended Thermodynamics Beyond the Monatomic Gas. Springer.CrossRefGoogle Scholar
Rykov, V.A. 1975 A model kinetic equation for a gas with rotational degrees of freedom. Fluid Dyn. 10 (6), 959966.CrossRefGoogle Scholar
Shakhov, E.M. 1968 Generalization of the Krook kinetic relaxation equation. Fluid Dyn. 3 (5), 9596.CrossRefGoogle Scholar
Singh, S.K. & Ansumali, S. 2015 Fokker–Planck model of hydrodynamics. Phys. Rev. E 91 (3), 033303.CrossRefGoogle ScholarPubMed
Singh, S.K., Thantanapally, C. & Ansumali, S. 2016 Gaseous microflow modeling using the Fokker–Planck equation. Phys. Rev. E 94 (6), 063307.CrossRefGoogle ScholarPubMed
Struchtrup, H. 2004 Stable transport equations for rarefied gases at high orders in the Knudsen number. Phys. Fluids 16 (11), 39213934.CrossRefGoogle Scholar
Succi, S. 2001 The Lattice Boltzmann Equation: for Fluid Dynamics and Beyond. Oxford University Press.Google Scholar
Tam, C.K.W. & Webb, J.C. 1993 Dispersion-relation-preserving finite difference schemes for computational acoustics. J. Comput. Phys. 107 (2), 262281.CrossRefGoogle Scholar
Tsutahara, M., Kataoka, T., Shikata, K. & Takada, N. 2008 New model and scheme for compressible fluids of the finite difference lattice Boltzmann method and direct simulations of aerodynamic sound. Comput. Fluids 37 (1), 7989.CrossRefGoogle Scholar
Wang, Z., Yan, H., Li, Q. & Xu, K. 2017 Unified gas-kinetic scheme for diatomic molecular flow with translational, rotational, and vibrational modes. J. Comput. Phys. 350, 237259.CrossRefGoogle Scholar
Wang-Chang, C.S. & Uhlenbeck, G.E. 1951 Transport phenomena in polyatomic gases. Research Rep. CM-681. University of Michigan Engineering.Google Scholar
Watari, M. 2007 Finite difference lattice Boltzmann method with arbitrary specific heat ratio applicable to supersonic flow simulations. Physica A 382 (2), 502522.CrossRefGoogle Scholar
Wu, L., White, C., Scanlon, T.J., Reese, J.M. & Zhang, Y. 2015 A kinetic model of the Boltzmann equation for non-vibrating polyatomic gases. J. Fluid Mech. 763, 2450.CrossRefGoogle Scholar
Figure 0

Table 1. Velocities and their corresponding weights for the RD3Q41 model with $\theta _0 = 0.2948964908710633$.

Figure 1

Figure 1. Specific heat ratio in simulating sound propagation in different gases. The line represents the reference value with $\delta$ the number of effective rotational degrees of freedom for various gases as listed in table 2.

Figure 2

Figure 2. Grid convergence study showing a second-order convergence for air.

Figure 3

Table 2. Specific heat ratios of real fluids (Green & Southard 2019) and their effective rotational degrees of freedom.

Figure 4

Figure 3. (a) Pressure perturbations versus Pr with (b) a zoomed plot. The width of the Gaussian wave increases with a decrease in ${Pr}$.

Figure 5

Table 3. Variation of ${La}$ with ${Pr}$ at ${Kn} = 10^{-3}$.

Figure 6

Figure 4. Comparison of the pressure fluctuation along the centreline at time $t^*$ from lattice Boltzmann simulation (points) and exact solution (line).

Figure 7

Figure 5. Transients in a planar Couette flow.

Figure 8

Figure 6. Temperature profiles at $y=0.25L$, $0.5L$ and $0.75L$ (a) and at $x=0.1L$, $0.2L$ and $0.5L$ (b) in a two-dimensional heated cavity.

Figure 9

Figure 7. Temperature profiles at steady state (symbols) compared with the analytical solution (lines) at varying Prandtl numbers.

Figure 10

Table 4. Velocities and their corresponding weights for the RD3Q167 model with $\theta _0 = 0.7374708021487686$.

Figure 11

Figure 8. Density, Mach number and pressure after $\textrm {time}=0.1$ for Sod's test problem.

Figure 12

Figure 9. Vorticity profile around a circular cylinder after 10 convection times.