2.1. Problem and assumptions

Consider a gas in a cavity with a square cross-section $-L/2 \lt X_{1} \lt L/2$ , $-L/2 \lt X_{2} \lt L/2$ , where $X_{i}$ are the Cartesian coordinates (see figure 1). The cavity is infinitely long in the $X_{3}$ direction. Its upper part ( $X_{2} = L/2$ ) is closed by a planar lid moving at a constant speed $U_{w}$ in the $X_{1}$ direction. The lid and cavity walls are kept at a common uniform temperature $T_{0}$ . The average density of the gas over the cross-section is $\rho _{0}$ . No external force acts on the gas. In comparison with Hattori (Reference Hattori2024), the geometry is identical, whereas the lid motion is in a different direction. We investigate the steady behaviour of the gas under the assumptions that (i) the behaviour of the gas is described by the Bhatnagar–Gross–Krook (BGK) equation (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954; Welander Reference Welander1954); (ii) the gas molecules are diffusely reflected on the surface of the lid and cavity walls; (iii) the lid speed $U_{w}$ is so small compared with the most probable speed $(2RT_{0})^{1/2}$ of gas molecules that the equation and boundary condition can be linearised around a uniform equilibrium state at rest with density $\rho _{0}$ and temperature $T_{0}$ , where $R$ denotes the specific gas constant.

Figure 1.

Sketch of the problem.

2.2. Basic equation and boundary condition

We denote by $X_{i}=Lx_{i}$ ( $i = 1, 2, 3$ ) the position, by $(2RT_{0})^{1/2}\boldsymbol{\zeta }$ (or $(2RT_{0})^{1/2}\zeta _{i}$ ) the molecular velocity, by $\rho _{0}(2RT_{0})^{-3/2}E(\zeta )[1+\phi (x_{1},x_{2},\boldsymbol{\zeta })]$ the velocity distribution function (VDF) of gas molecules and by $U_{w} = (2RT_{0})^{1/2}u_{w}$ the lid speed, where $E(\zeta ) = \pi ^{-3/2}{\rm e}^{-\zeta ^{2}}$ and $\zeta = |\boldsymbol{\zeta }|$ . Moreover, let us denote by $\rho _{0}[1+\omega (x_{1},x_{2})]$ the density of the gas, by $(2RT_{0})^{1/2}u_{i}(x_{1},x_{2})$ the flow velocity, by $T_{0}[1+\tau (x_{1},x_{2})]$ the temperature, by $p_{0}[1+P(x_{1},x_{2})]$ the pressure and by $p_{0}[\delta _{\textit{ij}}+P_{\textit{ij}}(x_{1},x_{2})]$ the stress tensor, respectively. Here, $p_{0} = \rho _{0}RT_{0}$ is the reference pressure and $\delta _{\textit{ij}}$ is the Kronecker delta.

From the assumptions (i)–(iii) in § 2.1, the behaviour of the gas is described by the boundary-value problem of the linearised BGK equation for the VDF $\phi (x_{1},x_{2},\boldsymbol{\zeta })$ :

(2.1a) \begin{align} \zeta _{1}\frac {\partial \phi }{\partial x_{1}} + \zeta _{2}\frac {\partial \phi }{\partial x_{2}} &= \frac {1}{k}\left [\omega + 2\zeta _{j}u_{j} + \tau \left (\zeta ^{2}-\frac {3}{2}\right ) - \phi \right ] \big [(x_{1}, x_{2}) \in \varOmega ,\ \boldsymbol{\zeta } \in \mathbb{R}^{3}\big ], \\[-10pt] \nonumber \end{align}

(2.1b) \begin{align} \omega &= \int \phi E \mathrm{d} \boldsymbol{\zeta }, \\[-10pt] \nonumber \end{align}

(2.1c) \begin{align} u_{j} &= \int \zeta _{j}\phi E \mathrm{d} \boldsymbol{\zeta }, \\[-10pt] \nonumber\end{align}

(2.1d) \begin{align} \tau &= \frac {2}{3}\int \left (\zeta ^{2}-\frac {3}{2}\right )\phi E \mathrm{d} \boldsymbol{\zeta }, \\[-10pt] \nonumber \end{align}

(2.1e) \begin{align} k &= \frac {\sqrt {\pi }}{2}\textit{Kn}, \quad \textit{Kn} = \frac {\ell _{0}}{L}, \quad \ell _{0} = \frac {(8RT_{0}/\pi )^{1/2}}{A_{\textit{c}}\rho _{0}}, \\[-10pt] \nonumber \end{align}

(2.1f) \begin{align} \text{b.c.} \quad \phi &= \sigma + 2\zeta _{1}u_{w}\chi _{\textit{U}} \quad [(x_{1},x_{2}) \in \partial \varOmega ,\ \zeta _{i}n_{i} \gt 0], \\[-10pt] \nonumber \end{align}

(2.1g) \begin{align} \sigma &= -2\sqrt {\pi }\int _{\zeta _{j}n_{j}\lt 0} \zeta _{i}n_{i}\phi E \mathrm{d} \boldsymbol{\zeta } \quad [(x_{1},x_{2}) \in \partial \varOmega ]. \\[10pt] \nonumber \end{align}

In the above, $\varOmega$ and $\partial \varOmega$ correspond to the cross-section of the cavity and its boundary, respectively, and $n_{i}$ is the unit vector normal to the boundary pointed to the gas:

(2.2a) \begin{align} \varOmega &= \left \{(x_{1},x_{2}) \ \left | -\frac {1}{2} \lt x_{1} \lt \frac {1}{2},\ -\frac {1}{2} \lt x_{2} \lt \frac {1}{2} \right. \right \}\!, \\[-10pt] \nonumber \end{align}

(2.2b) \begin{align} \partial \varOmega &= \left \{(x_{1},x_{2}) \ \left | \left (x_{1}=\pm \frac {1}{2},\ -\frac {1}{2} \lt x_{2} \lt \frac {1}{2}\right ) \text{ or } \left (-\frac {1}{2} \lt x_{1} \lt \frac {1}{2},\ x_{2}=\pm \frac {1}{2}\right ) \right. \right \}\!, \\[-10pt] \nonumber \end{align}

(2.2c) \begin{align} n_{i} &= \begin{cases} \displaystyle (\mp 1, 0, 0) & \displaystyle \quad \left (x_{1}=\pm \frac {1}{2},\ -\frac {1}{2} \lt x_{2} \lt \frac {1}{2}\right )\!, \\[9pt] \displaystyle (0, \mp 1, 0) & \displaystyle \quad \left (-\frac {1}{2} \lt x_{1} \lt \frac {1}{2},\ x_{2}=\pm \frac {1}{2}\right )\!. \end{cases} \\[10pt] \nonumber \end{align}

The quantities $A_{\textit{c}}\rho _{0}$ and $\ell _{0}$ in (2.1e ) are, respectively, the collision frequency and the mean free path of gas molecules for the BGK model in the equilibrium state at rest with density $\rho _{0}$ and temperature $T_{0}$ , where $A_{\textit{c}}$ is a constant. The Knudsen number is denoted by $\textit{Kn}$ . In the boundary condition (2.1f ), $\chi _{\textit{U}}$ is an indicator function defined on the boundary $\partial \varOmega$ , which is unity on the lid, $-1/2 \lt x_{1} \lt 1/2$ and $x_{2} = 1/2$ , and vanishes elsewhere. In (2.1b )–(2.1d ) and (2.3) appearing soon, the integration domain is the whole space of the molecular velocity $\boldsymbol{\zeta }$ .

Equations (2.1b )–(2.1d ) define the perturbations of the density, flow velocity and temperature. Similarly, the perturbations of the pressure and stress tensor are defined as the following moments of the VDF $\phi$ :

(2.3a) \begin{align} P &= \omega + \tau = \frac {2}{3}\int \zeta ^{2}\phi E \mathrm{d} \boldsymbol{\zeta }, \end{align}

(2.3b) \begin{align} P_{\textit{ij}} &= 2\int \zeta _{i}\zeta _{j}\phi E \mathrm{d} \boldsymbol{\zeta }. \end{align}

The first equality in (2.3a ) represents the linearised equation of state for an ideal gas. As evident from (2.3b ), $P_{\textit{ij}}$ is a symmetric tensor.

We can assume without loss of generality that the solution $\phi (x_{1},x_{2},\zeta _{1},\zeta _{2},\zeta _{3})$ of problem (2.1) is even in the molecular velocity component $\zeta _{3}$ . With this assumption, the flow velocity component $u_{3}$ , which appears on the right-hand side of (2.1a ), and the stress components $P_{13}$ and $P_{23}$ vanish, i.e.

(2.4) \begin{equation} u_{3}(x_{1},x_{2}) \equiv 0, \qquad P_{13}(x_{1},x_{2}) \equiv 0, \qquad P_{23}(x_{1},x_{2}) \equiv 0. \end{equation}

As expected from the geometry and direction of lid motion (see also figure 1), the solution $\phi$ can be considered to be symmetric with respect to $x_{1} = 0$ as

(2.5) \begin{equation} \phi (x_{1},x_{2},\zeta _{1},\zeta _{2},\zeta _{3}) = -\phi (-x_{1},x_{2},-\zeta _{1},\zeta _{2},\zeta _{3}). \end{equation}

Considering (2.5), it is seen that, among the macroscopic quantities, the density $\omega$ , flow velocity component $u_{2}$ , temperature $\tau$ , diagonal stress components $P_{11}$ , $P_{22}$ and $P_{33}$ and the function $\sigma$ in the boundary condition (2.1f ) are all odd in $x_{1}$ , whereas the flow velocity component $u_{1}$ and the shear stress component $P_{12}$ are even in $x_{1}$ :

(2.6a) \begin{align} &g(x_{1},x_{2}) = -g(-x_{1},x_{2}) \quad [g = \omega , u_{2}, \tau , P_{11}, P_{22}, P_{33};\ (x_{1}, x_{2}) \in \varOmega ], \end{align}

(2.6b) \begin{align} &\sigma (x_{1},x_{2}) = -\sigma (-x_{1},x_{2}) \quad [(x_{1}, x_{2}) \in \partial \varOmega ], \end{align}

(2.6c) \begin{align} &g(x_{1},x_{2}) = g(-x_{1},x_{2}) \quad [g = u_{1}, P_{12};\ (x_{1}, x_{2}) \in \varOmega ]. \end{align}

Before closing this subsection, we briefly comment on the difference from the previous work (Hattori Reference Hattori2024). In that work, the spanwise flow was induced, and thus, the non-zero macroscopic quantities were limited to its speed $u_{3}$ and the shear stress components $P_{13}$ and $P_{23}$ , which vanish in the present problem. In contrast, in the present problem, the flow develops within the cross-section, and the in-plane velocity components $u_{1,2}$ , pressure $P$ , normal stress components $P_{11,22,33}$ , density $\omega$ and temperature $\tau$ become non-zero.

2.3. Integral equations

In the present problem (2.1), the boundary velocity is discontinuous at the upper corners $(x_{1},x_{2}) = (\pm 1/2, 1/2)$ : it is $(u_{w},0)$ on the lid $(x_{1},1/2)$ and $(0,0)$ on the side walls $(\pm 1/2,x_{2})$ . Then, the boundary condition (2.1f ) is also discontinuous, and hence, the VDF $\phi (x_{1},x_{2},\zeta _{1},\zeta _{2},\zeta _{3})$ is discontinuous for $\zeta _{1} \lessgtr 0$ , $\zeta _{2} \lt 0$ at the upper corners $(\pm 1/2, 1/2)$ . This discontinuity propagates along the characteristic lines of (2.1a ) that pass through the upper corners, which complicates the analysis. Fortunately, for the case of the linearised BGK equation with the diffuse reflection condition, which is adopted here also, this difficulty can be avoided by reducing the problem to a system of linear integral equations for macroscopic quantities (see, e.g. Hasegawa & Sone Reference Hasegawa and Sone1988; Aoki, Sone & Yano Reference Aoki, Sone and Yano1989; Sone Reference Sone2007; Taguchi & Charrier Reference Taguchi and Charrier2008). First, integrating (2.1a ) along its characteristic lines under the boundary condition (2.1f ) allows us to express $\phi$ in terms of the density $\omega$ , flow velocity components $u_{1}$ and $u_{2}$ , temperature $\tau$ and the function $\sigma$ in (2.1f ). Here, note that $u_{3} \equiv 0$ by (2.4). Substituting this result into (2.1b ), (2.1c ) ( $j = 1,2$ ), (2.1d ) and the definition (2.1g ) of $\sigma$ to eliminate the VDF $\phi$ , and then carrying out the integration over the molecular velocity $\boldsymbol{\zeta }$ , we obtain the following integral equations for the macroscopic quantities $\omega (\boldsymbol{x})$ , $u_{1}(\boldsymbol{x})$ , $u_{2}(\boldsymbol{x})$ , $\tau (\boldsymbol{x})$ and $\sigma (\boldsymbol{x})$ :

(2.7a) \begin{align} h(\boldsymbol{x}) - A[h](\boldsymbol{x}) - A_{w}[\sigma ](\boldsymbol{x}) &= u_{w} a(\boldsymbol{x}) \quad (\boldsymbol{x} \in \varOmega ), \end{align}

(2.7b) \begin{align} \sigma (\boldsymbol{x}) - B[h](\boldsymbol{x}) - B_{w}[\sigma ](\boldsymbol{x}) &= u_{w} b(\boldsymbol{x}) \quad (\boldsymbol{x} \in \partial \varOmega ). \end{align}

Here, $\boldsymbol{x} = (x_{1}, x_{2})$ is the position and $h(\boldsymbol{x}) = [\omega (\boldsymbol{x}) \ u_{1}(\boldsymbol{x}) \ u_{2}(\boldsymbol{x}) \ \tau (\boldsymbol{x})]^{\mathrm{T}}$ collectively represents the macroscopic quantities, all defined in the gas domain $\varOmega$ , as the column vector. The function $\sigma (\boldsymbol{x})$ , unlike $h$ , is defined only on the boundary surface $\partial \varOmega$ and, hence, enters the equations independently, rather than being included in $h$ . Below, the definitions of the linear integral operators $A$ , $A_{w}$ , $B$ and $B_{w}$ and the source terms $u_{w} a(\boldsymbol{x})$ and $u_{w} b(\boldsymbol{x})$ are given. Since this involves rather long expressions, we first explain the notation. In addition to the position $\boldsymbol{x} = (x_{1}, x_{2})$ , we denote by $\bar {\boldsymbol{x}} = (-x_{1},x_{2})$ its mirror image, by $\boldsymbol{x}_{*} = (x_{1*}, x_{2*})$ the integration variable and by $\boldsymbol{x}_{*}^{(U)} = (x_{1*}, 1/2)$ the position on the lid; they are two-dimensional vectors in the $x_{1}x_{2}$ plane. We denote by $\mathrm{d} \boldsymbol{x}_{*} = \mathrm{d} x_{1*} \mathrm{d} x_{2*}$ a surface element and by $\mathrm{d} l(\boldsymbol{x}_{*})$ a line element along the boundary at $\boldsymbol{x}_{*}$ . The right halves of the cavity cross-section $\varOmega$ and the cavity boundary $\partial \varOmega$ are denoted respectively by $\varOmega _{\textit{R}}$ and $\partial \varOmega _{\textit{R}}$ : $\varOmega _{\textit{R}} = \{ (x_{1},x_{2}) \ | \ (x_{1},x_{2}) \in \varOmega ,\ 0 \lt x_{1} \lt 1/2 \}$ , $\partial \varOmega _{\textit{R}} = \{ (x_{1},x_{2}) \ | \ (x_{1},x_{2}) \in \partial \varOmega ,\ 0 \lt x_{1} \le 1/2 \}$ . With the introduction of the column vector $h$ , matrix and vector notation will be used, and we also use a symbol $\circ$ denoting the elementwise product of two matrices (or two vectors). Now, the operators $A$ , $A_{w}$ , $B$ and $B_{w}$ and the functions $a(\boldsymbol{x})$ and $b(\boldsymbol{x})$ are given by

(2.8a) \begin{align} A[h](\boldsymbol{x}) &= \frac {1}{k\pi }\int _{\varOmega } \frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}K(\boldsymbol{x}_{*}-\boldsymbol{x})h(\boldsymbol{x}_{*}) \mathrm{d} \boldsymbol{x}_{*} \nonumber \\ &= \frac {1}{k\pi }\int _{\varOmega _{\textit{R}}} \bigg \{\frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}K(\boldsymbol{x}_{*}-\boldsymbol{x})h(\boldsymbol{x}_{*})\nonumber \\ &\quad +\frac {1}{|\boldsymbol{x}_{*}-\bar {\boldsymbol{x}}|}\Big [S_{K} \circ K(\boldsymbol{x}_{*}-\bar {\boldsymbol{x}})\Big ]\Big [S_{h}\circ h(\boldsymbol{x}_{*})\Big ]\bigg \} \mathrm{d} \boldsymbol{x}_{*}, \end{align}

(2.8b) \begin{align} \notag\\ K(\boldsymbol{z}) &= [K_{\textit{ij}}(\boldsymbol{z})] \quad (i,j = 1, 2, 3, 4), \nonumber \\ K_{11}(\boldsymbol{z}) &= J_{0}(|\boldsymbol{z}|/k), \ \ K_{12}(\boldsymbol{z}) = -2r_{1}(\boldsymbol{z})J_{1}(|\boldsymbol{z}|/k), \nonumber \\ K_{13}(\boldsymbol{z}) &= -2r_{2}(\boldsymbol{z})J_{1}(|\boldsymbol{z}|/k), \ \ K_{14}(\boldsymbol{z}) = J_{2}(|\boldsymbol{z}|/k)-J_{0}(|\boldsymbol{z}|/k), \nonumber \\ K_{21}(\boldsymbol{z}) &= -r_{1}(\boldsymbol{z})J_{1}(|\boldsymbol{z}|/k), \ \ K_{22}(\boldsymbol{z}) = 2r_{1}^{2}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k), \nonumber \\ K_{23}(\boldsymbol{z}) &= 2r_{1}(\boldsymbol{z})r_{2}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k), \ \ K_{24}(\boldsymbol{z}) = -r_{1}(\boldsymbol{z})[J_{3}(|\boldsymbol{z}|/k)-J_{1}(|\boldsymbol{z}|/k)], \nonumber \\ K_{31}(\boldsymbol{z}) &= -r_{2}(\boldsymbol{z})J_{1}(|\boldsymbol{z}|/k), \ \ K_{32}(\boldsymbol{z}) = 2r_{1}(\boldsymbol{z})r_{2}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k), \nonumber \\ K_{33}(\boldsymbol{z}) &= 2r_{2}^{2}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k), \ \ K_{34}(\boldsymbol{z}) = -r_{2}(\boldsymbol{z})[J_{3}(|\boldsymbol{z}|/k)-J_{1}(|\boldsymbol{z}|/k)], \nonumber \\ K_{41}(\boldsymbol{z}) &= (2/3)[J_{2}(|\boldsymbol{z}|/k)-\!J_{0}(|\boldsymbol{z}|/k)],\!\! \ \ K_{42}(\boldsymbol{z}) = -(4/3)r_{1}(\boldsymbol{z})[J_{3}(|\boldsymbol{z}|/k)-\!J_{1}(|\boldsymbol{z}|/k)], \nonumber \\ K_{43}(\boldsymbol{z}) &= -(4/3)r_{2}(\boldsymbol{z})[J_{3}(|\boldsymbol{z}|/k)-J_{1}(|\boldsymbol{z}|/k)], \\[-10pt] \nonumber K_{44}(\boldsymbol{z}) &= (2/3)[J_{4}(|\boldsymbol{z}|/k)-2J_{2}(|\boldsymbol{z}|/k)+(3/2)J_{0}(|\boldsymbol{z}|/k)], \\[-10pt] \nonumber \end{align}

(2.8c) \begin{align} r_{i}(\boldsymbol{z}) &= \frac {z_{i}}{|\boldsymbol{z}|}, \\[-10pt] \nonumber \end{align}

(2.8d) \begin{align} J_{n}(\xi ) &= \int _{0}^{\infty }t^{n}\exp \left (-t^{2}-\frac {\xi }{t}\right ) \mathrm{d} t \quad (\xi \ge 0), \\[-10pt] \nonumber \end{align}

(2.8e) \begin{align} S_{K} &= \begin{bmatrix} 1 & -1 & 1 & 1 \\ -1 & 1 & -1 & -1 \\ 1 & -1 & 1 & 1 \\ 1 & -1 & 1 & 1 \end{bmatrix}, \quad S_{h} = \begin{bmatrix} -1 \\ 1 \\ -1 \\ -1 \end{bmatrix} , \\[-10pt] \nonumber \end{align}

(2.8f) \begin{align} A_{w}[\sigma ](\boldsymbol{x}) &= -\frac {1}{\pi }\int _{\partial \varOmega } \frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}M(\boldsymbol{x}_{*}-\boldsymbol{x})r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{i}(\boldsymbol{x}_{*})\sigma (\boldsymbol{x}_{*}) \mathrm{d} l(\boldsymbol{x}_{*}) \nonumber \\ & = -\frac {1}{\pi }\int _{\partial \varOmega _{\textit{R}}}\bigg [\frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}M(\boldsymbol{x}_{*}-\boldsymbol{x})r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{i}(\boldsymbol{x}_{*})\sigma (\boldsymbol{x}_{*}) \nonumber \\ &\quad -\frac {1}{|\boldsymbol{x}_{*}-\bar {\boldsymbol{x}}|}[S_{M}\circ M(\boldsymbol{x}_{*}-\bar {\boldsymbol{x}})]r_{i}(\boldsymbol{x}_{*}-\bar {\boldsymbol{x}})n_{i}(\boldsymbol{x}_{*})\sigma (\boldsymbol{x}_{*})\bigg ]\mathrm{d} l(\boldsymbol{x}_{*}), \\[-10pt] \nonumber \end{align}

(2.8g) \begin{align} M(\boldsymbol{z}) &= \begin{bmatrix} J_{1}(|\boldsymbol{z}|/k) \\ -r_{1}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k) \\ -r_{2}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k) \\ (2/3)[J_{3}(|\boldsymbol{z}|/k)-J_{1}(|\boldsymbol{z}|/k)] \end{bmatrix}, \\[-10pt] \nonumber \end{align}

(2.8h) \begin{align} S_{M} &= \begin{bmatrix} 1 \\ -1 \\ 1 \\ 1 \end{bmatrix}, \\[-10pt] \nonumber \end{align}

(2.8i) \begin{align} a(\boldsymbol{x}) &= \frac {1}{\pi }\int _{-1/2}^{1/2}\frac {1}{\left|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right|}r_{2}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right) \nonumber \\ &\quad \times \begin{bmatrix} -2r_{1}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)J_{2}\left(\left|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right|/k\right) \\[6pt] 2r_{1}^{2}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)J_{3}\left(\left|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right|/k\right) \\[6pt] 2r_{1}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)r_{2}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)J_{3}\left(\left|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right|/k\right) \\[6pt] -(4/3)r_{1}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)\left[J_{4}\left(\left|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right|/k\right)-J_{2}\left(\left|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right|/k\right)\right] \end{bmatrix} \mathrm{d} x_{1*}, \\[-10pt] \nonumber \end{align}

(2.8j) \begin{align} B[h](\boldsymbol{x}) &= \frac {2}{k\sqrt {\pi }}\int _{\varOmega } \frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{i}(\boldsymbol{x})W(\boldsymbol{x}_{*}-\boldsymbol{x})h(\boldsymbol{x}_{*}) \mathrm{d} \boldsymbol{x}_{*} \nonumber \\ & = \frac {2}{k\sqrt {\pi }}\int _{\varOmega _{\textit{R}}}\bigg \{\frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{i}(\boldsymbol{x})W(\boldsymbol{x}_{*}-\boldsymbol{x})h(\boldsymbol{x}_{*}) \nonumber \\ &\quad +\frac {1}{|\boldsymbol{x}_{*}-\bar {\boldsymbol{x}}|}r_{i}(\boldsymbol{x}_{*}-\bar {\boldsymbol{x}})n_{i}(\bar {\boldsymbol{x}})\big [S_{W}\circ W(\boldsymbol{x}_{*}-\bar {\boldsymbol{x}})\big ]\big [S_{h}\circ h(\boldsymbol{x}_{*})\big ]\bigg \}\mathrm{d} \boldsymbol{x}_{*}, \\[-10pt] \nonumber \end{align}

(2.8k) \begin{align} W(\boldsymbol{z}) & = \begin{bmatrix} J_{1}(|\boldsymbol{z}|/k) \\ -2r_{1}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k) \\ -2r_{2}(\boldsymbol{z})J_{2}(|\boldsymbol{z}|/k) \\ J_{3}(|\boldsymbol{z}|/k)-J_{1}(|\boldsymbol{z}|/k) \end{bmatrix}^{\mathrm{T}}, \\[-10pt] \nonumber \end{align}

(2.8l) \begin{align} S_{W} &= \begin{bmatrix} 1 & -1 & 1 & 1 \end{bmatrix}, \\[-10pt] \nonumber \end{align}

(2.8m) \begin{align} B_{w}[\sigma ](\boldsymbol{x}) &= -\frac {2}{\sqrt {\pi }}\int _{\partial \varOmega } \frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{i}(\boldsymbol{x}) \nonumber \\ &\quad \times r_{j}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{j}(\boldsymbol{x}_{*})J_{2}\left (\frac {|\boldsymbol{x}_{*}-\boldsymbol{x}|}{k}\right )\sigma (\boldsymbol{x}_{*}) \mathrm{d} l(\boldsymbol{x}_{*}) \nonumber \\ & = -\frac {2}{\sqrt {\pi }}\int _{\partial \varOmega _{\textit{R}}}\bigg [\frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|}r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{i}(\boldsymbol{x}) \nonumber \\ &\quad \times r_{j}(\boldsymbol{x}_{*}-\boldsymbol{x})n_{j}(\boldsymbol{x}_{*})J_{2}\left (\frac {|\boldsymbol{x}_{*}-\boldsymbol{x}|}{k}\right )\sigma (\boldsymbol{x}_{*}) \nonumber \\ &\quad -\frac {1}{|\boldsymbol{x}_{*}-\bar {\boldsymbol{x}}|}r_{i}(\boldsymbol{x}_{*}-\bar {\boldsymbol{x}})n_{i}(\bar {\boldsymbol{x}}) r_{j}(\boldsymbol{x}_{*}-\bar {\boldsymbol{x}})n_{j}(\boldsymbol{x}_{*})J_{2}\left (\frac {|\boldsymbol{x}_{*}-\bar {\boldsymbol{x}}|}{k}\right )\sigma (\boldsymbol{x}_{*})\bigg ] \mathrm{d} l(\boldsymbol{x}_{*}), \\[-10pt] \nonumber \end{align}

(2.8n) \begin{align} b(\boldsymbol{x}) &= -\frac {4}{\sqrt {\pi }}\int _{-1/2}^{1/2}\frac {1}{|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}|}r_{i}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)n_{i}(\boldsymbol{x}) \nonumber \\ &\quad \times r_{1}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)r_{2}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)J_{3}\Big (\frac {|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}|}{k}\Big ) \mathrm{d} x_{1*}. \end{align}

Some additional notes regarding the above equations may be in order. The most right-hand side of (2.8a ) is obtained by splitting the integral over the cavity into its right half, $\varOmega _{\textit{R}}$ , and left half, $\varOmega _{\textit{L}} = \{ (x_{1},x_{2}) \ | \ (x_{1},x_{2}) \in \varOmega , -1/2 \lt x_{1} \lt 0 \}$ , and then rewriting the integral over the latter by changing the integration variable from $\boldsymbol{x}_{*} = (x_{1*}, x_{2*}) \in \varOmega _{\textit{L}}$ to $\bar {\boldsymbol{x}}_{*} = (-x_{1*}, x_{2*}) \in \varOmega _{\textit{R}}$ and then transcribing $\bar {\boldsymbol{x}}_{*}$ as $\boldsymbol{x}_{*}$ . In this process, we use the parity (2.6a ) and (2.6c ) for $h(\boldsymbol{x}_{*}) = [\omega (\boldsymbol{x}_{*}) \ u_{1}(\boldsymbol{x}_{*}) \ u_{2}(\boldsymbol{x}_{*}) \ \tau (\boldsymbol{x}_{*})]^{\mathrm{T}}$ in $x_{1*}$ , which implies that $h(\boldsymbol{x}_{*}) = S_{h} \circ h(\bar {\boldsymbol{x}}_{*})$ for $\boldsymbol{x}_{*} \in \varOmega _{\textit{L}}$ , as well as the corresponding property of the kernel function $K$ , $K(\boldsymbol{x}_{*}-\boldsymbol{x}) = S_{K} \circ K(\bar {\boldsymbol{x}}_{*}-\bar {\boldsymbol{x}})$ for $\boldsymbol{x}_{*} \in \varOmega _{\textit{L}}$ . The expressions on the most right-hand sides of (2.8f ), (2.8j ) and (2.8m ) are obtained in the same way. The functions $J_{n}$ , appearing in various places, are the Abramowitz functions, defined by (2.8d ). The source terms $u_{w} a(\boldsymbol{x})$ and $u_{w} b(\boldsymbol{x})$ are given as the integrals over the lid surface $-1/2 \lt x_{1*} \lt 1/2$ , $x_{2*} = 1/2$ (see (2.8i ) and (2.8n )). Since the integration domains in $A$ , $B$ , $A_{w}$ and $B_{w}$ are restricted to the right half of the cavity, $\varOmega _{\textit{R}}$ and $\partial \varOmega _{\textit{R}}$ , the integral equations (2.7) can be solved within $\varOmega _{\textit{R}}$ and $\partial \varOmega _{\textit{R}}$ . In (2.7), the VDF and the molecular velocity $\boldsymbol{\zeta }$ are eliminated, so that the discontinuity of the VDF mentioned at the beginning of this subsection does not manifest itself. Hence, (2.7) is advantageous for accurate numerical analysis. In comparison with Hattori (Reference Hattori2024), where a single equation for the spanwise velocity component $u_{3}$ is obtained, the present formulation leads to coupled equations for $h(\boldsymbol{x})$ and $\sigma (\boldsymbol{x})$ , with additional quantities such as $r_{i}$ and $n_{i}$ defined in (2.8c ) and (2.2c ).

We remark that the functions $a(\boldsymbol{x})$ and $b(\boldsymbol{x})$ in (2.7) are direction dependent at the upper right corner $(x_{1},x_{2}) = (1/2, 1/2)$ . That is, the limiting values of $a$ and $b$ depend on the direction from which the corner is approached, which can be verified from their expressions. Then, as expected from the structure of (2.7), the solution $h(\boldsymbol{x}) = [\omega (\boldsymbol{x}) \ u_{1}(\boldsymbol{x}) \ u_{2}(\boldsymbol{x}) \ \tau (\boldsymbol{x})]^{\mathrm{T}}$ and $\sigma (\boldsymbol{x})$ can also exhibit direction-dependent limiting values at the corner. To address this property related to the corner appropriately, we introduce the following plane polar coordinates, $(r_{+},\varphi _{+})$ and $(r_{-},\varphi _{-})$ , centred at the corners (see figure 2):

(2.9) \begin{equation} x_{1} = \frac {1}{2} - r_{\pm }\cos \varphi _{\pm }, \quad x_{2} = \pm \frac {1}{2} \mp r_{\pm }\sin \varphi _{\pm } \quad \big(0 \le \varphi _{\pm } \le \pi /2\big). \end{equation}

Here, $r_{\pm }$ is the distance from the upper or lower right corner $(1/2,\pm 1/2) \equiv \boldsymbol{x}_{\textit{c}}^{\pm }$ to the point $\boldsymbol{x}$ , and $\varphi _{\pm }$ is the angle of the vector $\boldsymbol{x}-\boldsymbol{x}_{\textit{c}}^{\pm }$ measured from $x_{2} = \pm 1/2$ towards $x_{1} = 1/2$ . Using these coordinates, for example, the limiting values of $u_{i}(\boldsymbol{x})$ when the upper right corner is approached can be written as $\lim _{r_{+}\to 0}u_{i}(\boldsymbol{x}(r_{+},\varphi _{+}))$ with $\varphi _{+}$ being fixed. In this way, any potential direction dependence of macroscopic quantities is properly handled as their dependence on the azimuthal angle $\varphi _{\pm }$ .

Figure 2.

Plane polar coordinates with its origin at a corner: (a) ( $r_{+}$ , $\varphi _{+}$ ), and (b) ( $r_{-}$ , $\varphi _{-}$ ).

By substituting the expression of the VDF $\phi$ in terms of the density $\omega$ , flow velocity components $u_{1}$ and $u_{2}$ , temperature $\tau$ and the function $\sigma$ , mentioned before (2.7), into (2.3b ), we obtain the following expression for the stress components $P_{11}$ , $P_{22}$ and $P_{12}$ :

(2.10) \begin{align} &P_{\textit{ij}}(\boldsymbol{x}) \nonumber \\ &= \frac {2}{k\pi }\int _{\varOmega }\frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|} \big \{ r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})r_{j}(\boldsymbol{x}_{*}-\boldsymbol{x})J_{2}(|\boldsymbol{x}_{*}-\boldsymbol{x}|/k)\omega (\boldsymbol{x}_{*}) \nonumber \\ &\quad -2r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})r_{j}(\boldsymbol{x}_{*}-\boldsymbol{x})r_{\ell }(\boldsymbol{x}_{*}-\boldsymbol{x})J_{3}(|\boldsymbol{x}_{*}-\boldsymbol{x}|/k)u_{\ell }(\boldsymbol{x}_{*}) \nonumber \\ &\quad + r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})r_{j}(\boldsymbol{x}_{*}-\boldsymbol{x})\big [J_{4}(|\boldsymbol{x}_{*}-\boldsymbol{x}|/k)-J_{2}(|\boldsymbol{x}_{*}-\boldsymbol{x}|/k)\big ]\tau (\boldsymbol{x}_{*}) \big \}\mathrm{d} \boldsymbol{x}_{*} \nonumber \\ &\quad -\frac {2}{\pi }\int _{\partial \varOmega } \frac {1}{|\boldsymbol{x}_{*}-\boldsymbol{x}|} r_{i}(\boldsymbol{x}_{*}-\boldsymbol{x})r_{j}(\boldsymbol{x}_{*}-\boldsymbol{x})J_{3}(|\boldsymbol{x}_{*}-\boldsymbol{x}|/k) r_{\ell }(\boldsymbol{x}_{*}-\boldsymbol{x})n_{\ell }(\boldsymbol{x}_{*})\sigma (\boldsymbol{x}_{*}) \mathrm{d} l(\boldsymbol{x}_{*}) \nonumber \\ &\quad -\frac {4u_{w}}{\pi }\int _{-1/2}^{1/2} r_{i}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)r_{j}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)r_{1}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right)r_{2}\left(\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right) \nonumber \\ &\quad \times J_{4}\left(\left|\boldsymbol{x}_{*}^{(U)}-\boldsymbol{x}\right|/k\right) \mathrm{d} x_{1*}. \end{align}

Here, $(i,j) = (1,1)$ , $(2,2)$ and $(1,2)$ . Once $h(\boldsymbol{x}) = [\omega (\boldsymbol{x}) \ u_{1}(\boldsymbol{x}) \ u_{2}(\boldsymbol{x}) \ \tau (\boldsymbol{x})]^{\mathrm{T}}$ and $\sigma (\boldsymbol{x})$ are determined by solving (2.7), $P_{11}$ , $P_{22}$ and $P_{12}$ follow from (2.10). Note that the present formulation also allows the stress component $P_{33}$ and the heat-flow vector to be obtained. However, for brevity, their expressions and results are omitted in the present paper.

The integral equations (2.7) for the macroscopic quantities $h(\boldsymbol{x})$ and $\sigma (\boldsymbol{x})$ are solved by successive approximation (Hasegawa & Sone Reference Hasegawa and Sone1988; Aoki et al. Reference Aoki, Sone and Yano1989; Sone Reference Sone2007; Taguchi & Charrier Reference Taguchi and Charrier2008). Specifically, starting from some initial guesses $h^{0}(\boldsymbol{x})$ and $\sigma ^{0}(\boldsymbol{x})$ , sequences $\{h^{n}\}$ and $\{\sigma ^{n}\}$ are constructed by the recurrence relations

(2.11a) \begin{align} h^{n}(\boldsymbol{x}) &= A[h^{n-1}](\boldsymbol{x}) + A_{w}[\sigma ^{n-1}](\boldsymbol{x}) + u_{w} a(\boldsymbol{x}) \quad (\boldsymbol{x} \in \varOmega _{\textit{R}};\ n = 1, 2, \ldots ), \end{align}

(2.11b) \begin{align} \sigma ^{n}(\boldsymbol{x}) &= B[h^{n-1}](\boldsymbol{x}) + B_{w}[\sigma ^{n-1}](\boldsymbol{x}) + u_{w} b(\boldsymbol{x}) \quad (\boldsymbol{x} \in \partial \varOmega _{\textit{R}};\ n = 1, 2, \ldots ), \end{align}

which correspond to (2.7). The solution $h(\boldsymbol{x})$ and $\sigma (\boldsymbol{x})$ is obtained as the limit of the sequences. In actual numerical computations based on (2.11), the integrals $A[h^{n-1}]$ , $B[h^{n-1}]$ , $A_{w}[\sigma ^{n-1}]$ and $B_{w}[\sigma ^{n-1}]$ need to be evaluated numerically using the known data $h^{n-1}$ and $\sigma ^{n-1}$ at the grid points. This is done by extending the procedure described in Hattori (Reference Hattori2024). Information on the numerical method, grid system and computational accuracy is summarised in Appendix A.1.

4.1. Flow velocity, stress and pressure fields over whole cavity

Figure 3.

Plots of the flow velocity fields and the streamlines in the cavity: (a) $k = 0.01$ , (b) $k = 1$ and (c) $k = 10$ . The arrows indicate the flow velocity components $(u_{1}/u_{w},u_{2}/u_{w})$ at their starting points with their scale being shown on the left shoulder of each panel. At the four corners $(\pm 1/2, \pm 1/2)$ , two values, that approached from the upper or lower boundary and that approached from the side walls are shown although they are not visible for the lower corners. Its close-up around the upper right corner is shown on the top of each panel, where the two values are indicated by arrows with filled and hollow arrowheads, respectively, the corner by a bullet, and the lid and right wall by solid lines. The solid lines represent the streamlines originating from the points $(0,0.15)$ , $(0,0)$ , $(0.1,-0.1)$ , $(0.2,-0.2)$ , $(0.3,-0.3)$ , $(0.4,-0.4)$ , $(0.45,-0.45)$ , $(0.49,-0.49)$ and $(-0.49,-0.49)$ .

Figure 3 shows the vector fields of the flow velocity $(u_{1}/u_{w},u_{2}/u_{w})$ and the streamlines for $k=0.01$ , $1$ and $10$ over the whole cavity cross-section. The quantities in the left half of the cavity are recovered from the numerical results in the right half by using the parity relations (2.6a ) and (2.6c ). As defined in (2.1e ), $k = (\sqrt {\pi }/2)\textit{Kn}$ , and hereafter the parameter $k$ is used instead of the Knudsen number $\textit{Kn}$ . As the figure indicates, a flow is induced by the lid motion and a clockwise vortex appears. For smaller $k$ (figure 3 a), the flow speed is faster. As $k$ increases, the velocity slip on the lid ( $x_{2} = 1/2$ ), i.e. the difference between the lid speed and the gas flow velocity, becomes larger, and the overall flow speed decreases. The position of the vortex centre shifts slightly downward with increasing $k$ . For $k=1$ and $10$ , judging from the drawn streamlines, only a single vortex is present in the cavity (figures 3 b and 3 c). On the other hand, for $k = 0.01$ , in addition to the primary vortex, a second small counter-rotating vortex is observed near the lower corners $(\pm 1/2, -1/2)$ (figure 3 a). This vortex corresponds to the Moffatt’s corner vortex in fluid mechanics. Considering the possible direction-dependent limiting values, at each corner, two representative flow velocity vectors are shown in figure 3: one with the corner being approached along the upper or lower boundary and one with the corner being approached along the side walls. At the lower corners the arrows are not visible. At the upper right corner, as shown in the close-ups, the vector points to the right along the lid when approached from the lid side, whereas it points downward along the right wall when approached from the right wall. This indicates that the flow velocity at the upper corners is actually direction dependent for various Knudsen numbers. The corner vortices near the lower corners and the direction dependence of the macroscopic quantities are further examined later in §§ 4.2 and 4.3.

Figure 4.

Contour plots of the stress components $P_{12}/u_{w}$ , $P_{11}/u_{w}$ and $P_{22}/u_{w}$ and the pressure $P/u_{w}$ on the $x_{1}x_{2}$ plane for (a–d) $k = 0.01$ , (e–h) $k = 1$ , and (i–l) $k = 10$ . The isolines are drawn with intervals of $0.15$ in (b), (f) and (j), while they are drawn with intervals of $0.1$ in the other panels.

Figure 4 shows contour plots of the stress components $P_{12}/u_{w}$ , $P_{11}/u_{w}$ and $P_{22}/u_{w}$ and the pressure $P/u_{w}$ for $k=0.01$ , $1$ and $10$ . In all cases, their magnitudes are particularly large on the lid and side walls. For non-small Knudsen numbers ( $k=1$ and $10$ ), the variations of the macroscopic quantities are relatively mild and are observed over the whole cavity (figure 4 e–l). In contrast, when the Knudsen number is small ( $k=0.01$ ), they vary steeply near the upper corners (figure 4 a–d). The behaviour of the stress, pressure and flow velocity near the corners for small Knudsen numbers is further examined later in § 4.4. For $k=0.01$ , the magnitudes of the stress and pressure in the bulk away from the boundaries are smaller than those for $k=1$ and $10$ . This is consistent with the fact that their magnitudes, determined from the Stokes equation (3.1b ), Newton’s law (3.2) and a smooth $O(1)$ velocity distribution $u_{i}/u_{w}$ in the bulk, are $O(k)$ and are thus small.

Before closing this subsection, let us mention thermal characteristics of the gas briefly. In transition and highly rarefied regimes, the temperature $\tau /u_{w}$ is higher near the upper right corner and lower near the upper left corner. The heat flows from the low-temperature region to the high-temperature one, which is the so-called anti-Fourier heat transfer (see, e.g. John et al. Reference John, Gu and Emerson2010, Reference John, Gu and Emerson2011; Mohammadzadeh et al. Reference Mohammadzadeh, Roohi, Niazmand, Stefanov and Myong2012), and the present results confirm its existence even within the linearised regime. As the Knudsen number decreases, the temperature variations become increasingly confined near the upper corners. The heat-flow streamlines exhibit a more complex pattern, including the formation of a small vortex near the bottom wall. The behaviour of heat flow near the boundary is consistent with a shear-induced heat flow within the Knudsen layer. For details, including the visualised fields, the reader is referred to supplementary material S1 available at https://doi.org/10.1017/jfm.2026.11787.

4.2. Corner vortices

In conventional fluid mechanics, it is theoretically shown that corner vortices appear as an infinite series of vortices that rotate in alternating directions, with their size decreasing and their strength decaying rapidly (Moffatt vortices; see Moffatt Reference Moffatt1964). For the Stokes flow in a square cavity corresponding to the present setting, within the achievable computational accuracy, multiple corner vortices near the lower corners are actually confirmed numerically (see, e.g. Pan & Acrivos Reference Pan and Acrivos1967; Shankar Reference Shankar1993; Biswas & Kalita Reference Biswas and Kalita2018), which supports the theoretical statement. Here, we examine the effect of a finite Knudsen number on this property, specifically whether corner vortices exist for a given Knudsen number and, if they do, whether their number is finite or infinite. If figure 3 is interpreted naively, one would conclude that a single corner vortex exists for $k = 0.01$ , whereas no corner vortices are present for $k = 1$ and $10$ . However, since this argument may depend on the scale of the figure and the choice of the starting points of streamlines, we take an alternative viewpoint below.

Figure 5.

Sketch of streamlines near the lower right corner $(1/2, -1/2)$ with the presence of a single corner vortex being assumed. The points A, B, C and D chosen in the figure are specified by the same azimuthal angle $\varphi _{-}$ and by their respective distances $r_{-}$ from the corner (see also figure 2 b). The azimuthal velocity component $u_{\varphi _{-}}$ will be negative at points A and B in the primary vortex and at point C on the cavity-centre side of the corner vortex, while it will be positive at point D closer to the corner.

Figure 5 is a schematic illustration of streamlines near the lower right corner $(1/2, -1/2)$ under an assumption that a single corner vortex is present. As indicated in the figure, if there exist vortices rotating in the opposite direction, the sign of the azimuthal velocity component $u_{\varphi _{-}}$ is expected to reverse, near the corner-facing side of the vortex closer to the corner, e.g. at point D in the figure, relative to its sign at other locations. Here, $u_{\varphi _{-}} = u_{i}e_{i}^{(\varphi _{-})}$ with $(e_{1}^{(\varphi _{-})},e_{2}^{(\varphi _{-})}) = (\sin \varphi _{-},\cos \varphi _{-})$ being the unit vector in the $\varphi _{-}$ direction.

Figure 6.

Profiles of the azimuthal velocity component near the lower right corner $(1/2, -1/2)$ for $\varphi _{-} = \pi /4$ . The quantity $u_{\varphi _{-}}/u_{w}$ is plotted as a function of the distance $r_{-}$ from the lower right corner. Panels (a–d) show the profiles in intervals $0 \lt r_{-} \le 0.25$ , $0 \lt r_{-} \le 0.08$ , $0 \lt r_{-} \le 0.004$ and $0 \lt r_{-} \le 0.004$ , respectively. Results for different sets of $k$ are shown in each panel. Note that the vertical scale differs from panel to panel.

With this in mind, figure 6 shows plots of $u_{\varphi _{-}}$ for various $k$ as a function of the distance $r_{-}$ , with the azimuthal angle being fixed at $\varphi _{-} = \pi /4$ . First, for $k = 10$ and $1$ , $u_{\varphi _{-}}$ is clearly negative for all $r_{-}$ , including the corner $r_{-}\to 0$ (figure 6 a). This implies that, under the present setting of a square cavity closed by a moving lid, corner vortices are not observed in highly rarefied to intermediate regimes, in contrast to the fluid-dynamic regime. The feature, however, depends on geometry and setting; indeed, corner vortices in highly rarefied and free-molecular regimes are reported for gas flows in an open triangular cavity (Ben-Adva & Manela Reference Ben-Adva and Manela2026) and over a backward-facing step (Manela & Gibelli Reference Manela and Gibelli2020).

On the scale of figure 6(a), for $k \le 0.1$ , $u_{\varphi _{-}}$ approaches zero as $r_{-}$ decreases. Thus, we show a close-up in figure 6(b). While $u_{\varphi _{-}}$ remains negative over the entire range for $k = 0.1$ and $0.03$ , for $k = 0.02$ , there exists an interval over which $u_{\varphi _{-}}$ is positive. This implies, as depicted in figure 5, the emergence of a corner vortex at this Knudsen number. As $k$ decreases further, the profiles approach the Stokes solution. Finally, figures 6(c) and 6(d) show further close-ups. Figure 6(c) indicates that, for $k$ smaller than $k = 0.02$ and $0.015$ , corresponding to the transitional stage of corner-vortex growth, $u_{\varphi _{-}}$ at each position $r_{-}$ decreases monotonically as $k$ decreases. For $k = 0.003$ and $0.002$ , $u_{\varphi _{-}}$ is still positive over the entire plotted range $0 \lt r_{-} \le 0.004$ (figure 6 d). On the other hand, for $k \le 0.0015$ , there appears an interval of $r_{-}$ over which $u_{\varphi _{-}}$ becomes negative again. This implies the emergence of a second corner vortex at these Knudsen numbers.

Figure 7.

Streamlines near the lower right corner $(1/2, -1/2)$ for various $k$ . Panels (a–d) and panels (i–l) show the region $0.4 \le x_{1} \le 0.5$ and $-0.5 \le x_{2} \le -0.4$ , while panels (e–h) and panels (m–p) provide close-ups of the region $0.494 \le x_{1} \le 0.5$ and $-0.5 \le x_{2} \le -0.494$ . The values of $k$ are $0.03$ in (a,e), $0.02$ in (b,f), $0.01$ in (c,g), $0.002$ in (d,h), $0.0015$ in (i,m), $0.001$ in (j,n), $0.0006$ in (k,o) and $0.0003$ in (l,p). Panels (a–d) and panels (i–l) show streamlines originating from ten points $(m,-m)$ , where $m = 0.42$ , $0.4325$ , $0.445$ , $0.448$ , $0.468$ , $0.475$ , $0.485$ , $0.495$ , $0.4975$ and $0.4977$ . Panels (e–h) and panels (m–p) show streamlines originating from seven points $(m,-m)$ , where $m = 0.495$ , $0.4957$ , $0.4964$ , $0.4972$ , $0.4986$ , $0.4994$ and $0.4999$ .

The behaviour of the azimuthal velocity $u_{\varphi _{-}}$ as a function of the distance $r_{-}$ from the lower right corner has been examined in figure 6. We now show in figure 7 the streamlines near the corner for various small values of $k$ . As seen below, the observed features are consistent with those inferred from figure 6. For $k = 0.03$ , the streamlines are not closed within the displayed region, and no corner vortex is observed (figures 7 a and 7 e). On the other hand, for $k = 0.02$ , a first corner vortex is observed, as shown in figure 7(b). Figures 7(b) and 7(f) indicate that another small vortex can also be seen very close to the corner. However, this is regarded as a remnant of the primary clockwise vortex, which originally occupied the entire region, during the growth process of the first corner vortex, rather than as a fully developed second corner vortex. For $k = 0.01$ and $0.002$ , the first corner vortex further develops and increases in size (figures 7 c and 7 d); in the latter case, it reaches approximately the same size as that for smaller $k$ (compare figure 7 d with figure 7 i–l). In these cases, no second corner vortex is observed. For $k = 0.0015$ , a structure corresponding to a second corner vortex is just beginning to emerge (figure 7 m). For $k = 0.001$ , $0.0006$ and $0.0003$ , a second corner vortex is clearly observed and grows in size (figure 7 n–p).

Before closing this subsection, we summarise the results. To examine the presence and structure of corner vortices, we have focused on the distribution of the azimuthal velocity (figure 6) and the streamlines (figure 7). Figure 6(a–d) indicates that, for each fixed Knudsen number, the azimuthal flow velocity component retains a constant sign over a certain range close to the corner and does not exhibit sign changes with respect to the distance $r_{-}$ . Thus, unless one considers the limit $k \to 0$ , the number of corner vortices is finite within kinetic theory. This is in contrast to the Stokes solution, for which their number is theoretically known to be infinite. According to figures 6 and 7, zero, one or two corner vortices are observed over the range of Knudsen numbers examined ( $k \ge 0.0003$ ). The presence or absence of each corner vortex at a given Knudsen number can be understood in terms of the size relationship between the mean free path and the size of each corner vortex. When the former is not small compared with the latter, rarefaction effects, such as velocity slip, allow the gas to flow tangentially along the walls, and the flow tends not to separate from the walls, suppressing the formation of corner vortices. On the other hand, when the mean free path is sufficiently small compared with the corner-vortex size, such rarefaction effects are suppressed at that scale, and the corner vortex, which is a fluid-dynamic behaviour associated with viscosity, no-slip on the walls and the geometry with a sharp corner, can develop. Our numerical results support this interpretation. Actually, no corner vortices are observed for large Knudsen numbers down to $k = 0.03$ . For $k \le 0.02$ , a first corner vortex, which is of the size slightly less than $0.1$ (see, e.g. figure 7 l), appears. Moreover, for $k \le 0.001$ , a second corner vortex, which is of the size slightly less than $0.006$ (see, e.g. figure 7 p), is clearly observed. Note that, when the azimuthal angle is chosen as $\varphi _{-} = \pi /8$ and $3\pi /8$ , the same properties observed for $\varphi _{-} = \pi /4$ in figure 6(a–d) can be confirmed. Note that, unlike the other figures, figures 6(d) and 7(l–p) contain numerical errors associated with the scale of the plots, as they focus on quantities of extremely small magnitude and streamlines corresponding to very small velocities near the corner. To allow readers to assess the magnitude of these errors, corresponding figures obtained using a coarser grid are provided in supplementary material S2; compare figures 6 and 7 and figures S.3 and S.4 in the supplementary material. Although differences are visible, they are of a magnitude that does not affect the statements made in this subsection.

4.3. Direction dependence of limiting values of macroscopic quantities at the corners

Figure 8.

Profiles of flow velocity components and shear stress for the BGK and Stokes solutions with $k = 0.01$ : (a) $u_{1}/u_{w}$ , (b) $u_{2}/u_{w}$ , and (c) $P_{12}/u_{w}$ . The quantities are plotted as functions of the azimuthal angle $\varphi _{+}$ with $r_{+}$ , the distance from the upper right corner $(1/2,1/2)$ , being fixed. The solid, dashed, dash-dotted, dash-double-dotted, long-dashed and dotted lines indicate the results for $r_{+} = 0.2025$ , $0.02260$ , $2.771 \times 10^{-3}$ , $2.770 \times 10^{-4}$ , $2.778 \times 10^{-5}$ and the limiting values as $r_{+} \to 0$ , respectively.

To directly examine whether the macroscopic quantities at the upper corner exhibit direction-dependent limiting values, we use the polar coordinates shown in figure 2(a) and plot the quantities as functions of the azimuthal angle $\varphi _{+}$ with the radial distance $r_{+}$ being fixed. If a quantity is not constant with respect to $\varphi _{+}$ as $r_{+} \to 0$ , it is direction dependent. Figures 8(a) and 8(b) show such plots of the flow velocity components $u_{1}$ and $u_{2}$ , respectively, for both the BGK and Stokes solutions with $k = 0.01$ . In these figures, the positions on the lid and right wall correspond to the angles $\varphi _{+}/\pi = 0$ and $\varphi _{+}/\pi = 0.5$ , respectively. The fact that $u_{1}$ on the right wall and $u_{2}$ on the lid are identically zero for both solutions indicates that the flow velocity components normal to the boundaries vanish, which is consistent with the boundary conditions (3.1c ) and (2.1f ). For the Stokes solution, $u_{1}$ on the lid and $u_{2}$ on the right wall coincide with the boundary velocities $1$ and $0$ , respectively, owing to the no-slip condition (3.1c ). In contrast, for the BGK solution, $u_{1}$ on the lid and $u_{2}$ on the right wall deviate from the boundary velocities, namely, velocity slip appears. The magnitude of slip increases as the distance $r_{+}$ from the upper right corner decreases. Owing to the velocity slip, the BGK profiles are flatter than the Stokes ones. Nevertheless, they approach a non-uniform profile as the corner is approached ( $r_{+} \to 0$ ). This implies that the flow velocity exhibits direction-dependent limiting values not only in the Stokes solution but also in the kinetic solution. We additionally show here the shear stress $P_{12}$ in figure 8(c), which implies that the direction dependence in the BGK model is observed also for other macroscopic quantities. Note that the influence of the Knudsen number on the direction dependence at the upper right corner is relatively small. Indeed, for $k = 10$ and $k = 1$ , the profiles at the corner $r_{+} \to 0$ differ quantitatively only slightly from those in figure 8 for $k = 0.01$ , and they are qualitatively similar. For $k \le 0.1$ , the profiles also agree quantitatively on the scale of figure 8. (For brevity, the corresponding figures are omitted.) The direction dependence can be more clearly understood in terms of the VDF $\phi$ . The VDF at a given point is determined along straight characteristic lines of the governing equation, taking into account the boundary data. For molecular velocities directed from the upper right corner into the interior, $\zeta _{1} \lt 0$ and $\zeta _{2} \lt 0$ , the VDF tends to the boundary data at the corner as $r_{+} \to 0$ . The selected side of the corner, either $(1/2-0, 1/2)$ (approached along the lid) or $(1/2, 1/2-0)$ (approached along the right wall), depends on the azimuthal angle $\varphi _{+}$ measured from the corner, even when the corner is approached ( $r_{+} \to 0$ ). Since the boundary data on the lid and the right wall differ, in particular due to the source term $2\zeta _{1}u_{w}$ associated with the lid motion, the limiting values selected from the two sides are inherently different at the corner. As a result, the macroscopic quantities, which are moments of the VDF, exhibit direction-dependent limiting values. Furthermore, since the VDF for $\zeta _{1} \lt 0$ and $\zeta _{2} \lt 0$ is determined without explicit contribution from intermolecular collisions, this mechanism is essentially independent of the Knudsen number. See Appendix C for details.

A natural question that arises next is whether the macroscopic quantities are also direction dependent at the lower corner. To examine this, figure 9 shows their limiting values as $r_{-} \to 0$ , plotted as functions of the azimuthal angle $\varphi _{-}$ for various $k$ . In all panels, the profiles become flatter as $k$ decreases, and for $k = 0.01$ they are almost horizontal on the scale of the figure. Thus, for small Knudsen numbers, even if the BGK solution is direction dependent, its degree is very slight. In contrast, for larger $k$ , all quantities clearly vary with $\varphi _{-}$ , exhibiting significant direction dependence. Hence, the BGK solution exhibits direction-dependent limiting values also at the lower corner, and its degree becomes more significant as the Knudsen number increases. To understand this behaviour, we consider the underlying mechanism. At the lower corner, since both walls are at rest, unlike at the upper corner, there is no direct mechanism that introduces a discontinuity in the boundary data (2.1f ). Nevertheless, the function $\sigma (x_{1}, x_{2})$ , which is determined together with the VDF $\phi$ , can become discontinuous at the corner. It corresponds to the mass flux of molecules with $\zeta _{j}n_{j} \lt 0$ impinging on the boundary. Owing to the discontinuous change in the direction of the normal vector $n_{i}$ at the corner, molecules with the same velocity $\zeta _{i}$ contribute differently to the mass fluxes on the two sides of the lower right corner. As a result, $\sigma (1/2,-1/2+0)$ and $\sigma (1/2-0,-1/2)$ can differ. Then, molecules with velocities directed from the lower right corner into the interior, $\zeta _{1} \lt 0$ and $\zeta _{2} \gt 0$ , give rise to direction-dependent limiting values through the same mechanism as at the upper corner. This behaviour is observed for large $k$ . For small $k$ , the gas state near the lower right corner approaches an equilibrium state at rest with the cavity-wall temperature, corresponding to an isotropic distribution in $\zeta _{i}$ . Since the difference between $\sigma (1/2,-1/2+0)$ and $\sigma (1/2-0,-1/2)$ vanishes in this equilibrium state, it becomes small for small $k$ , and the resulting direction dependence is correspondingly weak. Further details are provided in Appendix C.

Figure 9.

Profiles of flow velocity components and shear stress for the BGK solution as the lower right corner is approached ( $r_{-} \to 0$ ): (a) $u_{1}/u_{w}$ , (b) $u_{2}/u_{w}$ , and (c) $P_{12}/u_{w}$ . The quantities are plotted as functions of the azimuthal angle $\varphi _{-}$ .

As shown in the previous two paragraphs, at both the upper and lower corners, the direction-dependent limiting values of macroscopic quantities originate from the discontinuity of the boundary data (2.1f ) for the VDF $\phi$ at the corners. In fact, in Hattori (Reference Hattori2024) the boundary condition at the lower corner is continuous, with $\phi =0$ on both walls, and the macroscopic quantities exhibit no direction dependence there.

4.4. Flow velocity, stress and pressure fields near the upper corner for small Knudsen numbers

Figure 10.

Flow speed and streamlines near the upper right corner $(1/2,1/2)$ . Here $k = 0.01$ . (a) Stokes and (b) BGK solutions. The solid lines represent the streamlines originating from the points $(0.46,0.46)$ , $(0.47,0.47)$ , $(0.48,0.48)$ , $(0.49,0.49)$ and $(0.496,0.496)$ . The contour plots with the dashed isolines drawn with intervals of $0.1$ represent the flow speed $|\boldsymbol{u}|/u_{w}$ .

Figure 10 shows contour plots of the flow speed and streamlines near the upper right corner $(1/2, 1/2)$ for $k = 0.01$ . For the Stokes solution, the isolines of the flow speed extend radially from the corner (figure 10 a). Since the speed is constant along the boundaries owing to the no-slip condition, no isolines intersect either the lid ( $x_{2} = 1/2$ ) or the right wall ( $x_{1} = 1/2$ ). In contrast, for the BGK solution, several isolines intersect the boundaries (figure 10 b). The streamlines originating from the same points also differ between the two solutions. In this way, although the chosen Knudsen number is fairly small ( $k = 0.01$ ), the kinetic (BGK) solution yields results that are clearly distinct from those of the Stokes solution near the corner.

Figure 11.

Contour plots of the stress components $P_{12}/u_{w}$ , $P_{11}/u_{w}$ and $P_{22}/u_{w}$ and the pressure $P/u_{w}$ near the upper right corner $(1/2,1/2)$ . Here $k = 0.01$ . (a–d) Stokes and (e–h) BGK solutions. The isolines are drawn with intervals of $0.1$ in (a) and (e), while they are drawn with intervals of $0.2$ in the other panels.

Figure 11(a–d) shows the profiles of the stress components $P_{12}/u_{w}$ , $P_{11}/u_{w}$ and $P_{22}/u_{w}$ and the pressure $P/u_{w}$ near the upper right corner for the Stokes solution. In figure 11(b–d) the isolines closest to the corner correspond to those with the maximum value among the drawn lines. In the regions bounded by these isolines and the boundaries, $P_{11}$ , $P_{22}$ and $P$ do not saturate but diverge as the corner is approached, as explained in the last paragraph of § 3. In figure 11(a) the isolines with the maximum value $0$ and minimum value $-0.8$ penetrate into the corner, and the region associated with intermediate values becomes narrower as the corner is approached, which is consistent with the divergence of $P_{12}$ .

In contrast, for the BGK solution shown in figure 11(e–h), the stress components and pressure remain finite. The magnitudes of $P_{12}/u_{w}$ , $P_{11}/u_{w}$ , $P_{22}/u_{w}$ and $P/u_{w}$ do not exceed $0.782$ , $2.280$ , $1.817$ and $2.013$ , respectively, over the whole domain including the corner. This implies that the problem of divergence of stress and pressure at the corner for the Stokes solution (the corner singularity) in the lid-driven cavity flow is resolved within kinetic theory. Since here the result for only one specific case $k = 0.01$ is discussed, we next examine the gas behaviour for smaller Knudsen numbers in more detail.

Figure 12.

Profiles of the tangential flow velocity component and shear stress on the lid and right wall. Here $k = 0.1$ , $0.01$ and $0.001$ . Results are shown for (a) $u_{1}/u_{w}$ on the lid ( $\varphi _{+} = 0$ ), (b) $u_{2}/u_{w}$ on the right wall ( $\varphi _{+} = \pi /2$ ), (c) $P_{12}/u_{w}$ on the lid ( $\varphi _{+} = 0$ ), and (d) $P_{12}/u_{w}$ on the right wall ( $\varphi _{+} = \pi /2$ ). The solid and dashed lines indicate the BGK and Stokes solutions, respectively. The quantities are plotted as functions of $r_{+}$ , the distance from the upper right corner $(1/2, 1/2)$ .

Figure 12 shows the profiles of the tangential flow velocity components and shear stress on the lid and right wall for $k = 0.1$ , $0.01$ and $0.001$ as functions of the distance $r_{+}$ , where these surfaces correspond to the directions $\varphi _{+} = 0$ and $\pi /2$ , respectively (see figure 2 a). From figure 12(a) (or figure 12 b), it is seen that $u_{1}/u_{w}$ on the lid (or $u_{2}/u_{w}$ on the right wall) approaches the Stokes solution values $1$ (or $0$ ) prescribed by the no-slip condition in the overall region as $k$ decreases. However, as the corner is approached ( $r_{+} \to 0$ ), the velocity approaches a different common value $0.35\ldots$ (or $-0.26\ldots$ ) for the cases with any $k$ in the figure. According to figures 12(c) and 12(d), the shear stress $P_{12}/u_{w}$ agrees well with the Stokes solution of $O(k)$ for large $r_{+}$ when $k = 0.01$ and $0.001$ , whereas the BGK solution approaches the values $-0.78\ldots$ for $\varphi _{+} = 0$ and $-0.38\ldots$ for $\varphi _{+} = \pi /2$ as $r_{+} \to 0$ for all cases shown in the figure. The results imply that both the velocity slip and the stress remain $O(1)$ no matter how small $k$ is. Since they are $O(k)$ away from the corner (see, e.g. the plots for $r_{+} \ge 0.3$ in each panel of figure 12), which is consistent with the usual magnitude of shear slip and with Newton’s law, the large velocity slip and large stress are features of the gas behaviour near the corner.

Figure 13.

Plots of the $x_{1}$ component of the flow velocity and the shear stress for $k = 0.03$ , $0.01$ , $0.003$ , $0.001$ and $0.0003$ . In (a) and (b), $u_{1}/u_{w}$ and $P_{12}/u_{w}$ on the lid ( $\varphi _{+} = 0$ ) are plotted, respectively, as functions of $r_{+}/k$ , the rescaled distance from the upper right corner $(1/2, 1/2)$ . In (c) and (d), contour plots of $u_{1}/u_{w}$ and $P_{12}/u_{w}$ vs $((1/2-x_{1})/k, (1/2-x_{2})/k)$ , the position relative to the corner rescaled by $k$ , are shown with the isolines being drawn with intervals $0.1$ and $0.05$ , respectively. The types of lines used are common in (a–d); see the line legend in (a). The Stokes solution is shown only in (a) and (b).

In figure 12 the cavity size is used as the length scale, and the variation of macroscopic quantities becomes steeper for smaller $k$ as the corner is approached ( $r_{+} \to 0$ ). This behaviour can be understood more clearly by introducing rescaled coordinates based on the mean free path. Accordingly, we plot the tangential flow velocity component $u_{1}/u_{w}$ and the shear stress $P_{12}/u_{w}$ on the lid ( $\varphi _{+}=0$ ) as functions of the rescaled distance $r_{+}/k$ in figures 13(a) and 13(b). Contour plots of the same quantities in a region near the corner as functions of the rescaled position $((1/2-x_{1})/k, (1/2-x_{2})/k)$ are further shown in figures 13(c) and 13(d). The results for different values of $k$ collapse onto a common profile in all panels. This provides evidence for the presence of a region near the corner in which the macroscopic quantities vary on the scale of the mean free path, indicating an asymptotic structure of the solution as $k \to 0$ and $r_{+} \to 0$ . Note that analogous scaling properties are also observed for the temperature and heat flow. For details, the reader is referred to the last paragraph and figure S.2 of supplementary material S1.

4.5. Total forces on lid and cavity walls and local stress profiles

Let $p_{0} L F_{i\textit{U}}$ , $p_{0} L F_{i\textit{B}}$ , $p_{0} L F_{i\textit{R}}$ and $p_{0} L F_{i\textit{L}}$ ( $i = 1, 2$ ) denote the $X_{i}$ components of the total forces per unit length in the $X_{3}$ direction exerted by the gas on the lid, bottom wall, right and left side walls, respectively. Then, $F_{i\textit{U}}$ , $F_{i\textit{B}}$ , $F_{i\textit{R}}$ and $F_{i\textit{L}}$ are given by

(4.1a) \begin{align} F_{1\textit{U}} &= \int _{-1/2}^{1/2} P_{12}\left (x_{1},\frac {1}{2}\right ) \mathrm{d} x_{1} = 2\int _{0}^{1/2} P_{12}\left (x_{1},\frac {1}{2}\right ) \mathrm{d} x_{1}, \end{align}

(4.1b) \begin{align} F_{2\textit{U}} &= \int _{-1/2}^{1/2} \left [1+P_{22}\left (x_{1},\frac {1}{2}\right )\right ] \mathrm{d} x_{1} = 1, \end{align}

(4.1c) \begin{align} F_{1\textit{B}} &= -\int _{-1/2}^{1/2} P_{12}\left (x_{1},-\frac {1}{2}\right ) \mathrm{d} x_{1} = -2\int _{0}^{1/2} P_{12}\left (x_{1},-\frac {1}{2}\right ) \mathrm{d} x_{1}, \end{align}

(4.1d) \begin{align} F_{2\textit{B}} &= -\int _{-1/2}^{1/2} \left [1+P_{22}\left (x_{1},-\frac {1}{2}\right )\right ] \mathrm{d} x_{1} = -1, \end{align}

(4.1e) \begin{align} F_{1\textit{R}} &= \int _{-1/2}^{1/2} \left [1+P_{11}\left (\frac {1}{2},x_{2}\right )\right ] \mathrm{d} x_{2} = 1 + \int _{-1/2}^{1/2} P_{11}\left (\frac {1}{2},x_{2}\right ) \mathrm{d} x_{2}, \end{align}

(4.1f) \begin{align} F_{2\textit{R}} &= \int _{-1/2}^{1/2} P_{12}\left (\frac {1}{2},x_{2}\right ) \mathrm{d} x_{2}, \end{align}

(4.1g) \begin{align} F_{1\textit{L}} &= -\int _{-1/2}^{1/2} \left [1+P_{11}\left (-\frac {1}{2},x_{2}\right )\right ] \mathrm{d} x_{2} = -1 + \int _{-1/2}^{1/2} P_{11}\left (\frac {1}{2},x_{2}\right ) \mathrm{d} x_{2} = F_{1\textit{R}}-2, \end{align}

(4.1h) \begin{align} F_{2\textit{L}} &= -\int _{-1/2}^{1/2} P_{12}\left (-\frac {1}{2},x_{2}\right ) \mathrm{d} x_{2} = -\int _{-1/2}^{1/2} P_{12}\left (\frac {1}{2},x_{2}\right ) \mathrm{d} x_{2} = -F_{2\textit{R}}. \end{align}

Here, in (4.1a ), (4.1c ) and (4.1h ) we have used the fact that the stress component $P_{12}$ is even in $x_{1}$ ; in (4.1b ) and (4.1d ), that $P_{22}$ is odd in $x_{1}$ ; and in (4.1g ), that $P_{11}$ is odd in $x_{1}$ (see (2.6c ) and (2.6a )). It follows from (4.1) that it suffices to obtain $F_{1\textit{U}}$ , $F_{1\textit{B}}$ , $F_{1\textit{R}}-1$ and $F_{2\textit{R}}$ to determine the total force acting on each surface.

Figure 14.

The forces acting on the lid and cavity walls. Plots related to $F_{1\textit{U}}$ , $F_{1\textit{B}}$ , $F_{1\textit{R}}-1$ and $F_{2\textit{R}}$ as functions of $k$ . Results are shown for (a) ( $-F_{1\textit{U}}/u_{w}$ , $-F_{1\textit{B}}/u_{w}$ , $(F_{1\textit{R}}-1)/u_{w}$ , $-F_{2\textit{R}}/u_{w}$ ) in the semilog plot, (b) ( $-F_{1\textit{U}}/u_{w}$ , $-F_{1\textit{B}}/u_{w}$ , $(F_{1\textit{R}}-1)/u_{w}$ , $-F_{2\textit{R}}/u_{w}$ ) in the log–log plot, and (c) ( $-F_{1\textit{U}}/u_{w} k$ , $-F_{1\textit{B}}/u_{w} k$ , $(F_{1\textit{R}}-1)/u_{w} k$ , $-F_{2\textit{R}}/u_{w} k$ ) in the semilog plot. Horizontal lines in (a) and (b) indicate the free-molecular solution, $(-F_{1\textit{U}}/u_{w}, -F_{1\textit{B}}/u_{w}, (F_{1\textit{R}}-1)/u_{w}, -F_{2\textit{R}}/u_{w})|_{k \to \infty } = (0.68533, 0.005235, 0.34528, 0.18723)$ .

Figure 14(a) shows $F_{1\textit{U}}/u_{w}$ , $F_{1\textit{B}}/u_{w}$ , $(F_{1\textit{R}}-1)/u_{w}$ and $F_{2\textit{R}}/u_{w}$ as functions of $k$ . The numerically obtained results for the free-molecular ( $k \to \infty$ ) solution are also shown. The free-molecular solution corresponds to the solution of an equation obtained from (2.7) by dropping $A[h]$ and $B[h]$ , which contain the factor $1/k$ , and by replacing the Abramowitz functions $J_{n}$ with $J_{n}(0)$ in $A_{w}[\sigma ]$ , $B_{w}[\sigma ]$ , $a$ and $b$ . While $-F_{1\textit{U}}$ , $(F_{1\textit{R}}-1)$ and $-F_{2\textit{R}}$ increase monotonically with $k$ (or the Knudsen number $\textit{Kn}$ ), $-F_{1\textit{B}}$ takes a maximum value around $k = 0.4$ .

To better illustrate the behaviour of the forces for small $k$ , the same quantities as in figure 14(a) are shown in figure 14(b) again but as a log–log plot. The total force $F_{1\textit{B}}$ agrees well with the Stokes solution, $F_{1\textit{B}}/u_{w} = -0.32k$ , which is superposed in the figure and is proportional to $k$ . The good agreement between the Stokes and BGK solutions for the force can also be confirmed at the level of the shear stress profile $P_{12}(x_{1}, -1/2)$ on the bottom wall (the corresponding figure is omitted for brevity). This implies that, for small $k$ , the influence of the corner on the gas behaviour is relatively weak on the bottom wall.

On the other hand, a closer look at figure 14(b) reveals that, unlike $F_{1\textit{B}}$ , the other forces $F_{1\textit{U}}$ , $(F_{1\textit{R}}-1)$ and $F_{2\textit{R}}$ are not simply proportional to $k$ . In fact, they behave as $k \ln k$ . To confirm this, we show a semi-logarithmic plot of $F_{1\textit{U}}/u_{w}$ , $F_{1\textit{B}}/u_{w}$ , $(F_{1\textit{R}}-1)/u_{w}$ and $F_{2\textit{R}}/u_{w}$ , each divided by $k$ , in figure 14(c). Such plots have been used in the literature to examine spatial singularities of macroscopic quantities (Takata & Taguchi Reference Takata and Taguchi2017). The quantities $-F_{1\textit{U}}/u_{w} k$ , $(F_{1\textit{R}}-1)/u_{w} k$ and $-F_{2\textit{R}}/u_{w} k$ increase as $k$ decreases, and the markers align along the straight lines in the plot. This demonstrates that $F_{1\textit{U}}$ , $(F_{1\textit{R}}-1)$ and $F_{2\textit{R}}$ actually behave as $k \ln k$ . The straight lines in the plot correspond to

(4.2a) \begin{align} \frac {F_{1\textit{U}}}{u_{w}} &\approx 4.226 k\ln k + 7.091k, \end{align}

(4.2b) \begin{align} \frac {F_{1\textit{R}}-1}{u_{w}} &\approx -2.113 k\ln k - 3.389k, \end{align}

(4.2c) \begin{align} \frac {F_{2\textit{R}}}{u_{w}} &\approx 1.338 k\ln k + 2.813k, \\[10pt] \nonumber \end{align}

which are determined by a least-squares fit to the numerical results for $k \le 0.01$ .

As discussed in Hattori (Reference Hattori2024), the $k \ln k$ factor in the total forces can be understood through the stress profiles along each boundary surface. Specifically, this factor arises because the stress continues to grow, in accordance with the fluid-dynamic description, inversely with the distance from the upper corner up to a few tens of mean free paths, after which it shifts to another bounded profile closer to the corner. Figures 15(a) and 15(b) show the shear stress profiles on the lid to confirm this explanation. In figure 15(a) the results for $k = 0.1$ , $0.01$ , $0.001$ and $0.0003$ are plotted as functions of the distance $r_{+}$ on a log–log scale. As $k$ decreases, the magnitude of the stress increases approximately inversely with $r_{+}$ over a wider interval from about $r_{+} = 0.2$ towards the corner. This corresponds to the fluid-dynamic behaviour of the stress explained in the last paragraph of § 3.

Figure 15.

Profiles of the shear stress and pressure on the lid. In (a), $P_{12}/u_{w}$ is shown as a function of $r_{+}$ , the distance from the upper right corner $(1/2, 1/2)$ . In (b) and (c), $P_{12}/u_{w}$ and $P/u_{w}$ are shown as functions of $r_{+}/k$ , the rescaled distance, respectively.

To clarify how close to the corner this growth persists, figure 15(b) shows the same data with the distance $r_{+}$ rescaled by $k$ . In this plot, the BGK profiles for different values of $k$ overlap well over a wider range of $r_{+}/k$ as $k$ decreases. Moreover, the Stokes solution agrees well with the BGK solutions up to $r_{+}/k = 30$ or $40$ .

Taken together, figures 15(a) and 15(b) imply that the BGK solution increases inversely with the distance from the upper corner from the far field down to a few tens of mean free paths (e.g. $r_{+} = c k$ with $c = 30$ or $40$ ). Considering the integral of $k/r_{+}$ over the interval deduced from above, $c k \lt r_{+} \lt x_{0}$ with $x_{0}$ a constant about $0.2$ , we have $\int _{ck}^{x_{0}}(k/r_{+})dr_{+} = - k \ln (ck) + k \ln (x_{0})$ , leading to the $k \ln k$ dependence.

Figure 15(b) provides a compact representation of the corner singularity in the Stokes solution and its modification within kinetic theory for a stress component. When referring to the singularity, attention is often paid to pressure rather than stress (Taylor Reference Taylor1962; Moffatt Reference Moffatt1964), presumably because pressure appears explicitly as a primary variable in the Stokes equations, although it is not directly linked to the force acting on a surface. With this in mind, to examine the case of pressure, the corresponding plot is shown in figure 15(c). Although the convergence of the curves with respect to $k$ is somewhat slower than in the case of shear stress, a similar inverse dependence on the distance from the corner followed by a suppression of pressure growth near the corner is confirmed. While the discussion here focuses on the stress and pressure on the lid, analogous behaviour can be confirmed also along the side walls (the corresponding figure is omitted for brevity).

4.6. Consideration on the validity of macroscopic fluid-dynamic descriptions

It is useful to reflect on the validity of macroscopic hydrodynamic descriptions in comparison with the present kinetic theory results. First, let us consider the case of the Stokes equations under the no-slip condition, also considered in the present paper, for sufficiently small $k$ . The Stokes solution agrees well with the kinetic solution, except for the macroscopic quantities within a distance of the order of the mean free path from the upper corners and the total forces over the surfaces including the upper corners (see, e.g. figures 6(a) and 6(b), 12 and 15). Furthermore, the corner vortices near the lower corners, predicted by the Stokes equations (Moffatt Reference Moffatt1964), are also observed in the kinetic solution as $k$ decreases. Thus, the present kinetic theory results demonstrate that, while the singularity at the upper corners in the Stokes solution is locally strong, it does not significantly propagate into the bulk region. This also supports the conventional belief in numerical simulations of the Stokes or Navier–Stokes equations that reliable results in the bulk can be obtained up to the grid resolution used, provided that the numerical scheme is designed not to excessively pick up the corner singularities. Such treatments have been achieved, for example, by using piecewise approximations or, in the case of spectral methods using global basis functions, by analytically subtracting the singular part of the solution near the corners to prevent pollution of the numerical results (see, e.g. Ghia et al. Reference Ghia, Ghia and Shin1982; Botella & Peyret Reference Botella and Peyret1998).

Next, let us consider the case of the Stokes equations under shear-slip boundary conditions and, more generally, continuum descriptions incorporating slip effects, in the slip-flow regime (e.g. $k=0.1$ ), which are usually employed to describe gas behaviour in this regime. At least two points would require careful consideration here, putting aside numerical computations using them. The first concerns the singular nature of the slip condition. For the present lid-driven cavity flow problem, with the present notation, the condition is given for the streamfunction $\psi$ as $\psi = 0,\ -n_{i}\partial _{x_{i}}\psi = u_{w}\chi _{\textit{U}} -k b_{\textit{s}}(n_{i}\partial _{x_{i}})(n_{j}\partial _{x_{j}}\psi ) \ [(x_{1}, x_{2}) \in \partial \varOmega ]$ . Here, $b_{\textit{s}}$ is a shear-slip coefficient, whose value depends on the kinetic equation and the boundary condition one adopts, and $(n_{i}\partial _{x_{i}})(n_{j}\partial _{x_{j}}\psi )$ corresponds to the normal gradient of the tangential flow velocity on the boundaries. Since $n_{i}\partial _{x_{i}} = r_{+}^{-1}\partial _{\varphi _{+}}$ on the lid and $n_{i}\partial _{x_{i}} = -r_{+}^{-1}\partial _{\varphi _{+}}$ on the right side wall, the shear-slip term is expressed as $-k b_{\textit{s}} r_{+}^{-2}\partial _{\varphi _{+}}^{2}\psi$ and, thus, includes a factor $r_{+}^{-2}$ , which would alter the gas behaviour significantly compared with the no-slip case, at least near the upper corners ( $0 \lt r_{+} \ll 1$ ). Whether such a slip term improves the solution relative to the no-slip case or conversely leads to inaccurate results due to the apparent singular factor $r_{+}^{-2}$ remains unclear.

We proceed under the assumption that the slip term improves the solution relative to the no-slip case, which exhibits a corner singularity. The second point concerns the $O(k)$ corrections in the bulk against the no-slip case, which are typically the primary interest in slip-flow analyses. Recall that, in flows over smooth surfaces, variations within the Knudsen layer influence the entire bulk at the same order through the slip boundary conditions in general, whose magnitudes are usually both $O(k)$ . In the present cavity flow, there is a region near the upper corners where the macroscopic quantities exhibit $O(1)$ variations in the scale of the mean free path. The fact that these $O(1)$ variations near the corners do not alter the behaviour in the bulk at the same magnitude $O(1)$ , as evidenced by the good agreement with the Stokes solution in the bulk that does not take into account such kinetic behaviour near the corners, is likely due to the geometric localisation of the corner, making its overall effect on the bulk region relatively small. In a two-dimensional plane, boundaries are lines, whereas corners are points. Nevertheless, at higher orders of $k$ , including $O(k)$ , which is important in the slip-flow regime, how the effects of the corner regions appear in the bulk remains unclear. In particular, neither the Stokes equation nor the commonly used shear-slip boundary condition written above explicitly account for the presence of corners, and therefore, their validity in capturing such effects is not evident a priori. Clarifying this point would require highly accurate numerical solutions of the corresponding macroscopic equations to directly assess the $O(k)$ contributions. Furthermore, it remains unclear whether the singular part of the solution for such equations can be analytically obtained or not, which is practically necessary for accurate numerical computation, as their asymptotic behaviour in approaching the corners may fundamentally differ from that of the Stokes equations under the no-slip condition. Thus, such an extensive comparative study lies beyond the scope of the present work.