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 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.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 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.