A theoretical derivation of slip boundary conditions based on the Cercignani–Lampis–Lord scattering model

Abstract

To characterize fluid flow in the slip regime, the use of Navier–Stokes–Fourier (NSF) equations with slip boundary conditions is prevalent. This trend underscores the necessity of developing reliable and accurate slip boundary conditions. According to kinetic theory, slip behaviours are intrinsically linked to the gas scattering processes at the surface. The widely used Maxwell scattering model, which employs a single accommodation coefficient to describe gas scattering processes, reveals its limitations when the difference between accommodation coefficients in the tangential and normal directions becomes significant. In this work, we provide a derivation of velocity slip and temperature jump boundary conditions based on the Cercignani–Lampis–Lord scattering model, which applies two independent accommodation coefficients to describe the gas scattering process. A Knudsen layer correction term is introduced to account for the impact of the surface on the velocity distribution function, which is associated with the scattering model. The governing equation of the correction term is established based on the linearized Boltzmann equation. Additionally, two moments are derived to capture the collision effect in the Knudsen layer: a conserving moment of collision invariants, and an approximate higher-order conserving moment. These moments are then employed to determine the coefficients in the correction term. We demonstrate that the derived slip coefficients align closely with numerical results obtained by solving the Boltzmann equation in the Knudsen layer. Besides, we apply the derived slip boundary conditions within the framework of the NSF equations, yielding numerical results that exhibit excellent consistency with those obtained through molecular-level simulations.

1. Introduction

The rarefied gas effect is a vital and foundational concern demanding attention in hypersonic flows (Ivanov & Gimelshein 1998; Candler 2019) and gas flows for microelectromechanical systems (MEMS) (Reese, Gallis & Lockerby 2003; Karniadakis, Beskok & Aluru 2006). In the case of hypersonic aircraft design, with its consideration of extremely high flow velocities, remarkably low atmospheric density, and sharp geometric shapes, the rarefied gas effect cannot be overlooked in the design process. Similarly, in the context of MEMS, as manufacturing processes and technology advance, the characteristic sizes of devices can reach the micron or even nanometre scale, which is comparable to the mean free path of gas molecules. This trend necessitates the consideration of the rarefied gas effect in the design process of MEMS.

The degree of rarefaction in a gas flow is typically characterized by the Knudsen number ($Kn$), which is defined as the ratio of the molecular mean free path to a specific characteristic length. Based on $Kn$, gas flow is conventionally classified into four regimes: the continuum regime ($Kn < 0.001$), the slip regime ($0.001 < Kn < 0.1$), the transition regime ($0.1 < Kn < 10$), and the free-molecule regime ($Kn > 10$).

In slip and early transition regimes, the application of the Navier–Stokes–Fourier (NSF) equations with no-slip boundary conditions becomes inappropriate (Struchtrup 2005; Sharipov 2015). On the other hand, the widely used direct simulation Monte Carlo (DSMC) method (Bird 1994), which is computationally efficient for simulating non-equilibrium flows in the transitional regime, becomes impractical for simulating complex three-dimensional gas flows in the near-continuum regime (Fei et al. 2020; Feng et al. 2023). As a result, there is a growing interest in developing a physically accurate and computationally efficient approach within the framework of NSF equations (Lockerby & Reese 2008; Lofthouse, Scalabrin & Boyd 2008; Gu & Emerson 2009). In the case of slightly rarefied flows, a widely adopted approach is applying suitable slip boundary conditions (Struchtrup 2005; Sharipov 2015).

Various slip boundary conditions have been proposed to address non-equilibrium effects that emerge near the wall within the slip regime, including first- and second-order boundary conditions, and others (Gökçen, MacCormack & Chapman 1987; Myong 2004; Le et al. 2012; Wang, Ou & Chen 2023); the corresponding results are well documented in the reviews (Cao et al. 2009; Zhang, Meng & Wei 2012; Akhlaghi, Roohi & Stefanov 2023). Classically, first-order boundary conditions encompass the Maxwell model (Maxwell 1879) for velocity slip, and the Smoluchowski model (Smoluchowski von Smolan 1898) for temperature jump. These models suggest that the slip velocity and temperature jump are proportional to the normal gradients of velocity and temperature at the surface, respectively. Towards increased accuracy and wider applicability of the NSF equations, a variety of second-order boundary conditions have also been developed (Hadjiconstantinou 2003; Lockerby et al. 2004; Radtke et al. 2012; Zeng et al. 2023), but these models include parameter values that remain a subject of debate for specific problems. When determining the slip coefficients in the boundary conditions mentioned above, accommodation coefficients (ACs) (Lord 1992; Arkilic, Breuer & Schmidt 2001) are used to quantify the extent to which molecules adapt to the wall. Typical slip boundary conditions assume that ACs are the same in all directions; however, experimental measurements (Yamamoto, Takeuchi & Hyakutake 2006; Liang, Li & Ye 2013) and molecular dynamics simulations (Spijker et al. 2010; Sipkens & Daun 2018) reveal substantial discrepancies between ACs in the tangential and normal directions. These discrepancies affect the accuracy of numerical simulations using the slip boundary conditions mentioned above. Consequently, it is necessary to develop slip boundary conditions that distinguish the tangential and normal ACs.

Numerous methods have been developed to establish slip boundary conditions, including the moment method (Grad 1949; Struchtrup & Weiss 2000; Torrilhon & Struchtrup 2008; To et al. 2015; Li & Yang 2023), half-flux method (Patterson 1956; Shen 2006; Zhang et al. 2021) and variational method (Loyalka 1968; Klinc & Kuščer 1972; McCormick 2005). The method of directly solving the linearized Boltzmann equation has also been developed to obtain precise slip coefficients (Sharipov 2003a, 2011; Siewert 2003; Basdanis, Valougeorgis & Sharipov 2023). Additionally, Sone, Ohwada & Aoki (1989) and Ohwada, Sone & Aoki (1989) used a finite difference method to numerically analyse the behaviour of the gas over a plane wall for hard-sphere molecules and diffuse reflection boundary conditions. Based on the numerical results, the Knudsen layer of a slightly rarefied gas flow past a body and the slip coefficients are derived. Aoki et al. (2017) systematically derived the slip boundary conditions for the compressible NSF equations on the basis of the Chapman–Enskog solution of the Boltzmann equation. The resulting formulae of the slip boundary conditions are summarized with explicit values of the slip coefficients for the Maxwell scattering model, contributing to the improved accuracy of previous slip boundary conditions through a rigorous theoretical derivation.

Historically, the first quantitative description of velocity slip was given by Maxwell (1879). Advancements upon the Maxwell scattering model showed that by considering the conservation of momentum and energy within the Knudsen layer, the expression for the slip coefficients could be determined. This groundbreaking insight gave rise to the now-famous Maxwell slip boundary condition. Loyalka (1971b) observed that Maxwell's results were not sufficiently accurate. Expanding on Maxwell's method, he incorporated the effect of collisions in the Knudsen layer. Specifically, he differentiated the distribution function of incident gas molecules near and far from the surface, and used two moments to determine them. This effort yielded values of slip coefficients that are closely aligned with numerical solutions obtained by solving the Boltzmann equation.

Notably, the slip boundary conditions derived using the aforementioned methods depend on the scattering model employed. Currently, most corresponding results are based on the Maxwell scattering model, which incorporates only one AC. Despite its widespread use, the Maxwell scattering model has notable limitations. For instance, it does not distinguish the ACs between normal and tangential directions, fails to accurately reproduce phenomena such as the lobular re-emission patterns observed when molecular beams hit the surface (Cercignani & Lampis 1971), struggles to describe backscattering (Basdanis, Tatsios & Valougeorgis 2022), and does not account for the temperature independence of ACs, among other shortcomings.

A more advanced scattering model was proposed by Cercignani & Lampis (1971) and later extended by Lord (1991, 1995), and it is thus denoted the Cercignani–Lampis–Lord (CLL) scattering model. This model includes two ACs, one for tangential momentum (TMAC) and the other for normal energy (NEAC). The CLL scattering model not only accurately reproduces feather-like structures around specular reflection lines in thermal beam scattering experiments (Cercignani & Lampis 1971), but also successfully predicts the thermal pressure difference index, a critical parameter for Knudsen pumps, to be less than 0.5 in the free-molecular limit of thermal transpiration (Sharipov 2003b). In addition, the CLL scattering model can effectively describe backscattering, which is a phenomenon that can occur on rough surfaces; these results agree exceptionally well with experimental measurements (Kalempa & Sharipov 2020). Recently, Sharipov & Volkov (2022) employed the CLL scattering model to evaluate the impact of ACs on the aerothermodynamics of a sphere based on ab initio interatomic potentials (Sharipov 2022), and they found that the TMAC and NEAC strongly affect all aerothermodynamic characteristics. The calculated drag and energy transfer coefficients for different ACs exhibit variations of up to $30\,\%$ and $200\,\%$, respectively.

Because the CLL scattering model performs better than the Maxwell scattering model, it is intriguing to derive slip boundary conditions based on the CLL scattering model. For instance, Struchtrup (2013) derived slip coefficients using the moment method. However, the resulting values deviate by $10\,\%$ from the numerical results obtained by solving the Boltzmann equation (Struchtrup 2005). Klinc & Kuščer (1972), as well as McCormick (2005), derived results using the variational method. Nevertheless, these derivations involve higher-order ACs, rendering them challenging for practical application. Recently, Basdanis et al. (2023) obtained the interpolation expression of slip coefficients for the CLL scattering model using symbolic regression. Zhang et al. (2021) obtained the CLL slip boundary condition via the half-flux method, yielding results consistent with those of Struchtrup (2013). Zhang et al. (2021) further refined their results by applying an empirical correction factor, akin to Loyalka's correction (Loyalka 1971b) for the Maxwell slip boundary condition. While the results of Basdanis et al. (2023) and Zhang et al. (2021) perform well, they rely on an empirical correction that lacks a rigorous derivation and clear physical foundation.

In this study, we employ a two-moment method to derive the slip boundary conditions based on the CLL scattering model. In our derivation, we consider the collision effect in the Knudsen layer, and incorporate the surface's impact on the distribution function by introducing a Knudsen layer correction term. We validate our newly derived slip coefficients by comparing our results with previous numerical results obtained by solving the Boltzmann equation. Furthermore, we apply the derived slip boundary conditions to computational fluid dynamics (CFD) calculations, comparing the results with those obtained using the DSMC method to verify their accuracy.

The remainder of the paper is organized as follows. Section 2 presents the theoretical derivation of our newly developed slip boundary conditions. In § 3, the accuracy of these boundary conditions is verified, including comparisons of slip coefficients and comparisons of computational results between the CFD and DSMC methods. Section 4 provides concluding remarks and some discussion.

2. Derivation of slip boundary conditions

2.1. Knudsen layer and scattering kernels in kinetic theory

The moment method proposed by Grad (1949) serves as a bridge between the mesoscale and macroscale, and it offers a promising tool to derive slip boundary conditions. To better elucidate the derivation process of velocity slip and temperature jump, certain fundamental aspects of kinetic theory necessary for the subsequent calculation are revisited.

For any gas flow near a solid surface, there exists a thin layer known as the Knudsen layer (Shan et al. 2022; Qian, Wu & Wang 2023). The Knudsen layer thickness is approximately that of a few mean free paths. In this region, collisions between gas molecules and the solid surface dominate, while collisions between gas molecules are rare enough that gas molecules remain in the thermodynamically non-equilibrium state. Therefore, the linear constitutive relationships assumed in the NSF equations no longer apply. To accurately characterize the dynamic behaviour of molecules in this region, kinetic methods must be employed, including direct calculations of the Boltzmann equation, and using the DSMC method.

A typical velocity profile within the Knudsen layer is depicted in figure 1. Here, the solid curve represents the average gas velocity obtained using kinetic methods. The actual velocity slip is denoted as $u( 0 )$. By extrapolating the velocity values from the outer edge of the Knudsen layer towards the surface, a ‘fictitious’ average velocity ${u_s}$ at the surface can be determined. If this fictitious velocity is used as a boundary condition for the NSF equations, then the velocity solution obtained outside the Knudsen layer will be identical to that obtained from the Boltzmann equation (Ivchenko, Loyalka & Tompson 2007). All discussions below are based on a reference frame relative to the surface. For a moving surface, all velocities can be transformed using the wall's velocity $- {u_w}$. The temperature jump follows the same trend.

Figure 1. Schematic of the velocity profile of the Knudsen layer near a surface. Here, the surface is assumed to be stationary, $\lambda$ is the molecular mean free path, ${u_w}$ is the velocity of the surface and equals zero in this situation, $u( 0 )$ is the actual velocity slip at the surface, and ${u_s}$ is the ‘fictitious’ slip velocity that would be needed to ensure that an NSF solution (dashed line) provides an accurate prediction (solid line) beyond the Knudsen layer.

According to kinetic theory, the state of gas molecules can be described using a distribution function $f( {{x_k},{c_k},t} )$, where ${x_k}$ denotes spatial coordinates, ${c_k}$ denotes molecular velocity, $t$ denotes time, and $f( {{x_k},{c_k},t} ) \,{\rm {d}} {\boldsymbol {c}} \,{\rm {d}}\kern0.06em {\boldsymbol {x}}$ represents the number of molecules with velocity ${\boldsymbol {c}}\sim {\boldsymbol {c}} + {\rm {d}}{\boldsymbol {c}}$ at position ${\boldsymbol {x}}\sim {\boldsymbol {x}} + {\rm {d}}\kern0.06em {\boldsymbol {x}}$ at time $t$. After the distribution function is determined, macroscopic properties can be obtained by taking moments, i.e.

(2.1)\begin{equation} \left. \begin{aligned} \rho & = \int {mf\,{\rm d}{\boldsymbol{c}}} ,\\ \rho {\bar{c} _k} & = \int {m{c_k}f\,{\rm d}{\boldsymbol{c}},} \\ \rho\,\frac{3}{2}\,RT & = \int {\frac{1}{2}\,mc_k'c_k'f\,{\rm d}{\boldsymbol{c}}} , \end{aligned} \right\} \end{equation}

where $\rho$ is mass density, ${\bar {c}_k}$ is macroscopic flow velocity, $T$ is temperature, $R = {k}/{m}$ is the gas constant, $m$ is the molecular mass, $k$ is the Boltzmann constant, and $c_k' = {c_k} - {\bar {c} _k}$ is the peculiar velocity of molecules. In a state of equilibrium, the distribution function of gas molecules is typically Maxwellian, i.e.

(2.2)\begin{equation} {f_M} = \frac{n}{{{{( {2{\rm \pi} RT} )}^{3/2}}}}\exp \left( { - \frac{{c_k'c_k'}}{{2RT}}} \right). \end{equation}

At the surface, interactions between gas molecules and the surface occur, resulting in differences between the distribution functions of incident and reflected molecules. We denote the velocities of incident molecules as ${{\boldsymbol {c}}_i}$, and the velocities of reflected molecules as ${{\boldsymbol {c}}_r}$. Here, only non-adsorbing surfaces are considered, which signifies that molecules immediately reflect upon colliding with the surface. The distribution function of gas molecules at the surface is denoted as

(2.3)\begin{equation} {f_w} = \left\{ \begin{array}{@{}ll} {f^ - }, & {c_{{n_i}}} \le 0,\\[2pt] {f^ + }, & {c_{{n_r}}} > 0, \end{array}\right. \end{equation}

where ${f^ - }$ is the distribution function of incident molecules (${c_{{n_i}}} = {{\boldsymbol {c}}_i} \boldsymbol {\cdot } {\boldsymbol {n}}$ represents the normal component of velocity), ${f^ + }$ is the distribution function of reflected molecules, and the distribution functions of incident and reflected molecules can be related through a scattering kernel (Cercignani & Cercignani 1988), i.e.

(2.4)\begin{equation} {c_{{n_r}}}\,{f^ + }( {{{\boldsymbol{c}}_r}} ) =- \int_{{{\boldsymbol{c}}_{{n_i}}} < 0} {{c_{{n_i}}}\,{f^ - }( {{{\boldsymbol{c}}_i}} )\,R( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} )\,{\rm d}{\boldsymbol{c}}'}. \end{equation}

The scattering kernel $R( {{{\boldsymbol {c}}_i} \to {{\boldsymbol {c}}_r}} )$ provides the probability of a gas molecule with velocity ${{\boldsymbol {c}}_i}\sim {{\boldsymbol {c}}_i} + \textrm {{d}}{{\boldsymbol {c}}_i}$ impacting the surface and being reflected with velocity ${{\boldsymbol {c}}_r}\sim {{\boldsymbol {c}}_r} + \textrm {{d}}{{\boldsymbol {c}}_r}$, and this probability is determined by the local temperature ${T_w}$ and velocity ${u_w}$ at the surface.

The specific form of the scattering kernel depends on the scattering model, and it must satisfy certain mathematical properties as provided below.

Nonnegativity:

(2.5)\begin{equation} R( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} ) \ge 0. \end{equation}

Normalization:

(2.6)\begin{equation} \int_{{c_n} > 0} {R( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} )} \,{\rm d}{\boldsymbol{c}} = 1. \end{equation}

Reciprocity:

(2.7)\begin{equation} {c_{{n_i}}}\exp \left( { - \frac{{{{\boldsymbol{c}}_i^2}}}{{2R{T_w}}}} \right)R( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} ) =- {c_{{n_r}}}\exp \left( { - \frac{{{{\boldsymbol{c}}_r^2}}}{{2R{T_w}}}} \right)R( { - {{\boldsymbol{c}}_r} \to - {{\boldsymbol{c}}_i}} ). \end{equation}

A schematic of various scattering models is shown in figure 2. For the Maxwell scattering model, the scattering kernel is typically given by

(2.8) \begin{equation} {R_M}( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} ) = \sigma\,\frac{{{c_n}}}{{2{\rm \pi} {{( {R{T_w}} )^2}}}}\exp \left( { - \frac{{{\boldsymbol{c}^2}}}{{2R{T_w}}}} \right) + ( {1 - \sigma } )\delta ( {{{\boldsymbol{c}}_i} - {{\boldsymbol{c}}_r} + 2{\boldsymbol{n}}{c_n}} ), \end{equation}

where the first part represents the diffusive reflection component, which signifies that molecules fully adapt to the surface's conditions, and the second part represents the specular reflection component, which signifies that molecules’ tangential velocities will remain unaltered while their normal velocities will be inverted. Here, $\sigma$ is the AC, and it takes values between 0 and 1, with the lower and upper limits corresponding to completely specular and diffuse reflection, respectively.

Figure 2. A schematic of various scattering models. Here, ${{\boldsymbol {c}}_i}$ denotes the velocity of the incident molecule, ${{\boldsymbol {c}}_r}$ denotes the velocity of the reflected molecule, and ${\boldsymbol {n}}$ is the normal vector at the surface. (a) Specular reflection model. (b) Diffuse reflection model. (c) Maxwell scattering model. (d) CLL scattering model.

Compared to the Maxwell scattering model, which includes only one coefficient, the CLL scattering model employs two ACs: the TMAC ${\sigma _t}$, and the NEAC ${\alpha _n}$. Its specific form is

(2.9)$$\begin{align} {R_{CLL}}( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} ) &= \frac{1}{{2{\rm \pi} }}\,\frac{1}{{{\alpha _n}{\sigma _t}( {2 - {\sigma _t}} )}}\,\frac{{{c_{{n_r}}}}}{{{{( {R{T_w}} )^2}}}}\,{{\rm I}_0}\left( {\frac{{\sqrt {1 - {\alpha _n}}\, {c_{{n_r}}}{c_{{n_i}}}}}{{{\alpha _n}R{T_w}}}} \right)\nonumber\\ &\times \exp \left[ { - \frac{{{c_{{n_r}}^2} + ( {1 - {\alpha _n}} ){c_{{n_i}}^2}}}{{2R{T_w}{\alpha _n}}} - \frac{{{{\left( {{{\boldsymbol{c}}_{{t_r}}} - (1 - {\sigma _t}){{\boldsymbol{c}}_{{t_i}}}} \right)^2}}}}{{2R{T_w}{\sigma _t}( {2 - {\sigma _t}} )}}} \right]. \end{align}$$

Here, ${{\boldsymbol {c}}_t}$ represents the tangential component of the velocity, which is a two-dimensional vector, and ${\textrm {I}_0}$ is the modified Bessel function of the first kind (order zero):

(2.10)\begin{equation} {{\rm I}_0}( x ) = \frac{1}{{2{\rm \pi} }}\int_0^{2{\rm \pi} } {\exp ( {x\cos \phi } )\,{\rm d}\phi }. \end{equation}

Note that ${\alpha _n}$ takes values between 0 and 1, while ${\sigma _t}$ takes values between 0 and 2, and ${\sigma _t} > 1$ represents the presence of molecular backscattering, which may occur on a rough surface. In extreme cases, when ${\sigma _t} = 2$, ${\alpha _n} = 0$, molecules scatter back completely; when ${\sigma _t} = 1$, ${\alpha _n} = 1$, molecules undergo complete diffuse reflection; and when ${\sigma _t} = 0$, ${\alpha _n} = 0$, molecules experience complete specular reflection. According to Sharipov & Moldover (2016), the typical values of the TMAC vary in the range $0.4 \le {\sigma _t} \le 1$, while the NEAC practically varies across the entire range, i.e. $0.01 \le {\alpha _n} \le 1$. Besides, in the CLL scattering model, the two tangential and one normal component of the scattering kernel are independent of each other, i.e.

(2.11)\begin{equation} {R_{CLL}}( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} ) = R( {{u_i},{u_r}} )\,R( {{w_i},{w_r}} )\,R( {{v_i},{v_r}} ). \end{equation}

Specifically, for the tangential velocity $u$, the scattering kernel is

(2.12)\begin{equation} R( {{u_i},{u_r}} ) = \frac{1}{{\sqrt {2{\rm \pi} RT_w{\sigma _t}( {2 - {\sigma _t}} )} }}\exp \left( { - \frac{{{{( {{u_r} - ( {1 - {\sigma _t}} ){u_i}} )^2}}}}{{2RT_w{\sigma _t}( {2 - {\sigma _t}} )}}} \right). \end{equation}

We can also introduce the tangential energy accommodation coefficient (TEAC) ${\alpha _t}$, which is related to the TMAC ${\sigma _t}$ by

(2.13)\begin{equation} {\alpha _t} = {\sigma _t}( {2 - {\sigma _t}} ). \end{equation}

Therefore, the scattering kernel for the tangential velocity $u$ can also be expressed as

(2.14)\begin{equation} R( {{u_i},{u_r}} ) = \frac{1}{{\sqrt {2{\rm \pi} RT_w{\alpha _t}} }}\exp \left( { - \frac{{{{\left( {{u_r} - \sqrt {1 - {\alpha _t}}\,{u_i}} \right)^2}}}}{{2RT_w{\alpha _t}}}} \right). \end{equation}

Note that in this case, ${\alpha _t}$ takes values between 0 and 1, and the case $1 < {\sigma _t} \le 2$ cannot be covered. The tangential velocity $w$ follows the same form. For normal velocity $v$, the scattering kernel is

(2.15)\begin{equation} R( {{v_i},{v_r}} ) = \frac{{{v_r}}}{{RT_w{\alpha _n}}}\,{{\rm I}_0}\left( {\frac{{\sqrt {1 - {\alpha _n}}\,{v_r}{v_i}}}{{RT_w{\alpha _n}}}} \right)\exp \left( { - \frac{{{v_r^2} + \left( {1 - {\alpha _n}} \right){v_i^2}}}{{2RT_w{\alpha _n}}}} \right). \end{equation}

In completing the corresponding calculations, partial integrals for the tangential and normal scattering kernels are provided in Appendix A.

2.2. The effect of the Knudsen layer on the distribution function

For a system that does not deviate significantly from the equilibrium state, an asymptotic solution of the Boltzmann equation can be obtained using the Chapman–Enskog expansion (Chapman & Cowling 1990). A concise overview of this procedure is provided in Appendix B. Truncated to first order, corresponding to the NSF equations, the distribution function can be written as

(2.16) \begin{equation} f = {f_0}\left[ {1 - \frac{{4\kappa {\beta ^2}}}{{5nk}}\left( {\frac{{\boldsymbol{c}'^2}}{{2RT}} - \frac{5}{2}} \right)c_i'\,\frac{{\partial \ln T}}{{\partial {x_i}}} - \frac{{4\mu {\beta ^4}}}{\rho }\left( {c_i'c_j' - \frac{1}{3}\,c_i'c_j'{\delta _{ij}}} \right)\frac{{\partial {{\bar{c} }_i}}}{{\partial {x_j}}}} \right], \end{equation}

where $\beta$ is the reciprocal of the most probable thermal speed, i.e. $\beta = {1}/{{\sqrt {2RT} }}$, $\mu$ is the molecular viscosity coefficient, and $\kappa$ is the thermal conductivity coefficient. According to Chapman & Cowling (1990), $\mu$ and $\kappa$ can be determined as

(2.17) \begin{equation} \left.\begin{aligned} \mu & = \frac{{2kT}}{{15}}\int_{ - \infty }^\infty B \left( {\frac{{\boldsymbol{c}'}}{{\sqrt {2RT} }}} \right)^4 {f_0}\,{\rm d}\boldsymbol{c}',\\ \kappa & = \frac{{2{k^2}T}}{{3m}}\int_{ - \infty }^\infty A \left( {\frac{{\boldsymbol{c}'}}{{\sqrt {2RT} }}} \right)^4 {f_0}\,{\rm d}\boldsymbol{c}', \end{aligned}\right\} \end{equation}

where $A$ and $B$ are defined in (B13)–(B16).

We now consider a semi-infinite expanse of a gas bounded by a flat plate located at $y = 0$, lying in the $x$–$z$ plane. At large distances from the plate, the distribution function can be described by the Chapman–Enskog expansion as shown above. However, in the vicinity of the plate, the presence of the Knudsen layer may invalidate the Chapman–Enskog solution. To address this, the influence of the Knudsen layer can be accounted for by introducing a Knudsen layer correction term $\varPhi ( {{\boldsymbol {c}},y} )$, which is a function of velocity and the spatial coordinate $y$. We further assume that the gas has a velocity $\bar {u} ( y )$ in the $x$-direction, and the gas maintains a constant gradient $\partial \bar {u}/\partial y$ and $\partial T/\partial y$ perpendicular to the plate. These gradients are sufficiently small. In this case, the distribution function can be further written as

(2.18) \begin{equation} f = {f_0}\left[ {1 - \frac{{4\kappa {\beta ^2}}}{{5nk}}\left( {\frac{{\boldsymbol{c}'^2}}{{2RT}} - \frac{5}{2}} \right)v'\,\frac{{\partial \ln T}}{{\partial y}} - \frac{{4\mu {\beta ^4}}}{\rho }\,u'v'\,\frac{{\partial \bar{u} }}{{\partial y}} + \varPhi } \right]. \end{equation}

The Knudsen layer appears very small from a macroscopic perspective but is sufficiently large when viewed from a microscopic perspective. For this reason, we denote the Knudsen layer correction term at the surface as $\varPhi ( {{\boldsymbol {c}},0} )$, and at the outer edge of the Knudsen layer as $\varPhi ( {{\boldsymbol {c}},\infty } )$. At the surface, the Knudsen layer correction terms of the incident distribution function ${\varPhi ^ - }( {{\boldsymbol {c}},0} )$ and reflected distribution function ${\varPhi ^ + }( {{\boldsymbol {c}},0} )$ differ due to molecules hitting the surface and then scattering out. The relationship between them can be determined by the scattering model, as mentioned previously in (2.4). Due to collisions between molecules, the Knudsen layer correction term becomes the same in a very short distance from the surface, so ${\varPhi ^ - }( {{\boldsymbol {c}},\infty } )$ and ${\varPhi ^ + }( {{\boldsymbol {c}},\infty } )$ are equal to each other. We denote the incident distribution function as ${f^ - }$ and the reflected distribution function as ${f^ + }$, so (2.18) can be written as

(2.19) \begin{equation} {f^ \pm } = {f_0}\left[ {1 - \frac{{4\kappa {\beta ^2}}}{{5nk}}\left( {\frac{{\boldsymbol{c}'^2}}{{2RT}} - \frac{5}{2}} \right)v'\,\frac{{\partial \ln T}}{{\partial y}} - \frac{{4\mu {\beta ^4}}}{\rho }\,u'v'\,\frac{{\partial \bar{u} }}{{\partial y}} + {\varPhi ^ \pm }} \right]. \end{equation}

Notably, ${f^ + }$ and ${f^ - }$ are different only at the surface; in the main region, they are identical.

Next, we will determine the governing equation of the Knudsen layer correction term $\varPhi$. Given that the system does not deviate significantly from equilibrium, the linearized Boltzmann equation can be used (Sharipov 2015), i.e.

(2.20)\begin{equation} \frac{{\partial h}}{{\partial t}} + {c_k}\,\frac{{\partial h}}{{\partial {x_k}}} + {F_k}\,\frac{{\partial h}}{{\partial {c_k}}} = L[ h ], \end{equation}

where the linearized collision operator is defined as

(2.21)\begin{equation} L[ h ] = \int_{ - \infty }^\infty {\int_0^{4{\rm \pi} } {\left( {{h^*} + {h^{1*}} - h - {h^1}} \right)f_0^1{c_r}\sigma \,{\rm d}\varOmega \,{\rm d}{\boldsymbol{c}_1}} }, \end{equation}

where $h$ is the deviation of the distribution function from equilibrium, which contains the solution of Chapman–Enskog expansion part $\phi$ as well as the added Knudsen layer correction terms, i.e. $\varPhi$. The relationship between them can be expressed as $f = {f_M}[ {1 + h} ] = {f_M}[ {1 + ( {\phi + \varPhi } )} ]$.

Following the principle of linear superposition, we subtract the Chapman–Enskog expansion part from (2.20), and the governing equation of the Knudsen layer correction term $\varPhi$ can be obtained. For the one-dimensional problem under consideration, the equation can be simplified as

(2.22)\begin{equation} v'\,\frac{{\partial \varPhi }}{{\partial y}} = L[ \varPhi ]. \end{equation}

Notably, in (2.22) and in the following contexts, we use a reference frame moving with the gas at the local average velocity. Therefore, we use peculiar velocities instead of instantaneous velocities.

2.3. The velocity slip problem

Now let us analyse the velocity slip problem. Considering the scenario depicted in figure 1, and neglecting changes in temperature, the distribution function can be approximated to first order as

(2.23)\begin{equation} {f^ \pm } = {f_0}\left[ {1 - \frac{{4\mu {\beta ^4}}}{\rho }\,u'v'\,\frac{{\partial \bar{u} }}{{\partial y}} + {\varPhi ^ \pm }} \right]. \end{equation}

According to the derivation of Maxwell, ${\varPhi ^ \pm }$ are assumed to take the form

(2.24)$$\begin{gather} {\varPhi ^ - }( {{\boldsymbol{c}},0} ) = {a_0}\,\frac{{u'}}{{\sqrt {2RT_0} }}, \end{gather}$$

(2.25)$$\begin{gather}{\varPhi ^ \pm }( {{\boldsymbol{c}},\infty } ) = {a_0}\,\frac{{u'}}{{\sqrt {2RT_0} }}, \end{gather}$$

where ${a_0}$ is a parameter to be determined, and ${T_0}$ is the equilibrium temperature, which is taken as the wall temperature ${T_w}$. Note that (2.24) and (2.25) rely on the strong assumption that the incident distribution function remains unchanged within the Knudsen layer, leading to significant uncertainties that cannot be ignored (Loyalka 1971b). To rectify this issue, we can assume that ${\varPhi ^ \pm }$ take the form

(2.26)$$\begin{gather} {\varPhi ^ - }\left( {{\boldsymbol{c}},0} \right) = {a_w}\frac{{u'}}{{\sqrt {2RT_0} }}, \end{gather}$$

(2.27)$$\begin{gather}{\varPhi ^ \pm }\left( {{\boldsymbol{c}},\infty } \right) = {a_0}\frac{{u'}}{{\sqrt {2RT_0} }}, \end{gather}$$

which uses a similar method to that of Loyalka (1971b). From a physical perspective, ${a_{0w}} \ne {a_0}$ indicates the inclusion of collision effects. If incident molecules in the Knudsen layer do not collide with other molecules, then the form of the distribution function is uniform throughout the Knudsen layer. As collisions between incident and reflected molecules are inevitable, discrepancies in distribution functions arise at the edge of the Knudsen layer and at the surface, leading to a clear distinction between ${a_{0w}}$ and ${a_0}$. Notably, compared to the macroscopic variables obtained from the original distribution function (2.16), the macroscopic variables such as $\rho,{\bar {c} _i},T$ determined by (2.23), (2.26) and (2.27) are unchanged according to (2.1), except for the tangential mean velocity $\bar {u}$, which corresponds to the velocity slip caused by the existence of the Knudsen layer.

Because ${\varPhi ^ - }( {{\boldsymbol {c}},0} )$ is known, the expression of ${\varPhi ^ + }( {{\boldsymbol {c}},0} )$ can be further determined using (2.4) with the CLL scattering model. After rearranging, ${\varPhi ^ + }( {{\boldsymbol {c}},0} )$ is written as

(2.28)\begin{align} {\varPhi ^ + }( {{\boldsymbol{c}},0} ) = [ {{A_0}( {1 - {\sigma _t}} ) + 1} ]\,\frac{\mu }{{pRT_0}}\,\frac{{\partial \bar{u} }}{{\partial y}}\,{u_r'}{v_r}' + {A_1}( {1 - {\sigma _t}} )\,\frac{\mu }{{pRT_0}}\,\frac{{\partial \bar{u} }}{{\partial y}}\,{u_r'} + ( {1 - {\sigma _t}} )\,\frac{{a{_{0w}}}}{{\sqrt {2RT_0} }}\,{u_r'}. \end{align}

Note that in the aforementioned calculation, integrals from Appendix A have been utilized, and ${A_0}$ and ${A_1}$ are two coefficients that are already determined.

The next step is to determine the two unknown parameters, i.e. ${a_{0w}}$ and ${a_0}$. This will use the moments of (2.22). Notably, following the principle of Chapman–Enskog expansion, the collision operator defined in (2.21) satisfies the relations

(2.29)\begin{gather} L[ \psi ] = 0, \end{gather}

(2.30)\begin{gather} L\left[ {B\,\frac{1}{{2RT}}\,u'v'} \right] = \frac{1}{{RT}}\,u'v', \end{gather}

(2.31)\begin{gather} L\left[ {A\,\frac{{v'}}{{\sqrt {2RT} }}} \right] = \frac{{v'}}{{\sqrt {2RT} }}\left( {\frac{{\boldsymbol{c}'^2}}{{2RT}} - \frac{5}{2}} \right), \end{gather}

where $A$, $B$ and $\psi$ are the same as in Appendix B.

It is useful to introduce a scalar product, which is defined as

(2.32)\begin{equation} ( {{\rho _1}( {x,{\boldsymbol{c}}} ),{\rho _2}( {x,{\boldsymbol{c}}} )} ) = \int {{\rho _1}( {x,{\boldsymbol{c}}} )\exp \left( { - \frac{{{\boldsymbol{c}}'^2}}{{2RT_0}}} \right){\rho _2}( {x,{\boldsymbol{c}}} )} \,{\rm d}{\boldsymbol{c}}'. \end{equation}

And the collision operator $L$ has the property (Sharipov 2015)

(2.33)\begin{equation} \left( {L[ \alpha ],\beta } \right) = \left( {L[ \beta ],\alpha } \right), \end{equation}

where $\alpha$ and $\beta$ are two arbitrary functions of velocity.

The determination of ${a_{0w}}$ and ${a_0}$ is divided into two steps.

1. (i) First, taking the scalar product of both sides of (2.22) with respect to $u'$, we have

   (2.34)\begin{gather} \text{left-hand side} = \left( {u',v'\,\frac{{\partial \varPhi }}{{\partial y}}} \right) = \frac{\partial }{{\partial y}}\left( {u'v',\varPhi } \right), \end{gather}
   (2.35)\begin{gather} \text{right-hand side} = \left( {u',L[ \varPhi ]} \right) = \left( {L[ {u'} ],\varPhi } \right) = 0. \end{gather}
   Note that $u'$ and $v'$ are peculiar velocities. Since $\partial \bar {u}/\partial y$ is sufficiently small, we can assume that the peculiar velocities are independent of the spatial coordinates. Therefore, in (2.34), $u'$ and $v'$ can enter inside the derivative sign with respect to $y$, allowing the derivative sign to be taken outside of the integration. Combining the above two equations, we obtain
   (2.36)\begin{equation} \frac{\partial }{{\partial y}}\left( {u'v',\varPhi } \right) = 0. \end{equation}
   Thus in the $y$-direction,
   (2.37)\begin{equation} \left( {u'v',\varPhi } \right) = {\rm const.} \end{equation}

2. (ii) Second, taking the scalar product of both sides of (2.22) with respect to $(B/ 2RT_0)u'v'$, we have

   (2.38) \begin{gather} \text{left-hand side} = \left( {B\,\frac{1}{{2RT_0}}\,u'v',v'\,\frac{{\partial \varPhi }}{{\partial y}}} \right) = \frac{\partial }{{\partial y}}\left( {B\, \frac{1}{{2RT_0}}\,u'v'^2,\varPhi } \right), \end{gather}
   (2.39) \begin{gather} \text{right-hand side} = \left( {B\,\frac{1}{{2RT_0}}\,u'v',L[ \varPhi ]} \right) = \left( {L\left[ {B\,\frac{1}{{2RT_0}}\,u'v'} \right],\varPhi } \right)\nonumber\\ = \left( {\frac{1}{{RT_0}}\,u'v',\varPhi } \right). \end{gather}
   Note that in (2.38), we also assume that the peculiar velocities and the equilibrium temperature ${T_0}$ are independent of the spatial coordinates. If we assume that $\varPhi$ takes the form as (2.26)–(2.28), then (2.39) equals zero. Combining the above two equations, we obtain
   (2.40)\begin{equation} \frac{\partial }{{\partial y}}\left( {B\,\frac{1}{{2RT_0}}\,u'v'^2,\varPhi } \right) = 0. \end{equation}
   Thus in the $y$-direction,
   (2.41)\begin{equation} \left( {B\,\frac{1}{{2RT_0}}\,u'v'^2,\varPhi } \right) = {\rm const.} \end{equation}
   The scalar function $B$ can be expanded in the Sonine polynomials as
   (2.42)\begin{equation} B = \sum_{r = 1}^\infty {{b_r}S_{5/2}^{( {r - 1} )}\left[ {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}}} \right]} . \end{equation}
   If we truncate to the first term, then the scalar function $B$ can be approximated as
   (2.43)\begin{equation} B \approx {b_1}, \end{equation}
   which is a constant independent of velocity. In this manner, we can extract $B /{{2RT_0}}$ outside of the scalar product (2.41). Thus in the $y$-direction,
   (2.44)\begin{equation} \left( {u'v{'^2},\varPhi } \right) \approx {\rm const.} \end{equation}

From (2.37) and (2.44), the expressions for ${a_{0w}}$ and ${a_0}$ can be determined. After simplifying, we obtain

(2.45)\begin{align} {a_{0w}} &= \frac{\mu }{p}\,\sqrt {\rm \pi}\,\frac{{\partial \bar{u} }}{{\partial y}}\, \frac{1}{{{\sigma _t}}}\left[ {( {1 - {\sigma _t}} )\left( {{A_0} + {A_1}\,\frac{{\sqrt 2 }}{{\sqrt {{\rm \pi} RT_0} }}} \right) + 1} \right], \end{align}

(2.46) $$\begin{align} a_0 &= \left\{ \frac{2 - \sigma _t}{2} \, \frac{( 1 - \sigma _t )}{{{\sigma _t}}}\left[\left( {A_0} + {A_1}\, \dfrac{{\sqrt 2 }}{{\sqrt {{\rm \pi} RT_0} }} \right) + 1\right]\right.\notag\\ &\quad + \left. \frac{1}{{\rm \pi} }\left[ {{( }{1 - {\sigma _t}} {)}\left( {2{A_0} + {A_1}\,\frac{{\sqrt {\rm \pi} }}{{\sqrt {2RT_0} }}} \right) + 2} \right] \vphantom{\left\{ \frac{2 - \sigma _t}{2} \, \frac{( 1 - \sigma _t )}{{{\sigma _t}}}\left[\left( {A_0} + {A_1}\, \dfrac{{\sqrt 2 }}{{\sqrt {{\rm \pi} RT_0} }} \right) + 1\right]\right.}\right\}\nonumber\\ &\quad \times \frac{\mu }{p}\,\sqrt {\rm \pi}\,\frac{{\partial \bar{u} }}{{\partial y}}. \end{align}$$

As per Appendix A, we determine ${A_0}$ as $1 - {\alpha _n}$ and ${A_1}$ as ${\sqrt{({{\rm \pi} RT_0} }}/{{\sqrt 2 }}){\alpha _n}$. Therefore, the expression for ${a_0}$ can be simplified as

(2.47)\begin{equation} {a_0} = \left\{ {\frac{{2 - {\sigma _t}}}{{{\sigma _t}}}\left( {1 + \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\sigma _t}} \right) - \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\alpha _n}( {1 - {\sigma _t}} )} \right\} \frac{\mu }{p}\,\sqrt {\rm \pi}\,\frac{{\partial \bar{u} }}{{\partial y}}. \end{equation}

Thus far, the distribution function at the outer edge of the Knudsen layer that considers the collision effect has been determined. Finally, the velocity slip can be obtained by taking the moment of the distribution function, i.e.

(2.48)\begin{align} {u_s} &= \frac{1}{\rho }\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {mu'{f^{( 0 )}}\left( {1 - \frac{{4\mu {\beta ^4}}}{\rho }\,u'v'\,\frac{{\partial \bar{u} }}{{\partial y}} + {a_0}\,\frac{{u'}}{{\sqrt {2RT_0} }}} \right)} } } \,{\rm d}{\boldsymbol{c}}'\nonumber\\ &= {a_0} \, \sqrt {\frac{{RT_0}}{2}} \nonumber\\ &= \left[ {\frac{{2 - {\sigma _t}}}{{{\sigma _t}}}\left( {1 + \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\sigma _t}} \right) - \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\alpha _n}( {1 - {\sigma _t}} )} \right] \frac{\mu }{p}\,\sqrt {\rm \pi}\,\frac{{\partial \bar{u} }}{{\partial y}}\,\sqrt {\frac{{RT_0}}{2}} . \end{align}

Notably, this slip coefficient formula is dependent not only on the TMAC but also on the NEAC, although the dependence is minimal.

2.4. The temperature jump problem

As shown in figure 3, the temperature profile within the Knudsen layer near a surface exhibits a similarity to the velocity profile within the Knudsen layer near a surface. Neglecting changes in velocity, the distribution function can be approximated to first order as

(2.49) \begin{equation} {f^ \pm } = {f_0}\left[ {1 - \frac{{4\kappa {\beta ^2}}}{{5nk}}\left( {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right)v'\,\frac{{\partial \ln T}}{{\partial y}} + {\varPhi ^ \pm }} \right]. \end{equation}

Figure 3. Schematic of the temperature profile of the Knudsen layer near a surface. Here, $\lambda$ is the molecular mean free path, ${T_w}$ is the temperature of the surface, $T( 0 )$ is the actual temperature of gas at the surface, and ${T_\infty }$ is the ‘fictitious’ temperature of gas that would be needed to ensure that an NSF solution (dashed line) provides an accurate prediction (solid line) beyond the Knudsen layer.

If Maxwell's assumption was adopted, then ${\varPhi ^ \pm }$ takes the form

(2.50)$$\begin{gather} {\varPhi ^ - }( {{\boldsymbol{c}},0} ) = {a_1}\left( {\frac{{{\boldsymbol{c}}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right), \end{gather}$$

(2.51)$$\begin{gather}{\varPhi ^ \pm }( {{\boldsymbol{c}},\infty } ) = {a_1}\left( {\frac{{{\boldsymbol{c}}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right). \end{gather}$$

Considering the influence of the Knudsen layer on the incident distribution function, two independent parameters are introduced to obtain a more accurate distribution function, i.e.

(2.52)$$\begin{gather} {\varPhi ^ - }( {{\boldsymbol{c}},0} ) = {a_{1w}}\left( {\frac{{{\boldsymbol{c}}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right), \end{gather}$$

(2.53)$$\begin{gather}{\varPhi ^ \pm }( {{\boldsymbol{c}},\infty } ) = {a_1}\left( {\frac{{{\boldsymbol{c}}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right). \end{gather}$$

Notably, compared to the macroscopic variables obtained from the original distribution function, the mean velocity ${\bar {c} _i}$ determined by (2.49), (2.52) and (2.53) is unchanged according to (2.1). However, the density and temperature are different from those in the main flow due to the existence of the Knudsen layer.

Next, substituting (2.52) into (2.4) and using the CLL scattering kernel, after carrying out straightforward calculations, an expression for ${\varPhi ^ + }( {{\boldsymbol {c}},0} )$ corresponding to the discontinuous term of reflected molecules at the surface can be obtained:

(2.54) \begin{align} {\varPhi ^ + }( {{\boldsymbol{c}},0} ) &=- \frac{\mu }{{2pRT_0}}\,\frac{{\partial \ln T}}{{\partial y}}\left[ {( {{A_0}( {1 - {\alpha _t}} ) + 1} )v'{\boldsymbol{c}}{'^2} + {A_1}( {1 - {\alpha _t}} ){\boldsymbol{c}}{'^2} - ( {{A_0}( {1 - {\alpha _t}} ) - {B_0}} )v{'^3}} \right.\nonumber\\ &\quad - \left.( {{A_1}( {1 - {\alpha _t}} ) - {B_1}} )v{'^2} + ( {{A_0}2RT_0{\alpha _t} - {A_0}5RT_0 + {B_2} - 5RT_0} )v' \right.\nonumber\\ &\quad + \left. {A_1}2RT_0{\alpha _t} - {A_1}5RT_0 + {B_3} \vphantom{v'{\boldsymbol{c}}{'^2}} \right]\nonumber\\ &\quad + \frac{{{a_{1w}}}}{{2RT_0}}\left[ {u{'^2}( {1 - {\alpha _t}} ) + w{'^2}( {1 - {\alpha _t}} ) + v{'^2}( {1 - {\alpha _n}} ) + 2RT_0\left( {{\alpha _t} + {\alpha _n} - \frac{5}{2}} \right)} \right]. \end{align}

In the corresponding calculation, integrals from Appendix A are used, with ${A_0}$, ${A_1}$ and ${B_0}\unicode{x2013} {B_3}$ incorporated as coefficients that are already determined.

Using the same process as in § 2.3, the determination of ${a_{1w}}$ and ${a_1}$ from (2.22) is divided into two steps.

1. (i) First, taking the scalar product of both sides of (2.22) with respect to $\boldsymbol {c}{'^2}$, we have

   (2.55)\begin{gather} \text{left-hand side} = \left( {{\boldsymbol{c}}{'^2},v'\,\frac{{\partial \varPhi }}{{\partial y}}} \right) = \frac{\partial }{{\partial y}}\left( {{\boldsymbol{c}}{'^2}v',\varPhi } \right), \end{gather}
   (2.56)\begin{gather} \text{right-hand side} = \left( {{\boldsymbol{c}}{'^2},L[ \varPhi ]} \right) = \left( {L[ {{\boldsymbol{c}}{'^2}} ],\varPhi } \right) = 0. \end{gather}
   Combining the above two equations, we obtain
   (2.57)\begin{equation} \frac{\partial }{{\partial y}}\left( {{\boldsymbol{c}}{'^2}v',\varPhi } \right) = 0. \end{equation}
   Thus in the $y$-direction,
   (2.58)\begin{equation} \left( {{\boldsymbol{c}}{'^2}v',\varPhi } \right) = {\rm const.} \end{equation}

2. (ii) Second, taking the scalar product of both sides of (2.22) with respect to $A {{v'}}/{{\sqrt {2RT_0} }}$, we have

   (2.59)\begin{gather} \text{left-hand side} = \left( {A \, \frac{{v'}}{{\sqrt {2RT_0} }},v'\,\frac{{\partial \varPhi }}{{\partial y}}} \right) = \frac{\partial }{{\partial y}}\left( {A \, \frac{1}{{\sqrt {2RT_0} }}\,v{'^2},\varPhi } \right), \end{gather}
   (2.60)$$\begin{gather} \text{right-hand side} = \left( {A \, \frac{{v'}}{{\sqrt {2RT_0} }},L[ \varPhi ]} \right) = \left( {L\left[ {A \, \frac{{v'}}{{\sqrt {2RT_0} }}} \right],\varPhi } \right) \nonumber\\ = \left( {\frac{{v'}}{{\sqrt {2RT_0} }}\left( {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right),\varPhi } \right). \end{gather}$$
   Note that in (2.55) and (2.59), the assumption that peculiar velocities and the equilibrium temperature ${T_0}$ are independent of spatial coordinates has been used. If we assume that $\varPhi$ takes the form as (2.52)–(2.54), then (2.60) equals zero. Combining the above two equations, we obtain
   (2.61)\begin{equation} \frac{\partial }{{\partial y}}\left( {A \, \frac{1}{{\sqrt {2RT_0} }}\,v{'^2},\varPhi } \right) = 0. \end{equation}
   Thus in the $y$-direction,
   (2.62)\begin{equation} \left( {A \, \frac{1}{{\sqrt {2RT_0} }}\,v{'^2},\varPhi } \right) = {\rm const.} \end{equation}
   The scalar function $A$ can be expanded in the Sonine polynomials:
   (2.63)\begin{equation} A = \sum_{r = 0}^\infty {{a_r}S_{3/2}^{( r )}\left[ {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}}} \right]} . \end{equation}
   If we truncate to first order, then the scalar function $A$ can be approximated as
   (2.64) \begin{align} A &\approx {a_1} \, S_{3/2}^{( 1 )}\left[ {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}}} \right]\nonumber\\ &= {a_1} \left[ { - \left( {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right)} \right]. \end{align}
   Substituting (2.64) into (2.62), we obtain that in the $y$-direction,
   (2.65)\begin{equation} \left( {v{'^2}\left( {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right),\varPhi } \right) \approx {\rm const.} \end{equation}

From (2.58) and (2.65), the expressions for ${a_{1w}}$ and ${a_1}$ can be determined. After simplifying, we obtain

(2.66) \begin{align} {a_{1w}} &= \frac{\mu }{p}\,\frac{{\partial \ln T}}{{\partial y}}\,\sqrt {2{\rm \pi} RT_0} \times \frac{1}{{8( {{\alpha _t} + {\alpha _n}} )}}\left[ ( {4{\alpha _t} + 11} ){A_0} + 8( {{\alpha _t} + 2} )\,\frac{{{A_1}}}{{\sqrt {2{\rm \pi} RT_0} }} - 21{B_0} \right. \notag\\ &\quad - \left. 24\,\frac{{{B_1}}}{{\sqrt {2{\rm \pi} RT_0} }} - 5\,\frac{{{B_2}}}{{RT_0}} - 8\,\frac{{{B_3}}}{{RT_0}}\,\frac{1}{{\sqrt {2{\rm \pi} RT_0} }} - 10 \right] , \end{align}

(2.67) \begin{align} {a_1} &= \sqrt {\rm \pi}\,\frac{\mu }{p}\,\frac{{\partial \ln T}}{{\partial y}}\,\sqrt {RT_0}\times\frac{{10 - 2{\alpha _t} - 3{\alpha _n}}}{{4( {{\alpha _t} + {\alpha _n}} )}} \nonumber\\ &\qquad \times \left[ \frac{{\begin{array}{@{}c@{}}( {4{\alpha _t} + 11} )\,\dfrac{{{A_0}}}{2} + 4( {{\alpha _t} + 2} )\,\dfrac{{{A_1}}}{{\sqrt {2{\rm \pi} RT_0} }} - \dfrac{{21}}{2}\,{B_0} - 12\,\dfrac{{{B_1}}}{{\sqrt {2{\rm \pi} RT_0} }} - \dfrac{5}{2}\,\dfrac{{{B_2}}}{{RT_0}} - 4\,\dfrac{{{B_3}}}{{RT_0}}\,\dfrac{1}{{\sqrt {2{\rm \pi} RT_0} }} - 5\end{array}}}{{10}} \right.\nonumber\\ &\qquad + \left. {\frac{{2( {{\alpha _t} + {\alpha _n}} )}}{{5{\rm \pi} }} \, \frac{{\begin{array}{@{}c@{}}( {4{\alpha _t} - 1} ){A_0} + ( {{\alpha _t} - 1} ){A_1}\,\dfrac{{\sqrt {2{\rm \pi} } }}{{\sqrt {RT_0} }} - 12{B_0} - \dfrac{3}{2}\,{B_1}\,\dfrac{{\sqrt {2{\rm \pi} } }}{{\sqrt {RT_0} }} - \dfrac{{{B_2}}}{{RT_0}} - 13\end{array}}}{{10 - 2{\alpha _t} - 3{\alpha _n}}}} \right] . \end{align}

As per Appendix A, it is proper to take ${A_0}$ as $1 - {\alpha _n}$, ${A_1}$ as $({{\sqrt {{\rm \pi} RT_0} }}/{{\sqrt 2 }}){\alpha _n}$, ${B_0}$ as $1 - {\alpha _n}$, ${B_3}$ as $3\sqrt {{{\rm \pi} }/{2}} \, {( {RT_0} )^{3/2}}{\alpha _n}$, and ${B_1} = {B_2} = 0$. As a result, the expression for ${a_1}$ can be transformed to

(2.68) \begin{align} {a_1} & = \sqrt {\rm \pi} \, \frac{\mu }{p}\,\frac{{\partial \ln T}}{{\partial y}}\,\sqrt {RT_0} \nonumber\\ &\quad \times \left\{\frac{{10 - 2{\alpha _t} - 3{\alpha _n}}}{{4( {{\alpha _t} + {\alpha _n}} )}}\left[ {\frac{{2{\alpha _t} + 3{\alpha _n} - 10}}{{10}} + \frac{2}{{5{\rm \pi} }}\,\frac{{( {{\alpha _t} + {\alpha _n}} )( {4{\alpha _t} + 9{\alpha _n} - 26} )}}{{10 - 2{\alpha _t} - 3{\alpha _n}}}} \right]\right.\notag\\ &\quad + \left. \frac{{4 - {\rm \pi}}}{{10{\rm \pi} }}\,{\alpha _n}( {1 - {\alpha _t}} ) \right\}. \end{align}

Thus far, the distribution function at the outer edge of the Knudsen layer that considers the collision effect has been determined. Finally, the temperature jump can be obtained by taking the moment of the distribution function, i.e.

(2.69) \begin{align} &\frac{{{T_\infty } - {T_w}}}{{{T_w}}}\nonumber\\ &\quad = \frac{2}{{3{n_0}}}\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {\left( {\frac{{\boldsymbol{c}{'^2}}}{{2R{T_0}}} - \frac{3}{2}} \right){f^{( 0 )}}\left( {\frac{\mu }{p}\,v'\left( {\frac{{\boldsymbol{c}{'^2}}}{{2R{T_0}}} - \frac{5}{2}} \right)\frac{{\partial \ln T}}{{\partial y}} + {a_1}\left( {\frac{{\boldsymbol{c}{'^2}}}{{2R{T_0}}} - \frac{5}{2}} \right)} \right)} } } \,{\rm d}{\boldsymbol{c}}\nonumber\\ &\quad = {a_1}\nonumber\\ &\quad = \left[ {\frac{{10 - 2{\alpha _t} - 3{\alpha _n}}}{{4( {{\alpha _t} + {\alpha _n}} )}}\left( {\frac{{2{\alpha _t} + 3{\alpha _n} - 10}}{{10}} + \frac{2}{{5{\rm \pi} }}\,\frac{{( {{\alpha _t} + {\alpha _n}} )( {4{\alpha _t} + 9{\alpha _n} - 26} )}}{{10 - 2{\alpha _t} - 3{\alpha _n}}}} \right) + \frac{{4 - {\rm \pi}}}{{10{\rm \pi} }}\,{\alpha _n}( {1 - {\alpha _t}} )} \right] \nonumber\\ &\qquad \times \sqrt {\rm \pi} \, \frac{\mu }{p}\,\frac{{\partial \ln T}}{{\partial y}}\,\sqrt {RT_0} . \end{align}

3. Results and validation

Equations (2.48) and (2.69) constitute the primary contributions of this paper. These two formulae explicitly show the dependence of slip coefficients on TEAC and NEAC with very concise forms, facilitating their use in practical applications. In this section, a comprehensive verification and evaluation of these results is reported. In § 3.1, we will compare the slip coefficients and jump coefficients obtained from (2.48) and (2.69) with the numerical solutions provided by Sharipov (2011). In § 3.2, we will implement these boundary conditions in CFD and compare the results with DSMC to further assess the accuracy of these formulae.

3.1. Comparison of the slip coefficients

To facilitate comparison, we express the boundary conditions for the velocity slip and temperature jump in a general form as

(3.1)$$\begin{gather} {u_s} = {C_m}\,\frac{\mu }{{\beta p}}\,\frac{{\partial \bar{u} }}{{\partial y}}, \end{gather}$$

(3.2)$$\begin{gather}T - {T_w} = {C_t}\,\frac{\mu }{{\beta p}}\,\frac{{\partial T}}{{\partial y}}, \end{gather}$$

where ${C_m}$ and ${C_t}$ are the velocity slip and temperature jump coefficients, respectively, $\mu$ is the molecular viscosity, and $p$ is the local pressure. The relationship between $\mu$ and $p$ can be expressed as $\mu = \tau p$, where $\tau$ is the mean collision time of molecules.

Transforming the derived theoretical results (2.48) and (2.69) into the above two forms, we can determine the slip coefficients as

(3.3)$$\begin{gather} {C_m} = \left[ {\frac{{2 - {\sigma _t}}}{{{\sigma _t}}}\left( {1 + \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\sigma _t}} \right) - \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\alpha _n}( {1 - {\sigma _t}} )} \right] \frac{{\sqrt {\rm \pi}}}{2}, \end{gather}$$

(3.4)$$\begin{gather} {C_t} = \frac{3}{8} \, \frac{{10 - 2{\alpha _t} - 3{\alpha _n}}}{{ {{\alpha _t} + {\alpha _n}}}}\left( {1 - \frac{{2{\alpha _t} + 3{\alpha _n}}}{{10}} + \frac{2}{{5{\rm \pi} }}\,\frac{{( {26 - 4{\alpha _t} - 9{\alpha _n}} )( {{\alpha _t} + {\alpha _n}} )}}{{10 - 2{\alpha _t} - 3{\alpha _n}}}} \right) \sqrt {\rm \pi}\nonumber\\ + \frac{3}{2} \, \frac{{{\rm \pi} - 4}}{{10{\rm \pi} }}\,{\alpha _n}( {1 - {\alpha _t}} )\,\sqrt {\rm \pi}. \end{gather}$$

In order to facilitate comparison and give a clear picture of the position of this work, we summarize the first-order velocity slip and temperature jump coefficients obtained from the existing literature, derived from both the Maxwell scattering model and the CLL scattering model, in tables 1 and 2. Notably, the results of Struchtrup (2013) for the CLL scattering model are not accurate enough, while the corrected results of Zhang et al. (2021) rely on an empirical correction. In contrast, our newly developed slip models are more physically grounded, and are accurate enough to be used in practical applications.

Table 1. Velocity slip boundary coefficients based on the Maxwell scattering model and CLL scattering model.

Table 2. Temperature jump boundary coefficients based on the Maxwell scattering model and CLL scattering model.

Interestingly, assuming that the ACs are identical in the tangential and normal directions, (3.3) and (3.4) will be simplified to

(3.5)$$\begin{gather} {C_m^*} = \frac{{2 - \alpha }}{\alpha }\,( {1 + 0.1366\alpha } ) \, \frac{{\sqrt {\rm \pi}}}{2} - 0.1211\alpha ( {1 - \alpha } ), \end{gather}$$

(3.6)$$\begin{gather}{C_t^*} = \frac{{2 - \alpha }}{\alpha }\,( {1 + 0.1621\alpha } ) \, \frac{{15}}{{16}}\,\sqrt {\rm \pi}- 0.073\alpha ( {1 - \alpha } ), \end{gather}$$

which are very close to the slip coefficients obtained by Loyalka (1971b, 1989) for the Maxwell scattering model:

(3.7)$$\begin{gather} C_m' = \frac{{2 - \alpha }}{\alpha }\,( {1 + 0.1366\alpha } ) \, \frac{{\sqrt {\rm \pi}}}{2}, \end{gather}$$

(3.8)$$\begin{gather}C_t' = \frac{{2 - \alpha }}{\alpha }\,( {1 + 0.1621\alpha } ) \, \frac{{15}}{{16}}\,\sqrt {\rm \pi}. \end{gather}$$

To assess the accuracy of the slip coefficients, we compared the results obtained from (3.3) and (3.4) with the numerical solutions obtained by Sharipov (2011) by solving the simplified Boltzmann equation, known as the SBGK model. Notably, Struchtrup (2013) also used the moment method to derive slip coefficients for the CLL scattering model. Therefore, we have also compared our results with those of Struchtrup (2013). Figure 4 shows that the velocity slip coefficient remains relatively constant with the NEAC, while the temperature jump coefficient varies significantly. Additionally, both the velocity slip and temperature jump coefficients decrease as the AC increases. Our theoretically derived slip coefficients closely match values from numerical solutions, while the results of Struchtrup (2013) show some deviation. Further analysis indicates that for the Struchtrup (2013) model, the maximum relative error in the velocity slip coefficient is approximately 13 %, and the relative error in the temperature jump coefficient is approximately 17 %. For the present model, the maximum relative error in the velocity slip coefficient is less than 1.1 %, and the relative error in the temperature jump coefficient does not exceed 4.4 %.

Figure 4. Comparison of the theoretical slip coefficients derived in this work with those of Struchtrup (2013) and the numerical results of Sharipov (2011).

3.2. Comparison between CFD and DSMC

The DSMC method is capable of capturing both near-equilibrium transport mechanisms (viscosity, thermal conductivity and molecular diffusion) and non-equilibrium effects (velocity slip, temperature jump, etc.); it has been employed extensively to obtain standard results for validating slip boundary conditions (Bird 1994). We now incorporate the slip boundary conditions into CFD simulations, and subsequently compare the obtained results with those from the DSMC method. This scenario presents a real test of the slip boundary conditions that we propose.

The DSMC simulations in this work are based on the open-source software SPARTA (Plimpton et al. 2019). The collision model used is the variable hard sphere (VHS) model (Bird 1994). Parameters related to the VHS model are determined through the equations

(3.9)\begin{gather} \mu = {\mu _{ref}}{\left( {\frac{T}{{{T_{ref}}}}} \right)^\omega }, \end{gather}

(3.10)\begin{gather} {\mu _{ref}} = \frac{{15\sqrt {{\rm \pi} mk{T_{ref}}} }}{{2{\rm \pi} d_{ref}^2( {5 - 2\omega } )( {7 - 2\omega } )}}, \end{gather}

where ${T_{ref}}$ and ${d_{ref}}$ represent the reference values for temperature and diameter in the VHS model, respectively, $\omega$ is the viscosity index, and $k$ is the Boltzmann constant. The CFD simulations were performed using the open-source software OpenFOAM (Jasak, Jemcov & Tukovic 2007), specifically version 10. And the density-based compressible solver rhoCentralFOAM (Greenshields et al. 2010) is chosen for use.

We extend the Maxwell & Smoluchowski slip boundary conditions provided in rhoCentralFOAM to incorporate the slip boundary conditions derived herein. To eliminate the interference of other factors, the calculation of viscosity coefficients in CFD utilizes the power-law model, which is defined in (3.9) and (3.10), with parameters aligned with the DSMC set-up. Argon (Ar) is chosen as the working gas in the simulations, and the relevant parameters are listed in table 3. The parameters corresponding to 1000 K refer to the settings of Lofthouse et al. (2008).

Table 3. Parameters used in CFD and DSMC.

We first examined the high-speed Couette flow, allowing a comparison of both the velocity slip and temperature jump results from the velocity and temperature fields. After that, we conducted simulations of a two-dimensional hypersonic cylinder flow. Notably, these benchmarks have been thoroughly investigated by Sharipov & Strapasson (2013), Sharipov (2015) and Sharipov & Volkov (2022), who used CLL scattering model as gas–surface interaction in combination with ab initio potential as an alternative intermolecular collision model. In our study, these examples are used to verify the accuracy and reliability of the slip boundary conditions that we derived.

3.2.1. Couette flow

In our simulation, the temperature of the plates, the initial temperature of Ar, and ${T_{ref}}$ are taken as 273 K. The velocities of the upper plate and the lower plate are set as $337.29\ \textrm {m}\ \textrm {s}^{-1}$ and $-337.29\ \textrm {m}\ \textrm {s}^{-1}$, respectively, which are equal to the most probable thermal speed of molecules at 273 K. The separation distance $L$ between the two plates is 0.001 m.

We consider a moderate rarefied case ($Kn =0.05$). The results are compared between CFD with no-slip and CLL slip boundary conditions derived by Struchtrup (2013), and the present model. Note that Struchtrup (2013) also utilized the moment method to derive the slip model. The DSMC results serve as a benchmark. For the DSMC and CLL slip simulations, two scenarios for the ACs are chosen: one with ${\alpha _t} = 1.0$, ${\alpha _n} = 0.5$, and another with ${\alpha _t} = 0.5$, ${\alpha _n} = 1.0$. As shown in figure 5, using a no-slip boundary condition would lead to significant deviation from the DSMC result. The velocity profiles from both CLL slip models closely align with the DSMC results, while for the temperature profile, the results using the present model are more accurate in the main region. This difference is particularly pronounced when ${\alpha _t}$ is large. In the vicinity of the wall, the present model seems to give a larger temperature jump than the actual value. This is not surprising, as all slip boundary conditions provide ‘fictitious’ values instead of ‘real’ values at the wall, in order to ensure that an NSF solution provides an accurate prediction beyond the Knudsen layer, as depicted in figure 3.

Figure 5. Velocity and temperature profiles for Couette flow with a global $Kn$ of 0.05. Comparison of the CFD results with no-slip, Maxwell slip and CLL slip boundary conditions derived herein. The DSMC results serve as a benchmark, for (a,c) ${\alpha _t} = 1.0$, ${\alpha _n} = 0.5$, and (b,d) ${\alpha _t} = 0.5$, ${\alpha _n} = 1.0$.

3.2.2. Two-dimensional hypersonic flow around a cylinder

For the two-dimensional hypersonic flow around a cylinder, we refer to the case of Lofthouse et al. (2008), and the relevant parameters are shown in table 4. The working gas is Ar, and the global $Kn$ is based on the cylinder diameter $D$. The incoming velocity is set as $2624.1\ \textrm {m}\ \textrm {s}^{-1}$, which corresponds to $Ma=10$. Here, ${T_\infty }$ and ${p_\infty }$ are the temperature and pressure at the initial state, respectively. The temperature of the cylinder surface is set as 500 K, and the reference temperature needed for the calculation of viscosity is chosen as 1000 K. Note that the curved wall effect and the thermal creep effect need to be considered in this case. Following Lockerby et al. (2004) and Loyalka (1971a), the velocity slip boundary conditions are modified as

(3.11)\begin{align} {u_s} = \left[ {\frac{{2 - {\sigma _t}}}{{{\sigma _t}}}\left( {1 + \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\sigma _t}} \right) - \frac{{4 - {\rm \pi}}}{{2{\rm \pi} }}\,{\alpha _n}( {1 - {\sigma _t}} )} \right] \frac{\mu }{p}\,\sqrt {\rm \pi}\left( {\frac{{\partial \bar{u} }}{{\partial y}} + \frac{{\partial \bar{v} }}{{\partial x}}} \right)\sqrt {\frac{{RT}}{2}} + \frac{3}{4}\,\frac{\mu }{{\rho T}}\,\frac{{\partial T}}{{\partial x}} \end{align}

for the present result, and

(3.12)\begin{equation} {u_s} = \frac{{2 - {\sigma _t}}}{{{\sigma _t}}} \, \frac{\mu }{p}\,\sqrt {\rm \pi}\left( {\frac{{\partial \bar{u} }}{{\partial y}} + \frac{{\partial \bar{v} }}{{\partial x}}} \right)\sqrt {\frac{{RT}}{2}} + \frac{3}{4}\,\frac{\mu }{{\rho T}}\,\frac{{\partial T}}{{\partial x}} \end{equation}

for the Struchtrup result.

Table 4. Simulation parameters for the case of two-dimensional hypersonic flow around a cylinder.

A comparison of temperature contours is presented in figure 6. The CFD results, using no-slip and CLL slip boundary conditions derived in this work, are depicted in the upper part of each plot, while the DSMC results are shown in the lower part. From the results, we can find that use of the CLL slip model performs well, especially in the wake regions of the cylinder. Note that discrepancies persist at the shock front and in the boundary layer. In these regions, the local $Kn$ exceeds 0.1, although the global $Kn$ is 0.01. Therefore, the inaccuracy of CFD results in these specific regions is caused mainly by the linear constitutive relations used in the NSF equations, and the selection of slip boundary conditions is only one contributing factor.

Figure 6. Comparison of temperature contours of CFD results for the circular cylinder flow using different boundary conditions with DSMC, with a global $Kn$ of 0.01, for (a,b) ${\alpha _t} = 1.0$, ${\alpha _n} = 0.5$, and (c,d) ${\alpha _t} = 0.5$, ${\alpha _n} = 1.0$, for (a,c) no-slip, (b,d) CLL slip.

At the stationary point, the corresponding results are presented in table 5, and show that using the CLL slip model derived in this work has quite high accuracy. For the case $Kn = 0.01$, the heat flux predicted using the CLL slip model proposed by Struchtrup (2013) has relative error $2.8\,\%$, while the relative errors of the results predicted using the present model do not exceed $1.6\,\%$. For the case $Kn = 0.05$, the heat flux predicted using the CLL slip model proposed by Struchtrup (2013) has relative error $9.3\,\%$, while the relative errors of the results predicted using the present model do not exceed $4.9\,\%$.

Table 5. Comparison of heat flux $\left ( {q/\left ( {\tfrac {1}{2}{\rho _\infty }V_\infty ^3} \right )} \right )$ at stationary points around a cylinder ($Ma = 10$) with DSMC results.

4. Conclusions and discussions

In this study, we presented a derivation of velocity slip and temperature jump boundary conditions based on the CLL scattering model. To incorporate the surface's impact on the distribution function, a Knudsen layer correction term was introduced. In addition, two moments were adopted to address the collision effect in the Knudsen layer. We demonstrated that the results obtained using the derived slip coefficients closely match the numerical results obtained by solving the Boltzmann equation. Moreover, we applied the slip boundary conditions to the framework of NSF equations for two benchmarks, including Couette flow and two-dimensional hypersonic flow around a cylinder, to compare the CFD results and DSMC results. The slip boundary conditions derived in this work perform well in simulating flows in the slip regime.

In contrast to the existing slip boundary conditions that utilize a single AC, the slip boundary conditions developed in this work incorporate two ACs: the TEAC and NEAC. It is noteworthy that the values of the TEAC and NEAC are dependent on many factors, including surface materials, surface temperature, gas temperature and gas velocity, among other influential variables, and they are typically unequal in most scenarios. For instance, Yamamoto et al. (2006) reported that in the case of nitrogen molecules impinging on a platinum surface contaminated with xenon molecules, the TEAC and NEAC are 0.52 and 0.61, respectively. The discrepancies between these values lead to a difference in temperature jump coefficients of approximately 13.2 % comparing values predicted by the Maxwell slip model (3.8) and the CLL slip model (3.4). Spijker et al. (2010) found that in the case of argon molecules impinging on a clean platinum surface at 300 K, the TEAC and NEAC are 0.28 and 0.46, respectively. The discrepancies between these values result in a significant difference in temperature jump coefficients predicted by these two models, reaching 41.0 %. We believe that the slip boundary conditions based on the CLL slip model are promising for providing more reasonable predictions than the Maxwell slip model, as long as the independent values of the TEAC and NEAC can be measured accurately.

Nowadays, while it is increasingly feasible to utilize the numerical solution of the Boltzmann equation for determining the slip coefficient, this method necessitates solving the Boltzmann equation again for each unique combination of TMAC and NEAC scenarios, demanding substantial computational resources. Conversely, the direct implementation of the proposed slip models in this study offers considerable convenience for engineering applications. On the other hand, while the flourishing development of machine learning technology enables the discovery of fitting formulae using existing data (Zhang & Ma 2020; Ma et al. 2024), such formulae often lack a strong physical background and may exhibit significant deviations when applied beyond the scope of the training set. In contrast, the methodology utilized in this study, while rooted in a classical approach, yields slip boundary condition models that possess clear physical meaning and robust applicability. Compared to the asymptotic theory developed by Sone et al. (2000) and Aoki et al. (2017), this approach is rather simple to implement. We demonstrated that by considering the collision effect in the Knudsen layer, sufficiently accurate slip boundary conditions could be achieved. Importantly, these slip boundary conditions are valid across the entire range of ACs. Among all slip models derived from the CLL scattering model using kinetic approaches thus far (Klinc & Kuščer 1972; McCormick 2005; Struchtrup 2013; Zhang et al. 2021; Basdanis et al. 2023), our model stands out for its clear physical foundation and high accuracy.

Several topics warrant further discussion.

First, in the course of our derivation, the CLL scattering model is employed. Despite its favourable characteristics, the CLL scattering model has certain flaws, and numerous researchers have dedicated significant efforts to improving the scattering model (Wu & Struchtrup 2017; Wang et al. 2021; Chen et al. 2023). The theoretical derivation method used in this work is not limited to the Maxwell or CLL scattering model; after a more advanced scattering kernel is expressed explicitly, the corresponding slip model can be derived.

Second, the scenario addressed in this study assumes that molecules experience complete scattering upon surface impact. Nonetheless, according to Myong's adsorption theory (Myong 2004; Myong et al. 2005), some molecules hitting the surface may not immediately scatter back but instead adsorb to it. Recently, Chen et al. (2023) introduced scattering kernels that take into account adsorption effects. Using a similar approach, slip models that account for both scattering and adsorption can be derived and deserve further investigation.

Third, we focus on only the first-order slip models in the present work (see (2.23) and (2.49)). On the other hand, Hadjiconstantinou (2003) and Hadjiconstantinou & Al-Mohssen (2005) argued that second-order slip models perform better in rarefied cases. If we take the second-order Chapman–Enskog expansion, then a second-order slip model corresponding to the CLL model can also be derived. This approach could be useful for exploring in future work.