1. Introduction
The rarefied gas effect is a vital and foundational concern demanding attention in hypersonic flows (Ivanov & Gimelshein Reference Ivanov and Gimelshein1998; Candler Reference Candler2019) and gas flows for microelectromechanical systems (MEMS) (Reese, Gallis & Lockerby Reference Reese, Gallis and Lockerby2003; Karniadakis, Beskok & Aluru Reference Karniadakis, Beskok and Aluru2006). 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 Reference Struchtrup2005; Sharipov Reference Sharipov2015). On the other hand, the widely used direct simulation Monte Carlo (DSMC) method (Bird Reference Bird1994), 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. Reference Fei, Zhang, Li and Liu2020; Feng et al. Reference Feng, Tian, Zhang, Fei and Wen2023). 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 Reference Lockerby and Reese2008; Lofthouse, Scalabrin & Boyd Reference Lofthouse, Scalabrin and Boyd2008; Gu & Emerson Reference Gu and Emerson2009). In the case of slightly rarefied flows, a widely adopted approach is applying suitable slip boundary conditions (Struchtrup Reference Struchtrup2005; Sharipov Reference Sharipov2015).
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 Reference Gökçen, MacCormack and Chapman1987; Myong Reference Myong2004; Le et al. Reference Le, White, Reese and Myong2012; Wang, Ou & Chen Reference Wang, Ou and Chen2023); the corresponding results are well documented in the reviews (Cao et al. Reference Cao, Sun, Chen and Guo2009; Zhang, Meng & Wei Reference Zhang, Meng and Wei2012; Akhlaghi, Roohi & Stefanov Reference Akhlaghi, Roohi and Stefanov2023). Classically, first-order boundary conditions encompass the Maxwell model (Maxwell Reference Maxwell1879) for velocity slip, and the Smoluchowski model (Smoluchowski von Smolan Reference Smoluchowski von Smolan1898) 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 Reference Hadjiconstantinou2003; Lockerby et al. Reference Lockerby, Reese, Emerson and Barber2004; Radtke et al. Reference Radtke, Hadjiconstantinou, Takata and Aoki2012; Zeng et al. Reference Zeng, Zhao, Jiang and Chen2023), 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 Reference Lord1992; Arkilic, Breuer & Schmidt Reference Arkilic, Breuer and Schmidt2001) 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 Reference Yamamoto, Takeuchi and Hyakutake2006; Liang, Li & Ye Reference Liang, Li and Ye2013) and molecular dynamics simulations (Spijker et al. Reference Spijker, Markvoort, Nedea and Hilbers2010; Sipkens & Daun Reference Sipkens and Daun2018) 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 Reference Grad1949; Struchtrup & Weiss Reference Struchtrup and Weiss2000; Torrilhon & Struchtrup Reference Torrilhon and Struchtrup2008; To et al. Reference To, Vu, Lauriat and Léonard2015; Li & Yang Reference Li and Yang2023), half-flux method (Patterson Reference Patterson1956; Shen Reference Shen2006; Zhang et al. Reference Zhang, Luan, Deng, Tian and Liang2021) and variational method (Loyalka Reference Loyalka1968; Klinc & Kuščer Reference Klinc and Kuščer1972; McCormick Reference McCormick2005). The method of directly solving the linearized Boltzmann equation has also been developed to obtain precise slip coefficients (Sharipov Reference Sharipov2003a, Reference Sharipov2011; Siewert Reference Siewert2003; Basdanis, Valougeorgis & Sharipov Reference Basdanis, Valougeorgis and Sharipov2023). Additionally, Sone, Ohwada & Aoki (Reference Sone, Ohwada and Aoki1989) and Ohwada, Sone & Aoki (Reference Ohwada, Sone and Aoki1989) 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. (Reference Aoki, Baranger, Hattori, Kosuge, Martalò, Mathiaud and Mieussens2017) 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 (Reference Maxwell1879). 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 (Reference Loyalka1971b) 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 Reference Cercignani and Lampis1971), struggles to describe backscattering (Basdanis, Tatsios & Valougeorgis Reference Basdanis, Tatsios and Valougeorgis2022), and does not account for the temperature independence of ACs, among other shortcomings.
A more advanced scattering model was proposed by Cercignani & Lampis (Reference Cercignani and Lampis1971) and later extended by Lord (Reference Lord1991, Reference Lord1995), 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 Reference Cercignani and Lampis1971), 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 Reference Sharipov2003b). 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 Reference Kalempa and Sharipov2020). Recently, Sharipov & Volkov (Reference Sharipov and Volkov2022) employed the CLL scattering model to evaluate the impact of ACs on the aerothermodynamics of a sphere based on ab initio interatomic potentials (Sharipov Reference Sharipov2022), 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 (Reference Struchtrup2013) 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 Reference Struchtrup2005). Klinc & Kuščer (Reference Klinc and Kuščer1972), as well as McCormick (Reference McCormick2005), derived results using the variational method. Nevertheless, these derivations involve higher-order ACs, rendering them challenging for practical application. Recently, Basdanis et al. (Reference Basdanis, Valougeorgis and Sharipov2023) obtained the interpolation expression of slip coefficients for the CLL scattering model using symbolic regression. Zhang et al. (Reference Zhang, Luan, Deng, Tian and Liang2021) obtained the CLL slip boundary condition via the half-flux method, yielding results consistent with those of Struchtrup (Reference Struchtrup2013). Zhang et al. (Reference Zhang, Luan, Deng, Tian and Liang2021) further refined their results by applying an empirical correction factor, akin to Loyalka's correction (Loyalka Reference Loyalka1971b) for the Maxwell slip boundary condition. While the results of Basdanis et al. (Reference Basdanis, Valougeorgis and Sharipov2023) and Zhang et al. (Reference Zhang, Luan, Deng, Tian and Liang2021) 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 (Reference Grad1949) 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. Reference Shan, Wang, Wang, Zhang and Guo2022; Qian, Wu & Wang Reference Qian, Wu and Wang2023). 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 Reference Ivchenko, Loyalka and Tompson2007). 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.
\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.
\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
\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 Reference Cercignani and Cercignani1988), i.e.
\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:
\begin{equation} R( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} ) \ge 0. \end{equation}Normalization:
\begin{equation} \int_{{c_n} > 0} {R( {{{\boldsymbol{c}}_i} \to {{\boldsymbol{c}}_r}} )} \,{\rm d}{\boldsymbol{c}} = 1. \end{equation}Reciprocity:
\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
\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
$$\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):
\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 (Reference Sharipov and Moldover2016), 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.
\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
\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
\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
\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
\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 Reference Chapman and Cowling1990). 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
\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 (Reference Chapman and Cowling1990), 
$\mu$ and 
$\kappa$ can be determined as
\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
\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
\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 Reference Sharipov2015), i.e.
\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
\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
\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
\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
$$\begin{gather} {\varPhi ^ - }( {{\boldsymbol{c}},0} ) = {a_0}\,\frac{{u'}}{{\sqrt {2RT_0} }}, \end{gather}$$
$$\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 Reference Loyalka1971b). To rectify this issue, we can assume that 
${\varPhi ^ \pm }$ take the form
$$\begin{gather} {\varPhi ^ - }\left( {{\boldsymbol{c}},0} \right) = {a_w}\frac{{u'}}{{\sqrt {2RT_0} }}, \end{gather}$$
$$\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 (Reference Loyalka1971b). 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
\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
\begin{gather} L[ \psi ] = 0, \end{gather}
\begin{gather} L\left[ {B\,\frac{1}{{2RT}}\,u'v'} \right] = \frac{1}{{RT}}\,u'v', \end{gather}
\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
\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 Reference Sharipov2015)
\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.
(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)Note that\begin{gather} \text{right-hand side} = \left( {u',L[ \varPhi ]} \right) = \left( {L[ {u'} ],\varPhi } \right) = 0. \end{gather}$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)Thus in the\begin{equation} \frac{\partial }{{\partial y}}\left( {u'v',\varPhi } \right) = 0. \end{equation}$y$-direction,
(2.37)\begin{equation} \left( {u'v',\varPhi } \right) = {\rm const.} \end{equation}(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)Note that in (2.38), we also assume that the peculiar velocities and the equilibrium temperature
\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}
${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)Thus in the\begin{equation} \frac{\partial }{{\partial y}}\left( {B\,\frac{1}{{2RT_0}}\,u'v'^2,\varPhi } \right) = 0. \end{equation}$y$-direction,
(2.41)The scalar function\begin{equation} \left( {B\,\frac{1}{{2RT_0}}\,u'v'^2,\varPhi } \right) = {\rm const.} \end{equation}$B$ can be expanded in the Sonine polynomials as
(2.42)If we truncate to the first term, then the scalar function\begin{equation} B = \sum_{r = 1}^\infty {{b_r}S_{5/2}^{( {r - 1} )}\left[ {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}}} \right]} . \end{equation}$B$ can be approximated as
(2.43)which is a constant independent of velocity. In this manner, we can extract\begin{equation} B \approx {b_1}, \end{equation}$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
\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}
$$\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
\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.
\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
\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
$$\begin{gather} {\varPhi ^ - }( {{\boldsymbol{c}},0} ) = {a_1}\left( {\frac{{{\boldsymbol{c}}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right), \end{gather}$$
$$\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.
$$\begin{gather} {\varPhi ^ - }( {{\boldsymbol{c}},0} ) = {a_{1w}}\left( {\frac{{{\boldsymbol{c}}{'^2}}}{{2RT_0}} - \frac{5}{2}} \right), \end{gather}$$
$$\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:
\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.
(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)Combining the above two equations, we obtain\begin{gather} \text{right-hand side} = \left( {{\boldsymbol{c}}{'^2},L[ \varPhi ]} \right) = \left( {L[ {{\boldsymbol{c}}{'^2}} ],\varPhi } \right) = 0. \end{gather}(2.57)Thus in the\begin{equation} \frac{\partial }{{\partial y}}\left( {{\boldsymbol{c}}{'^2}v',\varPhi } \right) = 0. \end{equation}$y$-direction,
(2.58)\begin{equation} \left( {{\boldsymbol{c}}{'^2}v',\varPhi } \right) = {\rm const.} \end{equation}(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)Note that in (2.55) and (2.59), the assumption that peculiar velocities and the equilibrium temperature$$\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}$$${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)Thus in the\begin{equation} \frac{\partial }{{\partial y}}\left( {A \, \frac{1}{{\sqrt {2RT_0} }}\,v{'^2},\varPhi } \right) = 0. \end{equation}$y$-direction,
(2.62)The scalar function\begin{equation} \left( {A \, \frac{1}{{\sqrt {2RT_0} }}\,v{'^2},\varPhi } \right) = {\rm const.} \end{equation}$A$ can be expanded in the Sonine polynomials:
(2.63)If we truncate to first order, then the scalar function\begin{equation} A = \sum_{r = 0}^\infty {{a_r}S_{3/2}^{( r )}\left[ {\frac{{\boldsymbol{c}{'^2}}}{{2RT_0}}} \right]} . \end{equation}$A$ can be approximated as
(2.64)Substituting (2.64) into (2.62), we obtain that in the
\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}
$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
\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}
\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
\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.
\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 (Reference Sharipov2011). 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
$$\begin{gather} {u_s} = {C_m}\,\frac{\mu }{{\beta p}}\,\frac{{\partial \bar{u} }}{{\partial y}}, \end{gather}$$
$$\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
$$\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}$$
$$\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 (Reference Struchtrup2013) for the CLL scattering model are not accurate enough, while the corrected results of Zhang et al. (Reference Zhang, Luan, Deng, Tian and Liang2021) 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
$$\begin{gather} {C_m^*} = \frac{{2 - \alpha }}{\alpha }\,( {1 + 0.1366\alpha } ) \, \frac{{\sqrt {\rm \pi}}}{2} - 0.1211\alpha ( {1 - \alpha } ), \end{gather}$$
$$\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 (Reference Loyalka1971b, Reference Loyalka1989) for the Maxwell scattering model:
$$\begin{gather} C_m' = \frac{{2 - \alpha }}{\alpha }\,( {1 + 0.1366\alpha } ) \, \frac{{\sqrt {\rm \pi}}}{2}, \end{gather}$$
$$\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 (Reference Sharipov2011) by solving the simplified Boltzmann equation, known as the SBGK model. Notably, Struchtrup (Reference Struchtrup2013) 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 (Reference Struchtrup2013). 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 (Reference Struchtrup2013) show some deviation. Further analysis indicates that for the Struchtrup (Reference Struchtrup2013) 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 (Reference Struchtrup2013) and the numerical results of Sharipov (Reference Sharipov2011).
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 Reference Bird1994). 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. Reference Plimpton, Moore, Borner, Stagg, Koehler, Torczynski and Gallis2019). The collision model used is the variable hard sphere (VHS) model (Bird Reference Bird1994). Parameters related to the VHS model are determined through the equations
\begin{gather} \mu = {\mu _{ref}}{\left( {\frac{T}{{{T_{ref}}}}} \right)^\omega }, \end{gather}
\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 Reference Jasak, Jemcov and Tukovic2007), specifically version 10. And the density-based compressible solver rhoCentralFOAM (Greenshields et al. Reference Greenshields, Weller, Gasparini and Reese2010) 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. (Reference Lofthouse, Scalabrin and Boyd2008).
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 (Reference Sharipov and Strapasson2013), Sharipov (Reference Sharipov2015) and Sharipov & Volkov (Reference Sharipov and Volkov2022), 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 (Reference Struchtrup2013), and the present model. Note that Struchtrup (Reference Struchtrup2013) 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. (Reference Lofthouse, Scalabrin and Boyd2008), 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. (Reference Lockerby, Reese, Emerson and Barber2004) and Loyalka (Reference Loyalka1971a), the velocity slip boundary conditions are modified as
\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
\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 (Reference Struchtrup2013) 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 (Reference Struchtrup2013) 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. (Reference Yamamoto, Takeuchi and Hyakutake2006) 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. (Reference Spijker, Markvoort, Nedea and Hilbers2010) 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 Reference Zhang and Ma2020; Ma et al. Reference Ma, Zhang, Feng, Xing and Wen2024), 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. (Reference Sone, Bardos, Golse and Sugimoto2000) and Aoki et al. (Reference Aoki, Baranger, Hattori, Kosuge, Martalò, Mathiaud and Mieussens2017), 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 Reference Klinc and Kuščer1972; McCormick Reference McCormick2005; Struchtrup Reference Struchtrup2013; Zhang et al. Reference Zhang, Luan, Deng, Tian and Liang2021; Basdanis et al. Reference Basdanis, Valougeorgis and Sharipov2023), 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 Reference Wu and Struchtrup2017; Wang et al. Reference Wang, Song, Qin and Luo2021; Chen et al. Reference Chen, Gibelli, Li and Borg2023). 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 Reference Myong2004; Myong et al. Reference Myong, Reese, Barber and Emerson2005), some molecules hitting the surface may not immediately scatter back but instead adsorb to it. Recently, Chen et al. (Reference Chen, Gibelli, Li and Borg2023) 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 (Reference Hadjiconstantinou2003) and Hadjiconstantinou & Al-Mohssen (Reference Hadjiconstantinou and Al-Mohssen2005) 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.
Funding
This work was supported by the National Natural Science Foundation of China (grant nos. 12272028 and 92371102).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Some integrations of scattering kernels
For the CLL scattering model, the tangential and normal scattering kernels are
\begin{gather} R( {{u_i},{u_r}} ) = \frac{1}{{\sqrt {2{\rm \pi} RT{\alpha _t}} }}\exp \left( { - \frac{{{{\left( {{u_r} - \sqrt {1 - {\alpha _t}}\,{u_i}} \right)^2}}}}{{2RT{\alpha _t}}}} \right) , \end{gather}
\begin{gather} R( {{w_i},{w_r}} ) = \frac{1}{{\sqrt {2{\rm \pi} RT{\alpha _t}} }}\exp \left( { - \frac{{{{\left( {{w_r} - \sqrt {1 - {\alpha _t}}\,{w_i}} \right)^2}}}}{{2RT{\alpha _t}}}} \right) , \end{gather}
\begin{gather} R( {{v_i},{v_r}} ) = \frac{{{v_r}}}{{RT{\alpha _n}}}\,{{\rm I}_0}\left( {\frac{{\sqrt {1 - {\alpha _n}}\,{v_r}{v_i}}}{{RT{\alpha _n}}}} \right)\exp \left( { - \frac{{{v_r^2} + ( {1 - {\alpha _n}} ){v_i^2}}}{{2RT{\alpha _n}}}} \right) , \end{gather}
where 
${\textrm {I}_0}$ is the modified Bessel function of the first kind (order zero), i.e.
\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}After straightforward calculations, the following integral results can be obtained for the tangential velocity scattering kernel:
\begin{align} &\int_{ - \infty }^\infty {R( {{u_i},{u_r}} )}\,R( {{w_i},{w_r}} )\exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){\rm d}{u_i}\,{\rm d}{w_i} \nonumber\\ &\qquad = \exp \left( { - \frac{{{u_r^2} + {v_i^2} + {w_r^2}}}{{2RT}}} \right),\end{align}
\begin{align} &\int_{ - \infty }^\infty {R( {{u_i},{u_r}} )}\,R( {{w_i},{w_r}} )\exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){u_i}\,{\rm d}{u_i}\,{\rm d}{w_i} \nonumber\\ &\qquad= \exp \left( { - \frac{{{u_r^2} + {v_i^2} + {w_r^2}}}{{2RT}}} \right){u_r}\sqrt {1 - {\alpha _t}} ,\end{align}
\begin{align} &\int_{ - \infty }^\infty {R( {{u_i},{u_r}} )}\,R( {{w_i},{w_r}} )\exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){u_i^2}\,{\rm d}{u_i}\,{\rm d}{w_i} \nonumber\\ &\qquad= \exp \left( { - \frac{{{u_r^2} + {v_i^2} + {w_r^2}}}{{2RT}}} \right)\left[ {{u_r^2}( {1 - {\alpha _t}} ) + RT{\alpha _t}} \right] , \end{align}
\begin{align} &\int_{ - \infty }^\infty {R( {{u_i},{u_r}} )} R( {{w_i},{w_r}} )\exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){u_i^3}\,{\rm d}{u_i}\,{\rm d}{w_i} \nonumber\\ &\qquad= \exp \left( { - \frac{{{u_r^2} + {v_i^2} + {w_r^2}}}{{2RT}}} \right){u_r}\sqrt {1 - {\alpha _t}} \left[ {{u_r^2}( {1 - {\alpha _t}} ) + 3RT{\alpha _t}} \right] . \end{align}For the normal velocity scattering kernel, due to the presence of Bessel functions, only a part of the integrals can be simplified into such a concise form:
\begin{align} &\int_{ 0 }^\infty {R( {{v_i},{v_r}} )} \exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){v_i}\,{\rm d}{v_i} \nonumber\\ &\qquad= \exp \left( { - \frac{{{u_i^2} + {v_r^2} + {w_i^2}}}{{2RT}}} \right){v_r},\end{align}
\begin{align} &\int_{ 0 }^\infty {R( {{v_i},{v_r}} )} \exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){v_i^3}\,{\rm d}{v_i} \nonumber\\ &\qquad= \exp \left( { - \frac{{{u_i^2} + {v_r^2} + {w_i^2}}}{{2RT}}} \right){v_r}\left[ {{v_r^2}( {1 - {\alpha _n}} ) + 2RT{\alpha _n}} \right] . \end{align}When the integrand contains even powers of normal velocity, the integral results still involve Bessel functions, making it difficult to obtain an analytical solution. Therefore, there is a need to approximate the integral results. For ease of calculation and analysis, we assume that the integral results can be approximated as follows:
\begin{align} &\int_{ 0 }^\infty {R( {{v_i},{v_r}} )} \exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){v_i^2}\,{\rm d}{v_i}\nonumber\\ &\qquad = \exp \left( { - \frac{{{u_i^2} + {v_r^2} + {w_i^2}}}{{2RT}}} \right)\left( {{A_0}{v_r^2} + {A_1}{v_r}} \right), \end{align}
\begin{align} &\int_{ 0 }^\infty {R( {{v_i},{v_r}} )} \exp \left( { - \frac{{{u_i^2} + {v_i^2} + {w_i^2}}}{{2RT}}} \right){v_i^4}\,{\rm d}{v_i}\nonumber\\ &\qquad = \exp \left( { - \frac{{{u_i^2} + {v_r^2} + {w_i^2}}}{{2RT}}} \right)\left( {{B_0}{v_r^4} + {B_1}{v_r^3} + {B_2}{v_r^2} + {B_3}{v_r}} \right), \end{align}
where 
${A_0},{A_1},{B_0},{B_1},{B_2},{B_3}$ are constants related to 
${\alpha _n}$.
These coefficients can be determined in the following way. In the case of 
${\alpha _n} = 1$, the Bessel function vanishes, so the integrations (A11) and (A12) have analytical solutions. If we assume that 
${A_0},{B_0},{B_1},{B_2}$ are functions of 
$1 - {\alpha _n}$, then according to the result of 
${\alpha _n} = 1$, the constants 
${A_1}$ and 
${B_3}$ can be determined as 
${\sqrt{({{\rm \pi} RT} }}/{{\sqrt 2 }}) {\alpha _n}$ and 
$3\sqrt {{{\rm \pi} }/{2}} \, {( {RT} )^{3/2}}{\alpha _n}$, respectively. For simplicity, we set 
${A_0}$ as 
$1 - {\alpha _n}$, 
${B_0}$ as 
$1 - {\alpha _n}$, and 
${B_1} = {B_2} = 0$. The comparison results between this approximation and the concise numerical solutions are shown in figure 7, which indicates that our approximation has quite high precision.
Appendix B. The Chapman–Enskog solution of the Boltzmann equation
The equilibrium solution of the Boltzmann equation is 
${f_0} = {f_M}$, i.e. the Maxwell distribution function. The Chapman–Enskog solution results from an expansion of the distribution function for small departures from 
${f_0}$. The original form of the Boltzmann equation is
\begin{equation} \frac{{\partial f}}{{\partial t}} + {c_i}\,\frac{{\partial f}}{{\partial {x_i}}} + {F_i}\,\frac{{\partial f}}{{\partial {c_i}}} = \int_{ - \infty }^\infty {\int_0^{4{\rm \pi} } {\left( {{f^*}\,{f^{1*}} - f\,{f^1}} \right)g\sigma \,{\rm d}\varOmega \,{\rm d}{\boldsymbol{c}_1}} }, \end{equation}
where 
$g$ denotes the relative velocity for the collision of two molecules, and 
$\sigma \,\textrm {d}\varOmega$ denotes the differential cross-section. Just as 
$f$ denotes the value of the distribution function 
$f$ at 
$\boldsymbol {c}$, 
${f^1}$ denotes the value of 
$f$ at 
${\boldsymbol {c}_1}$. Similarly, 
${f^*}$ and 
${f^{1*}}$ are used to denote the values of 
$f$ at 
${\boldsymbol {c}^*}$ and 
$\boldsymbol {c}_1^*$, respectively. And we will consider that only external forces 
${F_i}$ are independent of molecular velocity.
To develop the expansion for the Boltzmann equation, we will first write (B1) in a non-dimensional form. The non-dimensional variables are obtained by dividing 
${c_i}$ by a local reference velocity 
${c_r}$, 
${x_i}$ by a characteristic length 
$L$, 
$t$ by 
$L/{c_r}$, 
${F_i}$ by 
$c_r^2/L$, 
$f$ by 
${n_r}/c_r^3$, and 
$\sigma \,\textrm {d}\varOmega$ by 
${\nu _r}/( {{n_r}{c_r}} )$, where 
${n_r}$ is a reference density, and 
${\nu _r}$ is a reference collision frequency. If we denote the non-dimensional variables by a caret placed over the symbol (e.g. 
$\hat {f}$), we then obtain
\begin{equation} \xi \left[ {\frac{{\partial \hat{f}}}{{\partial \hat{t}}} + {{\hat{c}}_i}\,\frac{{\partial \hat{f}}}{{\partial {{\hat{x}}_i}}} + {{\hat{F}}_i}\,\frac{{\partial \hat{f}}}{{\partial {{\hat{c}}_i}}}} \right] = \int_{ - \infty }^\infty \int_0^{4{\rm \pi} } \left(\,{{{\hat{f}}^*}\,\{\,{\hat{f}}^{1*}\} - \hat{f}\,{{\hat{f}}^1}} \right)\hat{g}\hat{\sigma} \,{\rm d}\varOmega \,{\rm d}{{\hat{\boldsymbol{c}}}_1} , \end{equation}with
\begin{equation} \xi \equiv \frac{{{c_r}}}{{L{\nu _r}}}. \end{equation}
The parameter 
$\xi$ is a measure of the degree of departure from local equilibrium; for small values of 
$\xi$, we expand 
$f$ as the power series
\begin{equation} f = {f_0}( {1 + {\phi _1} + {\phi _2} + \cdots } ), \end{equation}and
\begin{equation} \hat{f} = {\hat{f}_0}\left( {1 + \xi {{\hat{\phi} }_1} + {\xi ^2}{{\hat{\phi} }_2} + \cdots } \right). \end{equation}We substitute the expansion (B5) into the non-dimensional Boltzmann equation (B2), and notice that
\begin{equation} {\hat{f}_0^*}\,\hat{f}_0^{1*} - \hat{f}_{0}\,\hat{f}_0^1 = 0. \end{equation}We then obtain
\begin{equation} \xi \left[ {\frac{{\partial {{\hat{f}_0}}}}{{\partial \hat{t}}} + {{\hat{c}}_i}\,\frac{{\partial {{\hat{f}}_0}}}{{\partial {{\hat{x}}_i}}} + {{\hat{F}}_i}\,\frac{{\partial {{\hat{f}_0}}}}{{\partial {{\hat{c}}_i}}}} \right] = \xi \int_{ - \infty }^\infty \int_0^{4{\rm \pi} } \left( \hat{\phi} _1^* + \hat{\phi}_1^{1*} - \hat{\phi} _1 - \hat{\phi}_1^1 \right)\hat{f}_{0}\,\hat{f}_0^1\,\hat{g}\hat{\sigma} \,{\rm d}\varOmega \,{\rm d}{{\hat{\boldsymbol{c}}}_1} \end{equation}
plus terms of order 
${\xi ^2}$. We now return to our original dimensional variables:
\begin{equation} \frac{{\partial {f_0}}}{{\partial t}} + {c_i}\,\frac{{\partial {f_0}}}{{\partial {x_i}}} + {F_i}\,\frac{{\partial {f_0}}}{{\partial {c_i}}} = \int_{ - \infty }^\infty {\int_0^{4{\rm \pi} } {\left( {{\phi _1^*} + {\phi _1^{1*}} - {\phi _1} - \phi _1^1} \right){f_0}\,f_0^1\,g\sigma \,{\rm d}\varOmega \,{\rm d}{\boldsymbol{c}_1}} }. \end{equation}
After explicit differentiation of 
${f_0}$ and eliminating the time derivatives by means of the conservation equations, the left-hand side of (B8) takes the form
\begin{align} \frac{{\partial
{f_0}}}{{\partial t}} + {c_i}\,\frac{{\partial
{f_0}}}{{\partial {x_i}}} + {F_i}\,\frac{{\partial
{f_0}}}{{\partial {c_i}}} = {f_0}\left[
{{c}_i'\left( {\frac{{\boldsymbol{c}{'^2}}}{{2RT}} -
\frac{5}{2}} \right)\frac{{\partial \ln T}}{{\partial
{x_i}}} + \frac{1}{{RT}}\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{align}
The right-hand side of (B9) is, to within a factor, the integral 
$I[ Q ]$ that is defined as
\begin{equation} I[ Q ] =- \int_{ - \infty }^\infty {\int_0^{4{\rm \pi} } {\left( {{Q^*} + {Q^{1*}} - Q - {Q^1}} \right){f_0}\,f_0^1\,{c_r}\sigma \,{\rm d}\varOmega \,{\rm d}{\boldsymbol{c}_1}} }. \end{equation}With (B9) and (B10), (B8) becomes
\begin{equation} {f_0}\left[ {{c}_i'\left( {\frac{{\boldsymbol{c}'{^2}}}{{2RT}} - \frac{5}{2}} \right)\frac{{\partial \ln T}}{{\partial {x_i}}} + \frac{1}{{RT}}\left( {c_i'c_j' - \frac{1}{3}\,c_i'c_j'{\delta _{ij}}} \right)\frac{{\partial {{\bar{c} }_i}}}{{\partial {x_j}}}} \right] =- I[ {{\phi _1}} ]. \end{equation}After some insightful observation, a partial solution of (B11) can then be written in the form
\begin{equation} {\phi _1} =- \frac{1}{n}\left[ {\sqrt {2RT}\,{A_j}\,\frac{{\partial \ln T}}{{\partial {x_j}}} + {B_{ij}}\,\frac{{\partial \overline{c_i}}}{{\partial {x_j}}} + \psi } \right], \end{equation}with
\begin{gather} {A_j} = A\left( {\frac{{\boldsymbol{c}'}}{{\sqrt {2RT} }},T} \right)\frac{{c_j'}}{{\sqrt {2RT} }}, \end{gather}
\begin{gather} {B_{ij}} = B\left( {\frac{{\boldsymbol{c}'}}{{\sqrt {2RT} }},T} \right) \frac{1}{{2RT}}\left( {c_i'c_j' - \frac{1}{3}\,c_i'c_j'{\delta _{ij}}} \right). \end{gather}
And the two scalar functions 
$A$ and 
$B$ must satisfy the relations
\begin{gather} I\left[ {A \, \frac{{c_j'}}{{\sqrt {2RT} }}} \right] = {f_0}\,\frac{{c_j'}}{{\sqrt {2RT} }} \left( {\frac{{\boldsymbol{c}{'^2}}}{{2RT}} - \frac{5}{2}} \right), \end{gather}
\begin{gather} I\left[ {B \, \frac{1}{{2RT}}\left( {c_i'c_j' - \frac{1}{3}\,c_i'c_j'{\delta _{ij}}} \right)} \right] = 2{f_0} \, \frac{1}{{2RT}}\left( {c_i'c_j' - \frac{1}{3}\,c_i'c_j'{\delta _{ij}}} \right). \end{gather}
As 
$\psi$ must satisfy the equation 
$I[ \psi ] = 0$, 
$\psi$ is a linear combination of the collision invariants and can be absorbed into the scalar function 
$A$.
These results can be combined and the partial solution summarized as follows:
\begin{gather} f \approx {f_0}( {1 + {\phi _1}} ), \end{gather}
\begin{gather} {\phi _1} =- \frac{1}{n}\left[ {\sqrt {2RT}\,{A_j}\,\frac{{\partial \ln T}}{{\partial {x_j}}} + {B_{ij}}\,\frac{{\partial \overline{c_i}}}{{\partial {x_j}}}} \right], \end{gather}
where 
${A_j}$ and 
${B_{ij}}$ take the forms in (B13) and (B14), and the two remaining unknown functions 
$A$ and 
$B$ satisfy (B15) and (B16).
