2. Single-particle distribution function

The Boltzmann equation provides the kinetic theory description of the evolution of thermodynamic systems

(2.1) \begin{align} \frac {\partial f}{\partial t}+ \frac {\partial }{\partial \boldsymbol{x}} \boldsymbol{\cdot }({\boldsymbol{c}} f) = {\mathscr{J}}(f,f), \end{align}

where $f = f (\boldsymbol{x}, \boldsymbol{c}, t)$ is the single-particle distribution function, ${\mathscr{J}}(f, f)$ is the binary collision operator and $\boldsymbol{x}$ , $\boldsymbol{c}$ , $t$ are the location vector, molecular velocity and time, respectively. The solution of the Boltzmann equation (2.1) is the single-particle distribution function which serves as the connecting link between the microscopic world and the macroscopic world. Once the form of the distribution function is known on solving the Boltzmann equation, the macroscopic field variables can be obtained by taking the moments of the distribution function (Chapman & Cowling Reference Chapman and Cowling1970; Cercignani Reference Cercignani1975; Sone Reference Sone2000; Zohar et al. Reference Zohar, Lee, Lee, Jiang and Tong2002; Agrawal Reference Agrawal2011; Agrawal et al. Reference Agrawal, Kushwaha and Jadhav2020). The general Onsager-consistent approach employed in the present work for evaluating an approximate form of the distribution function is discussed in detail in Mahendra & Singh (Reference Mahendra and Singh2013), Singh & Agrawal (Reference Singh and Agrawal2016), Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023). Here, we provide only the key steps for brevity in the present section.

With the help of Chapman–Enskog expansion, the zeroth-order Maxwellian ( $f_0$ ) and first-order ( $f_1$ ) distribution function can be employed to derive the Euler and N–S equations, respectively. The expressions for $f_0$ and $f_1$ are given as

(2.2) \begin{align} f_0(\boldsymbol{x}, \boldsymbol{c}, t) &= \frac {\rho }{m}\left (\frac {\beta }{\pi }\right )^{(3/2)} \exp { [-\beta ( \lvert \boldsymbol{c}-\boldsymbol{u} \rvert )^2]}, \end{align}

(2.3) \begin{align} f_{1} (\boldsymbol{x}, \boldsymbol{c}, t) &= f_0(\boldsymbol{x}, \boldsymbol{c}, t) - \underbrace {\sum _j \varUpsilon _j \odot X_j,}_{\bar {f}_{1}}, \end{align}

where $m$ is the molecular mass, $\beta = 1/(2 R T)$ , $R$ and $T(\boldsymbol{x}, t)$ are the specific gas constant and absolute temperature, respectively, and $\rho (\boldsymbol{x}, t)$ and $\boldsymbol{u}(\boldsymbol{x}, t)$ are the bulk density and velocity, respectively. The symbol $\odot$ denotes a full tensorial contraction for tensors of the same order. In (2.3), $f_{1}$ is obtained using the iterative refinement technique (Mahendra & Singh Reference Mahendra and Singh2013), and is expressed in terms of thermodynamic forces $X_j$ and the corresponding microscopic conjugate fluxes $\varUpsilon _j$ . Here, the subscript $j = [\tau , q]$ , where $\tau$ represents destabilising viscous processes and $q$ denotes thermal non-equilibrium processes. It is important to emphasise that these formulations of $X_j$ and $\varUpsilon _j$ are as per the OSP.

Similar to the first-order correction term ( $\bar {f}_{1}$ ) to the Maxwellian distribution function, the second-order correction term ( $\bar {f}_{2}$ ) has been derived using the iterative refinement technique in Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023, Reference Yadav, Jonnalagadda and Agrawal2024) as

(2.4) \begin{align} {f} &= f_0 +\bar {f}_{1} + \underbrace {\left ({\varUpsilon }_{\tau \tau }\odot {X}_{\tau }\right ) \odot X_\tau + ({\varUpsilon }_{qq}\odot {X}_{q})\odot X_q}_{{\bar {f}_{2,1}}} \nonumber \\ &\quad +\, \underbrace {{\varUpsilon }_{\tau }\odot \left ({X}_{\tau \tau }\odot {X}_{\tau }\right )+{\varUpsilon }_{q}\odot ({X}_{qq}\odot {X}_{q})}_{{\bar {f}_{2,2}}}, \end{align}

where the explicit expressions of ${\bar {f}_{2,1}}$ and ${\bar {f}_{2,2}}$ are documented in Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023) and Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2024), respectively, and which depend upon either the Euler or N–S equations for substituting the material derivatives in terms of spatial derivatives. For completeness and clarity, explicit expressions are presented in Appendix A, which serves as a direct reference for the formulations discussed in this section. However, readers interested in a comprehensive derivation and its relationship with ${\textit{Kn}}$ are encouraged to refer to the works of Mahendra & Singh (Reference Mahendra and Singh2013), Singh & Agrawal (Reference Singh and Agrawal2016), Singh et al. (Reference Singh, Jadhav and Agrawal2017) and Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023, Reference Yadav, Jonnalagadda and Agrawal2024), where these aspects are thoroughly detailed.

3. Generalised set of 13-moment equations

The generalised, three-dimensional, 13-moment transport equations are obtained after evaluating the moments of the Boltzmann equation. Explicitly, this evaluation, which involves computing $\langle \varPsi , ({\partial f}/{\partial t})+ ({\partial }/{\partial \boldsymbol{x}}) \boldsymbol{\cdot }({\boldsymbol{c}} f) \rangle = \langle \varPsi , {\mathscr{J}}(f,f) \rangle$ , where $\varPsi = \lbrace 1, c_i, ({C_i^2}/{2}), C_{\langle i}C_{j\rangle }, C_i({C_k^2}/{2})\rbrace$ , yields the following equations (Grad Reference Grad1949):

(3.1) \begin{align} \frac {\partial \rho }{\partial t}+\frac {\partial \rho u_k}{\partial x_k} &= 0, \end{align}

(3.2) \begin{align} \rho \frac {\partial u_i}{\partial t} +\rho u_k \frac {\partial u_i}{\partial x_k} +\frac {\partial p}{\partial x_i} &+\frac {\partial \sigma _{\textit{ik}}} {\partial x_k}= {\rho G_i}, \end{align}

(3.3) \begin{align} \rho \frac {\partial \epsilon }{\partial t} + \rho u_k \dfrac {\partial \epsilon }{\partial x_k}+\frac {\partial q_k}{\partial x_k} &+p\frac {\partial u_k}{\partial x_k} + \sigma _{\textit{ij}}\frac {\partial u_i}{\partial x_j}={0}, \end{align}

(3.4) \begin{align} \frac {\partial \sigma _{\textit{ij}}}{\partial t}+ u_k\frac {\partial \sigma _{\textit{ij}}}{\partial x_k}+ \frac {4}{5} \frac {\partial q_{\langle i}}{\partial x_{j\rangle }} + 2\sigma _{k\langle i}\frac {\partial u_{j \rangle }}{\partial x_k} &+ 2p\frac {\partial u_{\langle i}}{\partial x_{j \rangle }} + \sigma _{\textit{ij}}\frac {\partial u_k}{\partial x_k} \nonumber \\ +\, \frac {\partial }{\partial x_k} \underline {\left \langle C_{\langle i} C_j C_{k \rangle }, f \right \rangle } & \left \langle C_{\langle i} C_{j \rangle }, {\mathscr{J}}(f,f) \right \rangle \!, \end{align}

(3.5) \begin{align} \frac {\partial q_i}{\partial t} + u_k \frac {\partial q_i}{\partial x_k} + \frac {5}{2} \left (\frac {p}{\rho } \frac {\partial p }{\partial x_i} - \frac {p^2}{\rho ^2} \frac {\partial \rho }{\partial x_i} \right ) &+ \frac {\partial }{\partial x_j} \underline {\frac {1}{2} \big \langle {\lvert C \rvert }^2 C_{\langle i} C_{j \rangle }, f\big\rangle } -\frac {5}{2} \frac {p}{\rho } \frac {\partial \sigma _{\textit{ik}} }{\partial x_k} \nonumber \\ - \frac {\sigma _{\textit{ik}}}{\rho } \frac {\partial p} {\partial x_k} + \frac {1}{6} \frac {\partial }{\partial x_i} \underline {\big\langle { \lvert C \rvert }^4, ( f-f_0) \big \rangle } & -\, \frac {\sigma _{\textit{ij}}}{\rho } \frac {\partial \sigma _{jk}}{\partial x_k} +\frac {7}{5} q_k \frac {\partial u_i}{\partial x_k }\nonumber + \frac {7}{5} q_i \frac {\partial u_k}{\partial x_k} \\ +\, \frac {2}{5} q_k \frac {\partial u_k}{\partial x_i} + \frac {\partial u_j}{\partial x_k} \underline {\left \langle C_{\langle i} C_j C_{k \rangle }, f \right \rangle } =& \frac {1} {2} \big\langle {\lvert \boldsymbol{C} \rvert }^2 C_i, {\mathscr{J}}(f,f) \big\rangle \!. \end{align}

Here, $G_i$ represents body force, $q_i$ denotes the heat flux and $\sigma _{\textit{ij}}$ the stress tensor, and $C_{\langle i} C_j C_{k \rangle }$ is a third-order trace-free symmetric tensor, defined as (Struchtrup Reference Struchtrup2005; Agrawal et al. Reference Agrawal, Kushwaha and Jadhav2020)

(3.6) \begin{align} w_{\langle ijk \rangle } = w_{(ijk)} - \frac {1}{5} \left ( w_{(ill)} \delta _{jk} + w_{(jll)} \delta _{\textit{ik}} + w_{(kll)} \delta _{\textit{ij}} \right )\!, \end{align}

where round brackets indicate the symmetric part of a tensor. Note that only the underlined terms involving angular brackets with a comma denote moments, whereas angular brackets with indices denote traceless symmetric tensors. Here, $\epsilon = ({3}/{2}) R T$ is the specific internal energy per unit mass for a monatomic ideal gas, while $p=\rho R T$ is the thermodynamic pressure as represented by the ideal gas law. Equations (3.1)–(3.3) are the three standard hydrodynamic conservation equations. Note that the production terms in these equations, $\langle \lbrace 1, C _i, {C_i^2}/{2} \rbrace , {\mathscr{J}}(f,f) \rangle$ , vanish in accordance with the principles of conservation of mass, momentum and energy for particles undergoing elastic collisions. For (3.4) and (3.5), which describe the evolution of the stress tensor and heat flux vector, the production terms appear on the right-hand side. For Maxwell molecules, these production terms obtained using the BGK collision model are given as (Truesdell & Muncaster Reference Truesdell and Muncaster1980)

(3.7) \begin{align} \left \langle C_{\langle i} C_{j \rangle }, {\mathscr{J}}(f,f) \right \rangle & = -\frac {p}{\mu } \sigma _{\textit{ij}}, \end{align}

(3.8) \begin{align} \frac {1} {2} \big \langle {\lvert \boldsymbol{C} \rvert }^2 C_i, {\mathscr{J}}(f,f) \big \rangle & = -\frac {2}{3} \frac {p}{\mu } q_i, \end{align}

where $\sigma _{\textit{ij}}$ and $q_i$ contain first- and second-order contributions. The dynamic viscosity and thermal conductivity are temperature-dependent functions given by $\mu (T) = \mu _0(T/T_0)^\varphi$ and $\kappa (T) = \kappa _0(T/T_0)^\varphi$ , where $\mu _0$ and $\kappa _0$ correspond to their values at a reference temperature $T_0$ , and $\varphi$ represents the interaction potential between gaseous molecules. The value of $\varphi$ is unity for Maxwell molecules and $1/2$ for hard-sphere molecules (García-Colín et al. Reference García-Colín, Velasco and Uribe2008). Equations (3.4) and (3.5) contain three unknown underlined higher-order moments (closure relations) that need to be computed to obtain a closed set of 13 moment equations.

3.1. Closure relations

To achieve closure for (3.1)–(3.5), all three underlined closure relations have been evaluated by substituting the derived distribution function (2.4) in ${\langle C_{\langle i} C_j C_{k \rangle }, f \rangle }, ({{1}/{2}) \langle {\lvert C \rvert }^2 C_{\langle i} C_{j \rangle }, f\rangle }$ and $\langle { \lvert C \rvert }^4, ( f-f_0) \rangle$ , which requires integration with respect to peculiar velocity ( $C_i$ ). By following Jin & Slemrod (Reference Jin and Slemrod2001) and Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023), we, in the first step, make it more compact using the ideal gas equation, the N–S, and Fourier laws by replacing the N–S stress and heat flux terms present in the closure relations. In the second step, we neglect the fourth-order term in Knudsen number from the closure term to maintain the third-order accuracy of the governing set of SO13 equations since the remaining terms are of third-order in Knudsen number (Timokhin et al. Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017; Jadhav & Agrawal Reference Jadhav and Agrawal2021a ). Note that these fourth-order and third-order terms result from employing the N–S and Euler equations-based distribution function terms to evaluate the closure relations. After algebraic simplification, we obtain the following closure relationship:

(3.9) \begin{align} \kern-12pt \!\! \left \langle C_{\langle i}C_j C_{k \rangle }, f_2\right \rangle & = - 6 \frac {p R T t_{r (\tau )}^{2} }{\rho } \frac {\partial \rho }{\partial x_{\langle i}} \frac {\partial u_{j}}{\partial x_{k \rangle }} \!+\! t_{r (\tau )}\! \left (\!- {2{a_s} R T }\frac {\partial \sigma _{\langle ij} }{\partial x_{k \rangle }} \!+\! \frac {2 (8{a_s}\!-\!27) }{15} q_{\langle i} \frac {\partial u_{j}}{\partial x_{k \rangle }}\! \right )\!, \end{align}

(3.10) \begin{align} \kern-3.5pt \left \langle \frac {1}{2}C_k^2 C_{\langle i} C_{j \rangle }, f_2\right \rangle & = \frac {7 }{2}\frac { p }{\rho } \sigma _{i j}+ t_{r (\tau )} \left (- \frac {24{a_{\textit{qv}}} }{5} R T \frac {\partial q_{\langle i}}{\partial x_{j \rangle }} + \frac {75}{10}\frac {1}{4} \frac {R T}{\mu }\sigma _{k\langle i}\sigma _{j\rangle k } \right. \nonumber \\ & \left. \quad +\, \frac {145 }{10}\frac {1}{2} R T\sigma _{k\langle i} \frac {\partial u_{j} }{\partial x_{k \rangle }} - \frac {1}{2}\left [\frac {49}{3} - \frac {14 \varphi }{3}\right ] R T \sigma _{i j} \frac {\partial u_{l}}{\partial x_{l}} + \frac {91 }{15} \frac {1}{ \rho } q_{\langle i} \frac {\partial p}{\partial x_{j \rangle }} \right ) \nonumber \\ & \quad +\, \frac {56- 64{a_{\textit{qv}}} }{25}\frac {1}{ p} q_{\langle i} q_{j \rangle }, \end{align}

(3.11) \begin{align} \kern-3.5pt \left \langle {\lvert {C} \rvert }^4, \left (f_2-f_0\right )\right \rangle & = t_{r (\tau )} \left (\frac {24{a_{\textit{qs}}}}{5} R T \frac {\partial q_k}{\partial x_{k}} + 12\frac { q_{k} }{\rho } \frac {\partial p}{\partial x_{k}} \right ) + \frac {8\big (8{a_{\textit{qs}}} +35\big)}{25} \frac {1}{ p} q_k q_{k}, \end{align}

where parameter $t_{r(\tau )}$ represents the relaxation time for momentum. Here, $a_{s} = 1.61$ , $a_{\textit{qv}} = 1/2$ and $a_{\textit{qs}} = -5/2$ satisfy the compatibility condition, which represent the additive invariant property of kinetic theory (Balakrishnan et al. Reference Balakrishnan, Agarwal and Yun1999; Agarwal et al. Reference Agarwal, Yun and Balakrishnan2001) and are independent of the Knudsen number ( ${\textit{Kn}}$ ). For more in-depth understanding, interested readers are encouraged to refer to Balakrishnan et al. (Reference Balakrishnan, Agarwal and Yun1999), Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023) and Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2024). This leads to linear stability, OSP consistency, and compliance with the second law of thermodynamics, as demonstrated later in the manuscript. Note that, despite the above-mentioned modification, while deriving the final form of the closure reactions, the accuracy level of these equations remain intact. However, the modification alters the mathematical characteristics of (3.9)–(3.11), resulting in both linear and nonlinear forms of the equations needing an equivalent number of boundary conditions. Equations (3.1)–(3.5) along with the closure relations (3.9)–(3.11) form the closed set of HOT equations being proposed in this work.

4. Linear stability analysis

In this section, we establish the one-dimensional stability of the linearised equations presented in the previous section. For this purpose, we first represent relevant quantities in non-dimensional form as follows (Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2024):

(4.1) \begin{align} \rho = \rho _o(1+\bar {\rho }),\, T&=T_o(1+\bar {T}),\, p = p_o(1+\bar {p}),\, u_1=\sqrt {RT_o}\bar {u}_{1},\,\sigma _{11} = \rho _o R T_o \bar {\sigma }_{11},\nonumber \\ q_{1} &= \rho _o\left (\sqrt {R T_o}\right )^3\bar {q}_{1},\, x = H \bar {x}, \,t = \frac {H}{\sqrt {RT_o}}\bar {t}. \end{align}

The subscript ‘ $o$ ’ indicates the variable at equilibrium state, the subscript ‘ $1$ ’ indicates a relevant quantity in the $x$ -direction and $H$ represents the characteristic length. Meanwhile, the quantities with an overbar ( $\bar {\boldsymbol{\cdot }}$ ) represent small perturbations around the equilibrium or rest state. Assuming that the perturbations are small, the analysis focuses only on linear terms, significantly simplifying the calculations. As a result, a reduced set of linearised and non-dimensionalised equations starting from (3.1)–(3.5), with closure relations given by (3.9)–(3.10), are obtained as

(4.2a) \begin{align} \frac {\partial {{\bar {\rho }}}}{\partial \bar {t}} + \frac {\partial {{\bar {u}}_{1}}}{\partial \bar {x}} &= 0, \end{align}

(4.2b) \begin{align} {\frac {\partial {{\bar {u}}_{1}}}{\partial \bar {t}}}+\frac {\partial {{\bar {T}}}}{\partial \bar {x}} + \frac {\partial {{\bar {\rho }}}}{\partial \bar {x}} + \frac {\partial {{\bar {\sigma }}_{1 1 }}}{\partial \bar {x}} &= 0 , \end{align}

(4.2c) \begin{align} \frac {3 }{2} \frac {\partial {{\bar {T}}}}{\partial \bar {t}} + \frac {\partial {{\bar {q}}_{1}}}{\partial \bar {x}} + \frac {\partial {{\bar {u}}_{1}}}{\partial \bar {x}} &= 0 , \end{align}

(4.2d) \begin{align} \frac {\partial {{\bar {\sigma }}_{1 1 }}}{\partial \bar {t}} - \frac {483 }{250} Kn \frac {\partial ^{2} {{\bar {\sigma }}_{1 1 }}}{\partial \bar {x}^{2}} + \frac {8 }{15} \frac {\partial {{\bar {q}}_{1}}}{\partial \bar {x}} + \frac {4 }{3} \frac {\partial {{\bar {u}}_{1}}}{\partial \bar {x}} + \frac {{{\bar {\sigma }}_{1 1 }}}{Kn} &= 0 , \end{align}

(4.2e) \begin{align} \frac {\partial {{\bar {q}}_{1}} }{\partial \bar {t}} - \frac {18 }{5} Kn \frac {\partial ^{2} {{\bar {q}}_{1}}}{\partial \bar {x}^{2}} + \frac {5 }{2} \frac {\partial {{\bar {T}}}}{\partial \bar {x}} + \frac {\partial {{\bar {\sigma }}_{1 1 }}}{\partial \bar {x}} + \frac {2 {{\bar {q}}_{1}}}{3 Kn} &= 0 . \end{align}

Now, we employ the normal mode perturbation method and assume the solution of primary variables to be of the form

(4.3) \begin{align} \bar {\rho } &=\rho _A \exp {(\omega \bar {t}+i \kappa \bar {x})}, \,\bar {u}= u_{1A} \exp {(\omega \bar {t}+i \kappa \bar {x})},\, \bar {T}= T_A \exp {(\omega \bar {t}+i \kappa \bar {x})} ,\nonumber \\ \bar {\sigma }_{11} & = \bar {\sigma }_{11A}\exp {(\omega \bar {t}+i \kappa \bar {x})},\,\bar {q}_{1}= \bar {q}_{1A}\exp {(\omega \bar {t}+i \kappa \bar {x})}. \end{align}

Here, the variables $\kappa$ , $\omega$ and the subscript ‘ $A$ ’ represent the wavenumber, wave frequency and complex amplitude of the plane wave, respectively. By substituting these solutions (4.3) into (4.2a )–(4.2e ), a relationship between $\kappa$ and $\omega$ is obtained. This relationship, known as the dispersion relation, can be expressed as

(4.4) \begin{align} - \frac {3 \omega ^{5}}{2} &+ \omega ^{4} \left (- \frac {36 \kappa ^{2}}{5} - \frac {5}{2}\right ) + \omega ^{3} \left (- \frac {162 \kappa ^{4}}{25} - \frac {72 \kappa ^{2}}{5} - 1\right ) + \omega ^{2} \left (- \frac {111 \kappa ^{4}}{5} - 8 \kappa ^{2}\right ) \nonumber \\ & + \omega \left (- \frac {54 \kappa ^{6}}{5} - \frac {31 \kappa ^{4}}{2} - \frac {5 \kappa ^{2}}{3}\right ) +\left (- 3 \kappa ^{6} - \frac {5 \kappa ^{4}}{2} \right ) =0 .\end{align}

Upon obtaining (4.4), we examine the roots of this fifth-order polynomial to analyse the stability of the proposed equations. These roots, which can have both real and complex parts, provide valuable insights into the system’s behaviour when subjected to disturbances in space. As a result, we consider a real wavenumber and complex frequency represented by $\omega = \omega _r(\kappa ) + i\omega _i(\kappa )$ .

In order to ensure stability, it is essential to satisfy the condition $\omega _r(\kappa )= {\textit{Re}} (\omega ) \leqslant 0$ , indicating that the real part, ${\textit{Re}} (\omega )$ , of the frequency is negative. This ensures that the local amplitude of primary variables decreases over time in the presence of spatial disturbances. The stability analysis is reinforced by the solutions obtained from (4.4), which yield five roots for $\omega$ . These roots are shown in figure 1(a), clearly illustrating the stability criterion. Further, variation of the attenuation coefficient ( $\omega _r(\kappa )$ ) has been shown in figure 1(b), in which $ Kn\,( = \mu _o \sqrt {R T_o }/(p_o H)) \approx \kappa$ is adopted for convenience, as demonstrated in Balakrishnan et al. (Reference Balakrishnan, Agarwal and Yun1999). This assumption serves as an approximate representation of the validity and stability range of the proposed equations. It is important to note that the absence of this assumption does not affect the generality of the stability analysis of the present equations. Figure 1(b) shows that all five roots are negative for any value of the Knudsen number. Thus, the proposed SO13 moment equations are unconditionally stable for small spatial disturbances and are therefore free from the Bobylev instability, a well-known issue that plagues the Burnett and other HOT equations.

Figure 1. (a) Stability curve of the SO13 equations due to spatial perturbations. (b) Variation of attenuation coefficient with Knudsen number.

5. Compliance with Onsager symmetry principle

Here, we follow the procedure outlined by Romero & Velasco (Reference Romero and Velasco1995) and express the perturbed field variables ( ${\hat{\rho}}, {\hat{u}}_i, {\hat{\epsilon}}, {\hat{\sigma}}_{\textit{ij}}, {\hat{q}}_{i}$ ) in terms of thermodynamic forces in the corresponding hydrodynamic fields ( $X_\rho , X_{u,i}, X_\epsilon , X_{\sigma _{\textit{ij}}}, X_{q_i}$ ). The linearisation is performed about a local equilibrium state using the following perturbation structure: $ \rho = \rho _o + {\hat{\rho}}, \quad T = T_o + {\hat{T}}, \quad p = p_o + {\hat{p}}, \quad u_i = {\hat{u}}_i, \quad \sigma _{\textit{ij}} = {\hat{\sigma}}_{\textit{ij}}, \quad q_i = {\hat{q}}_i,$ where the subscript $o$ denotes equilibrium values and variables with a hat ( ${\hat{\,}}$ ) represent small deviations around the equilibrium state. Explicitly, these forces are given as (Singh & Agrawal Reference Singh and Agrawal2016)

(5.1) \begin{align} X_\rho = \frac {1 }{\rho _o T_o}{\hat{p}}-\frac {p_o}{\rho _o}X_\epsilon ,\, X_{u_i}= \frac {1}{T_o} {\hat{u}}_i, \, X_\epsilon = \frac {1}{T^2_o} {\hat{T}}, \, X_{\sigma _{\textit{ij}}} = \frac {1}{2p_oT_o}{{\hat{\sigma}}_{\textit{ij}}},\, X_{q_i} = \frac {2}{5 p_o R T^2_o}{\hat{q}}_i. \end{align}

We further substitute the perturbed field variables and the closure relations given by (3.9)–(3.10) in (3.1)–(3.5). Note that the ideal gas law is explicitly employed to express the pressure in terms of density. The governing equations are then linearised after substitution, which results in the following simplified SO13 equations containing both first- and second-order derivative terms:

(5.2) \begin{align} \frac {\partial {\hat{\rho}}}{\partial t} &= -\rho _o\frac {\partial {\hat{u}}_k}{\partial x_k}, \end{align}

(5.3) \begin{align} \rho _o\frac {\partial {\hat{u}}_i}{\partial t} &= -\frac {\partial {\hat{p}}}{\partial x_i} -\frac {\partial {\hat{\sigma}}_{\textit{ik}}}{\partial x_k}, \end{align}

(5.4) \begin{align} \rho _o \frac {\partial \hat{\epsilon }}{\partial t} &= -\frac {\partial {\hat{q}}_k}{\partial x_k} -p_o\frac {\partial {\hat{u}}_{k}}{\partial x_k},\\[-6pt] \nonumber \end{align}

(5.5) \begin{align} \frac {\partial {\hat{\sigma}}_{\textit{ij}}}{\partial t} = \frac {161R T_o \mu _{o} }{100p_{o}} \frac {\partial }{\partial x_{k}} \frac {\partial {{\hat{\sigma}}}_{\langle i j}}{\partial x_{k \rangle }} &-\frac {4}{5}\frac {\partial {\hat{q}}_{\langle i}}{\partial x_{ j \rangle }} -2p_o\frac {\partial {\hat{u}}_{\langle i}} {\partial x_{j \rangle }} - \frac {p_o}{\mu _o} {\hat{\sigma}}_{\textit{ij}}, \\[-6pt] \nonumber\end{align}

(5.6) \begin{align} \frac {\partial {\hat{q}}_{i}}{\partial t} = \frac {1}{2}\frac {R T_o \mu _{o} }{ p_{o}} \frac {\partial }{\partial x_{j}} \frac {\partial {{\hat{q}}}_{\langle i}}{\partial x_{j \rangle }} & -{\frac {5}{12}\frac {R T_o \mu _{o} }{ p_{o}} \frac {\partial ^{2}{{\hat{q}}}_{ k}}{\partial x_{i}\partial x_{k }} }- \frac {5}{2}R p_o\frac {\partial {\hat{T}}}{\partial x_{ i }} - \frac {p_o}{\rho _o}\frac {\partial {\hat{\sigma}}_{\textit{ik}}} {\partial x_{k}} - \frac {2}{3} \frac {p_o}{\mu _o} {\hat{q}}_{i}, \end{align}

where $\mu _o$ denotes equilibrium values. Finally, (5.1) is used to replace the perturbed quantities present in the linearised moment equations (5.2)–(5.6), which yields

(5.7) \begin{align} \frac {\partial {\hat{\rho}}}{\partial t}& = - T_{o} \rho _{o} \frac {\partial X_{u_ i}}{\partial x_{i}} , \end{align}

(5.8) \begin{align} \frac {\partial {\hat{u}}_i}{\partial t} = - T_{o} p_{o} \frac {\partial X_{\epsilon }}{\partial x_{i}} & - 2 T_{o} p_{o} \frac {\partial X_{\sigma _{ i j}}}{\partial x_{j}} - T_{o} \rho _{o} \frac {\partial X_{\rho }}{\partial x_{i}}, \end{align}

(5.9) \begin{align} \frac {\partial \hat{\epsilon }}{\partial t}= - \frac {5 R T_{o}^{2} p_{o} }{2} \frac {\partial X_{q_i}}{\partial x_{i}} & - T_{o} p_{o} \frac {\partial X_{u_i}}{\partial x_{i}} , \end{align}

(5.10) \begin{align} \frac {\partial {{\hat{\sigma}}_{\textit{ij}}}}{\partial t}= \frac {161 }{50} { R T_{o} \mu _{o} } \frac {\partial }{\partial x_{k}} \frac {\partial X_{\sigma _{ \langle i j} }}{\partial x_{k \rangle }} & - 2 R T_{o}^{2} p_{o} \frac {\partial X_{q _{\langle i}}}{\partial x_{j \rangle }} - \frac {2 T_{o} p_{o}^{2}}{\mu _{o}}X_{\sigma _{ i j}} - 2 T_{o} p_{o} \frac {\partial X_{u_{\langle i}}}{\partial x_{j \rangle }} , \end{align}

(5.11) \begin{align} \frac {\partial {\hat{q}}_i}{\partial t} = \frac {5 }{24} { R^{2} T_{o}^{2} \mu _{o} } \frac {\partial ^{2} X_{q \langle i}}{\partial x_{j}\partial x_{j \rangle }} & - \frac {5 R T_{o}^{2} p_{o}^{2}}{3 \mu _{o}} X_{q_i}- \frac {5 R T_{o}^{2} p_{o} }{2} \frac {\partial X_{\epsilon }}{\partial x_{i}} - 2 R T_{o}^{2} p_{o} \frac {\partial X_{\sigma _{i k}}}{\partial x_{k}}. \end{align}

Further, (5.7)–(5.11) can be recast into a compact form as

(5.12) \begin{align} \frac {\partial \alpha _i }{\partial t} = - \sum _{j} L_{\textit{ij}}X_j, \end{align}

where $X_j$ is provided in (5.1), and $\alpha _i$ represents all the primary variables present in the proposed equations and the matrix of phenomenological coefficients, $L_{\textit{ij}}$ , is given as

(5.13) \begin{align} {L_{\textit{ij}} =}{\left [\begin{matrix}0 & - T_{o} \rho _{o} \boldsymbol{\nabla }& 0 & 0 & 0\\- T_{o} \rho _{o} \boldsymbol{\nabla }& 0 & - T_{o} p_{o} \boldsymbol{\nabla }& - 2 T_{o} p_{o} \boldsymbol{\nabla }& 0\\0 & - T_{o} p_{o} \boldsymbol{\nabla }& 0 & 0 & - \frac {5 R T_{o}^{2} p_{o} \boldsymbol{\nabla }}{2}\\0 & - 2 T_{o} p_{o} \boldsymbol{\nabla }& 0 & \frac {161 R T_{o} \mu _{o} \boldsymbol{\nabla} ^{2}}{50} - \frac {2 T_{o} p_{o}^{2}}{\mu _{o}} & - 2 R T_{o}^{2} p_{o} \boldsymbol{\nabla }\\0 & 0 & - \frac {5 R T_{o}^{2} p_{o} \boldsymbol{\nabla }}{2} & - 2 R T_{o}^{2} p_{o} \boldsymbol{\nabla }& \frac {5 R^{2} T_{o}^{2} \mu _{o} \boldsymbol{\nabla} ^{2}}{24} - \frac {5 R T_{o}^{2} p_{o}^{2}}{3 \mu _{o}}\end{matrix}\right ]}. \end{align}

Here, $\boldsymbol{\nabla}$ denotes the first-order spatial derivative operator (i.e. $\partial / \partial x_i$ ), and $\boldsymbol{\nabla} ^2$ represents the Laplacian operator. It is important to note that the matrix $L_{\textit{ij}}$ does not contain any terms involving the operator $({\partial ^2}/{\partial x_i \partial x_k})$ .

A system of equations of the form of (5.12) is compliant with the OSP if the Hermitian conjugate of the matrix of phenomenological coefficients, $L^{\dagger }$ , satisfies the relation $L^{\dagger } = DLD$ (McLennan Reference McLennan1974; Romero & Velasco Reference Romero and Velasco1995). Here, $D$ is a diagonal matrix with elements $D_{\textit{ii}} = \mp 1$ based on the odd/even parity of the primary variables under time reversal. For (5.7)–(5.11), $D$ is explicitly evaluated to be

(5.14) \begin{equation} D= \left [ \begin{matrix} & 1 & 0 & 0 & 0 & 0 &\\ & 0 & -1 & 0 & 0 & 0 &\\ & 0 & 0 & 1 & 0 & 0 &\\ & 0 & 0 & 0 & 1 & 0 &\\ & 0 & 0 & 0 & 0 &-1 & \end{matrix} \right ]\!. \end{equation}

For matrix $L$ given in (5.13), $L^\dagger$ is given as

(5.15) \begin{equation} {\left [\begin{matrix}0 & T_{o} \rho _{o} \boldsymbol{\nabla }& 0 & 0 & 0\\T_{o} \rho _{o} \boldsymbol{\nabla }& 0 & T_{o} p_{o} \boldsymbol{\nabla }& 2 T_{o} p_{o} \boldsymbol{\nabla }& 0\\0 & T_{o} p_{o} \boldsymbol{\nabla }& 0 & 0 & \frac {5 R T_{o}^{2} p_{o} \boldsymbol{\nabla }}{2}\\0 & 2 T_{o} p_{o} \boldsymbol{\nabla }& 0 & \frac {161 R T_{o} \mu _{o} \boldsymbol{\nabla} ^{2}}{50} - \frac {2 T_{o} p_{o}^{2}}{\mu _{o}} & 2 R T_{o}^{2} p_{o} \boldsymbol{\nabla }\\0 & 0 & \frac {5 R T_{o}^{2} p_{o} \boldsymbol{\nabla }}{2} & 2 R T_{o}^{2} p_{o} \boldsymbol{\nabla }& \frac {5 R^{2} T_{o}^{2} \mu _{o} \boldsymbol{\nabla} ^{2}}{24} - \frac {5 R T_{o}^{2} p_{o}^{2}}{3 \mu _{o}}\end{matrix}\right ]}, \end{equation}

which can be trivially shown to satisfy $L^\dagger = DLD$ , thus proving that (5.7)–(5.11) satisfy the OSP. Most available HOT equations do not satisfy the OSP, and thus, the compliance of the proposed model with the OSP highlights a significant theoretical contribution.

6. Compliance with second law of thermodynamics

From Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2007), the bulk entropy generation rate for the Grad moment-based transport equations for the Maxwellian molecule is presented as

(6.1) \begin{align} \varSigma _b = \frac {1}{2Kn} \sigma _{\textit{ij}}\sigma _{\textit{ij}} + \frac {2}{5} \frac {Pr}{Kn} q_iq_i - \frac {1}{2} A_{\textit{ijk}} \frac {\partial \sigma _{\textit{ij}}}{\partial x_k} - \frac {9 }{20} Pr^2 B_{\textit{ik}} \frac {\partial q_i}{\partial x_k}, \end{align}

where subscript $b$ in $\varSigma _b$ denotes the bulk variation of the entropy generation rate and ${\textit{Pr}}$ is the Prandtl number. Here, the last two terms arise from the additional terms in the closure expression compared with the G13 equations. Note that the tensors $A_{\textit{ijk}}$ and $B_{\textit{ik}}$ represent only the linear terms in the closure expression. As a result, from (3.9)–(3.11), the expression of $A_{\textit{ijk}}$ and $B_{\textit{ij}}$ are given as

(6.2) \begin{align} A_{\textit{ijk}} = -2 Kn a_s \frac {\partial \sigma _{\langle ij}}{\partial x_{k\rangle }}, \quad B_{\textit{ij}} = -\frac {24}{5} a_{\textit{qv}} Kn \frac {\partial q_{\langle i}}{\partial x_{j\rangle }} + \frac {8}{5} Kn a_{\textit{qs}} \frac {\partial q_k}{\partial x_k}\delta _{\textit{ij}}, \end{align}

where only linear terms are considered for simplicity and tensor $A_{\textit{ijk}}$ represents the linear terms from 3.9, while $B_{\textit{ij}}$ corresponds to the sum of the linear terms from (3.10)–(3.11). Here, the coefficients $a_\alpha$ ( $\alpha \in \{s, qv, qs\}$ ) are already defined in § 3, where the first two coefficients are positive while the last one is negative. For the Maxwell molecule, these coefficients were obtained naturally from kinetic theory while following the compatibility conditions. Unlike our closure, the standard Grad closure does not include terms associated with $ a_\alpha$ ( $\alpha \in \{s, qv, qs\}$ ), effectively neglecting certain dissipative effects. Hence, setting $ a_\alpha = 0$ in our closure equations does not directly recover the Grad 13-moment closure due to the presence of additional nonlinear terms.

Here, $\varSigma _b$ must remain non-negative, $\varSigma _b \geqslant 0$ , for all variables such as $\rho , u_i, {T}, \sigma _{\textit{ij}}, q_i$ . This condition can be readily shown by substituting (6.2) in (6.1) and then by performing algebraic operations. Consequently, we obtain the following simplified expression for $\varSigma _b$  as:

(6.3) \begin{align} \varSigma _b &= \frac {1}{2Kn} \sigma _{\textit{ij}}\sigma _{\textit{ij}} + \frac {2}{5} \frac {Pr}{Kn} q_iq_i+ {1.61} Kn \frac {\partial \sigma _{\langle ij}}{\partial x_{k\rangle }} \frac {\partial \sigma _{\textit{ij}}}{\partial x_k} \nonumber \\ &\quad +\, \frac {9}{20} Pr^2 Kn \left ( \frac {12}{5} \frac {\partial q_{\langle i}}{\partial x_{j\rangle }} + {4} \frac {\partial q_k}{\partial x_k}\delta _{\textit{ij}}\right ) \frac {\partial q_i}{\partial x_j} .\end{align}

It is important to observe that cross-coupling terms are absent in (6.3) owing to the restriction of (6.1) to the linearised form of the governing equations, as documented in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2007). Since the last two terms of (6.3) appear with the same sign but with different coefficients, as in the case of the R13 equations presented in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2007), this ensures that $\varSigma _b$ remains positive under all conditions.

7. Poiseuille flow

This section aims to analytically solve the two-dimensional compressible plane Poiseuille flow problem driven by a streamwise external force ( $G$ ) by employing the semi-linearisation method (Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023; Yadav & Agrawal Reference Yadav and Agrawal2024). Despite its simplicity, planar Poiseuille flow serves as an ideal test case for validating the proposed equations due to its distinct non-equilibrium behaviours. Extensive studies have shown that, under dilute conditions, these flows exhibit phenomena beyond the predictive scope of the Navier–Stokes–Fourier framework, including a flow-aligned heat flux without a temperature gradient (Uribe & Garcia Reference Uribe and Garcia1999; Aoki, Takata & Nakanishi Reference Aoki, Takata and Nakanishi2002), a non-uniform pressure profile (Tij & Santos Reference Tij and Santos1994; Aoki et al. Reference Aoki, Takata and Nakanishi2002) and a characteristic temperature depression (Aoki et al. Reference Aoki, Takata and Nakanishi2002; Zheng, Garcia & Alder Reference Zheng, Garcia and Alder2002; Xu Reference Xu2003). This will be followed by comparing our analytical solution with different existing models and DSMC simulation data to validate the analytical solution.

7.1. Analytical solution for Poiseuille flow

The schematic of the force-driven plane Poiseuille flow problem is presented in figure 2, which has two parallel plates separated by a distance $H$ , and the origin is located at the centre of the channel’s entrance. The flow is driven by a body force (shown as ‘ $G$ ’ in the figure). To derive the solution for this problem, we simplify the problem depicted in figure 2 by assuming that the flow remains steady and that all primary variables depend only on the normal direction of the flow. In this analysis, we further assume that the viscosity ( $\mu$ ) remains constant and independent of the temperature. Additionally, it should be noted that the cross-stream velocity ( $v$ ) is zero due to the presence of impermeable bounding walls. By employing the assumptions mentioned above and simplifications for this flow problem, we finally obtain the following variables:

Figure 2. The diagram shows the compressible plane Poiseuille flow problem, driven by an external force ( $G$ ). The upper and lower plate temperatures are assumed to have the same temperature $T_w$ .

(7.1) \begin{align} { u_i = \begin{bmatrix} u_1(y) \\ 0 \\ 0 \end{bmatrix}},\, {\sigma _{\textit{ij}} = \begin{bmatrix} \sigma _{11}(y) & \sigma _{12}(y) & 0 \\ \sigma _{12}(y) & \sigma _{22}(y) & 0 \\ 0 & 0 & -\sigma _{11}(y) - \sigma _{22}(y) \end{bmatrix}},\, {{q}_i = \begin{bmatrix} q_1(y) \\ q_2(y) \\ 0 \end{bmatrix}}. \end{align}

We introduce a set of dimensionless variables in (7.2), which allows further simplification of (3.1)–(3.5)

(7.2) \begin{align} \bar {y} & =\frac {y}{H}, \quad \bar {\rho }= \frac {\rho -\rho _o}{\rho _o}, \quad \bar {T}= \frac {T-T_o}{T_o}, \quad \bar {G} = \frac {H}{RT_o}G , \quad \bar {u}_1= \frac {u_1}{\sqrt {RT_o}}, \quad \bar {\sigma }_{\textit{ij}}= \frac {\sigma _{\textit{ij}}}{\rho _o {R T_o}}, \nonumber\\ \bar {q}_1 &= \frac {q_1}{\rho _o\sqrt {{RT_o}}^3}. \end{align}

The substitution of (7.2) in (3.1)–(3.5) and considering the assumptions for the present analysis, we obtain

(7.3a) \begin{align} x\hbox{-}\text{Momentum equation}&:- \bar {G} + \frac {{\rm d} \bar {\sigma }_{1 2}}{{\rm d} \bar {y}} = 0 , \end{align}

(7.3b) \begin{align} y\hbox{-}\text{Momentum equation}&:\frac {{\rm d} {{\bar {T}}}}{{\rm d} \bar {y}} + \frac {{\rm d} {{\bar {\rho }}}}{{\rm d} \bar {y}} + \frac {{\rm d} {{\bar {\sigma }}_{2 2 }}}{{\rm d} \bar {y}} = 0, \end{align}

(7.4) \begin{align} \text{Energy equation}&: {\bar {\sigma }}_{1 2 } \frac {{\rm d} {{\bar {u}}_{1}} }{{\rm d} \bar {y}} + \frac {{\rm d} {{\bar {q}}_{2}}}{{\rm d} \bar {y}} = 0, \end{align}

(7.5a) \begin{align} \text{Stress equation for}\, \bar {\sigma }_{11}&: \frac {1412 Kn^{2} }{1875} \frac {{\rm d} {{\bar {\sigma }}_{1 2 }}}{{\rm d} \bar {y}} \frac {{\rm d}^{2} {{\bar {u}}_{1}}}{{\rm d} \bar {y}^{2}} - \frac {161 Kn }{150}\frac {{\rm d}^{2} {{\bar {\sigma }}_{1 1 }}}{{\rm d} \bar {y}^{2}} + \frac {161 Kn }{375} \frac {{\rm d}^{2} {{\bar {\sigma }}_{2 2 }}}{{\rm d} \bar {y}^{2}} \nonumber \\ & + \frac {4 }{3} {{\bar {\sigma }}_{1 2 }} \frac {{\rm d} {{\bar {u}}_{1}}}{{\rm d} \bar {y}} - \frac {4 }{15} \frac {{\rm d} {{\bar {q}}_{2}}}{{\rm d} \bar {y}} +\frac {1}{Kn} {{\bar {\sigma }}_{1 1 }}= 0, \end{align}

(7.5b) \begin{align} \text{Stress equation for}\, \bar {\sigma }_{12} &: - \frac {644 Kn }{375} \frac {{\rm d}^{2} {{\bar {\sigma }}_{1 2 }} }{{\rm d} \bar {y}^{2}} + \frac {2 }{5} \frac {{\rm d} {{\bar {q}}_{1}} }{{\rm d} \bar {y}} + \frac {{\rm d} {{\bar {u}}_{1}} }{{\rm d} \bar {y}} + \frac {1}{Kn} {{\bar {\sigma }}_{1 2 }} = 0, \end{align}

(7.5c) \begin{align} \text{Stress equation for}\, \bar {\sigma }_{22}&: - \frac {353 Kn^{2} }{625} \frac {{\rm d} {{\bar {\sigma }}_{1 2 }}} {{\rm d} \bar {y}} \frac {{\rm d}^{2} {{\bar {u}}_{1}} }{{\rm d} \bar {y}^{2}} - \frac {483 Kn }{250}\frac {{\rm d}^{2} {{\bar {\sigma }}_{2 2 }} }{{\rm d} \bar {y}^{2}} - \frac {2 }{3}{{\bar {\sigma }}_{1 2 }} \frac {{\rm d} {{\bar {u}}_{1}}}{{\rm d} \bar {y}} \nonumber \\ & + \frac {8 }{15} \frac {{\rm d} {{\bar {q}}_{2}}}{{\rm d} \bar {y}} + \frac {1}{Kn}{{\bar {\sigma }}_{2 2 }} = 0 ,\end{align}

(7.6a) \begin{align} \text{Heat flux equation for}\, \bar {q}_1&: - \frac {6 }{5}Kn \frac {{\rm d}^{2}{{\bar {q}}_{1}}}{{\rm d} \bar {y}^{2}} + \frac {1}{2} \frac {{\rm d} {{\bar {\sigma }}_{1 2 }}}{{\rm d} \bar {y}} + \frac {1}{3 Kn} {{\bar {q}}_{1}} = 0, \end{align}

(7.6b) \begin{align} \text{Heat flux equation for}\, \bar {q}_2&: - \frac {29 Kn }{12} {{\bar {\sigma }}_{1 2 }} \frac {{\rm d}^{2} {{\bar {u}}_{1}} }{{\rm d} \bar {y}^{2}} - \frac {18 Kn }{5}\frac {{\rm d}^{2} {{\bar {q}}_{2}}}{{\rm d} \bar {y}^{2}} - \frac {2067 Kn }{500} \frac {{\rm d} {{\bar {\sigma }}_{1 2 }}}{{\rm d} \bar {y}} \frac {{\rm d} {{\bar {u}}_{1}}}{{\rm d} \bar {y}} \nonumber \\ &+ \frac {2 }{5}{{\bar {q}}_{1}} \frac {{\rm d} {{\bar {u}}_{1}} }{{\rm d} \bar {y}} + \frac {1}{4}{{\bar {\sigma }}_{1 2 }} \frac {{\rm d} {{\bar {\sigma }}_{1 2 }} }{{\rm d} \bar {y}} + \frac {5 }{2} \frac {{\rm d} {{\bar {T}}}}{{\rm d} \bar {y}} + \frac {{\rm d} {{\bar {\sigma }}_{2 2 }}}{{\rm d} \bar {y}} + \frac {2 }{3 Kn} {{\bar {q}}_{2}} = 0 .\end{align}

The above set of equations is expressed in terms of the Knudsen number ( ${\textit{Kn}} = \mu _0 \sqrt {RT_0}/(p_0 H)$ ), which characterises the degree of rarefaction. Note that the Knudsen number is determined based on the equilibrium state $\rho _0, T_0$ . In (7.3a )–(7.6b ), linear terms are retained along with key nonlinear contributions. Following Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023), the equations include nonlinear viscous heating terms, $\bar {\sigma }_{21}({{\rm d} \bar {u}_1}/{{\rm d}\bar {y}})$ and $\bar {q}_{1}({{\rm d} \bar {u}_1}/{{\rm d}\bar {y}})$ , as well as their combinations, which appear in the SO13 equation and higher-order closure relations. These terms account for key phenomena, including viscous heating and non-hydrodynamic effects. Specifically, the first term in (7.6b ), which is first-order in Knudsen number, arises from this consideration.

Integration of (7.3a ) results in

(7.7) \begin{align} \bar {\sigma }_{21} = \bar {G}\bar {y}+C_{\sigma _{21}}, \end{align}

where $C_{\sigma _{21}}$ represents an integration constant (appropriate subscripts have been added to the constants to indicate their origin). Subsequently, by solving the coupled equations (7.6a ) and (7.5b ) while incorporating (7.7), we obtain the following simultaneous solution:

(7.8) \begin{align} \bar {q}_1 &= C_{{q}_{1 1}} \sinh \left ( A \right )+ C_{{q}_{1 2}} \cosh \left ( A \right ) - {\frac {3 \bar {G} Kn}{2}} , \end{align}

(7.9) \begin{align} {{\bar {u}}_{1}} &= {C_{u_{1}}} - {\frac {C_{\sigma _{21}} }{Kn}\bar {y}} - \frac {2 C_{q_{11}} }{5} \sinh \left ( A \right ) - \frac {2 C_{q_{12}} }{5} \cosh \left ( A \right ) + \frac {3 \bar {G} Kn}{5} - {\frac {\bar {G} }{2 Kn}\bar {y}^{2}}, \end{align}

where $A= { ({\sqrt {5} \bar {y}}/{3 Kn} )}$ , the integration constants are represented by $C_{q_{1 1}}$ , $C_{q_{1 2}}$ and $C_{u_{1}}$ . Upon incorporating (7.7)–(7.9) into (7.4), we obtain the subsequent expression

(7.10) \begin{align} {{\bar {q}}_{2}} & = C_{q_{21}} + \frac {2 C_{{q}_{11}} C_{\sigma _{21}} }{5} \sinh \left ( A \right ) - \frac {6 \sqrt {5} C_{{q}_{11}} \bar {G} Kn }{25} \cosh \left ( A \right ) + \frac {2 C_{{q}_{11}} \bar {G} \bar {y} }{5} \sinh \left ( A \right ) \nonumber \\ &\quad +\, \frac {2 C_{{q}_{12}} C_{\sigma _{21}} }{5} \cosh \left ( A \right ) - \frac {6 \sqrt {5} C_{{q}_{12}} \bar {G} Kn }{25}\sinh \left ( A \right ) + \frac {2 C_{{q}_{12}} \bar {G} \bar {y} }{5} \cosh \left ( A \right ) + \frac {C_{\sigma _{21}}^{2} }{Kn}\bar {y} \nonumber \\ &\quad +\, \frac {C_{\sigma _{21}} \bar {G}}{Kn} \bar {y}^{2} + \frac {\bar {G}^{2} }{3 Kn} \bar {y}^{3}, \end{align}

where $C_{q_{21}}$ is an integration constant. Thereafter, integrating (7.5a ) and (7.5c ) individually yields the solution for the normal stresses ( ${\bar {\sigma }}_{1 1 }$ , and ${\bar {\sigma }}_{22}$ ) as follows:

(7.11a) \begin{align} {{\bar {\sigma }}_{1 1 }} & = C_{\sigma _{111}} \sinh \left ( C \right ) + C_{\sigma _{112}} \cosh \left ( C \right ) + \frac {8 C_{\sigma _{21}}^{2}}{5} \nonumber \\ & \quad- \frac {912 \sqrt {5} C_{{q}_{11}} C_{\sigma _{21}} }{1199}\cosh \left ( A \right ) + \frac {10580600056 C_{{q}_{11}} \bar {G} Kn }{107820075} \sinh \left ( A \right ) \nonumber \\ &\quad -\frac {912 \sqrt {5} C_{{q}_{11}} \bar {G} } {1199} \bar {y} \cosh \left ( A \right ) - \frac {912 \sqrt {5} C_{{q}_{12}} C_{\sigma _{21}}}{1199} \sinh \left ( A \right ) \nonumber \\ &\quad + \frac {10580600056 C_{{q}_{12}} \bar {G} Kn }{107820075} \cosh \left ( A \right ) - \frac {912 \sqrt {5} C_{{q}_{12}} \bar {G} }{1199} \bar {y} \sinh \left ( A \right ) \nonumber \\ &\quad + \frac {16 C_{\sigma _{21}} \bar {G} }{5} \bar {y}- \frac {C_{\sigma _{221}} }{2} \sinh {\left (B \right )} - \frac {C_{\sigma _{222}} }{2} \cosh {\left (B \right )} + \frac {9784 \bar {G}^{2} Kn^{2}}{1875} + \frac {8 \bar {G}^{2} }{5}\bar {y}^{2}, \end{align}

(7.11b) \begin{align} {{\bar {\sigma }}_{2 2 }} & = C_{\sigma _{221}} \sinh \left ( B \right ) + C_{\sigma _{222}} \cosh \left ( B \right ) - \frac {3251 \bar {G}^{2} Kn^{2}}{625} + \frac {24 \sqrt {5} C_{{q}_{11}} C_{\sigma _{21}} }{11} \cosh \left ( A \right ) \nonumber \\ & \quad- \frac {1723268 C_{{q}_{11}} \bar {G} Kn }{9075} \sinh \left ( A \right ) + \frac {24 \sqrt {5} C_{{q}_{11}} \bar {G} }{11} \bar {y} \cosh \left ( A \right ) + \frac {24 \sqrt {5} C_{{q}_{12}} C_{\sigma _{21}} }{11}\sinh \left ( A \right ) \nonumber \\ &\quad - \frac {1723268 C_{{q}_{12}} \bar {G} Kn }{9075} \cosh \left ( A \right ) + \frac {24 \sqrt {5} C_{{q}_{12}} \bar {G} }{11} \bar {y} \sinh \left ( A \right ) - \frac {6 C_{\sigma _{21}}^{2}}{5} - \frac {12 C_{\sigma _{21}} \bar {G} }{5} \bar {y} \nonumber \\ & \quad- \frac {6 \bar {G}^{2} }{5}\bar {y}^{2}, \\[9pt] \nonumber \end{align}

where $B ={({5 \sqrt {4830} \bar {y}}/{483 Kn} )}$ , $C = {({5 \sqrt {966} \bar {y}}/{161 Kn} )}$ , while $ C_{\sigma _{ 111}}$ , $ C_{\sigma _{ 112}}$ , $ C_{\sigma _{ 221}}$ and $ C_{\sigma _{ 222}}$ are integration constants. Finally, thermodynamic variables ( $T$ , and $p$ ) are obtained after integrating (7.6b ) and (7.3b ) individually as

(7.12a) \begin{align} \bar {T} &= C_{T} + \frac {4 C_{{q}_{11}}^{2} }{125}\cosh ^{2}\left ( A \right ) + \frac {8 C_{{q}_{11}} C_{{q}_{12}} }{125} \sinh \left ( A \right ) \cosh \left ( A \right ) - \frac {9623 \sqrt {5} C_{{q}_{11}} C_{\sigma _{21}} }{12375} \cosh \left ( A \right ) \nonumber \\ &\quad + \frac {85851704 C_{{q}_{11}} \bar {G} Kn }{1134375} \sinh ( A ) - \frac {9623 \sqrt {5} C_{{q}_{11}} \bar {G} }{12375} \bar {y} \cosh ( A ) + \frac {4 C_{{q}_{12}}^{2} }{125} \cosh ^{2} ( A ) \nonumber \\ &\quad - \frac {9623 \sqrt {5} C_{{q}_{12}} C_{\sigma _{21}} }{12375}\! \sinh ( A ) \!+\! \frac {85851704 C_{{q}_{12}} \bar {G} Kn }{1134375}\! \cosh ( A ) \!-\! \frac {9623 \sqrt {5} C_{{q}_{12}}\! \bar {G} }{12375} \bar {y} \sinh ( A ) \nonumber \\ &\quad - \frac {2 C_{\sigma _{21}}^{2} }{15 Kn^{2}} \bar {y}^{2} + \frac {3299 C_{\sigma _{21}} \bar {G} }{3750} \bar {y} - \frac {4 C_{\sigma _{21}} \bar {G} }{45 Kn^{2}} \bar {y}^{3} - \frac {2 C_{\sigma _{221}} }{5} \sinh ( B ) - \frac {2 C_{\sigma _{222}} }{5} \cosh ( B ) \nonumber \\ &\quad - \frac {4 C_{q 21} }{15 Kn} \bar {y} + \frac {3299 \bar {G}^{2} }{7500} \bar {y}^{2} - \frac {\bar {G}^{2} }{45 Kn^{2}}\bar {y}^{4}, \end{align}

(7.12b) \begin{align} \bar {p} & = C_{p} \!-\! \frac {24 \sqrt {5} C_{{q}_{11}} C_{\sigma _{21}} }{11}\cosh \left ( A \right ) \!+\! \frac {1723268 C_{{q}_{11}} \bar {G} Kn }{9075} \sinh \left ( A \right ) \!-\! \frac {24 \sqrt {5} C_{{q}_{11}} \bar {G} }{11}\bar {y} \cosh \left ( A \right ) \nonumber \\ &\quad - \frac {24 \sqrt {5} C_{{q}_{12}} C_{\sigma _{21}} }{11} \sinh \left ( A \right ) + \frac {1723268 C_{{q}_{12}} \bar {G} Kn }{9075} \cosh \left ( A \right ) - \frac {24 \sqrt {5} C_{{q}_{12}} \bar {G}}{11} \bar {y} \sinh \left ( A \right ) \nonumber \\ & \quad + \frac {6 C_{\sigma _{21}}^{2}}{5} + \frac {12 C_{\sigma _{21}} \bar {G} }{5} \bar {y}- C_{\sigma _{221}} \sinh \left ( B \right ) - C_{\sigma _{222}} \cosh \left ( B \right ) + \frac {3251 \bar {G}^{2} Kn^{2}}{625} + \frac {6 \bar {G}^{2} }{5}\bar {y}^{2}, \end{align}

where $C_T$ and $C_p$ are integration constants.

It is important to highlight that the N–S equations predict zero streamwise heat flux, while the G13 equations predict a constant streamwise heat flux. Note that the velocity ( $\bar {u}_1$ ), heat flux ( $\bar {q}_1,\,\bar {q}_2$ ) and shear stress component ( $\bar {\sigma }_{21}$ ) explicitly depend on the forcing term ( $\bar {G}$ ), whereas other physical quantities ( $\bar {T},\,\bar {p}$ ) and the normal stress components ( $\bar {\sigma }_{11},\,\bar {\sigma }_{22}$ ) become non-zero due to their coupling with the previously mentioned parameters.

Note that the expression for density ( $\bar {\rho }$ ) is derived from pressure and temperature, as it is a dependent variable. Hence, it is not presented here for brevity. Within this context, it is crucial to note that the integration constants $C_{\sigma _{21}}$ , $C_{q_{11}}$ , $C_{\sigma _{ 221}}$ and $C_{q_ {2 1}}$ , which are present in (7.7)–(7.10), become zero due to the symmetry displayed by the streamwise velocity and temperature. Incorporating these results into (7.10)–(7.12b ) results in the following analytical solution of various quantities for the problem:

(7.13a) \begin{align} \bar {\sigma }_{21} &= \underline {\bar {G}\bar {y}}, \end{align}

(7.13b) \begin{align} \bar {q}_1 &= C_{q_{12}} \cosh \left ( A \right ) - {\frac {3 \bar {G} Kn}{2}}, \end{align}

(7.13c) \begin{align} {{\bar {u}}_{1}} &= \underline {C_{u_{1}} - \frac {\bar {G} \bar {y}^{2}}{2 Kn}} - \frac {2 C_{q_{12}} }{5}\cosh \left ( A \right ) + \frac {3 \bar {G} Kn}{5}, \end{align}

(7.13d) \begin{align} {{\bar {q}}_{2}} &= - \frac {6 \sqrt {5} C_{q_{12}} \bar {G} Kn }{25} \sinh \left ( A \right ) + \frac {2 C_{q_{12}} \bar {G} }{5} \bar {y} \cosh \left ( A \right ) + \underline {\frac {\bar {G}^{2} }{3 Kn}\bar {y}^{3}}, \end{align}

(7.13e) \begin{align} {{\bar {\sigma }}_{1 1 }} &= C_{\sigma _{111}} \sinh {\left (C \right )} + C_{\sigma _{112}} \cosh {\left (C \right )} - \frac {C_{\sigma _{221}} }{2} \sinh {\left (B \right )} - \frac {912 \sqrt {5} C_{{q}_{12}} \bar {G} }{1199} \bar {y} \sinh {\left (A \right )} \nonumber \\&\quad + \frac {10580600056 C_{{q}_{12}} \bar {G} Kn }{107820075} \cosh {\left (A \right )} - \frac {C_{\sigma _{222}} }{2} \cosh {\left (B \right )} + \frac {9784 \bar {G}^{2} Kn^{2}}{1875} + \frac {8 \bar {G}^{2} }{5}\bar {y}^{2}, \end{align}

(7.13f) \begin{align} {{\bar {\sigma }}_{2 2 }} &=- \frac {1723268 C_{{q}_{12}} \bar {G} Kn }{9075}\cosh ( A ) + \frac {24 \sqrt {5} C_{{q}_{12}} \bar {G} }{11} \bar {y} \sinh ( A ) \nonumber \\&\quad + C_{\sigma _{221}} \sinh ( B ) + C_{\sigma _{222}} \cosh ( B ) - \frac {3251 \bar {G}^{2} Kn^{2}}{625} \!-\! \frac {6 \bar {G}^{2} }{5}\bar {y}^{2}, \end{align}

(7.13g) \begin{align} \bar {T} &= \underline { C_{T} \!-\! \frac {\bar {G}^{2} }{45 Kn^{2}} \bar {y}^{4} } \!+\! \frac {3299 \bar {G}^{2} }{7500} \bar {y}^{2} \!+\! \frac {4 C_{{q}_{12}}^{2} }{125}\cosh ^{2} ( A ) + \frac {85851704 C_{{q}_{12}} \bar {G} Kn }{1134375} \cosh ( A ) \nonumber \\&\quad - \frac {9623 \sqrt {5} C_{{q}_{12}} \bar {G} }{12375} \bar {y} \sinh ( A ) - \frac {4 C_{q 21} }{15 Kn} \bar {y} - \frac {2 C_{\sigma _{221}} }{5}\sinh ( B ) - \frac {2 C_{\sigma _{222}} }{5} \cosh ( B ) ,\end{align}

(7.13h) \begin{align} \bar {p} &= \underline {Cp} + \frac {1723268 C_{{q}_{12}} \bar {G} Kn }{9075} \cosh \left ( A \right ) - \frac {24 \sqrt {5} C_{{q}_{12}} \bar {G} }{11}\bar {y} \sinh \left ( A \right ) - C_{\sigma _{221}} \sinh \left ( B \right ) \nonumber \\&\quad - C_{\sigma _{222}} \cosh \left ( B \right ) + \frac {3251 \bar {G}^{2} Kn^{2}}{625} + \frac {6 \bar {G}^{2} }{5} \bar {y}^{2} .\end{align}

The underlined terms in (7.13) represent the specific components that arise from the solution derived by employing the N–S equations for the present problem. We also compare our solution with that of the G13 equations in § 8. All the other terms, including the hyperbolic cosine and sine, are the result of the Knudsen layer responsible for rarefaction effects. The solutions obtained for the force-driven plane Poiseuille flow problem in (7.13) will be compared and validated against other model solutions and DSMC data in the subsequent subsection.

7.2. Validation

In this subsection, we compare the solution obtained from the simplified form of SO13 equations (7.13) with the N–S, G13, R13 equations and DSMC data reported in Zheng et al. (Reference Zheng, Garcia and Alder2002, Reference Zheng, Garcia and Alder2003). This comparison requires evaluating the integration constants in the solution using accurate boundary conditions, which are crucial for determining these constants. However, note that extended models beyond the Navier–Stokes–Fourier framework require additional boundary conditions: the Burnett equations necessitate higher-order gradients, while Grad-type methods involve prescribing conditions for higher-order moments. Since ill-posed boundary conditions can lead to unphysical results, their careful formulation is essential. However, deriving precise boundary conditions remains an active research area for most HOT equations. It should be noted that deriving boundary conditions for the proposed equations is beyond the scope of this work. While we plan to address this issue in future work, as done in Rana & Struchtrup (Reference Rana and Struchtrup2016) and Rana et al. (Reference Rana, Gupta, Sprittles and Torrilhon2021), we choose to proceed without being hindered by the lack of boundary conditions. Therefore, we adopt the approach suggested in Uribe & Garcia (Reference Uribe and Garcia1999), García-Colín et al. (Reference García-Colín, Velasco and Uribe2008), Yadav & Agrawal (Reference Yadav and Agrawal2024) and Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2024). Accordingly, we use data from DSMC simulations to evaluate the integration constants and thereby obtain a complete solution for all the variables.

In the simplified solution of (7.13), five integration constants require five boundary conditions for the corresponding variables at one wall, as listed in table 1, while symmetry determines those at the other. It implies that the integration constants are evaluated by matching the analytical solutions with DSMC data at the wall. Note that in classical hydrodynamics, where Knudsen boundary layers are absent, integration constants linked to it vanish.

Table 1. Integration coefficients for various models based on DSMC data at the wall.

As discussed above, we utilise the DSMC data reported in Zheng et al. (Reference Zheng, Garcia and Alder2002, Reference Zheng, Garcia and Alder2003) to evaluate the boundary conditions at $\bar {y} = \pm 0.5$ for all the primary variables. In a similar manner, we also apply this step to every result presented in the current work for other equations, namely, the N–S, G13 and R13 equations, to make the comparisons on a common ground. Note that the DSMC simulation was conducted for a non-dimensional force of $\bar {G} = 0.2355$ and a Knudsen number of $\textit{Kn}=0.072$ , with the Knudsen number adjusted to match the definition used in the present study. In addition to these non-dimensional parameters, the boundary values of the primary variables obtained from the DSMC data are employed to evaluate the integration constants and their values are summarised in table 1. These values are subsequently used to obtain their profiles across the channel in the following section.

Figure 3 presents a comprehensive comparison of the results from the derived analytical solution of the proposed SO13, other (N–S, G13, R13) equations and DSMC data for the problem considered under the above-mentioned conditions. As discussed above, the remaining integration constants ( $C_{{q}_{12}}$ and $C_{u_{1}}$ ), which arise after applying the symmetry conditions in (7.7)–(7.13d ), are evaluated using the available DSMC data.

Figure 3. Cross-stream variation of (a) shear stress ( $\bar {\sigma }_{21}$ ), (b) streamwise heat flux ( $\bar {q}_1$ ), (c) streamwise velocity ( $\bar {u}_1$ ) and (d) cross-stream heat flux ( $\bar {q}_2$ ). The solution is compared with the corresponding results from the N–S, G13 and R13 equations and the DSMC data (reported by Zheng et al. Reference Zheng, Garcia and Alder2002, Reference Zheng, Garcia and Alder2003) for ${\textit{Kn}} = 0.072$ and $\bar {G} = 0.2355$ .

Figure 4. Cross-stream variation of (a) normal stress ( $\bar {\sigma }_{22}$ ), (b) pressure ( $\bar {p}$ ), (c) temperature ( $\bar {T}$ ) and (d) density ( $\rho$ ). The solution is also compared with the corresponding results from the N–S, G13 and R13 equations and the DSMC data (reported by Zheng et al. Reference Zheng, Garcia and Alder2002, Reference Zheng, Garcia and Alder2003) for ${\textit{Kn}} = 0.072$ and $\bar {G} = 0.2355$ .

Figure 3(a) focuses on the variation of shear stress ( $\bar {\sigma }_{21}$ ), which exhibits a linear profile. Notably, the profiles derived from the SO13 and other equations are indistinguishable and show good agreement with the DSMC data. Similarly, the streamwise heat fluxes ( $\bar {q}_1$ ) obtained from the SO13 and R13 equations overlap (figure 3 b) and align well with the DSMC data, in which the profile transitions from negative in the bulk region to positive near the wall boundaries. The change in sign indicates that the tangential heat flux is in the opposite direction to the flow direction in the bulk region. Also, we observe steep variation near the walls while the heat flux remains constant in the bulk region. The solution of the N–S equations gives zero tangential heat flux, while the G13 equations yield a constant value, which matches the solution of the SO13 equations at the centre of the microchannel. For streamwise velocity ( $\bar {u}_1$ ) and cross-stream heat flux ( $\bar {q}_2$ ), as shown in figures 3(c) and 3(d), respectively, the results of all equations set are indistinguishable and consistent with the DSMC data.

In figure 4, we compare the analytical solutions for normal stress ( $\bar {\sigma }_{22}$ ), pressure ( $\bar {p}$ ), temperature ( $\bar {T}$ ) and density ( $\bar {\rho }$ ) as obtained using the SO13 equations with those of the N–S, G13, R13 equations and benchmark with the DSMC data. The variation of $\bar {\sigma }_{22}$ across the microchannel (figure 4 a) is better captured with the SO13 equations as compared with the N–S, G13 and R13 solutions. In particular, the N–S equations essentially give zero stresses, while the G13 equations fail to capture the bi-modal profile. The non-uniform pressure profile across the microchannel is accurately captured with the SO13 equations (figure 4 b), with the N–S equations predicting a constant pressure profile.

The results for temperature are particularly interesting, with the DSMC simulations showing a characteristic temperature dip at the centre. As shown in figure 4(c), the N–S equations fail to capture this dip while the results of the SO13 equations are in better agreement with the DSMC data as compared with those of the G13 and R13 equations. Figure 4(d) compares the normalised density across the channel height. The results of the SO13 equations match closely with the DSMC data, thereby again establishing the superiority of the SO13 equations.