2.1. The R26 constitutive relations

The system of 26 moment equations (2.3)–(2.7) and (2.11)–(2.13) contains the additional moments $\varPhi _{ijkl}$, $\varPsi _{ijk}$ and $\varOmega _{i}$, i.e.

(2.15)$$\begin{gather} \varPhi_{ijkl} = \int_{\mathbb{R}^{3}} C_{\langle i}C_{j}C_{k}C_{l\rangle }f\, \mathrm{d}\boldsymbol{c} , \end{gather}$$

(2.16)$$\begin{gather}\varPsi_{ijk} = \int_{\mathbb{R}^{3}}C^{2}C_{\langle i}C_{j}C_{k\rangle }f\, \mathrm{d}\boldsymbol{c} -9 m_{ijk}, \end{gather}$$

(2.17)$$\begin{gather}\varOmega_{i} = \int_{\mathbb{R}^{3}} C^{4}C_{i}\,f\, \mathrm{d}\boldsymbol{c} - 28 q_{i} , \end{gather}$$

which need to be specified in order to close the system of the 26 moment equations. The R26 constitutive relations, which provides a closure for the 26 moment equations, follow from a CE-like expansion on the evolution equations for the moments $\varPhi _{ijkl}$, $\varPsi _{ijk}$ and $\varOmega _{i}$ in a 45-moment theory (see Gu & Emerson Reference Gu and Emerson2009). The R26 constitutive relations read as (Gu & Emerson Reference Gu and Emerson2009)

(2.18)\begin{equation} \left. \begin{aligned} \varPhi_{ijkl} & ={-}4\frac{{Kn}}{{Pr}_{\varPhi}}\frac{\partial m_{\langle ijk}}{\partial x_{l\rangle}}, \\ \varPsi_{ijk} & ={-}\frac{27}{7}\frac{{Kn}}{{Pr}_{\varPsi}}\frac{\partial R_{\langle ij}}{\partial x_{k\rangle}}, \\ \varOmega_{i} & ={-}\frac{7}{3}\frac{{Kn}}{{Pr}_{\varOmega}} \left(\frac{\partial \varDelta}{\partial x_i} + \frac{12}{7}\frac{\partial R_{ij}}{\partial x_j}\right), \end{aligned} \right\} \end{equation}

where the values of the parameters are ${Pr}_{\varPhi }=2.097$, ${Pr}_{\varPsi }=1.698$ and ${Pr}_{\varOmega }=1$ for MM (Gu & Emerson Reference Gu and Emerson2009). The system of 26 moment equations (2.3)–(2.7) and (2.11)–(2.13) along with the constitutive relations (2.18) form the system of the LR26 equations.

4.1. Boundary conditions for the scalar fluxes

The scalar fluxes in (4.5) and their corresponding forces are collected in $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ as

(4.7a,b) \begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{J}} = \begin{bmatrix} \mathcal{V}_n \\ q_{n} \\ m_{nnn} \\ \varPsi_{nnn} \\ \varOmega_{n} \end{bmatrix} \quad\text{and}\quad \boldsymbol{\mathcal{F}} = \begin{bmatrix} \mathcal{P}+\sigma_{nn} \\ \mathcal{T}+\frac{2}{5}\sigma_{nn}+\frac{1}{5}R_{nn} +\frac{1}{15}\varDelta \\ \frac{3}{4}\sigma_{nn}+\frac{3}{28}R_{nn} +\frac{5}{12}\varPhi_{nnnn} \\ \frac{3}{56} R_{nn} \\ \frac{1}{120} \varDelta + \frac{1}{70} R_{nn} \end{bmatrix}. \end{equation}

Substituting $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ from (4.7a,b) in (4.6), we obtain the PBC for the scalar fluxes. Nevertheless, the matrix $\boldsymbol {\mathcal {L}}$ is still unknown. The elements of the matrix $\boldsymbol {\mathcal {L}}$ can be obtained by comparing the coefficients of $\mathcal {P}$, $\mathcal {T}$, $\sigma _{nn}$, $R_{nn}$ and $\varDelta$ in the corresponding terms of the PBC and MBC for the scalar fluxes. The MBC for the scalar fluxes (obtained through the Maxwell accommodation model) read as (Rana et al. Reference Rana, Lockerby and Sprittles2018b)

(4.8)\begin{equation} \left. \begin{aligned} \mathcal{V}_n & ={-} \varsigma \left(\mathcal{P} - \frac{1}{2} \mathcal{T} + \frac{1}{2} \sigma_{nn} - \frac{1}{28} R_{nn} - \frac{1}{120} \varDelta +b_1 \varPhi_{nnnn}\right), \\ q_n & ={-}\frac{\mathcal{V}_n}{2} - \varkappa \left(2 \mathcal{T} + \frac{1}{2} \sigma_{nn} + \frac{5}{28} R_{nn} + \frac{1}{15} \varDelta +b_2 \varPhi_{nnnn} \right), \\ m_{nnn} & ={-}\frac{2}{5} \mathcal{V}_n + \varkappa \left(\frac{2}{5} \mathcal{T} - \frac{7}{5} \sigma_{nn} - \frac{1}{14} R_{nn} + \frac{1}{75} \varDelta +b_3 \varPhi_{nnnn}\right), \\ \varPsi_{nnn} & = \frac{6}{5} \mathcal{V}_n + \varkappa \left(\frac{6}{5} \mathcal{T} + \frac{9}{5} \sigma_{nn} - \frac{93}{70} R_{nn} + \frac{1}{5} \varDelta +b_4 \varPhi_{nnnn}\right), \\ \varOmega_n & = 3 \mathcal{V}_n + \varkappa \left(8 \mathcal{T} + 2 \sigma_{nn} - R_{nn} - \frac{4}{3} \varDelta +b_5 \varPhi_{nnnn}\right), \end{aligned} \right\} \end{equation}

where

(4.9a–e)\begin{gather} b_1 ={-}\tfrac{1}{24}, \quad b_2 ={-} \tfrac{1}{12}, \quad b_3 ={-} \tfrac{13}{30}, \quad b_4 ={-} \tfrac{1}{2}, \quad b_5 = 0, \end{gather}

(4.10a,b)\begin{gather} \varsigma = \frac{\vartheta}{2-\vartheta} \sqrt{\frac{2}{\rm \pi}} \quad\text{and} \quad \varkappa = \frac{\vartheta + \chi (1-\vartheta)}{2-\vartheta - \chi (1-\vartheta)} \sqrt{\frac{2}{\rm \pi}}, \end{gather}

with $\vartheta$ being the evaporation/condensation coefficient and $\chi$ being the accommodation coefficient. Clearly, for $\vartheta =0$, $\mathcal {V}_n=0$ from (4.8); hence, we obtain the boundary conditions for a non-evaporating boundary (Gu & Emerson Reference Gu and Emerson2009). It is worthwhile noting that the MBC (4.8) with (4.9a–e) cannot be expressed in the force-flux formalism (4.6), and we shall show later that this leads to violation of the Onsager reciprocity relations and/or the second law at the interface.

Comparison of the coefficients of $\mathcal {P}$, $\mathcal {T}$, $\sigma _{nn}$, $R_{nn}$ and $\varDelta$ in the PBC for the scalar fluxes (obtained on substituting $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ from (4.7a,b) in (4.6)) and MBC (4.8) yields

(4.11)\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{L}} = \varsigma \begin{bmatrix} 1 & -\frac{1}{2} & -\frac{2}{5} & \frac{6}{5} & 3 \\ -\frac{1}{2} & \frac{1}{4} & \frac{1}{5} & -\frac{3}{5} & -\frac{3}{2} \\ -\frac{2}{5} & \frac{1}{5} & \frac{4}{25} & -\frac{12}{25} & -\frac{6}{5} \\ \frac{6}{5} & -\frac{3}{5} & -\frac{12}{25} & \frac{36}{25} & \frac{18}{5} \\ 3 & -\frac{3}{2} & -\frac{6}{5} & \frac{18}{5} & 9 \end{bmatrix} +\varkappa \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & -\frac{2}{5} & -\frac{6}{5} & -8 \\ 0 & -\frac{2}{5} & \frac{52}{25} & -\frac{44}{25} & \frac{8}{5} \\ 0 & -\frac{6}{5} & -\frac{44}{25} & \frac{916}{25} & -\frac{72}{5} \\ 0 & -8 & \frac{8}{5} & -\frac{72}{5} & 224 \end{bmatrix}. \end{equation}

Both the matrices on the right-hand side of (4.11) are symmetric and positive semidefinite. Therefore, the matrix $\boldsymbol {\mathcal {L}}$ in (4.11) is also symmetric and positive semidefinite as $\varsigma$ and $\varkappa$ are non-negative. Thus, the PBC for the scalar fluxes are obtained from (4.6) with $\boldsymbol {\mathcal {J}}$, $\boldsymbol {\mathcal {F}}$ and $\boldsymbol {\mathcal {L}}$ being given by (4.7a,b) and (4.11). On simplification, the PBC for the scalar fluxes are given by (4.8) with

(4.12a–e)\begin{equation} b_1 ={-}\tfrac{1}{6}, \quad b_2 ={-} \tfrac{1}{6}, \quad b_3 ={-} \tfrac{13}{15}, \quad b_4 = \tfrac{11}{15}, \quad b_5 ={-} \tfrac{2}{3}. \end{equation}

4.2. Boundary conditions for the vector fluxes

The vector fluxes in (4.5) and their corresponding forces are collected in $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ as

(4.13a,b)\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{J}} = \begin{bmatrix} \bar{\sigma}_{ni} \\ \bar{R}_{ni} \\ \bar{\varPhi}_{nnni} \end{bmatrix} \quad\text{and}\quad \boldsymbol{\mathcal{F}} = \begin{bmatrix} \bar{\mathscr{V}}_{i}+\frac{2}{5}\bar{q}_{i} +\bar{m}_{nni} \\ \frac{1}{5}\bar{q}_{i} + \frac{1}{7} \bar{m}_{nni} +\frac{1}{14} \bar{\varPsi}_{nni} +\frac{1}{70} \bar{\varOmega}_{i} \\ \frac{5}{8}\bar{m}_{nni} \end{bmatrix}. \end{equation}

Substituting $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ from (4.13a,b) in (4.6), we obtain the PBC for the vector fluxes. Nevertheless, the matrix $\boldsymbol {\mathcal {L}}$ is still unknown. The elements of the matrix $\boldsymbol {\mathcal {L}}$ can be obtained by comparing the coefficients of $\bar {\mathscr {V}}_{i}$, $\bar {q}_{i}$ and $\bar {m}_{nni}$ in the corresponding terms of the PBC and MBC. The MBC for the vector fluxes (obtained through the Maxwell accommodation model) read as (Gu & Emerson Reference Gu and Emerson2009; Struchtrup et al. Reference Struchtrup, Beckmann, Rana and Frezzotti2017)

(4.14)\begin{equation} \left. \begin{aligned} \bar{\sigma}_{ni} & ={-} \varkappa \left( \bar{\mathscr{V}}_{i} +\tfrac{1}{5}\bar{q}_{i} +\tfrac{1}{2}\bar{m}_{nni} + b_6 \bar{\varPsi}_{nni} + b_7 \bar{\varOmega}_{i}\right), \\ \bar{R}_{ni} & ={-} \varkappa \left( -\bar{\mathscr{V}}_{i} + \tfrac{11}{5} \bar{q}_{i} + \tfrac{1}{2} \bar{m}_{nni} + b_8 \bar{\varPsi}_{nni} + b_9 \bar{\varOmega}_{i}\right), \\ \bar{\varPhi}_{nnni} & ={-} \varkappa \left( -\tfrac{4}{7}\bar{\mathscr{V}}_{i} -\tfrac{12}{35} \bar{q}_{i}+\tfrac{9}{7} \bar{m}_{nni} + b_{10} \bar{\varPsi}_{nni} + b_{11} \bar{\varOmega}_{i}\right), \end{aligned} \right\} \end{equation}

where

(4.15a–f)\begin{align} b_6 ={-}\tfrac{1}{36}, \quad b_7 ={-}\tfrac{1}{280}, \quad b_8 = \tfrac{5}{12}, \quad b_9 = \tfrac{5}{56}, \quad b_{10} = \tfrac{1}{14}, \quad b_{11} ={-} \tfrac{3}{490}. \end{align}

Comparison of the coefficients of $\bar {\mathscr {V}}_{i}, \bar {q}_{i}, \bar {m}_{nni}$ in the PBC for the vector fluxes (obtained on substituting $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ from (4.13a,b) in (4.6)) and MBC (4.14) yields

(4.16)\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{L}} = \varkappa \begin{bmatrix} 1 & -1 & -\frac{4}{7} \\ -1 & 13 & -\frac{4}{7} \\ -\frac{4}{7} & -\frac{4}{7} & \frac{152}{49} \end{bmatrix}. \end{equation}

The matrix $\boldsymbol {\mathcal {L}}$ in (4.16) is symmetric and positive semidefinite as $\varkappa$ is non-negative. Thus, the PBC for the vector fluxes are obtained from (4.6) with $\boldsymbol {\mathcal {J}}$, $\boldsymbol {\mathcal {F}}$ and $\boldsymbol {\mathcal {L}}$ being given by (4.13a,b) and (4.16). On simplification, the PBC for the vector fluxes are given by (4.14) with

(4.17a–f)\begin{equation} b_6 ={-}\tfrac{1}{14}, \quad b_7 ={-}\tfrac{1}{70}, \quad b_8 = \tfrac{13}{14}, \quad b_9 = \tfrac{13}{70}, \quad b_{10} ={-} \tfrac{2}{49}, \quad b_{11} ={-} \tfrac{2}{245}. \end{equation}

4.3. Boundary conditions for the tensor fluxes

4.3.1. Boundary conditions for the rank-2 tensor fluxes

The rank-2 tensor fluxes in (4.5) and their corresponding forces are collected in $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ as

(4.18a,b)\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{J}} = \begin{bmatrix} \tilde{m}_{n ij} \\ \tilde{\varPsi}_{n ij} \end{bmatrix} \quad\text{and}\quad \boldsymbol{\mathcal{F}} = \begin{bmatrix} \frac{1}{2} \tilde{\sigma}_{ij} + \frac{1}{14}\tilde{R}_{ij} + \frac{1}{2} \tilde{\varPhi}_{nn ij} \\ \frac{1}{28} \tilde{R}_{ij} \end{bmatrix}. \end{equation}

Substituting $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ from (4.18a,b) in (4.6), we obtain the PBC for the rank-2 tensor fluxes. Nevertheless, the matrix $\boldsymbol {\mathcal {L}}$ is still unknown. The elements of the matrix $\boldsymbol {\mathcal {L}}$ can be obtained by comparing the coefficients of $\tilde {\sigma }_{ij}$ and $\tilde {R}_{ij}$ in the corresponding terms of the PBC and MBC for the rank-2 tensor fluxes.

The MBC for the rank-2 tensor fluxes (obtained through the Maxwell accommodation model) read as

(4.19)\begin{equation} \left. \begin{aligned} \tilde{m}_{nij} & ={-} \varkappa \left( \tilde{\sigma}_{ij} + \tfrac{1}{14}\tilde{R}_{ij} + b_{12} \tilde{\varPhi}_{nn ij} \right), \\ \tilde{\varPsi}_{nij} & = \varkappa \left( \tilde{\sigma}_{ij} - \tfrac{15}{14}\tilde{R}_{ij} + b_{13} \tilde{\varPhi}_{nn ij} \right), \end{aligned} \right\} \end{equation}

where

(4.20a,b)\begin{equation} b_{12} = \tfrac{1}{2} \quad\text{and}\quad b_{13} ={-} \tfrac{1}{2}. \end{equation}

Comparison of the coefficients of $\tilde {\sigma }_{ij}$ and $\tilde {R}_{ij}$ in the PBC for the rank-2 tensor fluxes (obtained on substituting $\boldsymbol {\mathcal {J}}$ and $\boldsymbol {\mathcal {F}}$ from (4.18a,b) in (4.6)) and MBC (4.19) yields

(4.21)\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{L}} = \varkappa \begin{bmatrix} 2 & -2 \\ -2 & 34 \end{bmatrix}. \end{equation}

The matrix $\boldsymbol {\mathcal {L}}$ in (4.16) is symmetric and positive semidefinite as $\varkappa$ is non-negative. Thus, the PBC for the rank-2 tensor fluxes are obtained from (4.6) with $\boldsymbol {\mathcal {J}}$, $\boldsymbol {\mathcal {F}}$ and $\boldsymbol {\mathcal {L}}$ being given by (4.18a,b) and (4.21). On simplification, the PBC for the rank-2 tensor fluxes are given by (4.19) with

(4.22a,b)\begin{equation} b_{12} = 1 \quad\text{and}\quad b_{13} = 1. \end{equation}

4.3.2. Boundary conditions for the rank-3 tensor fluxes

The rank-3 tensor flux in (4.5) and its corresponding forces are

(4.23a,b)\begin{equation} J_r = \check{\varPhi}_{n ijk} \quad\text{and}\quad \mathcal{F}_s = \tfrac{1}{6} \check{m}_{ijk}. \end{equation}

Substituting $\mathcal {J}_r$ and $\mathcal {F}_s$ from (4.23a,b) in (4.6), we obtain the PBC for the rank-3 tensor flux. Nevertheless, the coefficient $\mathcal {L}_{rs}$ in (4.6) is still unknown but can be obtained by comparing the coefficients of $\check {m}_{ijk}$ in the PBC and MBC for the rank-3 tensor flux.

The MBC for rank-3 tensor flux (obtained through the Maxwell accommodation model) read as

(4.24)\begin{equation} \check{\varPhi}_{nijk} ={-}\varkappa \left(\check{m}_{ijk} + b_{14} \check{\varPsi}_{ijk}\right), \end{equation}

where

(4.25)\begin{equation} b_{14} = \tfrac{1}{18}. \end{equation}

Comparison of the coefficients of $\check {m}_{ijk}$ in the PBC for rank-3 tensor fluxes (obtained on substituting $\mathcal {J}_r$ and $\mathcal {F}_s$ from (4.23a,b) in (4.6)) and MBC (4.24) yields

(4.26)\begin{equation} \mathcal{L}_{rs} = 6 \varkappa. \end{equation}

Thus, the PBC for the rank-3 tensor flux is obtained from (4.6) with $\mathcal {J}_r$, $\mathcal {F}_s$ and $\mathcal {L}_{rs}$ being given by (4.23a,b) and (4.26). On simplification, the PBC for the rank-3 tensor flux is given by (4.24) with

(4.27)\begin{equation} b_{14} = 0. \end{equation}

5.1. Poiseuille flow and thermal transpiration flow

Let us first consider the steady state [$\partial _t(\cdot )=0$] of a gas confined between two stationary ($v^{w}_{i}=0$), isothermal ($T^{w}=1$) and fully diffusive ($\chi =1$) walls of a channel located at $y=\mp 1/2$; see figure 1. The flow is assumed to be fully developed and driven either by a uniform external force $F$ (Poiseuille flow) in the positive $x$-direction or by a constant temperature gradient parallel to the walls $\tau$ (thermal transpiration flow). The normals for the bottom and top walls are pointing in $\pm y$-directions, respectively.

Figure 1.

Schematic of the Poiseuille flow ($F=1$, $\tau =0$) and thermal transpiration flow ($F=0$, $\tau =1$) in a channel.

The solution for the velocity and heat flux, for the problem under consideration, from the LR26 equations read as

(5.1a)$$\begin{gather} v_{x} = C_{1}+\underline{\frac{F}{2{Kn}}\left( \frac{1}{4}-y^{2}\right)} -\frac{2}{5} \sum_{i=1}^{4}\left( C_{i}^{K}-\frac{5}{2}\beta_{i}\right)\exp \left( \frac{y\gamma _{i}}{{Kn}}\right), \end{gather}$$

(5.1b)$$\begin{gather}q_{x} = \underline{-\frac{5{Kn}}{2\Pr }\tau} -\frac{{Kn}}{\Pr } F+\sum_{i=1}^{4}C_{i}^{K}\exp \left( \frac{y\gamma _{i}}{{Kn}}\right), \end{gather}$$

where $C_{1}$ and $C^{K}_{1}$, $C^{K}_{2}$, $C^{K}_{3}$ and $C^{K}_{4}$ are the constants of integration. Moreover, due to symmetry about the $x$-axis, $C^{K}_{2}=C^{K}_{1}$, $C^{K}_{4}=C^{K}_{3}$, $\beta _2 = \beta _1$ and $\beta _4 = \beta _3$. The unspecified three constants of integration, i.e. $C_1$, $C^{K}_{1}$ and $C^{K}_{3}$, are to be evaluated from the three (vector) boundary conditions at $y=-1/2$.

The underlined terms in (5.1) represent the classical solution of the NSF equations without velocity slip for the Poiseuille flow. The constant term $C_{1}$ in (5.1) accounts for the velocity slip and terms with exponential functions describe the Knudsen layer contributions to the velocity and heat flux. In classical hydrodynamics, the velocity profile for the Poiseuille flow is parabolic; the LR26 equations add velocity slip $C_{1}$ and Knudsen layer contributions. Similarly, the Fourier's contribution to the heat flux is given by the underlined term in (5.1b). The coefficients $\gamma _{i}$'s ($i\in \{1,2,3,4\}$) in (5.1) determine the thickness of the Knudsen layer; hence, the LR26 equations describe the Knudsen layer as a superposition of four exponential functions. The coefficients $\beta _1$, $\beta _3$ and $\gamma _{i}$'s ($i\in \{1,2,3,4\}$) for MM read as $\beta _1 = -0.24697$, $\beta _3 = 2.25494$, $\gamma _1 = -0.510199$, $\gamma _2 = 0.510199$, $\gamma _3 = -1.26497$ and $\gamma _4 = 1.26497$ (Gu et al. Reference Gu, Emerson and Tang2010; Gu & Emerson Reference Gu and Emerson2014).

The total mass flow and heat flow rates of the gas are, respectively, defined as

(5.2a,b)\begin{equation} M = \int_{{-}1/2}^{1/2}v_{x}\,\mathrm{d}y \quad\text{and}\quad Q = \int_{{-}1/2}^{1/2}q_{x}\,\mathrm{d}y. \end{equation}

We shall use the superscripts ‘$P$’ and ‘$T$’ for the Poiseuille ($F=1$, $\tau =0$) and transpiration ($F=0$, $\tau =1$) flows, respectively.

Figure 2(a) exhibits the mass flow rate, $M^{P}$, in the Poiseuille flow plotted over the Knudsen number. The results computed from the MBC applied to the LR26 equations with the MM collision model are compared with the LBE data from Wu et al. (Reference Wu, Zhang, Liu, Zhang and Reese2017). The Knudsen minimum is observed in the Poiseuille flow, in which the mass flow rate $M^{P}$ first decreases with the Knudsen number, attains a minimum value around a critical Knudsen number and then increases; see figure 2(a). The LR26 equations with PBC closely follow the mass flow rate profile from the LBE up to ${Kn}\lesssim 5$.

Figure 2.

Poiseuille flow and thermal transpiration flow: (a) the total mass flux $M^{P}$ in the Poiseuille flow and (b) violation of the Onsager reciprocity relation, $M^{T}=Q^{P}$, for the MBC in the case of the LR26 equations (main panel) and in the case of the LR13 equations (inset). The triangles in the left panel denote the results obtained from the LBE in Wu et al. (Reference Wu, Zhang, Liu, Zhang and Reese2017).

Due to the microscopic reversibility, the Onsager reciprocity relations should hold, which in this case state that $M^{T}=Q^{P}$ (Wu et al. Reference Wu, Zhang, Liu, Zhang and Reese2017). It is evident from figure 2(b) that the MBC do not comply with the Onsager reciprocity relations for the LR26 (and also for the LR13; shown in the inset) equations, while the PBC do.

Figure 3 illustrates the entropy generation rates at the interface computed with the LR26 equations (main panels) and the LR13 equations (insets) in the case of (a) the Poiseuille flow and (b) the transpiration flow. It is clear from figure 3(a) that the MBC for the LR26 equations (as well as for the LR13 equations) always violate the second law of thermodynamics (i.e. the entropy production $\varSigma$ becomes negative) for large Knudsen numbers. For the transpiration flow case (figure 3b), the violation of the second law of thermodynamics by the MBC for the LR26 equations occurs already at a relatively small Knudsen number (${Kn}<1$), which is within the transition regime. Surprisingly, the LR13 equations with MBC in the transpiration case do not contravene with the second law of thermodynamics whilst they do in the Poiseuille flow case around ${Kn}=1$; see the insets in figure 3. Nevertheless, the PBC always predict a positive entropy production in both the cases for both models in contrast to the conventional MBC.

Figure 3.

The entropy generation rates at the boundary for (a) the Poiseuille flow and (b) the thermal transpiration flow. The main panels and insets exhibit the results computed with the LR26 and LR13 equations, respectively.

5.2. Evaporation of a spherical liquid droplet

We consider a liquid droplet of (fixed) initial radius $\hat {R}_0$ with a given interface temperature $\hat {T}_l$ and the corresponding saturation pressure $\hat {p}_{{sat}}$, immersed in its own vapour; see figure 4. The far-field (dimensional) temperature $\hat {T}_{\infty }$, pressure $\hat {p}_{\infty }$ and the initial radius of the droplet $\hat {R}_0$ are used for the purpose of non-dimensionalisation.

Figure 4.

Schematic of an evaporating liquid droplet surrounded by its own vapour.

We consider two cases: (i) the flow is driven by a unit pressure difference while the temperature of the liquid is equal to the far-field temperature [i.e. $\varDelta p := (\hat {p}_{{sat}}-\hat {p}_{\infty })/\hat {p}_{\infty }=1$ and $\varDelta T := (\hat {T}_{l}-\hat {T}_{\infty })/\hat {T}_{\infty }=0$], and (ii) the flow is driven by a unit temperature difference while the saturation pressure is the same as the far-field pressure (i.e. $\varDelta p =0$ and $\varDelta T =1$). The problem has been investigated in detail by Rana et al. (Reference Rana, Lockerby and Sprittles2018b). We skip the details here for the sake of succinctness and focus mainly on looking at the Onsager reciprocity relations; for computational (and other) details, the reader is referred to Rana et al. (Reference Rana, Lockerby and Sprittles2018b). The Onsager reciprocity relations for this problem lead to $Q^{\varDelta p}=M^{\varDelta T}$ (Chernyak & Margilevskiy Reference Chernyak and Margilevskiy1989), where $Q^{\varDelta p}$ is the (dimensionless) heat flux in the pressure-driven process, $M^{\varDelta T}$ is the (dimensionless) mass flux in the temperature-driven process, and the superscripts ‘$\varDelta p$’ and ‘$\varDelta T$’ are used to denote the pressure- and temperature-driven processes, respectively.

The mass and heat fluxes in the pressure- and temperature-driven processes computed from the R13 and R26 theories using the PBC and MBC are illustrated in figure 5 for different values of the Knudsen number. The top and bottom rows of figure 5 exhibit the results for the pressure- and temperature-driven processes, respectively, while the left and right columns display the mass and heat fluxes, respectively. For comparison purpose, the results obtained by solving the LBE (Sone, Takata & Wakabayashi Reference Sone, Takata and Wakabayashi1994) are also included and they are depicted by symbols in the figure. In general, the results from both the R13 and R26 theories agree with those from the LBE for ${Kn} \lesssim 0.5$. However, for the mass flux in the pressure-driven case (figure 5a), the R26 equations on using the PBC (blue solid lines) give an excellent match with the LBE results (symbols) for ${Kn} \lesssim 0.5$ whereas the results from the other theories are comparatively higher than the LBE results even for small Knudsen numbers. Let us focus on panels (b,c) of figure 5, which illustrate the cross-effects, i.e. the heat flux in the pressure-driven case and mass flux in the temperature-driven case, respectively, and account for the verification of the Onsager reciprocity relations. The LBE predicts the same values the heat flux in the pressure-driven case and mass flux in the temperature-driven case, thanks to the Onsager reciprocity relations (Chernyak & Margilevskiy Reference Chernyak and Margilevskiy1989; Sone et al. Reference Sone, Takata and Wakabayashi1994). Evidently from panels (b,c) of figure 5, the heat flux in the pressure-driven case computed with the LR26/LR13 theory using the PBC (blue solid line for the LR26 equations and green solid line for the LR13 equations) seems to be the same as the mass flux in the temperature-driven case computed with the LR26/LR13 theory using the PBC (blue solid line for the LR26 equations and green solid line for the LR13 equations), whereas the heat flux in the pressure-driven case computed with the LR26/LR13 theory using the MBC (red dashed line for the LR26 equations and magenta dashed line for the LR13 equations) seems to be different from the mass flux in the temperature-driven case computed with the LR26/LR13 theory using the MBC (red dashed line for the LR26 equations and magenta dashed line for the LR13 equations). In order to have a clear picture of the Onsager reciprocity relations, the difference of $Q^{\varDelta p}$ and $M^{\varDelta T}$ is plotted in figure 6, which clearly shows that both the LR13 and LR26 equations on using the MBC for this problem violate the Onsager reciprocity relations while both the LR13 and LR26 equations on using the PBC respect the Onsager reciprocity relations for all Knudsen numbers.

Figure 5.

Profiles of the mass and heat fluxes plotted over the Knudsen number for the problem of evaporation of a spherical liquid droplet: (a) mass flux in the pressure-driven case, (b) heat flux in the pressure-driven case, (c) mass flux in the temperature-driven case and (d) heat flux in the temperature-driven case. The LBE data (symbols) are taken from Sone et al. (Reference Sone, Takata and Wakabayashi1994).

Figure 6.

Violation of the Onsager reciprocity relations for an evaporating droplet.

We have also computed the entropy generation rate at the interface for this problem and found (but not shown here for the sake of conciseness) that it remains positive for the R13 and R26 equations with the PBC in both pressure- and temperature-driven cases, also remains positive for the R13 equations with the MBC in the pressure-driven case but does not remain positive for the R13 equations with the MBC in the temperature-driven case for all Knudsen numbers. The R26 equations with the MBC yield a positive entropy generation rate for this problem.

6.1. Analytic solution of the LR26 equations

The analytic solutions of the LR26 equations (6.5)–(6.12) for the aforesaid problem can be obtained by converting these partial differential equations to ordinary differential equations via assuming an explicit dependency on the azimuthal angle $\theta$. To this end, the vectorial and tensorial components with an odd number of indices in $\theta$ are taken to be proportional to $\sin {\theta }$ while scalars and tensorial components with an even number of indices in $\theta$ are taken to be proportional to $\cos {\theta }$ (Torrilhon Reference Torrilhon2010). The components for vector fields, namely the velocity $\boldsymbol {v}(r,\theta )$ and the heat flux $\boldsymbol {q}(r,\theta )$, are taken to be

(6.13)\begin{equation} \left. \begin{aligned} v_r(r, \theta) & = u_{\infty} [1+a_{1}(r)] \cos{\theta}, \quad v_{\theta}(r, \theta) ={-}u_{\infty}[1+a_{2}(r)] \sin{\theta}, \\ q_r(r, \theta) & = \alpha_1(r) \cos{\theta}, \quad q_\theta(r, \theta) = \alpha_2(r) \sin{\theta}. \end{aligned} \right\} \end{equation}

The scalars, namely $T(r,\theta )$, $p(r,\theta )$ and $\varDelta (r,\theta )$, are taken to be

(6.14a–c)\begin{equation} T(r,\theta) = c(r) \cos{\theta}, \quad p(r,\theta) = d(r) \cos{\theta}, \quad \varDelta(r,\theta) = \delta(r)\cos{\theta}. \end{equation}

The components of the second- and third-rank tensors are taken to be

(6.15)\begin{equation} \left. \begin{aligned} \sigma_{rr}(r, \theta) & = b_1(r)\cos{\theta}, \quad \sigma_{r\theta}(r, \theta) = b_2(r)\sin{\theta}, \\ R_{rr}(r, \theta) & = \beta_1(r)\cos{\theta}, \quad R_{r\theta}(r, \theta) = \beta_2(r) \sin{\theta}, \\ m_{rrr}(r, \theta) & = \kappa_1(r)\cos{\theta}, \quad m_{rr\theta}(r, \theta) = \kappa_2(r) \sin{\theta}. \end{aligned} \right\} \end{equation}

In (6.13)–(6.15), $a_{1}(r)$, $a_{2}(r)$, $\alpha _{1}(r)$, $\alpha _{2}(r)$, $c(r)$, $d(r)$, $\delta (r)$, $b_{1}(r)$, $b_{2}(r)$, $\beta _{1}(r)$, $\beta _{2}(r)$, $\kappa _{1}(r)$ and $\kappa _{2}(r)$ are thirteen unknowns in total, which need to be determined from thirteen equations (6.5)–(6.10).

Substitution of constitutive relations (6.11)–(6.12) along with solution ansatz (6.13)–(6.15) in (6.5)–(6.10) results into a set of eight first-order ordinary differential equations (emanating from (6.5)–(6.7)) and five second-order ordinary differential equations (from (6.8)–(6.10)). The latter include second derivatives of $\delta (r)$, $\beta _1(r)$, $\beta _2(r)$, $\kappa _1(r)$ and $\kappa _2(r)$ due to regularised constitutive relations (6.11)–(6.12). With this, we obtain an analytic solution for the unknowns $a_{1}(r)$, $a_{2}(r)$, $\alpha _{1}(r)$, $\alpha _{2}(r)$, $c(r)$, $d(r)$, $\delta (r)$, $b_{1}(r)$, $b_{2}(r)$, $\beta _{1}(r)$, $\beta _{2}(r)$, $\kappa _{1}(r)$ and $\kappa _{2}(r)$, which has been relegated to Appendix E for better readability. Expressions (E1)–(E13) contain eight integration constants $C_{1,2,3}$ and $K_{1,2,3,4, 5}$ that are determined from boundary conditions (4.8)–(4.14).

Let us investigate a few qualitative features of these solutions. These solutions are the sums of the bulk contributions (which are polynomials in $1/r$) and the Knudsen layer contributions (terms containing the exponential functions $\varTheta _i(r)$); see Appendix E. The bulk solutions contain three constants of integration $C_{1,2,3}$ while Knudsen layers introduce five additional constants – $K_{1,2}$ for vector variables and $K_{3,4,5}$ for scalar variables. The bulk part can further be divided into the classical NSF contribution and second-order contribution (in the Knudsen number) that is neglected in the classical NSF theory. In the hydrodynamic limit (i.e. Stokes flow without velocity slip and temperature jump), $C_1=-3$ and $C_2=1/2$ while $C_3=0$ (Lamb Reference Lamb1932). The last relation ($C_3=0$) implies an isothermal flow condition. Hence, the non-zero bulk contributions to the temperature and heat flux stem from the second-order terms, i.e. the terms with $C_3$ in (E3a) and (E1c,d). The NSF equations with a velocity-slip boundary condition predict $C_1$ and $C_2$ as a function of the Knudsen number ${Kn}$ and the accommodation coefficient $\chi$ (Torrilhon Reference Torrilhon2010).

For the Knudsen layer contributions in (E1)–(E13), the LR13 theory (Torrilhon Reference Torrilhon2010) approximates Knudsen layers with two exponential functions whilst the LR26 theory gives superposition of five exponential functions to describe the Knudsen layers. On the contrary, the NSF theory does not predict Knudsen layers at all. In the following, we shall show that the LR26 theory gives a remarkably improved flow description with added details to the Knudsen layers.

6.2. Gas flow past a rigid sphere

First, let us consider the classical problem of a gas flow past a rigid sphere. The same problem was investigated by Torrilhon (Reference Torrilhon2010) using the LR13 and linear NSF equations with velocity slip.

Figure 8 illustrates the profiles of the radial and polar components of the velocity, pressure and temperature as functions of the radial distance $r$ starting from the surface of the sphere $r=1$. The exact solutions for the problem under consideration obtained with the LR26 equations in the present work are compared with those obtained with the LR13 and linear NSF solutions in Torrilhon (Reference Torrilhon2010) and with the numerical solutions of the LBE for three values of the Knudsen number ${\mathsf{Kn}} = 0.1, 0.4, 1$ given in Takata et al. (Reference Takata, Sone and Aoki1993). Note that the Knudsen number ${\mathsf{Kn}}$ used in Takata et al. (Reference Takata, Sone and Aoki1993) relates to the Knudsen number ${Kn}$ used in the present work via ${Kn}= 1.27\,{\mathsf{Kn}}/\sqrt {2}$. Therefore, the Knudsen number values ${\mathsf{Kn}} = 0.1, 0.4, 1$ in Takata et al. (Reference Takata, Sone and Aoki1993) correspond to the Knudsen number values ${Kn} \approx 0.09, 0.36, 0.9$, respectively, used in the present work. As expected, for small values of the Knudsen number [red curves and symbols (${Kn} = 0.09$)], all theories, including the NSF theory, give a good match with the LBE data for the velocities. However, for moderate values of the Knudsen number [blue curves and symbols (${Kn} = 0.36$)], the results from the NSF theory show deviations from those from kinetic theory (LBE results) while the LR13 theory still gives reasonably good results. The mismatch between the solutions from the NSF equations and LBE is attributed to the lack of Knudsen layer description by the former. For even larger values of the Knudsen number [green curves and symbols (${Kn} = 0.9$)], even the results from the LR13 theory start deviating from the LBE results. Notably, the LR13 theory in this case predicts $u_\theta / (u_\infty \sin {\theta }) >1$ near the surface of the sphere, which means that the slip velocity on the top of the sphere ($\theta ={\rm \pi} /2$) is larger than the far-field velocity – a result which is clearly unphysical and contradicted by more accurate theories (LBE and LR26). Hence, one may conclude that for large values of the Knudsen number, the results from LR13 theory cannot be trusted. On the other hand, the LR26 theory offers an excellent match upto ${Kn} \sim 1$. Interestingly, the LR26 theory also predicts $u_\theta / (u_\infty \sin {\theta }) > 1$ for ${Kn}\gtrsim 5$, which deduces that the LR26 theory can predict reasonably good results for the Knudsen numbers $0<{Kn}\lesssim 5$.

Figure 8.

Profiles of the (a) radial velocity, (b) polar velocity (c) pressure and (d) temperature of the gas past a rigid sphere as functions of the radial distance from the surface of the sphere for various Knudsen numbers. Results obtained from the LR26 theory (solid lines) are compared with those obtained from the LR13 (dashed lines) and NSF (dotted lines) theories (Torrilhon Reference Torrilhon2010) and with the numerical solution of the LBE (symbols) from Takata, Sone & Aoki (Reference Takata, Sone and Aoki1993).

Another remarkable feature in figure 8(d) is that the NSF theory predicts a non-negative temperature profile whereas all other theories, including the LBE, give a non-positive temperature profile. To demonstrate the consequences of this, the heat flux lines and temperature contours obtained from LR26, LR13 and NSF theories are illustrated in figure 9. It is evident from the figure that all the theories, including the NSF theory, do predict the thermal polarisation – an effect that, for the problem under consideration, shows different temperatures in the front and back sides of the sphere in the absence of nonlinear dissipations and external temperature differences (Aoki & Sone Reference Aoki and Sone1987; Takata et al. Reference Takata, Sone and Aoki1993; Torrilhon Reference Torrilhon2010). Nevertheless, the LR13 and LR26 theories show a hot region in the front side of the sphere and a cold region in the back side of the sphere while the temperature prediction by the NSF theory is the other way round. The temperature predictions by the LR13 and LR26 theories are in agreement with the findings of Takata et al. (Reference Takata, Sone and Aoki1993) based on the LBE and of Torrilhon (Reference Torrilhon2010) based on the LR13 equations, and also make sense because the temperature is likely to be more on the front side of the sphere due to compression of the gas at the front side of the sphere and is likely to be comparatively lesser on the back side of the sphere due to expansion of the gas at the back side of the sphere. Furthermore, the LR13 and LR26 theories predict that the heat flows from the cold region (back side of the sphere) to the hot region (front side of the sphere), which is a non-Fourier effect. The NSF theory cannot capture this effect and shows the heat flowing from the hot region to the cold region. Such a non-Fourier effect is well established in lid-driven cavity problems, (see, e.g. Rana, Torrilhon & Struchtrup Reference Rana, Torrilhon and Struchtrup2013, and references therein). By comparing the top and bottom rows of figure 9, it is also clear that the non-uniformity in temperature (i.e. the temperature difference between the temperatures at the front and back sides of the sphere) is more for a larger Knudsen number. This is also in agreement with the findings of Takata et al. (Reference Takata, Sone and Aoki1993).

Figure 9.

Heat flux lines drawn on top of the temperature contours for ${Kn} = 0.09$ (a–c) and ${Kn} = 0.36$ (d–f) from the LR26, LR13 and NSF theories.

Figure 10 illustrates the drag force normalised with the Stokes drag, $F_{{Stokes}} = 6{\rm \pi} {Kn}u_{\infty }$, plotted over the Knudsen number. The experimental data from Goldberg (Reference Goldberg1954) (fitted by Bailey et al. Reference Bailey, Barber, Emerson, Lockerby and Reese2005) is included for comparison. As expected, all theories agree with the experimental data for small values of the Knudsen number. As the Knudsen number increases, the drag decreases due to slip, which in not embedded in the Stokes formula and, hence, the Stokes formula gives vanishing drag force for all values of the Knudsen number. Notably, the drag force from the NSF theory with velocity-slip boundary conditions attains a non-zero positive value as ${Kn} \to \infty$ while the other theories (LR13 and LR26) predict a vanishing drag in this limit – an observation tantamount to experimental results. Furthermore, the LR26 theory provides an excellent prediction for the drag force for values of the Knudsen number beyond unity (upto ${Kn} \approx 10$). On the contrary, the LR13 theory underpredicts the drag force significantly for ${Kn} \gtrsim 0.2$.

Figure 10.

Normalised drag force from the NSF (dot-dashed red line), R13 (dashed green line) and R26 (solid blue line) theories as a function of the Knudsen number. The circles denote the experimental data of Goldberg (Reference Goldberg1954) (fitted by Bailey et al. Reference Bailey, Barber, Emerson, Lockerby and Reese2005). The results are normalised with the Stokes drag $F_{{Stokes}}=6{\rm \pi} {Kn} u_\infty$.

6.3. Gas flow past an evaporating spherical droplet

Finally, we consider the case of an evaporating droplet placed in a saturated vapour ($p_{{sat}} = p_\infty = 0$) flowing past the droplet with velocity $u_\infty$, i.e. phase change is driven by the flow field, not by the far field. The assumptions made in the previous section remain true here, i.e. the shape and size of the droplet remain fixed and the internal motion of the droplet is ignored. In addition, the evaporation coefficient $\vartheta$ is assumed to be unity.

Figures 11 again illustrates the profiles of the velocity components, pressure and temperature as functions of the radial direction $r$ starting from the surface of the sphere $r=1$. The exact solutions for the problem under consideration obtained with the LR26 equations (solid lines) and LR13 equations (dashed lines) are compared with those obtained with the numerical solutions of the LBE for three values of the Knudsen number ${Kn} = 0.1, 0.4, 1$ given in Sone et al. (Reference Sone, Takata and Wakabayashi1994). Notice that the Knudsen number ${Kn}$ is the same as that used in Takata et al. (Reference Takata, Sone and Aoki1993). Therefore, similarly to the above, these Knudsen number values correspond to ${Kn} = 0.09, 0.26, 0.9$, respectively.

Figure 11.

Profiles of the (a) radial velocity, (b) polar velocity (c) pressure and (d) temperature of the gas past an evaporating spherical droplet as functions of the radial distance from the surface of the droplet for various Knudsen numbers. Results obtained from the LR26 theory (solid lines) are compared with those obtained from the LR13 theory (dashed lines) and with the numerical solution of the LBE (symbols) from Sone et al. (Reference Sone, Takata and Wakabayashi1994).

As expected, when the Knudsen number is small [red lines and symbols (${Kn} = 0.09$)], both the LR13 and LR26 theories (dashed and solid lines) agree well with the results from the LBE (symbols). However, as the Knudsen number increases, the differences between the results from the two theories become noticeable – in particular, for the radial velocity and temperature. At higher Knudsen numbers (${Kn}\gtrsim 0.5$), the Knudsen layer description by the LR13 theory is inaccurate whilst the LR26 theory remains valid even for ${Kn}= 1$.

The normalised drag force in the case of a gas flow past an evaporating droplet is illustrated in figure 12. The solutions obtained from the NSF, LR13 and LR26 theories are compared against the numerical solutions of the LBE given in Sone et al. (Reference Sone, Takata and Wakabayashi1994). In contrast to the non-evaporating droplet case, all theories in this case predict vanishing drag as ${Kn}\to \infty$. For large values of the Knudsen number, the NSF theory overpredicts the drag while the LR13 theory underpredicts it. Again, the LR26 theory gives accurate predictions for the drag force even for large values of the Knudsen number (${Kn} \lesssim 10$).

Figure 12.

Normalised drag force from the LR13 (dashed green line) and LR26 (solid blue line) theories and from the numerical solution of the LBE (circles) from Sone et al. (Reference Sone, Takata and Wakabayashi1994) as a function of the Knudsen number. The results are normalised with the Stokes drag $F_{{Stokes}}=6{\rm \pi} {Kn} u_\infty$.