Reduced kinetic model of polyatomic gases

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.

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 1970; Cercignani 1988). Unlike the continuum Navier–Stokes–Fourier hydrodynamics equation, the Boltzmann equation is a valid description even at highly non-equilibrium states (Mott-Smith 1951; Liepmann, Narasimha & Chahine 1962; Cercignani 1988), encountered in the presence of strong shock waves at high Mach number (ratio of flow speed to sound speed) and in a highly rarefied flow characterized by a large Knudsen number (ratio of the mean free path to characteristic length scale) (Oh, Oran & Sinkovits 1997; Struchtrup 2004; Ansumali et al. 2007b). However, any analysis of the integro-differential Boltzmann equation is a formidable task even for the simplest problems. Thus, one often models the Boltzmann dynamics via a simplified collision term that converts the evolution equation to a partial differential equation (Bhatnagar, Gross & Krook 1954; Lebowitz, Frisch & Helfand 1960; Holway 1966; Shakhov 1968; Gorban & Karlin 1994; Levermore 1996; Andries et al. 2000; Ansumali et al. 2007a; Singh & Ansumali 2015; Agrawal, Singh & Ansumali 2020). An important example is the Bhatnagar–Gross–Krook (BGK) model (Bhatnagar et al. 1954), 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 1966; Andries et al. 2000), the quasi-equilibrium models (Gorban & Karlin 1994; Levermore 1996), the Shakhov model (Shakhov 1968) and the Fokker–Planck model (Singh & Ansumali 2015; Singh, Thantanapally & Ansumali 2016) are used as kinetic models with tunable Prandtl number. The ES–BGK model (Holway 1966; Andries et al. 2000) is an elegant but simple improvement over the BGK model. This model assumes that the distribution function relaxes to an anisotropic Gaussian distribution within mean free time $\tau$. The anisotropic Gaussian in itself evolves towards the MB distribution with a second time scale. The presence of a second time scale as free parameter ensures that the time scales related to momentum and thermal diffusivity are independent, 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 2008; Wu et al. 2015). 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 1951; Wu et al. 2015; Chapman & Cowling 1970). However, it is often not feasible to do any analysis on the Boltzmann-type equation for polyatomic gases. Therefore, several simplifications to model polyatomic gases have also been proposed. They essentially incorporate the rotational kinetic energy by decomposing the two-particle distribution function into two independent single-particle distribution functions (Morse 1964; Andries et al. 2000; Kataoka & Tsutahara 2004; Watari 2007; Nie, Shan & Chen 2008; Tsutahara et al. 2008; Larina & Rykov 2010; Wu et al. 2015; Wang et al. 2017; Bernard, Iollo & Puppo 2019). Furthermore, a thermodynamic framework and extensions thereof were developed for modelling highly nonequilibrium phenomena in dense and rarefied polyatomic gases where the Navier–Stokes–Fourier theory is no longer valid (Arima et al. 2012; Müller & Ruggeri 2013; Ruggeri & Sugiyama 2015). A few BGK-like models have also been proposed for polyatomic gases which accept the Prandtl number as a tunable parameter (Andries et al. 2000; Brull & Schneider 2009).

Hydrodynamic simulations for a realistic system require the development of reduced-order models to account for rotational degrees of freedom ideally without increasing the phase-space dimensionality. Indeed, the standard approach is to demonstrate that the two-particle distribution function describing the translational and rotational 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. 2000). 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 2020a). This model supplemented the standard Boltzmann BGK equation with an advection–relaxation equation for the evolution of rotational energy. It preserved the correct conservation laws for diatomic gases in the hydrodynamic limit and satisfied the $H$ theorem. However, the model was restricted to diatomic gases and a Prandtl number $7/5$, limiting its application for heat transfer problems.

We propose a kinetic model of polyatomic gases to tune Prandtl number, specific heat ratio and bulk viscosity in a physically transparent fashion. To do so, we write a new collision kernel which is a linear combination of the ES–BGK and BGK kernels that are locally relaxing to different temperatures at different 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 2001; Ansumali et al. 2007b; Atif et al. 2017; Atif, Kolluru & Ansumali 2022). 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 2004; Nie et al. 2008; Chen et al. 2010; Huang, Lan & Li 2020). 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.2a–c)\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 1970; Cercignani 1988; McQuarrie 2000). The nonlinear integro-differential Boltzmann collision kernel is often replaced by simpler models that should recover the following essential features such as the conservation laws, the $H$ theorem, and that the MB distribution is the zero of collision (Cercignani 1988).

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. 1954) 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 1970):

(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 1965) also describes the collision as simple relaxation process but, unlike the BGK model, it overcomes the restriction on the Prandtl number. The extra ingredient for 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. 2000). 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 1964; Rykov 1975; Kuščer 1989). 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 1975; Andries et al. 2000). 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 1978; Kuščer 1989). 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 2009)

(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. (2000) 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 1958). 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. (2020a) 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 1978). 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 2009):

   (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. 2000), 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 2001; Atif, Namburi & Ansumali 2018). In particular, we validate the proposed kinetic model using a 41-velocity crystallographic lattice from Kolluru et al. (2020b), 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 2016). 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.4a–c)\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 2003), 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 2019) 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 1987)

(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 2005)

(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 1993)

(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 1965). 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 2002) 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 1997):

(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 2007)

(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. 2015)

(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. 2021) 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.