2.1. Problem statement

Figure 1.

Schematic of the problem. A square cavity of side length $L^*$ is filled with a perfect monatomic gas. The horizontal and vertical walls are maintained at temperatures $T_{h}^*$ and $T_{c}^*$ , respectively, where $T_{h}^*\gt T_{c}^*$ . Each wall is labelled a number, and the indicated angles $\theta _i$ ( $i=1,\ldots ,4$ ) correspond to the angles introduced in (3.11).

Figure 1. Long description

A schematic of a square cavity with side length L*. The cavity is filled with a perfect monatomic gas. The horizontal walls are maintained at temperature T_h, and the vertical walls are maintained at temperature T_c, where T_h is greater than T_c. Each wall is labeled with a number: the top wall is labeled 1, the bottom wall is labeled 2, the right wall is labeled 3, and the left wall is labeled 4. The angles θ1, θ2, θ3, and θ4 are indicated at the intersection of the dashed lines, which represent the directions of gas flow due to temperature differences. The x* and y* axes are shown at the center of the cavity, indicating the coordinate system used for analysis.

A schematic of the problem is shown in figure 1. We consider a perfect monatomic gas enclosed in a two-dimensional square cavity of side length $L^*$ (hereafter, asterisks denote dimensional quantities). The total gas mass contained in the cavity per unit length (in its spanwise direction) is $m^*_{\textit{tot}}$ . No external body forces, including gravity, act on the gas. The top and bottom cavity walls are maintained at a constant temperature $T_{h}^*$ , and are labelled as walls 1 and 3, respectively, and the left and right boundaries are fixed at temperature $T_{c}^*$ and labelled wall numbers 2 and 4, respectively. Without loss of generality, we assume $T_{h}^* \gt T_{c}^*$ . The origin of the $(x^*,y^*)$ coordinate system is located at the cavity centre point, and we seek to analyse the steady flow field generated in the cavity at arbitrary gas rarefaction rates.

In the framework of gas kinetic theory and the present two-dimensional steady-flow setup, the gas state is governed by the probability density function $f^* = f^*(x^*,y^*,\boldsymbol{\xi }^*)$ of finding a molecule at a position $\boldsymbol{r^*}=(x^*,y^*)$ with molecular velocity $\boldsymbol{\xi }^*$ . To analyse the velocity space, a cylindrical coordinate system is adopted, where $\boldsymbol{\xi }^* = (\xi _r^*,\theta _\xi ,\xi _z^*)$ . Here, $\xi _r^* \in [0,\infty )$ denotes the molecular speed in the plane of motion, $\theta _\xi$ marks the in-plane counterclockwise velocity direction and $\xi _z^* \in (-\infty ,\infty )$ is the velocity component normal to the cavity. To model gas–surface interactions between the gas particles and the solid walls, we apply the diffuse Maxwell boundary condition (Sone Reference Sone2007),

(2.1) \begin{equation} f^*({\boldsymbol{\widetilde {r}}_{i}^*},\boldsymbol{\xi }^*\boldsymbol{\cdot }\hat{\boldsymbol{n}}\gt 0)=\frac {\rho ^*_i(\boldsymbol{\widetilde {r}}_i^*)}{(\pi 2\mathcal{R}^*_gT^*_i)^{3/2}}\exp {\left [-\frac {{\xi ^*}^2}{2\mathcal{R}^*_gT^*_i}\right ]}, \end{equation}

stating that the reflected particles achieve thermal equilibrium with their reflecting boundaries. In (2.1), $\boldsymbol{\widetilde {r}}_i^*$ marks some point along the $i$ th wall $(i=1,\ldots ,4)$ , $\hat{\boldsymbol{n}}$ is a unit vector directed normal to the wall and into the gas domain and $\rho _i^*(\boldsymbol{\widetilde {r}}_i^*)$ are yet unknown functions associated with the mass flux of particles reflected from the $i$ th wall. In addition, $U_{{mp}_{i}}=(2\mathcal{R}_g^* T_i^*)^{1/2}$ denotes the most probable speed of gas particles at the $i$ th wall temperature, where $\mathcal{R}_g^*$ is the specific gas constant and $T_i^*$ is the temperature of the $i$ th wall. In line with figure 1, the temperatures of the walls are $T_{1}^*=T_{3}^* = T_{h}^*$ and $T_{2}^*=T_{4}^* = T_{c}^*$ .

To render the problem dimensionless, we scale the position by the cavity side length $L^*$ , the velocity by the $T_h^*$ -based wall temperature most probable molecular speed, $U_{{mp}_{h}}^*=(2\mathcal{R}^*_g T^*_{h})^{1/2}$ , and the temperature by $T^*_{h}$ . The density is normalised by the average density based on the total mass of the gas confined in the cavity, $\rho _0^*=m^*_{\textit{tot}}/L^{*2}$ , and the pressure is non-dimensionalised by $\rho _0^*\mathcal{R}^*_gT_h^*$ . The system non-dimensional description is then governed by the cold-to-hot wall temperature ratio and the gas mean Knudsen number,

(2.2) \begin{equation} R_T=T_c^*/T_h^* \quad \mathrm{and} \quad \textit{Kn} = \lambda ^*/L^*, \end{equation}

respectively, where $\lambda ^*$ denotes the mean free path of a gas molecule. Assuming an HS interaction gas model, $\lambda ^* = m^*/(\pi \sqrt {2}\rho _0^* d^{*2})$ , where $m^*$ and $d^*$ denote the mass and diameter of a gas particle, respectively (Sone Reference Sone2007). Substituting $\lambda ^*$ into (2.2), we obtain

(2.3) \begin{equation} \textit{Kn} = \frac {m^*}{\pi \sqrt {2}\rho ^*_0 d^{*2}L^*} \ . \end{equation}

The problem statement is completed by imposing a normalisation condition for the gas density $\rho$ , fixing the total gas mass in the cavity,

(2.4) \begin{equation} \int _{-1/2}^{1/2}\int _{-1/2}^{1/2}\rho (x,y)\,\mathrm{d}x\,\mathrm{d}y=1 . \end{equation}

In what follows, we introduce the kinetic models to be used to formulate the collision operator in the Boltzmann equation.

2.2. The model equations

Adopting the scaling introduced in § 2, and within the two-dimensional steady setup considered, the scaled velocity distribution function

(2.5) \begin{equation} f = \frac {(2\mathcal{R}_g^* T_{h}^*)^{3/2}}{\rho _0^*} f^* \end{equation}

is governed by the non-dimensional Boltzmann-type equation

(2.6) \begin{equation} \boldsymbol{\xi }\boldsymbol{\cdot }\boldsymbol{\nabla }\!f= \frac {1}{\textit{Kn}}\mathcal{J}(f), \end{equation}

where $\mathcal{J}(f)$ denotes the nonlinear integral collision operator, accounting for the effect of binary molecular collisions on the evolution of the distribution function. For HS gas, the corresponding scaled collision operator is given by (Kogan 1969).

(2.7) \begin{equation} \mathcal{J}_{\textit{HS}}(f) = \frac {1}{2\sqrt 2\pi }\int \left (f(\boldsymbol{r},\boldsymbol{\xi '})f(\boldsymbol{r},\boldsymbol{\xi '_1})-f(\boldsymbol{r},\boldsymbol{\xi })f(\boldsymbol{r},\boldsymbol{\xi _1})\right ) |\boldsymbol{g}\boldsymbol{\cdot }\hat {\boldsymbol{k}}|\,\mathrm{d}\varOmega \,\mathrm{d\boldsymbol{\xi _1}} , \end{equation}

where $\mathrm{d}\varOmega$ denotes the solid-angle element associated with its unit vector $\hat {\boldsymbol{k}}$ , $\boldsymbol{\xi }_1$ is the integration variable corresponding to the molecular velocity, and $\boldsymbol{g} = \boldsymbol{\xi }_1 - \boldsymbol{\xi }$ is the relative velocity between particles prior to collision. The post-collision molecular velocities are

(2.8) \begin{equation} \boldsymbol{\xi }' = \boldsymbol{\xi } + \hat{\boldsymbol{k}} (\hat{\boldsymbol{k}} \boldsymbol{\cdot }\boldsymbol{g}) \ \ \ \mathrm{and} \ \ \ \boldsymbol{\xi }_1' = \boldsymbol{\xi }_1 - \hat{\boldsymbol{k}} (\hat{\boldsymbol{k}} \boldsymbol{\cdot }\boldsymbol{g}). \end{equation}

In addition to the above, we employ the BGK (Bhatnagar et al. Reference Bhatnagar, Gross and Krook1954) and Shakhov (Reference Shakhov1968) collision models, which simplify the analysis. In the former,

(2.9) \begin{equation} \mathcal{J}_{\textit{BGK}}(f) = \nu (x,y)\,(f_M - f), \end{equation}

whereas in the latter

(2.10) \begin{equation} \mathcal{J}_{\textit{S}}(f) = \nu (x,y) \left [ f_M \left ( 1 + \frac {4}{15 \rho T^2} \, \boldsymbol{q} \boldsymbol{\cdot }(\boldsymbol{\xi } - \boldsymbol{u}) \left ( \frac {(\boldsymbol{\xi }-\boldsymbol{u})^2}{T} - \frac {5}{2} \right ) \right ) - f \right ]. \end{equation}

In (2.9) and (2.10),

(2.11) \begin{equation} f_M = \frac {\rho }{(\pi T)^{3/2}}\exp \left [-\frac {(\boldsymbol{\xi }-\boldsymbol{u})^2}{T}\right ] \ \ \ \mathrm{and} \ \ \ \nu (x,y) = \frac {\rho \sqrt {\pi T}}{2} \end{equation}

denote the local Maxwellian distribution and the non-dimensional local variation of the molecular collision frequency, respectively. In addition, $\boldsymbol{u}=(u_x,u_y)$ and $\boldsymbol{q}=(q_x,q_y)$ mark the macroscopic two-dimensional velocity and heat-flux vectors, respectively. In terms of velocity-space quadratures,

(2.12) \begin{eqnarray} &\rho & = \int f\,\mathrm{d}\boldsymbol{\xi } \ \ ,\ \ u_x = \frac {1}{\rho }\int \xi _r \cos {\theta _\xi } f \,\mathrm{d}\boldsymbol{\xi } \ \ , \ \ u_y = \frac {1}{\rho }\int \xi _r\sin {\theta _\xi } f\, \mathrm{d}\boldsymbol{\xi } \ ,\nonumber \\ &T& =\frac {2}{3\rho } \left (\int (\xi _r ^2+\xi _z^2)f\,\mathrm{d}\boldsymbol{\xi }-\rho (u_x^2+u_y^2)\right ) , \nonumber \\ &q_x&=\int \left [(\xi _r \cos {\theta _\xi }-u_x)^2+(\xi _r\sin {\theta _\xi }-u_y)^2+\xi _z^2\right ](\xi _r \cos {\theta _\xi }-u_x)f\,\mathrm{d}\boldsymbol{\xi } \nonumber \\ &\mathrm{and}& \nonumber \\ &q_y&=\int \left [(\xi _r \cos {\theta _\xi }-u_x)^2+(\xi _r\sin {\theta _\xi }-u_y)^2+\xi _z^2\right ](\xi _r \sin {\theta _\xi }-u_y)f\,\mathrm{d}\boldsymbol{\xi }, \end{eqnarray}

mark the hydrodynamic density, velocity components, temperature and heat flux components, respectively, where $\mathrm{d}\boldsymbol{\xi }=\xi _r\,\mathrm{d}\xi _r\,\mathrm{d}\theta _\xi \,\mathrm{d}\xi _z$ . Importantly, the Shakhov model surpasses the BGK kernel by yielding the correct continuum-limit $2/3$ Prandtl number for a monatomic gas, rather than the unity value obtained in the latter. It is therefore considered more reliable than the BGK model, particularly in problems coupling thermal and dynamical effects, as in the present problem. We further discuss the differences between the models’ predictions in § 6.

3.1. The free-molecular limit

The free-molecular distribution function $f^{(0)}$ satisfies (3.4), indicating that it remains constant along free-flight trajectories. Consequently, each particle carries a distribution function fixed by its last wall collision,

(3.7) \begin{equation} f^{(0)}(\boldsymbol{r},\boldsymbol{\xi })=\frac {\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)}{(\pi T_i)^{3/2}}\exp {\left [-\frac {\xi ^2}{T_i}\right ]} , \end{equation}

where $\boldsymbol{\widetilde {r}}_i$ follows the definition after (3.6). To obtain $\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)$ , we impose the no-penetration condition along the cavity walls. Applying (2.12), the free-molecular mass flux vector at the boundaries is

(3.8) \begin{eqnarray} (\rho ^{(0)} u_x^{(0)},\rho ^{(0)} u_y^{(0)})= \int \frac {\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)\sqrt {T_i}}{4\sqrt {\pi }}(\cos \theta _\xi ,\sin \theta _\xi )\,\mathrm{d}\theta _\xi . \end{eqnarray}

Taking

(3.9) \begin{equation} \rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i) = {A}/{\sqrt {T_i}} \ , \end{equation}

the no-penetration condition is satisfied at each point on the walls. In (3.9), $A$ is a normalisation constant determined by the density normalisation condition (cf. (2.4)). At free-molecular conditions, we denote this constant by $A^{(0)}$ , and its value is fixed through

(3.10) \begin{equation} \int _{-1/2}^{1/2}\int _{-1/2}^{1/2}\rho ^{(0)}(x,y)\,\mathrm{d}x\,\mathrm{d}y=1 , \end{equation}

such that the total mass inside the cavity is conserved. To obtain an explicit closed-form expression for $f^{(0)}$ , we introduce the angles $\theta _i$ ( $i=1,\ldots ,4$ ) separating particles arriving at each point from the different walls, as illustrated in figure 1. Following geometrical considerations, these are given by

(3.11) \begin{eqnarray} \nonumber \theta _1 &=\arctan \left (\frac {1/2-y}{1/2-x} \right ) ,\quad \theta _2=\arctan \left (\frac {1/2-y}{1/2+x} \right ) , \\ \theta _3 &=\arctan \left (\frac {1/2+y}{1/2+x} \right ) \quad \mathrm{and} \quad \theta _4=\arctan \left (\frac {1/2+y}{1/2-x} \right ). \end{eqnarray}

Combining (3.11) with (3.7) and (3.9), we obtain an explicit piecewise expression for $f^{(0)}$ ,

(3.12) \begin{equation} f^{(0)}(\boldsymbol{r},\boldsymbol{\xi })=\frac {A^{(0)}}{\pi ^{3/2} T_i^{2}}\exp {\left [-\frac {\xi ^2}{T_i}\right ]} , \quad i=\begin{cases} 1, &\pi +\theta _1\lt \theta _\xi \lt 2\pi -\theta _2 \\ 2,&2\pi -\theta _2\lt \theta _\xi \lt 2\pi +\theta _3 \\ 3,&2\pi +\theta _3\lt \theta _\xi \lt 3\pi -\theta _4 \\ 4,&3\pi -\theta _4\lt \theta _\xi \lt 3\pi +\theta _1. \end{cases} \end{equation}

The free-molecular hydrodynamic fields may now be calculated using (2.12), yielding

(3.13) \begin{eqnarray} \nonumber &\rho ^{(0)}(x,y)& \,= \frac {A^{(0)}}{2\pi }\left (\frac {1}{\sqrt {R_T}}-1\right )\left (\theta _1+\theta _2+\theta _3+\theta _4\right )+A^{(0)} ,\quad u_x^{(0)} = 0 , \quad u_y^{(0)} = 0, \\ \nonumber &T^{(0)}(x,y)& \,= \sqrt {R_T}\left ( \frac {2\pi (1+\sqrt {R_T})}{(1-\sqrt {R_T})(\theta _1+\theta _2+\theta _3+\theta _4)+2\pi \sqrt {R_T}}-1\right ), \\ \nonumber &q_x^{(0)}(x,y)&\,= \frac {A^{(0)}(1-R_T)}{2\sqrt {\pi }}\left ( \sin \theta _1-\sin \theta _2-\sin \theta _3+\sin \theta _4\right ) \quad \mathrm{and} \\ &q_y^{(0)}(x,y)& \, = \frac {A^{(0)}(1-R_T)}{2\sqrt {\pi }}\left (-\cos \theta _1-\cos \theta _2+\cos \theta _3+\cos \theta _4\right ). \end{eqnarray}

In accordance with the gas scaled equation of state, the pressure is given by $p^{(0)} = \rho ^{(0)} T^{(0)}$ . Using the normalisation condition in (3.10) along with (3.11) and (3.13), we obtain

(3.14) \begin{equation} A^{(0)}= {2\sqrt {R_T}}/{\left (1+\sqrt {R_T}\right )} , \end{equation}

which fixes the free-molecular solution.

3.2. The near-free-molecular correction

The solution for (3.5), satisfying the $O(\textit{Kn}^{-1})$ part of the boundary condition in (3.6), is

(3.15) \begin{equation} f^{(1)}(x,y,\boldsymbol{\xi })=\frac {\rho _i^{(1)}(\boldsymbol{\widetilde {r}}_i)}{(\pi T_i)^{3/2}}\exp {\left [-\frac {\xi ^2}{T_i}\right ]}+\frac {1}{\xi _r}\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\mathcal{J}(f^{(0)})\,\mathrm{d}l. \end{equation}

Here, $\rho _i^{(1)}(\boldsymbol{\widetilde {r}}_i)$ represent the corrections for the wall flux functions due to the effect of molecular collisions, to be determined through the impermeability condition. The straight-line integral in (3.15) is taken along the corresponding free-flight trajectory. It starts at $\boldsymbol{\widetilde {r}}_i= (x_i,y_i)$ , the location of recent wall reflection, and terminates at $\boldsymbol{r}=(x,y)$ . Notably, while $f^{(1)}$ becomes singular as $\xi _r\rightarrow 0$ , reflecting the breakdown of near-free-molecular conditions for slow particles, its hydrodynamic moments remain finite. Specifically, in the present two-dimensional setting, the velocity-space integrals are convergent due to the inverse-square-root form of the singularity, $\xi _r^{-1} = (\xi _x^2 + \xi _y^2)^{-1/2}$ (Smolderen Reference Smolderen1965; Sone Reference Sone2007). Logarithmic divergences in $\textit{Kn}$ have been reported in one-dimensional or unbounded domains (infinite or semi-infinite), yet such configurations are not considered here. More specifically, the convergence of the hydrodynamic fields in the present two-dimensional configuration may be justified using the cylindrical-space representation of the molecular velocity vector. Here, the integration measure

(3.16) \begin{equation} \mathrm{d}\boldsymbol{\xi } = \xi _r\,\mathrm{d}\xi _r\,\mathrm{d}\theta _\xi\, \mathrm{d}\xi _z \end{equation}

cancels out the $1/\xi _r$ singularity appearing in $f^{(1)}$ when evaluating the hydrodynamic density. Moreover, higher-order moments involve additional powers of $\xi _r$ , which further weaken the singular behaviour and ensure that all quadratures converge.

In principle, the asymptotic procedure may be followed to derive higher-order corrections. The second correction, $f^{(2)}$ , is expected to exhibit a stronger singularity, scaled as $\xi _r^{-2}$ , which is no longer integrable for the lower-order hydrodynamic moments, such as the density (Smolderen 1965). Consequently, the asymptotic expansion ceases to be valid for these quantities. Higher-order moments, however, involve higher powers of $\xi _r$ in the velocity-space integrals, and therefore remain convergent. Accordingly, the inclusion of $f^{(2)}$ may provide additional information on quantities such as the mass flux at moderately large Knudsen numbers, although its derivation and implementation are considerably more involved and are not followed here.

Parametrising the integral trajectory in (3.15) in terms of $\theta _\xi$ , we obtain

(3.17) \begin{align} (\widetilde {x}_1,\widetilde {y}_1) &= (x + (1/2 - y)\cot \theta _\xi , 1/2) , \quad &(\widetilde {x}_2,\widetilde {y}_2) &= (-1/2, y - (1/2 + x)\tan \theta _\xi ), \nonumber \\ (\widetilde {x}_3,\widetilde {y}_3) &= (x - (1/2 + y)\cot \theta _\xi , -1/2) , \quad &(\widetilde {x}_4,\widetilde {y}_4) &= (1/2, y + (1/2 - x)\tan \theta _\xi ), \end{align}

where

(3.18) \begin{equation} i= \begin{cases} 1,\quad &\pi +\theta _1\lt \theta _\xi \lt 2\pi -\theta _2, \\ 2,&2\pi -\theta _2\lt \theta _\xi \lt 2\pi +\theta _3, \\ 3,&2\pi +\theta _3\lt \theta _\xi \lt 3\pi -\theta _4, \\ 4,&3\pi -\theta _4\lt \theta _\xi \lt 3\pi +\theta _1. \end{cases} \end{equation}

The quantity $\mathcal{J}(f^{(0)})$ appearing in (3.15) is obtained by substituting the free-molecular solution from (3.7) and (3.13) into the pertinent collision kernel, namely the HS (2.7), BGK (2.9) or Shakhov (2.10) models. By symmetry, we have

(3.19) \begin{equation} \rho _1^{(1)}(x) \equiv \rho _3^{(1)}(x) \quad \mathrm{and} \quad \rho _2^{(1)}(y) \equiv \rho _4^{(1)}(y), \end{equation}

so that only two integral equations are required to determine the wall flux correction functions. The impermeability condition on the cavity walls takes the form

(3.20) \begin{equation} \int _{\boldsymbol{\xi }\boldsymbol{\cdot }\hat{\boldsymbol{n}}\gt 0}\boldsymbol{\xi }\boldsymbol{\cdot }\hat{\boldsymbol{n}}f^{(1)}(\boldsymbol{\widetilde {r}}_i,\boldsymbol{\xi })\,\mathrm{d}\boldsymbol{\xi }+ \int _{\boldsymbol{\xi }\boldsymbol{\cdot }\hat{\boldsymbol{n}}\lt 0}\boldsymbol{\xi }\boldsymbol{\cdot }\hat{\boldsymbol{n}}f^{(1)}(\boldsymbol{\widetilde {r}}_i,\boldsymbol{\xi })\,\mathrm{d}\boldsymbol{\xi }=0, \end{equation}

where the first and second integrals express the separate contributions of the outgoing and incoming particles to the macroscopic gas velocity normal to the surface, respectively. Considering the impermeability balance over the $y=1/2$ surface and starting with the contribution of reflected particles, we substitute (3.15) into the first integral in (3.20), to obtain

(3.21) \begin{equation} \int _{\boldsymbol{\xi }\boldsymbol{\cdot }\hat{\boldsymbol{n}}\gt 0}\xi _r\sin {\theta _\xi }f^{(1)}(x,1/2,\boldsymbol{\xi })\,\mathrm{d}\boldsymbol{\xi }=-\frac {\rho _1^{(1)}(x)}{2\sqrt {\pi }} . \end{equation}

Then, the contribution of the incoming particles to the flux balance, accounted for by the second integral in (3.20), is

(3.22) \begin{align} &\int _{\boldsymbol{\xi }\boldsymbol{\cdot }\hat{\boldsymbol{n}}\lt 0}\xi _r\sin {\theta _\xi }f^{(1)}(x,1/2,\boldsymbol{\xi })\,\mathrm{d}\boldsymbol{\xi } = \int _{0}^{\theta _3} \frac {\rho ^{(1)}_2(\widetilde {y}_2)\sqrt {R_T}}{4\sqrt {\pi }}\sin {\theta _\xi }\,\mathrm{d}\theta _\xi \nonumber \\ &\quad+\int _{\theta _3}^{\pi -\theta _4} \frac {\rho ^{(1)}_3(\widetilde {x}_3)}{4\sqrt {\pi }}\sin {\theta _\xi }\,\mathrm{d}\theta _\xi +\int _{\pi -\theta _4}^{\pi } \frac {\rho ^{(1)}_4(\widetilde {y}_4)\sqrt {R_T}}{4\sqrt {\pi }}\sin {\theta _\xi }\,\mathrm{d}\theta _\xi \nonumber \\ &\quad+\int _{-\infty }^\infty \int _0^\infty \int _{0}^{\pi } \sin {\theta _\xi }\left [\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\mathcal{J}(f^{(0)})\,\mathrm{d}l\right ]\,\mathrm{d}\boldsymbol{\xi }, \end{align}

where the integrations over $\xi _z$ and $\xi _r$ in the first three terms on the right-hand side are evaluated analytically. Combining (3.21) and (3.22), and using the symmetry relation (3.19), we obtain the no-penetration condition over wall 1,

(3.23) \begin{align} &\frac {\rho _1^{(1)}(x)}{2\sqrt {\pi }}-\int _{0}^{\theta _3} \frac {\rho ^{(1)}_2(\widetilde {y}_2)\sqrt {R_T}}{4\sqrt {\pi }}\sin {\theta _\xi }\,\mathrm{d}\theta _\xi -\int _{\theta _3}^{\pi -\theta _4} \frac {\rho ^{(1)}_1(\widetilde {x}_3)}{4\sqrt {\pi }}\sin {\theta _\xi }\,\mathrm{d}\theta _\xi \nonumber \\&\quad-\int _{\pi -\theta _4}^{\pi } \frac {\rho ^{(1)}_2(\widetilde {y}_4)\sqrt {R_T}}{4\sqrt {\pi }}\sin {\theta _\xi }\,\mathrm{d}\theta _\xi =\int _{-\infty }^\infty \int _0^\infty \int _{0}^{\pi } \sin {\theta _\xi }\left [\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\mathcal{J}(f^{(0)})\,\mathrm{d}l\right ]\,\mathrm{d}\boldsymbol{\xi }. \end{align}

Following a similar procedure to obtain the impermeability balance over wall 4, we find

(3.24) \begin{align} &\frac {\rho _2^{(1)}(y)\sqrt {R_T}}{2\sqrt {\pi }}-\int _{3\pi /2}^{2\pi -\theta _2} \frac {\rho ^{(1)}_1(\widetilde {x}_1)}{4\sqrt {\pi }}\cos {\theta _\xi }\,\mathrm{d}\theta _\xi -\int _{2\pi -\theta _2}^{2\pi +\theta _3} \frac {\rho ^{(1)}_2(\widetilde {y}_2)\sqrt {R_T}}{4\sqrt {\pi }}\cos {\theta _\xi }\,\mathrm{d}\theta _\xi \nonumber \\&\quad-\int _{2\pi +\theta _3}^{5\pi /2} \frac {\rho ^{(1)}_1(\widetilde {x}_3)}{4\sqrt {\pi }}\cos {\theta _\xi }\,\mathrm{d}\theta _\xi =\int _{-\infty }^\infty \int _0^\infty \int _{3\pi /2}^{5\pi /2} \cos {\theta _\xi }\left [\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\mathcal{J}(f^{(0)})\,\mathrm{d}l\right ]\,\mathrm{d}\boldsymbol{\xi }. \end{align}

In (3.23) and (3.24), some details on the evaluation of the right-hand side terms are in order. For the BGK and Shakhov models, the calculation of $\mathcal{J}_{\textit{BGK}}(f^{(0)})$ and $\mathcal{J}_{\textit{S}}(f^{(0)})$ is straightforward (see (2.9)–(2.11) and (3.13)). In contrast, for the HS model, the computation of $\mathcal{J}_{\textit{HS}}(f^{(0)})$ involves a five-dimensional quadrature (see (2.7)). Specifically, $\mathrm{d}\varOmega =\sin \phi _k\,\mathrm{d}\theta _k\,\mathrm{d}\phi _k$ , and, using spherical coordinates in the $\boldsymbol{\xi }_1$ -space,

(3.25) \begin{equation} \mathrm{d}\boldsymbol{\xi _1}=\xi _1^2\sin \phi _{\xi _1}\,\mathrm{d}\xi _1\,\mathrm{d}\theta _{\xi _1}\, \mathrm{d}\phi _{\xi _1}, \end{equation}

where $\theta _{\xi _1}$ and $\theta _k$ denote the in-plane azimuthal angles. These integrals are evaluated numerically using Simpson’s rule, implemented in vectorised form in MATLAB. We find that discretisations of $N_{\theta _k}=N_{\phi _k}=N_{\xi _1}=N_{\phi _{\xi _1}}\approx 20$ and $N_{\theta _{\xi _1}}\approx 60$ points are sufficient for obtaining converged results. Proceeding to evaluate the right-hand side quadratures, the $\xi _z$ and $\xi _r$ integrals may be carried out analytically for the BGK and Shakhov models, or numerically for the HS model. The evaluation of the right-hand side terms is completed by taking the path and $\theta _\xi$ quadratures, which are calculated numerically for all three models.

Equations (3.23) and (3.24) form a system of coupled non-homogeneous integral equations for the wall-flux correction functions. The equations are solved by discretising the fluxes at $N_w$ equally spaced points along the cavity walls, yielding a system of $2N_w$ linear non-homogeneous algebraic equations. The system is inverted using a MATLAB subroutine. We find that $N_w \approx 50$ points suffice to obtain converged results.

Notably, for the BGK and Shakhov models, the procedure requires negligible computational effort compared with the molecular simulation calculations described in § 5, and is completed within seconds on a standard Intel Core i7 processor. In contrast, the HS model calculation is more demanding, requiring several days on a 32-core workstation. Note that, for the HS case, the asymptotic solution needs to be computed only once for each value of $R_T$ , after which the solution may be readily obtained for any Knudsen number, as explained in the following subsection.

3.3. The hydrodynamic fields

Using (2.12) together with (3.1) and introducing

(3.26) \begin{equation} I^{(1)}_{i,j}=\int _{-\infty }^{\infty }\int _{0}^{\infty }f^{(1)}\xi _r^i\xi _z^j\,\mathrm{d}\xi _r\,\mathrm{d}\xi _z, \end{equation}

the near-free-molecular density, temperature and mass-flux components are given by the sum of their free-molecular and first-order corrections,

(3.27) \begin{equation}\begin{gathered} \rho =\rho ^{(0)}+\frac {1}{\textit{Kn}}\int _0^{2\pi }I^{(1)}_{1,0}\,\textrm{d}\theta _\xi , \quad T = \frac {2}{3\rho }\left (\frac {3}{2}p^{(0)}+\frac {1}{\textit{Kn}}\int _0^{2\pi }(I^{(1)}_{3,0}+I^{(1)}_{1,2})\textrm{d}\theta _\xi -\rho (u_x^2+u_y^2)\right ) , \\ \rho u_x=\frac {1}{\textit{Kn}}\int _0^{2\pi }I^{(1)}_{2,0}\cos \theta _\xi\, \textrm{d}\theta _\xi \quad \mathrm{and} \quad \rho u_y=\frac {1}{\textit{Kn}}\int _0^{2\pi }I^{(1)}_{2,0}\sin \theta _\xi \,\textrm{d}\theta _\xi , \end{gathered}\end{equation}

respectively. To fix the density field $\rho$ , the normalisation constant $A$ (see (3.9)) must be recalculated, since the first-order correction modifies the density by an $O(\textit{Kn}^{-1})$ addition. Consequently, the normalisation constant is given by

(3.28) \begin{equation} A=A^{(0)}+\textit{Kn}^{-1}A^{(1)}. \end{equation}

The corrected value of $A$ is obtained numerically by imposing the normalisation condition in (2.4), which ensures the conservation of the mass in the cavity.

Advantageously, and different from the numerical schemes outlined in § 5, once the first-order correction fields have been computed, the solution at any Knudsen number (for a given $R_T$ ) may be obtained through mere substitution of its value in (3.27). For the BGK and Shakhov models, the $\xi _r$ and $\xi _z$ integrals in (3.26) may be evaluated in a closed form. To this end, we substitute (3.15) together with (2.9) or (2.10) into (3.26), to obtain

(3.29) \begin{eqnarray} &I^{(1),\textit{BGK}}_{1,0}&\,= \frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)}{2\pi }+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{4}\left (\frac {\rho ^{(0)}}{\sqrt {T^{(0)}}}-\frac {\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)}{\sqrt {T_i}}\right )\,\mathrm{d}l, \nonumber \\ &I^{(1),\textit{BGK}}_{2,0}&\,= \frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)\sqrt {T_i}}{4\sqrt {\pi }}+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{4\sqrt {\pi }}\left (\rho ^{(0)}-\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)\right )\,\mathrm{d}l, \nonumber \\ &I^{(1),\textit{BGK}}_{3,0}&\, =\frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)T_i}{2\pi }+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{8}\left (\rho ^{(0)}\sqrt {T^{(0)}}-\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)\sqrt {T_i}\right )\,\mathrm{d}l, \nonumber \\ &I^{(1),\textit{BGK}}_{1,2}& \,= \frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)T_i}{4\pi }+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{4}\left (\rho ^{(0)}\sqrt {T^{(0)}}-\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)\sqrt {T_i}\right )\,\mathrm{d}l, \end{eqnarray}

for the BGK model, and

(3.30) \begin{eqnarray} &I^{(1),S}_{1,0}& \,= \frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)}{2\pi }+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{4}\left (\frac {\rho ^{(0)}}{\sqrt {T^{(0)}}}-\frac {\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)}{\sqrt {T_i}}-4\frac {q_x^{(0)}\cos {\theta _\xi }+q_y^{(0)}\sin {\theta _\xi }}{15(T^{(0)})^2\sqrt {\pi }}\right )\,\mathrm{d}l, \nonumber \\ &I^{(1),S}_{2,0}& \,= \frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)\sqrt {T_i}}{4\sqrt {\pi }}+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{4\sqrt {\pi }}\left (\rho ^{(0)}-\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)-\frac {q_x^{(0)}\cos {\theta _\xi }+q_y^{(0)}\sin {\theta _\xi }}{15(T^{(0)})^{3/2}}\right )\,\mathrm{d}l, \nonumber \\ &I^{(1),S}_{3,0}& \,=\frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)T_i}{2\pi }+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{8}\left (\rho ^{(0)}\sqrt {T^{(0)}}-\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)\sqrt {T_i}\right )\,\mathrm{d}l, \nonumber \\ &I^{(1),S}_{1,2}&\,= \frac {\rho ^{(1)}_i(\boldsymbol{\widetilde {r}}_i)T_i}{4\pi }+\int _{\boldsymbol{\widetilde {r}}_i}^{\boldsymbol{r}}\frac {\rho ^{(0)}\sqrt {T^{(0)}}}{4}\left (\rho ^{(0)}\sqrt {T^{(0)}}-\rho _i^{(0)}(\boldsymbol{\widetilde {r}}_i)\sqrt {T_i}\right )\,\mathrm{d}l, \end{eqnarray}

for the Shakhov model. These expressions will be used in § 6 to rationalise some of our findings. For the HS interaction model, the corresponding $I^{(1),\textit{HS}}_{i,j}$ quantities are computed numerically.

5.1. Numerical implementation of the particle solvers

To validate the asymptotic predictions and investigate the flow behaviour at intermediate rarefaction rates, particle-based numerical calculations were performed for the HS, BGK and Shakhov descriptions. The three solvers differ only in the way the cell-wise collision step is modelled, while sharing the same particle transport, wall-reflection, cell-sorting and macroscopic sampling procedures. The calculations were carried out in dimensional variables and were subsequently non-dimensionalised using the scales introduced in § 2.

The simulations were performed for an argon gas in a square cavity of side length $L^*=10^{-3}\,\mathrm{m}$ . The horizontal walls were maintained at $T_h^*=1000\,\mathrm{K}$ , while the vertical walls were kept at $T_c^*=200\,\mathrm{K}$ or $T_c^*=500\,\mathrm{K}$ corresponding to $R_T=0.2$ or $R_T=0.5$ , respectively. The molecular mass and hard-sphere diameter were taken as $m^*=6.63\times 10^{-26}\,\mathrm{kg}$ and $d^*=4.17\times 10^{-10}\,\mathrm{m}$ , respectively, yielding the hard-sphere collision cross-section $\sigma ^*=\pi d^{*2}=5.46\times 10^{-19}\,\mathrm{m^2}$ . The most probable molecular speed based on the hot-wall temperature was then ${U^*_{{mp}_{h}}} \approx 645\,\mathrm{m\,s}^{-1}$ . For each Knudsen number, the reference number density was evaluated based on the hard-sphere mean-free-path relation (cf. (2.3)),

(5.1) \begin{equation} n_0^* =\left ({\sqrt {2}\sigma ^* \textit{Kn} L^*}\right )^{-1}, \end{equation}

and the simulator-particle weight was then set to

(5.2) \begin{equation} F_N={n_0^* L^{*2}}/{N_p}, \end{equation}

where $N_p$ is the total number of simulator particles in the two-dimensional computational domain per unit span.

Particle calculations used a uniform Cartesian mesh with $200$ cells in each direction. The fixed time step was taken

(5.3) \begin{equation} \Delta t^*=0.15\,\frac {\Delta x^*}{{U^*_{{mp}_{h}}}} =1.162\times 10^{-9}\,\mathrm{s}. \end{equation}

The total number of simulator particles varied between $5.0\times 10^6$ and $2.2\times 10^7$ , depending on the model and Knudsen number, and corresponding to approximately $125$ – $550$ particles per cell on the $200\times 200$ grid. Most calculations used up to $2.2\times 10^7$ particles. Sampling was started after an initial transient interval, typically after $10^4$ – $6\times 10^4$ time steps, and the total number of time steps ranged from $2.0\times 10^6$ to $5.4\times 10^6$ . Selected high-statistics cases exploited the symmetry of the problem by computing only one quarter of the domain and imposing specular-reflection boundary conditions along the symmetry lines, $x=0$ and $y=0$ . These calculations used a $100\times 100$ simulated-quadrant grid (equivalent to $200\times 200$ over the full cavity), $N_{p,q}=8.0\times 10^7$ simulator particles in the quadrant ( $8000$ particles per simulated cell) and $\Delta t^*=1.162\times 10^{-9}\,\mathrm{s}$ . The corresponding particle weight was $F_N=n_0^*L^{*2}/(4N_{p,q})$ . Sampling commenced after $10^5$ time-steps, and was performed every two steps, with a total of $3\times 10^6$ time steps for these runs.

The numerical implementation was GPU-accelerated. CuPy was used for GPU-resident array operations, including cell sorting and particle indexing, while Numba-CUDA kernels were used for the main particle-level operations, including particle initialisation, ballistic transport, diffuse wall reflection and the BGK and Shakhov relaxation steps. The hard-sphere solver was also GPU accelerated, applying cell-wise NTC candidate-pair generation, the collision acceptance/rejection algorithm and post-collision velocity updates on the GPU.

The differences between the numerical schemes employed lie entirely in the implementation of the cell-wise collision step, as described in the following. This common computational framework makes the comparison between the three kinetic descriptions especially transparent, since differences in the results may be attributed to the collision model rather than to differences in particle transport, wall reflection, grid resolution or sampling procedure.

Prior to applying the numerical scheme to the current setup, we verified our calculations via solving the planar (one-dimensional) Fourier heat-conduction problem (Rader et al. Reference Rader, Gallis, Torczynski and Wagner2006). The problem was solved at $\textit{Kn} = 0.02$ for stationary and fully diffuse walls. The walls were positioned at $x=0$ and $x=L=1.0\times 10^{-3}\,\mathrm{m}$ and maintained at temperatures $1020\,\mathrm{K}$ and $980\,\mathrm{K}$ , respectively. The same molecular parameters have been used as in the cavity calculation. Three statistically independent realisations, each containing two million particles, have been simulated, yielding ensemble-averaged values for the thermal conductivities of the HS, BGK and Shakhov models. In the HS case, the calculation yielded $\kappa _{\textit{HS}}=0.02382\ \mathrm{W\ m}^{-1}\ \mathrm{K}^{-1}$ , which differs by only $0.57\,\,\%$ from the first Chapman–Enskog prediction (Wang et al. Reference Wang, Ramírez-Hinestrosa and Frenkel2020). In addition, the calculated ratios of gas heat conductivities, $\kappa _{HS}/\kappa _{\textit{BGK}} = 1.497$ and $\kappa _{\textit{Shakhov}}/\kappa _{\textit{HS}} = 1.033$ , were found within 1.7 % and 1.4 % deviations from their theoretical values, $1.473$ and $1.019$ , respectively. These results confirm that the implementations correctly distinguish the expected thermal-transport behaviour of the three collision models.

5.2. HS model: the DSMC method

To simulate the HS gas system, the DSMC method, first introduced by Bird (Reference Bird1994), was applied. Intermolecular collisions were modelled as discrete binary hard-sphere encounters using the standard no-time-counter (NTC) procedure. For each cell, the code estimated the number of candidate collision pairs from the instantaneous local particle population, the simulator-particle weight, the cell volume, the time step and an upper bound for the product of collision cross-section and relative speed. For a cell containing $N_c$ particles, the number of candidate collision pairs was evaluated in proportion to

(5.4) \begin{equation} \frac {1}{2}N_c(N_c-1) \frac {F_N(\sigma ^* g^*)_{max }\Delta t^*}{V_c^*}, \end{equation}

where $g^*$ is the relative speed between particles and $V_c^*$ is the cell volume per unit span. Candidate pairs were selected randomly and accepted with probability

(5.5) \begin{equation} P_{\textit{coll}}=\frac {\sigma ^* g^*}{(\sigma ^* g^*)_{max }}. \end{equation}

For accepted collisions, the post-collision relative velocity direction was sampled isotropically in the centre-of-mass frame.

5.3. BGK particle scheme

In the BGK mode, the collision stage does not simulate binary molecular encounters explicitly. Instead, the detailed collision process is replaced by stochastic relaxation toward a local Maxwellian equilibrium inside each cell. This particle implementation follows the relaxation strategy of Gallis & Torczynski (Reference Gallis and Torczynski2000), as summarised in the following.

At each time step, the local hydrodynamic moments were first computed from the particles residing in each cell. These included the cell-average velocity, temperature, number density and pressure. The local number density was evaluated as

(5.6) \begin{equation} n^*=\frac {N_cF_N}{V_c^*}, \end{equation}

and the pressure was obtained from the ideal-gas relation

(5.7) \begin{equation} p^*=n^*k_B T^*. \end{equation}

The local BGK collision frequency was then evaluated as

(5.8) \begin{equation} \nu ^*=\frac {p^*}{\mu ^*}, \end{equation}

where $\mu ^*$ is the hard-sphere viscosity evaluated at the local cell temperature. The collision frequency $\nu ^*=p^*/\mu ^*$ used in the particle solvers was calibrated to yield the same dimensionless collision frequency $\nu =\rho \sqrt {\pi T}/2$ (stated in (2.11)) and the same hard-sphere mean free path defining the Knudsen number in (2.3), so that the asymptotic and numerical results correspond to an identical rarefaction level.

The probability that a particle relaxes during one time step was set to

(5.9) \begin{equation} P_{\textit{relax}}=1-\exp (-\nu ^*\Delta t^*). \end{equation}

For each particle selected for relaxation, its velocity was resampled from the local Maxwellian distribution defined by the cell-averaged hydrodynamic velocity and temperature. Specifically, denoting by $\boldsymbol{u}_{\textit{cell}}^* =(u_{x,\textit{cell}}^*,u_{y,\textit{cell}}^*,u_{z,\textit{cell}}^*)$ the dimensional hydrodynamic velocity obtained from the particles in the cell, the post-relaxation velocity components were sampled as

(5.10) \begin{equation} \xi _\alpha ^* = u_{\alpha ,\textit{cell}}^* + \left (\frac {k_B T_{\textit{cell}}^*}{m^*}\right )^{1/2}G_\alpha , \quad \alpha =x,y,z, \end{equation}

where $T_{\textit{cell}}^*$ is the local cell temperature and $G_x,G_y,G_z$ are independent standard Gaussian random variables. Although the physical flow field is two-dimensional, the molecular velocity space is three-dimensional; hence, the $z$ -component is retained in the Maxwellian sampling. In the present configuration $u_{z,\textit{cell}}^*=0$ by symmetry, up to statistical fluctuations.

The BGK simulations used the same transport, wall-reflection, cell-sorting, grid and sampling procedures as the DSMC calculations.

5.4. Shakhov particle scheme

The Shakhov particle solver followed the same general relaxation structure as the BGK solver, but modified the post-relaxation sampling to account for the leading heat-flux correction and recover the correct monatomic-gas Prandtl number. As in the BGK calculation, the code first computed the local density, velocity, temperature, pressure and relaxation probability in each cell. In addition, the local heat-flux vector was estimated using the peculiar velocity vector,

(5.11) \begin{equation} \boldsymbol{c}^*=\boldsymbol{\xi }^*-\boldsymbol{u}^*, \end{equation}

as

(5.12) \begin{equation} \boldsymbol{q}^*=\frac {1}{2}m^* n^* \langle c^{*2}\boldsymbol{c}^*\rangle . \end{equation}

The relaxation probability was taken the same as in the BGK solver,

(5.13) \begin{equation} P_{\textit{relax}}=1-\exp (-\nu ^*\Delta t^*), \end{equation}

with $\nu ^*=p^*/\mu ^*$ . For each particle selected for relaxation, a candidate velocity was first sampled from the local Maxwellian distribution. The Maxwellian proposals were reweighted through an acceptance–rejection procedure, using the Shakhov correction weight

(5.14) \begin{equation} w= 1+(1-\textit{Pr}) \frac {2\boldsymbol{c}^*\boldsymbol{\cdot }\boldsymbol{q}^*}{5p^*\theta ^*} \left ( \frac {c^{*2}}{2\theta ^*}-\frac {5}{2} \right ), \quad \mathrm{where} \ \theta ^*=\frac {k_{\mathrm B}T^*}{m^*} \ \mathrm{and} \ \textit{Pr}=\frac {2}{3}, \end{equation}

such that the accepted post-relaxation velocities followed the Shakhov equilibrium distribution. In the implementation, the correction was regularised as

(5.15) \begin{equation} w_{\textit{lim}} = \min \left \{2, \max (0, w)\right \}, \end{equation}

with proposal envelope $A=2$ and acceptance probability $w_{\textit{lim}}/A$ . Prior to the finite-trial fallback, accepted samples followed the normalised regularised target, proportional to $f_M w_{\textit{lim}}$ , rather than the unmodified Shakhov target, in cells where clipping occurred. At most four proposals were attempted, and failure led to a Maxwellian fallback. The largest above-limit candidate-event rate was $0.0225\,\%$ of selected relaxations, and the aggregate fallback rate was $3.81\times 10^{-9}\,\%$ . These diagnostics indicate that the regularisation was rare, while the planar Fourier-flow validation independently tested the resulting transport coefficients.