1. Introduction
The gas most commonly encountered in everyday life is air. It is a mixture of polyatomic gases – primarily nitrogen, which constitutes approximately 78 % of air, and oxygen at approximately 21 % – along with trace amounts of constituents such as argon, carbon dioxide, ozone and water vapour; most of these are polyatomic gases. Air also contains tiny particles called aerosols, such as fog, dust, smoke and particulate air pollutants, so-called respirable ultrafine particulates (UFPs). Owing to their nanoscale size, UFPs can penetrate deep into the lungs and blood stream and can have severe health implications, including heart attacks and respiratory diseases (De Falco et al. Reference De Falco2017). Aerosols can also be produced by coughing or sneezing, and airborne viruses can pass from person to person through these nanoscale droplets; therefore, scientists have paid much attention to the dynamics of aerosols that waft through the air.
Modelling of micro/nano gas flows is especially intriguing because classical gas dynamics models fail to capture the flow characteristics (such as drag force on aerosols, mass and heat flux on the evaporation meniscus) in rarefied conditions. As the representative physical length scale (
$L$) of the flow becomes comparable to the mean free path (
$\lambda$) in the gas, i.e. the Knudsen number Kn (
$=\lambda /L$) 
$\approx 1$, the Navier–Stokes–Fourier (NSF) equations fail to predict thermal and fluid flow fields to sufficient accuracy. The Boltzmann equation (Kremer Reference Kremer2010) in such a situation offers a more accurate description of the micro/nano gas flows, but unfortunately is computationally expensive. As a work around over the last few years, an impressive body of research has been devoted to the development of extended hydrodynamic models for polyatomic gases (Zhdanov Reference Zhdanov1968; Mallinger Reference Mallinger1998; Kremer Reference Kremer2010; Cai & Li Reference Cai and Li2014; Rahimi & Struchtrup Reference Rahimi and Struchtrup2014), which may provide accurate and computationally efficient simulation-for-design capabilities down to nanoscales.
The notion of extended hydrodynamic models goes back to Grad's pioneering work on approximating the Boltzmann equation via moment equations (Grad Reference Grad1958). The 13 moment equations obtained via Grad's method yield unphysical sub-shocks and are also prone to losing their hyperbolicity in strong non-equilibrium conditions. Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003) and Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2004) developed a linearly stable system of regularized 13 moment (R13) equations for monoatomic gases using the order-of-magnitude method (Struchtrup Reference Struchtrup2005) by combining the ideas of the Chapman–Enskog (Chapman & Cowling Reference Chapman and Cowling1970) and Grad methods. The regularization process extends the validity of the Grad closures through the addition of higher-order derivatives and provides smooth shock structures. Rahimi & Struchtrup (Reference Rahimi and Struchtrup2014), using the same method, obtained a hierarchy of moment equations for polyatomic gases at different levels of accuracy. Cai & Li (Reference Cai and Li2014) proposed a numerical regularized moment method of arbitrary order (NRxx) for polyatomic gases with the ellipsoidal statistical Bhatnagar–Gross–Krook (BGK) collision operator.
Unfortunately, the extended hydrodynamic equations obtained by Grad's method, the Chapman–Enskog method or the order-of-magnitude method are not accompanied by a proper entropy balance with strictly non-negative production. This might lead to breakdown of solutions or unphysical results (Struchtrup & Nadler Reference Struchtrup and Nadler2020). Another approach towards approximation of the Boltzmann equation is the rational extended thermodynamics (RET) framework (Müller & Ruggeri Reference Müller and Ruggeri2013; Ruggeri & Sugiyama Reference Ruggeri and Sugiyama2015), which, however, has the same limitations as those for the moment equations (only the linearized equations have a proper entropy law). Another popular family of moment closure methods is the maximum-entropy hierarchy (Junk & Unterreiter Reference Junk and Unterreiter2002; Torrilhon Reference Torrilhon2016). The maximum-entropy closures ensure entropy inequality and hyperbolicity, but they suffer from high computational cost incurred to compute the closure. Furthermore, it is impossible to obtain a closed form solution from the higher-order maximum-entropy closures.
Recently, a phenomenological procedure for monoatomic gases was proposed (Rana, Gupta & Struchtrup Reference Rana, Gupta and Struchtrup2018) in which the entropy flux contains additional nonlinear contributions in stress and heat flux – an ansatz consistent with the RET. The resulting coupled constitutive relations (CCR) are accompanied by an entropy inequality, in which the entropy remains the equilibrium entropy as integrated from the equilibrium Gibbs equation, but entropy flux and entropy generation exhibit higher-order correction terms. These additional terms in the entropy flux add a few second-order correction terms to the NSF system. The CCR system allows us to formulate thermodynamically consistent boundary conditions and their efficacy in capturing rarefaction effects at small scales was shown by Rana et al. (Reference Rana, Gupta and Struchtrup2018, Reference Rana, Saini, Chakraborty, Lockerby and Sprittles2021b). In this work, we extend the theory of the CCR to polyatomic gases. The resulting model (polyatomic-CCR) is benchmarked with experimental data for the slow rarefied gas flow over different shapes; the model's prediction for the drag force is in excellent agreement with the experiments. The normal shocks structure in nitrogen is also computed for the Mach number 
${Ma} \lesssim 4$; again, there is a good match between the polyatomic-CCR model and the experiments. In large measure, the polyatomic-CCR model enables the study of small-scale physics that is not accessible by the conventional fluid dynamics at very little to no added computational cost.
The remainder of this paper is organized as follows. In § 2, we introduce the governing equations and classical irreversible thermodynamics (CIT) framework. In § 2.3, derivation of the polyatomic-CCR model is presented. In § 3, we introduce linearized and dimensionless equations followed by their solutions for slow rarefied gas flow over a spherical and doublet particle in §§ 4 and 6, respectively. A brief derivation of the Green's functions is also given in § 5. In § 6, the normal stock structure in nitrogen gas is computed. We conclude and discuss future directions in § 7.
2. The governing equations
2.1. Conservation laws for polyatomic gases
The conservation laws are the balance equations for the mass density 
$\rho$, the fluid velocity 
$\boldsymbol {v}$ and internal energy 
$u$. For a compressible, polyatomic gas these are formulated as (Kremer Reference Kremer2010)
$$\begin{gather} \frac{\partial \rho }{\partial t}+\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \left( \rho \boldsymbol{v}\right) =0, \end{gather}$$
$$\begin{gather}\frac{\partial \left( \rho \boldsymbol{v}\right) }{\partial t}+\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \left[ \rho \boldsymbol{v}\otimes \boldsymbol{v}+ \left( p+\varPi \right) \boldsymbol{\mathsf{I}}+\boldsymbol{\sigma} \right] =0, \end{gather}$$
$$\begin{gather}\frac{\partial \left( \rho u+\frac{1}{2}\rho v^{2}\right) }{\partial t} +\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \left[ \left( \rho u+\frac{1}{2}\rho v^{2}+p+\varPi \right) \boldsymbol{v}+\boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{v}+\boldsymbol{q}\right] =0, \end{gather}$$
where 
$\boldsymbol{\mathsf{I}}$ is the identity tensor, 
$p$ is the pressure, 
$\varPi$ is the dynamic pressure, 
$\boldsymbol {\sigma }$ is the symmetric and trace-free part of the pressure tensor and 
$\boldsymbol {q}$ is the heat flux. The internal energy 
$u$, of a polyatomic ideal gas is related to the thermodynamic temperature 
$T$ through the caloric equation 
${\rm d}u=(({3+\delta })/{2})\,{\rm d}\theta :={{\rm d}\theta }/({\gamma -1})$, i.e. the constant volume specific heat is given by 
$c_v= (({3+\delta })/{2}){R}$. Here, 
${R}$ is the ideal gas constant and 
$\gamma$ 
$(=c_p/c_v)$ is the specific heat ratio.
For convenience, we also define the temperature and the dynamic temperature in energy units as 
$\theta := {R}T{}$ and 
$\vartheta :=\varPi /\rho$. A value 
$\delta$ can be chosen such that the proper caloric equation of gas is recovered, i.e.
\begin{equation} \delta = \frac{5-3 \gamma}{\gamma-1}. \end{equation}
For example, for the nitrogen gas or for the dry air 
$\gamma \approx 1.4$, which gives 
$\delta =2$.
In molecular gases, 
$c_v$ depends on the temperature, which implies 
$\delta$ should also be a function of temperature. However, for processes where changes in the temperature are small, one can assume 
$c_v$ to be constant, hence, 
$\delta$ can be assumed constant. Throughout this article we consider such processes, and assume 
$\delta$ to be constant. Furthermore, the pressure 
$p$ in (2.1) is assumed to given by the ideal gas law, 
$p=\rho \theta$.
The conservation laws (2.1) need to be closed by constitutive relations between the field variables (
$\rho$, 
$\boldsymbol {v}$ and 
$\theta$) and their fluxes (
$\varPi$, 
$\boldsymbol {\sigma }$ and 
$\boldsymbol {q}$). These constitutive relations can be obtained by enforcing the second law inequality, as shown below.
2.2. The NSF constitutive relations from the CIT framework
The CIT framework is based on the hypothesis of local equilibrium, according to which, the Gibbs relation is assumed to be locally valid, i.e.
\begin{equation} \theta \,{\rm d}\eta ={\rm d}u-\frac{p}{\rho ^{2}}{\rm d}\rho , \end{equation}
where 
$\eta =s/{R}$ denotes the dimensionless entropy (
$s$ being the dimensional entropy). Multiplying the last equation with 
$\rho$, and replacing 
${\rm D}\rho /{\rm D}t$ and 
${\rm D}u/{\rm D}t$ by means of the conservation laws (2.1), one obtains the entropy balance equation as
\begin{equation} \frac{\partial \left( \rho \eta \right)} {\partial t}+\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \left[ \rho \eta \boldsymbol{v}+\frac{\boldsymbol{q}}{\theta }\right] ={-}\frac{1}{\theta } \varPi \boldsymbol{\bigtriangledown}\boldsymbol{\cdot} \boldsymbol{v}-\frac{1}{\theta }\boldsymbol{\sigma} \boldsymbol{:} \left\langle \boldsymbol{\bigtriangledown} \boldsymbol{v}\right\rangle -\frac{\boldsymbol{q}}{ \theta ^{2}}\boldsymbol{\cdot} \boldsymbol{\bigtriangledown} \theta . \end{equation}
Here, the angular parentheses around a matrix denote the symmetric and trace-free part, i.e. 
$\langle A\rangle =(A+A^{T})/2 - \boldsymbol{\mathsf{I}} \,{\rm Tr}(A)/3$; 
$\boldsymbol{\mathsf{I}}$ being the identity matrix and 
${\rm Tr}(A)$ the trace of 
$A$.
Comparing (2.4) with a general balance equation of the form
\begin{equation} \partial \left( \rho \eta \right) /\partial t+\boldsymbol{\bigtriangledown}\boldsymbol{\cdot} \left( \rho \eta \boldsymbol{v}+\boldsymbol{\varPsi }\right) =\varSigma , \end{equation}
it follows that the non-convective entropy flux 
$\boldsymbol {\varPsi }=\boldsymbol {q}/\theta$ and the entropy production
\begin{equation} \varSigma ={-}\frac{1}{\theta }\varPi \underline{\boldsymbol{\bigtriangledown}\boldsymbol{\cdot} \boldsymbol{v}}- \frac{1}{\theta }\boldsymbol{\sigma}\boldsymbol{:}\underline{\left\langle \boldsymbol{\bigtriangledown} \boldsymbol{v}\right\rangle }-\frac{\boldsymbol{q}}{\theta ^{2}}\boldsymbol{\cdot} \underline{\boldsymbol{\bigtriangledown} \theta }. \end{equation}
The expression for 
$\varSigma$ forms a bi-linear relationship in the thermodynamic fluxes (
$\varPi$, 
$\boldsymbol {\sigma }$, 
$\boldsymbol {q}$) and the underlined thermodynamic forces. To guarantee the positiveness of the entropy production 
$\varSigma$, it is sufficient to assume linear flux-force relations of the form
\begin{equation} \boldsymbol{\sigma } ={-}2\mu \left\langle \boldsymbol{\bigtriangledown} \boldsymbol{v}\right\rangle {, }\quad \varPi ={-}\mu _{b}\boldsymbol{\bigtriangledown}\boldsymbol{\cdot}\boldsymbol{v}\quad \text{and}\quad \boldsymbol{q}={-}\kappa \boldsymbol{\bigtriangledown} \theta , \end{equation}
where 
$\mu$, 
$\mu _{b}$ and 
$\kappa / {R}$ are positive phenomenological coefficients, identified as the viscosity, the bulk viscosity and the thermal conductivity of the gas, respectively.
2.3. Coupled constitutive relations from the Gibbs relation
To obtain second-order CCR, the entropy flux 
$\varPsi$ entails all the quadratic vector terms in 
$\boldsymbol {q}$, 
$\boldsymbol { \sigma }$ and 
$\varPi$, i.e.
\begin{equation} \boldsymbol{\varPsi }=\frac{\boldsymbol{q}}{\theta }-\alpha _{0} \frac{\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{q}}{p\theta }-\beta _{0}\frac{\varPi \boldsymbol{q}}{ p\theta }, \end{equation}
where 
$\alpha _{0}$, 
$\beta _{0}$ are phenomenological dimensionless coefficients. The ansatz (2.8) is motivated from the entropy flux of 14 moments theory (Pavić, Ruggeri & Simić Reference Pavić, Ruggeri and Simić2013), in which 
$\alpha _0=\beta _0 = 2/(5+\delta )$.
The internal degrees of freedom of such polyatomic matter have been addressed in the literature as additional energy variables, see for example, Pavić et al. (Reference Pavić, Ruggeri and Simić2013). Literature regarding such approaches appeared as early as 1967 (Kogan Reference Kogan1967) and has been used effectively thereafter in the macroscopic modelling of polyatomic gases (Liu & Müller Reference Liu and Müller1983; Pavić & Simić Reference Pavić and Simić2014; Takashi Arima & Sugiyama Reference Takashi Arima and Sugiyama2018), producing the correct caloric equation of state for polyatomic gases. Of course the exact values for the phenomenological coefficients (
$\alpha _0$, 
$\beta _0$) appearing in (2.9) depend on the intricate nature of the energy transfer mechanism among internal and transnational degrees of freedom. Nevertheless, these can also be obtained through asymptotic reduction of the kinetic equation, which we shall briefly discuss after establishing the second law inequality.
Substitution of the extended entropy flux 
$\boldsymbol {\varPsi }$ from the last equation in (2.5), and balancing the terms give an extended second law
\begin{equation} \frac{\partial \left( \rho \eta \right) }{\partial t}+\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \left[ \rho \eta \boldsymbol{v}+\frac{\boldsymbol{q}}{\theta}-\frac{\alpha _{0}}{ p\theta }\boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{q}-\frac{\beta_{0}}{p\theta }\varPi \boldsymbol{q} \right] =\varSigma , \end{equation}where the entropy production term now reads
\begin{align} \varSigma &={-}\frac{\mathbf{\sigma }}{\theta }:\left\langle \boldsymbol{\bigtriangledown} \boldsymbol{v}+\frac{\alpha _{0}}{p}\left\{\boldsymbol{\bigtriangledown} \boldsymbol{q}-\frac{\boldsymbol{q} \boldsymbol{\bigtriangledown} \left(p\theta \right)}{2 p \theta} \right\}\right\rangle -\frac{\varPi }{\theta }\left[ \boldsymbol{\bigtriangledown}\boldsymbol{\cdot} \boldsymbol{v}+\frac{\beta _{0}}{p}\left\{\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \boldsymbol{q}-\frac{\boldsymbol{{q}} \boldsymbol{\cdot}\boldsymbol{\bigtriangledown} \left(p\theta \right)}{2 p \theta} \right\}\right] \nonumber\\ &\quad -\frac{{\boldsymbol{q}}}{\theta ^{2}}\boldsymbol{\cdot} \left[ \boldsymbol{\bigtriangledown} \theta +\frac{ \alpha _{0}}{\rho }\left\{ \boldsymbol{\bigtriangledown}\boldsymbol{\cdot} \boldsymbol{\sigma }-\frac{\boldsymbol{ \sigma }\boldsymbol{\cdot} \boldsymbol{\bigtriangledown} \left(p\theta \right)}{2 p \theta}\right\} +\frac{\beta _{0}}{\rho }\left\{ \boldsymbol{\bigtriangledown} \varPi -\frac{\varPi \boldsymbol{\bigtriangledown} \left(p \theta\right)}{2 p \theta}\right\} \right] . \end{align}Again, the entropy production forms a bilinear form as a sum of products of the thermodynamic fluxes and generalized thermodynamic forces. Once again, in order to guarantee the positivity of the entropy production, a linear flux-force relationship is assumed, which leads to the following CCR:
\begin{equation} \left.\begin{gathered} \boldsymbol{\sigma } ={-}2\mu \left\langle \underline{\boldsymbol{\bigtriangledown} \boldsymbol{v}}+\frac{\alpha _{0}}{p}\left\{\boldsymbol{\bigtriangledown} \boldsymbol{q}-\frac{\boldsymbol{q} \boldsymbol{\bigtriangledown} \left(p\theta \right)}{2 p \theta} \right\}\right\rangle {, }\\ \varPi ={-}\mu _{b}\left[ \underline{\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \boldsymbol{v}}+\frac{\beta _{0}}{p}\left\{\boldsymbol{\bigtriangledown} \boldsymbol{\cdot} \boldsymbol{q}-\frac{\boldsymbol{q} \boldsymbol{\cdot}\boldsymbol{\bigtriangledown} \left(p\theta \right)}{2 p \theta} \right\}\right] \end{gathered}\right\}, \end{equation}and
\begin{equation} \boldsymbol{q} ={-}\kappa \left[ \underline{\boldsymbol{\bigtriangledown} \theta }+\frac{ \alpha _{0}}{\rho }\left\{ \boldsymbol{\bigtriangledown}\boldsymbol{\cdot} \boldsymbol{\sigma } -\frac{\boldsymbol{\sigma }\boldsymbol{\cdot}\boldsymbol{\bigtriangledown} \left(p\theta \right)}{2p\theta}\right\} +\frac{\beta _{0}}{\rho }\left\{ \boldsymbol{\bigtriangledown} \varPi -\varPi \frac{\boldsymbol{\bigtriangledown} \left(p \theta \right)}{2p\theta}\right\} \right] . \end{equation}The stress (2.11) given by the CCR theory entails the hydrodynamic stress (underlined term) and the non-Newtonian second-order corrections. Similarly, the heat flux equation (2.12) contains Fourier's law (underlined term), along with the second-order corrections to the heat flux.
In order to study non-equilibrium effects in rarefied polyatomic gases, Rahimi & Struchtrup (Reference Rahimi and Struchtrup2014) introduced a set of 36 moment equations and further utilized the order of magnitude method to identify the leading-order terms of all the moments appearing in the 36 moment theory. In their work the balance equations for the dynamic temperature 
$\vartheta =\varPi /\rho$ reads (35(d) in Rahimi & Struchtrup Reference Rahimi and Struchtrup2014)
\begin{align} &\underline{\rho \frac{{\rm D}\vartheta }{{\rm D}t}+\frac{\delta }{3}\frac{2}{3+\delta }( \varPi \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}+\sigma :\boldsymbol{\nabla} \boldsymbol{v})+\frac{\delta }{3}\frac{2}{5+\delta }\boldsymbol{\nabla}\boldsymbol{\cdot} \varGamma }+ \frac{\delta }{3}\frac{2}{3+\delta }\left( p\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{v}+\frac{2}{ 5+\delta }\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{q}\right)\nonumber\\ &\quad ={-}\frac{1}{\tau _{int}}\varPi, \end{align}
where 
$\textrm {D}/\textrm {D}t$ denotes the material time derivative, 
$\tau _{int}$ is the relaxation time for the dynamic temperature and one defines 
$\varGamma = \boldsymbol {q}^{tr}-({5}/{\delta })\boldsymbol {q}^{in}$; 
$\boldsymbol {q}^{tr}$ and 
$\boldsymbol {q}^{in}$ being the translational and the internal heat fluxes, respectively. Furthermore, using the Chapman–Enskog expansion (in 
$\tau _{in}$), is was shown in the same article that the underlined terms in last equation are higher order compared with the others, which can be assumed small in first-order approximations. Neglecting these higher-order terms, one gets the leading-order constitutive relation for 
$\varPi$, as
\begin{equation} \varPi ^{(1)}={-}\tau _{int}\frac{\delta }{3}\frac{2}{3+\delta }\left( p\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}+\frac{2}{5+\delta }\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}\right) . \end{equation}
Comparing the last equation with (2.11
$_2$), one immediately recognizes
\begin{equation} \tau _{int}=\frac{3}{\delta }\frac{3+\delta }{2}\frac{\mu _{b}}{p}\quad\text{and}\quad\beta _{0}=\frac{2}{5+\delta }. \end{equation}
Similarly, the Chapman–Enskog expansion – via introducing a smallness parameter 
$\varepsilon$, and assuming 
$\epsilon \sim O( \varepsilon ^{1/2})$, 
$\tau _{tr}\sim O( \varepsilon ^{\alpha /2})$, 
$\tau _{in}\sim O( \varepsilon ^{1/2})$ and 
$0<\alpha <1$ – of the balance equation for heat flux follows (from 35( f) in Rahimi & Struchtrup Reference Rahimi and Struchtrup2014), one again finds 
$\beta _{0}={2}/({5+\delta })$. Similar, arguments based upon the Chapman–Enskog expansion can be made to justify 
$\alpha _{0}={2 }/({5+\delta })$. Nevertheless, here, we shall point out that the values for these phenomenological coefficients may vary with different collision models; consideration of which is beyond the scope of this article and throughout this article we shall assume
\begin{equation} \alpha _{0}=\frac{2}{5+\delta }{,\quad }\beta _{0}=\frac{2}{5+\delta }{. } \end{equation}Henceforth, we shall refer to (2.11)–(2.12) as the polyatomic-CCR model.
3. Linearized and dimensionless equations
In this section, the equations we consider are dimensionless and linearized by introducing small perturbations from their values in a reference rest state defined by a constant pressure 
$p_{0}$ and a constant temperature 
$\theta _{0}$. The relations between the field variables and their dimensionless deviations (denoted with hat symbols) from their reference rest state are given as
\begin{equation}
\left.\begin{gathered} p
=p_{0}(1+\hat{p}) {, }\quad\theta =\theta _{0}(1+\hat{ \theta}) {, }\quad\boldsymbol{v}=\sqrt{\theta
_{0}}\hat{\boldsymbol{v}}{, } \quad\boldsymbol{\sigma }
=p_{0}\hat{\boldsymbol{\sigma}}{, } \quad\varPi =
p_{0}\hat{\varPi} {, }\\ \boldsymbol{q}=p_{0}\sqrt{\theta
_{0}}\hat{\boldsymbol{v}}\quad \text{and}\quad\boldsymbol{x
}=L\hat{\boldsymbol{x}}, \end{gathered}\right\}
\end{equation}
where 
$L$ is a characteristic length scale. The linearized and stationary (
$\partial /\partial t\equiv 0$) conservation laws are given by
\begin{equation}
\boldsymbol{\bigtriangledown} \boldsymbol{\cdot}\, \hat{\boldsymbol{v}} =0,\quad \boldsymbol{\bigtriangledown}
\boldsymbol{\cdot} [ (
\hat{p}+\hat{\varPi}) \boldsymbol{\mathsf{I}} + \hat{\boldsymbol{\sigma}}]
=0\quad\text{and}\quad
\boldsymbol{\bigtriangledown}\boldsymbol{\cdot}
\hat{\boldsymbol{q}}=0, \end{equation}
and the linearized stress, dynamic pressure and heat flux are specified as
\begin{align} \hat{\boldsymbol{\sigma}} ={-}2\,{Kn}\left\langle \boldsymbol{\bigtriangledown} \hat{\boldsymbol{v}}+\frac{2}{\delta +5}\boldsymbol{\bigtriangledown} \hat{\boldsymbol{q}}\right\rangle ,\quad\hat{\varPi}=0\quad\text{and}\quad \boldsymbol{q} ={-}\frac{{Kn}\,c_{p}}{{Pr}}\left[ \boldsymbol{\bigtriangledown} \hat{\theta}+\frac{2}{\delta +5}\boldsymbol{\bigtriangledown}\boldsymbol{\cdot} \hat{\boldsymbol{\sigma}} \right] . \end{align}
In (3.3a–c), the Knudsen number appears as the scaled viscosity 
${Kn}=\mu _{0}\sqrt {\theta _{0}}/(p_{0}L)$, and 
${Pr}= c_{p}\mu _{0}/\kappa$ is the Prandtl number. Here, it is assumed that all the field variables are measured relative to their equilibrium values at infinity, hence all perturbations must vanish as 
$|\hat {\boldsymbol {x}}|\rightarrow \infty$.
3.1. Mathematical similarities between the CCR and NSF equations
It is worth noting that the steady, linearized CCR equations (3.2a–c)–(3.3a–c) can be converted to the Stokes–Fourier-like system by introducing the following auxiliary variables:
\begin{equation} \hat{\boldsymbol{V}} :=\hat{\boldsymbol{v}}+\frac{2}{\delta +5}\boldsymbol{\bigtriangledown}\hat{\boldsymbol{q}}\quad\text{and}\quad \hat{\varTheta} :=\hat{\theta}-\frac{2}{\delta +5}\hat{p}. \end{equation}As a result, the linearized conservation laws (3.2a–c) and (3.3a–c) turn into
\begin{equation} \boldsymbol{\bigtriangledown}\boldsymbol{\cdot}\,\hat{\boldsymbol{V}}=0, \quad \boldsymbol{\bigtriangledown} \hat{p}= {Kn}\,\Delta \hat{\boldsymbol{V}}\quad\text{and}\quad\Delta \hat{\varTheta}=0, \end{equation}
which are mathematically the same as the Stokes equations in 
$\{\hat {\boldsymbol {V}}, \hat {p}\}$ and the Laplace equation in 
$\hat {\varTheta }$; however, now the Stokes equations are coupled with the Laplace equation via auxiliary variables 
$\hat {\boldsymbol {V}}$ and 
$\hat {\varTheta }$. Nevertheless, from a mathematical point of view, the classical analytic solutions and techniques from potential theory (Batchelor Reference Batchelor2000) for the Stokes and Laplace equations can be carried over to the slow, steady-state polyatomic-CCR model, as we show in next section.
In order to assess the accuracy of the newly derived CCR model, we shall apply these equations to two classical boundary value problems: (i) slow rarefied gas flow over a spherical particle, for which a closed form solution can be obtained, and (ii) slow rarefied gas flow over a doublet, which is solved numerically using the method of fundamental solution (MFS) (Lockerby & Collyer Reference Lockerby and Collyer2016; Rana et al. Reference Rana, Saini, Chakraborty, Lockerby and Sprittles2021b).
4. Slow rarefied gas flow over a spherical particle
First, we consider the case of a slow rarefied gas flow around a spherical particle of radius 
$L$. For this problem, solutions of the NSF equations are available (Batchelor Reference Batchelor2000), giving the pressure 
$\hat {p}$, the velocity in radial direction 
$\hat {u}_{r}$ and the azimuthal directions 
$\hat {u}_{\phi }$ as
\begin{equation} \left.\begin{gathered} \hat{p}={Kn}\frac{a_{1}}{\hat{r}^{2}}\cos \phi{, } \quad \hat{u}_{r} =\hat{U}_{\infty }\left[ 1+\frac{a_{1}}{\hat{r}}+\frac{a_{2}}{ \hat{r}^{3}}\right] \cos \phi\quad \text{and}\\ \hat{u}_{\varphi } ={-}\hat{U}_{\infty }\left[ 1+\frac{a_{1}}{2\hat{r}}- \frac{a_{2}}{2\hat{r}^{3}}\right] \sin \phi, \end{gathered}\right\} \end{equation}
in the spherical coordinate system. Here, 
$\hat {U} _{\infty }$(
$=U_{\infty }/\sqrt {\theta _{0}}$) is the far-field velocity of the gas and 
$\hat {r}$ (
$=r/L$) is the dimensionless radial distance, 
$\hat {r}$. The temperature and the heat flux are obtained as
\begin{equation} \left.\begin{gathered} \hat{\theta}= \frac{c_{1}}{\hat{r}^{2}}\cos \phi {, } \quad \hat{q}_{r} = \frac{2 c_p\, {Kn}}{{Pr}}\frac{c_1-\alpha_0 a_1\,{Kn}}{\hat{r}^3}\cos \phi\quad \text{and}\\ \hat{q}_{\phi} = \frac{c_p \,{Kn}}{{Pr}}\frac{c_1-\alpha_0 a_1 \,{Kn}}{\hat{r}^3}\sin \phi . \end{gathered}\right\} \end{equation}
The integration constants (
$a_{1}$, 
$a_{2}$ and 
$c_{1}$) are obtained by applying appropriate boundary conditions. The Stokes formula for the drag force exerted on a spherical particle due to the gas flow reads 
$\hat {F}_{S}=-6{\rm \pi} \,{Kn}\hat {U}_{\infty }$, which is obtained by applying the classical no-slip and no-jump boundary conditions at 
$\hat {r}=1$, i.e. 
$\hat {u}_{r}=\hat {u}_{\phi }=0$ and 
$\hat { \theta }=0$ to the NSF solutions, i.e. (4.1a–c)–(4.2a–c) with 
$\alpha _0=0$. The Stokes’ drag formula is valid for 
${Kn}\rightarrow 0$. Likewise, the NSF solutions with the standard first-order slip and temperature jump boundary conditions, i.e.
\begin{equation} \hat{u}_{r}=0,\quad \hat{u} _{\phi }={-}\sqrt{\frac{\rm \pi}{2} }\hat{\sigma}_{r\phi}\quad\text{and}\quad\hat{\theta}={-}\frac{2}{(4+\delta)}\sqrt{\frac{\rm \pi}{2} }\hat{q}_{r}, \end{equation}
yield drag force 
$\hat {F}_{D} = \hat {F}_{S} (2{Kn}+\sqrt {2/{\rm \pi} })/(3 \,{Kn}+\sqrt {2/{\rm \pi} })$. In (4.3a–c), the accommodation coefficients are taken to be unity (Sharipov Reference Sharipov2011).
The second-order boundary conditions which shall be used for the polyatomic-CCR solutions (4.1a–c)–(4.2a–c), read
\begin{equation} \hat{u}_{r}=0{, }\quad \hat{u} _{\phi }={-}\sqrt{\frac{\rm \pi}{2} }\hat{\sigma}_{r\phi }-\underline{\frac{3}{4}}\hat{q} _{\phi}\quad\text{and}\quad\hat{\theta}={-}\underline{\frac{2}{(4+\delta)}\sqrt{\frac{\rm \pi}{2}} }\hat{q}_{r}-\frac{3}{4}\hat{\sigma}_{rr}.\end{equation}
Here, the first boundary condition (4.4a–c
$_1$) represents the no-penetration condition on a solid wall. The slip boundary condition (4.4a–c
$_2$) defines the velocity slip at the wall, which relates the tangential gas velocity slip to the tangential shear stress and heat flux. The slip boundary condition, when compared with the second-order boundary conditions from the literature (Deissler Reference Deissler1964; Lockerby et al. Reference Lockerby, Reese, Emerson and Barber2004; Zhang, Meng & Wei Reference Zhang, Meng and Wei2012) gives the first-order and second-order slip coefficients as 
$A_1=1$ and 
$A_2 =({3}/{2{\rm \pi} })({1}/{Pr}) =0.682$; see Appendix A for details. The third boundary condition (4.4a–c
$_3$) is the temperature jump boundary condition, where the temperature jump coefficient 
$\zeta _{T}=({\sqrt {{\rm \pi} }}/{2})({2}/{(4+\delta )})({c_{p}}/{Pr})$ (Sharipov Reference Sharipov2011) and the accommodation coefficient is equal to one. An additional term in (4.4a–c
$_3$) appears as a consequence of the Onsager symmetry (Rana et al. Reference Rana, Gupta and Struchtrup2018, Reference Rana, Gupta, Sprittles and Torrilhon2021a).
It should be noted here that values for the slip and jump coefficients vary in the literature (Sharipov Reference Sharipov2011) and one may perform asymptotic expansions – taking into account the so-called Knudsen layer corrections (Rana & Struchtrup Reference Rana and Struchtrup2016) – to fit theory/experimental results. Indeed, in (4.4a–c
$_2$) the coefficient 
$3/4$ (which is equivalent to taking an appropriate value for the second-order slip coefficient 
$A_2$) was chosen so that it gives a good fit for the drag coefficient for a sphere; see Appendix A for further details along with the sensitivity analysis of the results on slip/jump coefficients. In particular, the considered choice of 
$A_2$ significantly overestimates the value of the thermal slip coefficient (Sharipov Reference Sharipov2011). Nevertheless, for the problems considered in this paper, the effects of the thermal slip coefficient are negligible.
Substitution of the CCR solutions (4.1a–c)–(4.2a–c) in (4.4a–c) gives integration constants (
$a_{1}$, 
$a_{2}$ and 
$c_{1}$), which in turn provide a closed form expression for the drag as
\begin{equation} \hat{F}_{D} = \hat{F}_{S} \left(1+kn\frac{1+\left(\dfrac{3kn}{2{\rm \pi} }\dfrac{1+4\dfrac{3\delta +19}{\delta +5}kn}{Pr +2 \dfrac{\delta +5}{\delta +4}kn}\right)}{1+kn\left( 2+\dfrac{9kn}{4{\rm \pi} }\dfrac{3-\delta }{\Pr + 2\dfrac{\delta +5}{\delta +4}kn}\right) }\right)^{{-}1}, \end{equation}
where, for brevity, we have introduced a modified Knudsen number 
$kn = \sqrt {{{\rm \pi} }/{2}}\,{Kn}$. Clearly, 
$\hat {F}_{D} \rightarrow \hat {F}_{S}$ as 
${Kn}\rightarrow 0$ and 
$\hat {F}_{D}/\hat {F}_{S}\rightarrow 0$ as 
${Kn}\rightarrow \infty$, hence the continuum and the free-molecular regimes are recovered, qualitatively.
In figure 1(a) the drag force 
$\hat {F}_D$ normalized with Stokes’ drag 
$\hat {F}_S$ is shown as a function of the Knudsen number. The results obtained from the polyatomic-CCR model (solid blue), first-order NSF (dashed red) and Stokes’ formula (thin green) are compared with the experimental results (black dot-dashed lines and symbols) given by Allen & Raabe (Reference Allen and Raabe1982).
Figure 1. (a) The normalized drag force (with Stokes’ drag) on a sphere against Knudsen number. The results computed from the NSF with first-order slip and the polyatomic-CCR model are compared with the experimental data by Allen & Raabe (Reference Allen and Raabe1982). (b) The drag force on a doublet against Knudsen number. Two different flow conditions are assumed, i.e. flow parallel to the line of centres (
$\varphi = 0$) and flow perpendicular to the line of centres (
$\varphi = {\rm \pi}/2$). The experimental data are taken from Cheng et al. (Reference Cheng, Allen, Gallegos, Yeh and Peterson1988).
For larger values of the Knudsen number, the force decreases. Notably, for the NSF with the first-order boundary conditions the drag force reaches a finite value (
$2/3$) as 
${Kn} \rightarrow \infty$, which contradicts experimental observations. The polyatomic-CCR results match the experiment quite remarkably, even for large values of the Knudsen number. In the above computations, transport properties are taken to be for dry air, i.e. 
$Pr =0.7$ and 
$\delta =2$, which gives the ratio of the specific heats 
$\gamma = 1.4$ – a value consistent with air/nitrogen at moderate temperatures.
5. Slow rarefied gas flow over a doublet
In this section we shall consider the non-trivial case of flow over a doublet for which we shall apply the MFS technique. However, before that, we first formulate the Green functions associated with the polyatomic-CCR model (3.2a–c)–(3.3a–c).
5.1. Green's functions and the MFS
Green's functions for the Stokes equations and the Laplace equations are obtained as a response of the fluid with regard to a point body force 
$\hat {\boldsymbol {f}}$ in the momentum equation (stokeslet) and a heat source of strength 
$\hat {g}$ (thermlet) in the energy balance equation (Lockerby & Collyer Reference Lockerby and Collyer2016; Rana et al. Reference Rana, Saini, Chakraborty, Lockerby and Sprittles2021b). Green's functions associated with the polyatomic-CCR equations follow from the stokeslet and thermlet, as
$$\begin{gather} \hat{\boldsymbol{v}}{\left( \hat{\boldsymbol{r}}\right) } =\frac{1}{8{\rm \pi} }\left(\frac{\hat{\boldsymbol{r}}\hat{\boldsymbol{r}}}{|\hat{\boldsymbol{r}}|^{3}}+\frac{\boldsymbol{\mathsf{I}}}{|\hat{\boldsymbol{r}}|}\right) \boldsymbol{\cdot}\hat{\boldsymbol{f}}+\underline{\frac{ 3\alpha_{0}^{2}c_{p}\,{Kn}^{2}}{4{\rm \pi} \,Pr |\hat{\boldsymbol{r}}|^{5}}\langle\hat{\boldsymbol{r}}\hat{\boldsymbol{r}}\rangle\boldsymbol{\cdot} \hat{\boldsymbol{f}}}{, } \end{gather}$$
$$\begin{gather}\hat{p}{\left( \hat{\boldsymbol{r}}\right) } =\frac{{Kn}}{4{\rm \pi} }\frac{\boldsymbol{ \hat{f}{\cdot} \hat{r}}}{|\boldsymbol{r}|^{3}}\quad\text{and }\quad\hat{\boldsymbol{\sigma}}{\left( \hat{\boldsymbol{r}}\right) } =\frac{3\,{Kn}}{ 4{\rm \pi} |\hat{\boldsymbol{r}}|^{5}}(\hat{\boldsymbol{f}}\boldsymbol{\cdot}\hat{\boldsymbol{r}}+\underline{2\alpha _{0}\,{Kn} \,\hat{g}}) \langle\hat{\boldsymbol{r}}\hat{\boldsymbol{r}}\rangle, \end{gather}$$
where 
$\hat {\boldsymbol {r}}= \hat {\boldsymbol {x}}-\hat {\boldsymbol {x}}^{s}$ is the position of an arbitrary point 
$\hat {\boldsymbol {x}}$ from the computational domain relative to the singularity point located at 
$\hat {\boldsymbol {x}}^{s}$ (a point outside of the computational domain).
The temperature 
$\hat {T}$ and heat flux 
$\hat {\boldsymbol {q}}$ are obtained as
\begin{equation} \hat{T}{\left( \hat{\boldsymbol{r}}\right) } =\frac{Pr }{4{\rm \pi} c_{p}}\frac{\hat{g}}{| \hat{\boldsymbol{r}}|}\quad\text{and} \quad \hat{\boldsymbol{q}}{\left( \hat{\boldsymbol{r}}\right) } =\frac{{Kn}}{4{\rm \pi} } \frac{\hat{g}}{|\hat{\boldsymbol{r}}|^{3}}\hat{\boldsymbol{r}}-\underline{\frac{ 3\alpha _{0}c_{p}\,{Kn}^{2}}{4{\rm \pi} \,Pr |\hat{\boldsymbol{r}}|^{5}}\langle\hat{\boldsymbol{r}}\hat{\boldsymbol{r}}\rangle\boldsymbol{\cdot} \hat{\boldsymbol{f}}}. \end{equation}
In the polyatomic-CCR Green's functions (5.1)–(5.2a,b), the coupling between the stokeslet and thermlet occurs due to the underlined term. Obviously, for the classic NSF equations (
$\alpha _{0}=0$), such coupling is absent and one obtains the classical stokeslet and thermlet.
5.2. Drag on a doublet
In figure 1b), we show the drag force on a doublet using the MFS technique developed by Lockerby & Collyer (Reference Lockerby and Collyer2016), and Rana et al. (Reference Rana, Saini, Chakraborty, Lockerby and Sprittles2021b) for the case of monoatomic gases. The essential idea behind MFS is to decompose the solutions of the linearized partial differential equations as a superposition of the Green's functions, where the weights are determined such that the boundary conditions are satisfied at some boundary nodes (collocation points).
In particular, we study a doublet moving with an angle 
$\varphi$ between the flow direction and the line joining the centres of the doublet. This is an interesting case to consider because (i) it is a full three-dimensional problem and no closed form analytic solution exists and (ii) experimental results are available for this problem in the literature (Cheng et al. Reference Cheng, Allen, Gallegos, Yeh and Peterson1988) for 
$\varphi =0$ and 
$\varphi ={\rm \pi} /2$ (Cheng et al. Reference Cheng, Allen, Gallegos, Yeh and Peterson1988). The overall agreement of the polyatomic-CCR model's drag force with the experimental data is excellent in both cases (
$\varphi = 0,{\rm \pi} /2$) and, as before, the NSF theory over-predicts the drag.
So far, we have considered steady-state slow rarefied flow conditions for which the polyatomic-CCR model can be transformed into the NSF-like equations and analytic/semi-analytic approaches can be utilized. For the considered cases, the dynamic temperature identically vanishes and the polyatomic nature of the gas only appears through 
${Pr}$, 
$\mu$, 
$c_p$, 
$c_v$, 
$\alpha _0$ and 
$\beta _0$, with the last two values being 
$2/5$ and 
$0$ for monoatomic gases, respectively. For nonlinear (for example, shocks in hyper-sonic flows) or transient processes (for example, wave propagation and Rayleigh–Brillouin scattering in polyatomic gases), other polyatomic effects might become manifest.
6. Normal shock structure in nitrogen
Finally, in order to assess the effects of nonlinear terms, we shall now use the polyatomic-CCR model to compute the shock structure in nitrogen gas. The normal shock structure problem is modelled as a one-dimensional problem so that the field variables (velocity, stress tensor and heat flux vector, etc.) have only a non-zero 
$x$-component. In the frame of reference of the moving shock (in the 
$x$-direction), the gas upstream and downstream is assumed to be in equilibrium state with 
$(\rho, v_x, \theta ) = (\rho _0, U_0, \theta _0)$ and 
$(\rho, v_x, \theta ) = (\rho _\infty, v_\infty, \theta _\infty )$, respectively. The upstream and downstream conditions are linked to each other by the Rankine–Hugoniot relations (Anderson Reference Anderson2003)
\begin{align} \frac{\rho_{\infty}}{\rho_0} = \frac{(\gamma +1) {Ma}^2}{2+{Ma}^2\left(\gamma -1\right)},\quad \frac{\rho_{\infty}v_{\infty}}{\rho_0\sqrt{\theta_0}} = \sqrt{\gamma}\,{Ma} \quad {\rm and}\quad \frac{\rho_{\infty}\theta_{\infty}}{\rho_0\theta_0} = \frac{1-\gamma +2 \gamma \,{Ma}^2}{1+\gamma}, \end{align}
where 
${Ma} = U_0/\sqrt {\gamma \theta _0}$ is the Mach number defined in terms of the upstream condition.
The mean free path based upon the upstream conditions can be defined as 
$\lambda _0 = (8 \sqrt {2}\mu _0)/(5\rho _0 \sqrt {{\rm \pi} \theta _0})$, where 
$\mu _0$ is the viscosity in the upstream state. For the shock problem, the polyatomic-CCR model reduces to
$$\begin{gather} \tau _{xx} := \varPi+\sigma_{xx} ={-}\frac{4\mu }{3}\left( 1+\frac{\delta}{2\left (3+\delta\right)}Z \right) \left( \frac{{\rm d}v_{x}}{{\rm d}\kern0.7pt x}+\frac{\alpha _{0}}{p}\left\{\frac{{\rm d}q_{x}}{{\rm d}\kern0.7pt x}-\frac{q_x}{2p\theta}\frac{{\rm d}\left(p \theta\right)}{{\rm d}\kern0.7pt x}\right\} \right){, } \end{gather}$$
$$\begin{gather}q_{x} ={-}\frac{\mu c_{p}}{Pr }\left[ \frac{{\rm d}\theta }{{\rm d}\kern0.7pt x}+\frac{\alpha _{0} }{\rho }\left\{ \frac{{\rm d}\left(\varPi+\sigma_{xx}\right)}{{\rm d}\kern0.7pt x}-\frac{\varPi+\sigma_{xx}}{2 p\theta}\frac{{\rm d} \left(p\theta\right) }{{\rm d}\kern0.7pt x} \right\} \right]. \end{gather}$$
Here, for convenience, we have defined 
$\tau _{xx} = \varPi +\sigma _{xx}$ and introduced 
$Z$ as the dimensionless rotational collision number. The viscosity is obtained from Sutherland's law and 
$Z$ from Parker's formula (Parker Reference Parker1959) as
\begin{align} \mu =\mu _0 \left(\frac{T}{T_0}\right)^{3/2}\frac{T_0+T_S}{T+T_S}\quad \text{and}\quad Z:=\frac{3(3+\delta)\mu_b}{2\delta\mu } = \frac{Z_\infty}{1+\dfrac{{\rm \pi} ^{3/2}}{2} \sqrt{\dfrac{T^*}{T}}+\left(\dfrac{{\rm \pi} ^2}{4}+{\rm \pi} \right)\dfrac{T^*}{T}}, \end{align}
where 
$T_s = 107\ \mathrm {K}$ is Sutherland's constant, 
$Z_\infty = 18.2$ and 
$T^{*} = 91.5\ \mathrm {K}$. The reduced system of ordinary differential equations is numerically solved using a second-order central finite difference scheme with boundary conditions (6.1a–c).
Figure 2 shows the profiles of reduced density, 
$\hat {\rho }:=(\rho -\rho _0)/(\rho _\infty -\rho _0)$ and reduced temperature, 
$\hat {\theta }:=(\theta -\theta _0)/(\theta _\infty -\theta _0)$ vs 
$x/\lambda _0$ for three Mach numbers, 
${Ma} = 1.7$ (2a), 
${Ma} =3.2$ (2b) and 
${Ma} =3.8$ (2c). For a smaller Mach number (
${Ma} = 1.7$), the macroscopic theories and the direct simulation Monte Carlo (DSMC) give a good agreement. At higher Mach numbers (2b,c), the upstream region of the shock is resolved accurately by the polyatomic-CCR model whereas some deviations are observed in the downstream region. Interestingly, for the NSF theory we observe an opposite behaviour. For the temperature shock profiles, the DSMC predictions appear to be slight narrower than the macroscopic theories. The NSF theory slightly over-predicts the temperature while the polyatomic-CCR theory under-predicts it. The main goal of the present study was to derive the polyatomic-CCR model and demonstrate its usefulness for some practical applications; a thorough test of these equations for the shock wave flow problems is beyond the scope of this article. Nevertheless, for the considered range of the Mach number, the overall results from the polyatomic-CCR model are in good agreement with the experimental data.
Figure 2. Comparison of reduced density and temperature profiles obtained using the polyatomic-CCR models with the experimental measurements of Alsmeyer (Reference Alsmeyer1976) for Mach numbers: (a) 
${Ma} = 1.7$ (b) 
${Ma} = 3.2$ and (c) 
${Ma} = 3.8$.
This is noteworthy given that the entropy current density we assume has been derived in Pavić et al. (Reference Pavić, Ruggeri and Simić2013). It employs the exponential distribution function in kinetic theory and the entropy of this distribution is analytically obtained by integrating over internal degrees of freedom and a particular velocity. This procedure crucially depends on the concept of local temperature which has been used in the exponential distribution. Despite the difficulties in this concept, it turns out that local temperature may thus be fruitfully applied to rarefied flows of polyatomic gases as well.
7. Conclusion and future directions
The CIT framework applied to rarefied flows requires justification for the local equilibrium hypothesis which is taken to hold in the case of rarefied gases. There have been several objections to this assumption, especially regarding the concept of local temperature and entropy. However, under the RET assumption of entropy current density subject to a different constitutive law, CCR modelling of rarefied gas flows extends the applicability beyond the classical NSF equation whilst maintaining the same mathematical and thermodynamical structure. Good agreement of CCR has been previously obtained with experiments in the case of monoatomic gases. In this article, we have proposed a polyatomic-CCR model, and obtained its solutions for some classical flow problems, namely slow flow past an aerosol with different shapes. Further, these have been benchmarked with experimental data from the literature. The polyatomic-CCR theory provides an excellent match with experimental observations over the rage of all Knudsen numbers. We have also computed the normal shock structure using the polyatomic-CCR model. A very good agreement was observed with the kinetic theory, as well as with the experimental data.
We have formulated Green's functions and a MFS technique for the polyatomic-CCR model, which, in principle, allows for simulations of slow polyatomic gas flow over any arbitrary shapes and sizes. Thus, from a modelling perspective, a simultaneously computationally efficient and accurate (to nano-scales) method down was proposed. For a future direction, a machine learning algorithm trained using the data from the polyatomic-CCR framework is planned, which can provide a simulation-for-design tool to study the multi-scale process of dispersion and the dynamics of UFPs in air.
The proposed model is found to possess three important limitations which could be worked upon in the future. They involve thermodynamic issues related to local equilibrium, temperature dependent specific heats and asymptotic consistency of higher-order boundary conditions. For example, the last issue assumes importance when one deals with flows with thermal transpiration (e.g. Knudsen pump). We comment upon these below. The proposed model assumes the validity of local equilibrium. This assumption cannot be justified when the bulk-to-shear-viscosity ratio is very large; for example, carbon dioxide gas has a bulk-to-shear-viscosity ratio 
${\sim }2000$. Conceptually, in addition, the definition of local thermodynamic quantities has been questioned. Temperature and entropy have been emphasized the most in this regard. The proportionality of internal energy and temperature has been especially a crucial assumption which requires addressing. In a scenario like this, an interesting approach could be to use an extended Gibbs equation in order to get a refined NSF equations or the corresponding CCR (Ruggeri Reference Ruggeri2015).
The polyatomic-CCR theory developed here assumes constant specific heats (hence 
$\delta$ to be constant). In light of the results in the § 6, it is worthwhile discussing the validity of this assumption. For the range of Mach number considered in this article (
$1.7 \leqslant {Ma} \leqslant 4$) the maximum temperature variations – from (6.1a–c) – are 
$1.46 \leqslant {\theta _{\infty }}/{\theta _0} \leqslant 4.0s$, which are seemingly high. However, for nitrogen gas, the change in 
$\delta$ between 200 and 800 K is approximately 
${\lesssim }7.5\,\%$ (Borgnakke & Sonntag Reference Borgnakke and Sonntag2009), therefore the assumption of 
$\delta$ being constant was justifiable. Nevertheless, for large enough temperature variations, this assumption may not hold true and an extension of the CCR model for temperature dependent specific heats (i.e. assuming 
$\delta = \delta (\theta )$) will be sought.
Furthermore, the boundary conditions (4.4a–c) employed in §§ 4–5 were obtained in a way to fit the drag coefficient on the sphere (by choosing an appropriate value of the second-order slip coefficient). The asymptotic analysis (see also Appendix A) suggests that these boundary conditions may lead to erroneous predictions in situations when thermal transpiration/thermal stress play an important role. A detailed analysis of such flows along with a consistent set of boundary conditions for the polyatomic-CCR theory is also planned for the future.
Funding
This work has been financially supported by ‘SRG’ project SRG/2021/000790 funded by the Science & Engineering Research Board, India.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Asymptotic analysis of the boundary conditions
Here, we shall carry out a comparison of boundary conditions (4.4a–c) with those from the literature and discuss qualitative effects of different terms appearing in the boundary conditions. The (linearized and dimensionless) generalized slip boundary condition relates the tangential gas velocity slip 
$\hat {V}_{i}$ to the tangential shear stress and heat flux, as (Lockerby et al. Reference Lockerby, Reese, Emerson and Barber2004)
\begin{align} \hat{V}_{i}t_{i} &=A_{1}\sqrt{\frac{{\rm \pi} }{2}}{Kn}\left( \frac{\partial \hat{v} _{i}}{\partial \hat{x}_{j}}+\frac{\partial \hat{v}_{j}}{\partial \hat{x}_{i}} \right) n_{i}t_{j}-A_{2}\left( \sqrt{\frac{{\rm \pi} }{2}}{Kn}\right) ^{2}\frac{ \partial }{\partial \hat{x}_{k}}\left( \frac{\partial \hat{v}_{i}}{\partial \hat{x}_{k}}+\frac{\partial \hat{v}_{k}}{\partial \hat{x}_{i}}\right) t_{i}\nonumber\\ &\quad +A_{3}{Kn}\frac{\partial \hat{\theta}}{\partial \hat{x}_{i}}t_{i}-A_{4}\sqrt{ \frac{{\rm \pi} }{2}}{Kn}^{2}\frac{\partial ^{2}\hat{\theta}}{\partial \hat{x}_{i} \hat{x}_{j}}n_{i}t_{j}, \end{align}
where 
$n_{i}$ is unit normal pointing from the boundary into the gas and 
$t_{i}$ is a tangential direction. In the last equation, 
$A_{1}$ and 
$A_{2}$ are the first-order and second-order slip coefficients, respectively. Furthermore, 
$A_{3}$ is the thermal slip coefficient (Sharipov Reference Sharipov2011) (in Sharipov's notation 
$A_{3}=\sigma _{T}$), which gives rise to thermal creep, and the last term (with 
$A_{4}$) is due to the thermal-stress slip flow (Lockerby et al. Reference Lockerby, Reese, Emerson and Barber2004) (in Lockerby's notation the mean free path is defined as 
$\lambda =\sqrt {{\rm \pi} /2}\,{Kn}\,L$, and in Sharipov's notation 
$A_{1}=2\sigma _{P}/\sqrt {{\rm \pi} }$ where 
$\sigma _{P}$ is defined as the viscous slip coefficient in Sharipov (Reference Sharipov2011, table 1), and 
$\ell$ is defined by as 
$\ell =\sqrt {2}\,{Kn}\,L$).
The boundary conditions (4.4a–c), when written in the Cartesian coordinates, read
$$\begin{gather} \hat{\sigma}_{ij}n_{i}t_{j} ={-}\sqrt{\frac{2}{{\rm \pi} }}\frac{\chi }{2-\chi } \frac{1}{\eta _{{VS}}}(\hat{V}_{i}+\eta ^{1}\hat{q}_{i}) t_{i}, \end{gather}$$
$$\begin{gather}\hat{q}_{k}n_{k} ={-}\sqrt{\frac{2}{{\rm \pi} }}\frac{\chi }{2-\chi }\frac{ 4+\delta }{2}\frac{1}{\eta _{{TJ}}}(\hat{\mathcal{T}}+\eta ^{1} \hat{\sigma}_{ij}n_{i}n_{j}). \end{gather}$$
Here, we have further introduced the accommodation coefficient 
$\chi$, the velocity slip coefficient 
$\eta _{{VS}}$ and the temperature jump coefficient 
$\eta _{{TJ}}$. For the calculations performed in §§ 4–5, we had taken 
$\chi =\eta _{{VS}}=\eta _{{TJ}}=1$ and 
$\eta ^{1}=3/4$; below, we provide a rationale for this choice.
Substituting the CCR (3.3a–c) into (A2) and using the conservation laws (3.2a–c), one gets
\begin{align} \hat{V}_{i}t_{i} &=\sqrt{\frac{{\rm \pi} }{2}}Kn\frac{2-\chi }{\chi }\eta _{ {VS}}\left( \frac{\partial \hat{v}_{i}}{\partial \hat{x}_{j}}+\frac{ \partial \hat{v}_{j}}{\partial \hat{x}_{i}}\right) n_{i}t_{j}+\eta ^{1}\frac{ Kn}{Pr }\frac{\partial \hat{\sigma}_{ik}}{\partial \hat{x}_{k}}t_{i} \nonumber\\ &\quad +\frac{\eta ^{1}c_{p}}{Pr}Kn\frac{\partial \hat{\theta}}{\partial \hat{x} _{i}}t_{i}+\sqrt{\frac{{\rm \pi} }{2}}\frac{2-\chi }{\chi }\eta _{{VS}}Kn \frac{2}{5+\delta }\left( \frac{\partial \hat{q}_{i}}{\partial \hat{x}_{j}}+ \frac{\partial \hat{q}_{j}}{\partial \hat{x}_{i}}\right) n_{i}t_{j}. \end{align}
In order to compare the last equation with (A1), 
$\hat {\sigma }_{ik}$ and 
$\hat {q}_{i}$ need be replaced by the corresponding NSF equations, i.e. by 
$\hat {\sigma }_{ik}^{{NSF}}$ and 
$\hat {q}_{i}^{{NSF}}$, to get
\begin{align} \hat{V}_{i}t_{i} &=\sqrt{\frac{{\rm \pi} }{2}}Kn\frac{2-\chi }{\chi }\eta _{ {VS}}\left( \frac{\partial \hat{v}_{i}}{\partial \hat{x}_{j}}+\frac{ \partial \hat{v}_{j}}{\partial \hat{x}_{i}}\right) n_{i}t_{j}-\eta ^{1}\frac{ Kn^{2}}{Pr }\frac{\partial }{\partial \hat{x}_{k}}\left( \frac{\partial \hat{v}_{i}}{\partial \hat{x}_{k}}+\frac{\partial \hat{v}_{k}}{\partial \hat{x}_{i}}\right)t_{i}\nonumber\\ &\quad +\frac{\eta ^{1}c_{p}}{Pr }Kn\frac{\partial \hat{\theta}}{\partial \hat{x} _{i}}t_{i}+\sqrt{\frac{{\rm \pi} }{2}}Kn\frac{2-\chi }{\chi }\eta _{{VS}} \frac{2}{5+\delta }\left( \frac{\partial \hat{q}^{{NSF}}_{i}}{\partial \hat{x}_{j}}+ \frac{\partial \hat{q}^{{NSF}}_{j}}{\partial \hat{x}_{i}}\right) n_{i}t_{j}. \end{align}Comparing the last equation with a generalized slip boundary condition (A1), one identifies
\begin{equation} A_{1}=\frac{2-\chi }{\chi }\eta _{{VS}}{, }\quad A_{2}=\frac{2}{{\rm \pi} } \frac{\eta ^{1}}{Pr}\quad \text{and}\quad A_{3}=\frac{\eta ^{1}c_{p}}{Pr }=\sigma _{T}. \end{equation} There exists a vast literature documenting values of the first-order and second-order slip coefficients (Zhang et al. Reference Zhang, Meng and Wei2012); a few are tabulated in table 1. As we can see from the table, the value for 
$A_{1}$ is around unity, however, no conclusion can be reached on the correct value of the coefficient 
$A_{2}$; the value for 
$A_{2}(=({3}/{2{\rm \pi} })({1}/{ Pr}))$ was chosen so that it gives a good match for the drag coefficient for sphere. However, for the chosen values of 
$A_{2}$, the presented boundary conditions yield 
$A_{3}=\frac {3}{4}{c_{p}}/{Pr}$, which is approximately 3.5 times the reported values; in Sharipov (Reference Sharipov2011), the value of 
$A_{3}$ ranges between 0.94 and 1.12. Therefore, one would expect that the presented model/boundary conditions will not provide adequate results for problems involving thermal creep/thermal stress. In order to tackle such problems, one might choose 
$\eta ^{1}$ appropriately (instead of 
$3/4$) to get the correct value for 
$A_{3}$; however, this is beyond the scope of this article.
Table 1. Values of the first-order and second-order slip coefficients (Zhang et al. Reference Zhang, Meng and Wei2012).
Similarly, substituting the CCR (3.3a–c) into (A3) and using the conservation laws, one obtains
\begin{equation} \hat{\mathcal{T}}=\zeta _{T}\sqrt{2}\,Kn\frac{\partial \hat{\theta}}{\partial \hat{x}_{k}}n_{k}+\sqrt{\frac{{\rm \pi} }{2}}\frac{2-\chi }{\chi }\frac{2\eta _{{TJ}}}{4+\delta }\frac{Kn}{Pr }\frac{\partial \hat{\sigma}_{kr}}{ \partial \hat{x}_{r}}n_{k}-\eta ^{1}\hat{\sigma}_{ij}n_{i}n_{j}, \end{equation}where
\begin{equation} \zeta _{T}=\frac{\sqrt{{\rm \pi} }}{2}\frac{2-\chi }{\chi }\frac{2\eta _{{TJ }}}{4+\delta }\frac{c_{p}}{Pr }. \end{equation}
Here, 
$\zeta _{T}$ is the temperature jump coefficient (Sharipov Reference Sharipov2011). The reported value for 
$\zeta _{T}$ by Sharipov (Reference Sharipov2011) is 
$\zeta _{T}=({\sqrt {{\rm \pi} }}/{2}) ({2}/({4+\delta }))({c_{p}}/{Pr })(({2-\chi })/{\chi } +0.17)$; which is approximately 
$17\,\%$ higher than considered value.
In order to quantify the change in drag coefficient with respect to the slight variations in 
$A_2$ (the second-order slip coefficient), 
$\eta _{{VS}}$, 
$\eta _{{TJ}}$, in figure 3(a) we plot the drag on a sphere with different values of 
$A_2$, while taking 
$\eta _{{VS}}=\eta _{{TJ}}=1$. The second-order slip coefficient has a notable effect on the drag, in particular at large values of the Knudsen number. On the other hand, the small variations in values for 
$\eta _{{VS}}$ 
$\eta _{{TJ}}$ lead to slight changes in the drag coefficients – as seen in figure 3(b) – of less than 10 % within the reported values of 
$\eta _{{VS}}$ and 
$\eta _{{TJ}}$.
Figure 3. (a) The normalized drag force on a sphere against Knudsen number. The results predicted by the CCR theory with different values of 
$A_2$, while taking 
$\eta _{{VS}}=\eta _{{TJ}}=1$. (b) Relative percentage error in drag force predicted by the CCR theory for 
$A_{2}=({3}/{2{\rm \pi} })({1}/{ Pr })$ for different values of 
$\eta _{{VS}}$ and 
$\eta _{{TJ}}$.
