Knudsen Layer Behaviour and Momentum Accommodation from Surface Roughness Modelling

This work analyses the formation of the Knudsen layer in micro/nanoscale flows by linking a rough wall collision model to a continuum flow model via asymptotic matching. Expressions for the accommodation coefficients in terms of the surface characteristics are derived, allowing for boundary layer analysis of rarefied flows without the use of prior determined accommodation coefficients. This derived model, through use of the Lennard–Jones parameters for a nanoscale system, allows for a prediction of the the effective Tangential Momentum Accommodation Coefficient (TMAC) in flows against ordered nanoscale surfaces.


Introduction
In micro/nanoscale flows, physical characteristic length scales, L, become comparable to, or smaller than, the molecular mean free path (MFP), λ. In such cases the Knudsen number, K n = λ/L, may no longer be smaller than order unity, indicating that rarefaction effects may be present. In such flows the physical behaviour of the fluid deviates from the continuum case and boundary effects become more pronounced [1]. A notable change in the boundary behaviour in micro/nanoscale flows is the increased effect of the Knudsen layer (KL), a region of thermodynamic non-equilibrium near a solid boundary with width of order λ. In the regime of approximately 0.01 < K n < 0.1 these boundary effects become significant, and the no-slip condition used throughout most continuum fluid mechanics breaks down and instead a slip condition must be used [1]. This regime is known as the slip regime. The surface effects in the slip regime are generally characterised by set of parameters called the accommodation coefficients, which model the transfer of some property from the fluid to the wall [2]. In general an accommodation coefficient of some property, q, may be represented as where the subscripts i and r represent the value of the property evaluated for incident and reflected particles, w represents the value at the wall and the angled brackets are the ensemble average of the property. Generally one or two accommodation coefficients are used to model the boundary behaviour, with the Tangential Momentum Accommodation Coefficient (TMAC) and the Energy Accommodation Coefficient (EAC) being commonly used [2,3]. The TMAC is generally found to be approximately unity for surfaces with microscale roughness [4] and as such, most analyses have investigated the case of α ∼ 1. Despite this, it is notable that atomically smooth surfaces, such as nanotubes or graphene, exhibit very low accommodation values. Predictions of this value have been calculated via simulation to be of the order 0.001 for light gases such as hydrogen and methane [5]. It is important to note that few methods of predicting the TMAC for ordered nanoscale systems currently exist, with most values being generated via molecular dynamics simulations [5,6].
In opposition to foundational works such as by Maxwell [7] and Welander [8], where a strained continuum model with given accommodation coefficients is considered, more recent works have considered the loss of momentum near a boundary by the interaction of rarefied particles with a wall with some geometric roughness. Works by Nicholson and Bhatia [9,10] have investigated the collisional behaviour of hard sphere molecules interacting with a hard wall, modelling a defect-free atomic surface such as graphene or a metallic lattice. These works investigated the boundary in a deterministic sense, where the wall shape is determined prior to analysis, and rebounds can be investigated geometrically. An opposing approach is the stochastic methodology, initially presented by Tsuji et al. [11], as the 'virtual wall model', where an incident molecule is reflected by a virtual wall angle governed by a probability density function (PDF). This work was later extended by Sommerfeld and Huber [12], introducing the shadow effect model (SEM), which was developed to model the accessibility of angled surfaces for a given incident angle of the fluid molecule. The effective wall PDF in this model is written as P S H (θ, γ ) = sin(γ + θ) sin γ P wall (θ )H (θ, γ ) (2) where P wall is the non-adjusted wall angle PDF, determined experimentally and H is a normalisation function for a given wall roughness and incident particle angle. In this work, H , takes into account the 'shadow effect' in which low-angled incident particles cannot reach certain wall angles, and is modelled as a Heaviside step function. An experimental component of this work measured the surface roughness of some common engineering materials. It was found that the non-adjusted PDF for these surfaces could be modelled well by normal distribution. Extensions to this model was continued in several works [13][14][15] which extend this model by considering the effect of multiple collisions and more accurate shadowing/screening functions. These models progressively become more complex and admit less analytic information. These past models have concerned the application to particle-laden multiphase flow, where the colliding particles are of the microscale but do not consider the particle-particle interactions which are expected at the outer limit of the boundary layer. As technology has progressed and nanoscale engineering has become a more mature field, works aimed at modelling of fluid mechanical systems of the nanoscale have been performed. A series of models aimed at modelling the additional effects in nanoscale gassurface interactions have emerged, such as the hard-cube [16] and soft-cube [17] models and their extensions as the washboard model by Tully [18]. These models all account for the atomic motions of both the fluid molecules and wall atoms, while including the effect of long-range potential forces, adsorption and heat transfer to varying degrees. These works have often investigated the common test of molecular beam scattering of (generally non-adsorbing) rarefied gases from either a Pt(100) or Pt(111) surface, with the Ar-Pt system being common [6]. It is noteworthy that few theoretical models for the prediction of the TMAC of these nanoscale systems have been developed. To the authors' knowledge, no simple expression estimating this parameter for arbitrary systems currently exist in the literature.
With this motivation, this work aims to connect the outer continuum flow and the free-molecular gas-surface model via boundary layer matching, connecting the rough-wall statistical methods as in previous works [11,12,19] with a formal statistical mechanics treatment in the outer layer as akin to that of Coron [20] and later Aoki et al. [21]. The paper is organised as follows. Section 2 outlines some mathematical preliminaries and introduces the governing equations and expressions of the boundary operators. In Sect. 3 the equations in the outer layer are derived. In Sect. 4 the Knudsen layer expansion is performed using both the Maxwell and Cerciginani-Lampis boundary models and the matching procedure to the outer layer is shown. In Sect. 5, the roughness sub-layer is introduced and the corrected Knudsen layer equation is derived. In Sect. 6, numerical analysis is performed on the equations and approximate relations are presented. Finally, in Sect. 7, this analysis is applied to an example boundary as in the works of Nicholson and Bhatia [9,10] to model nanoscale wall collisions, presenting general expressions for the prediction of the TMAC in these systems.

Preliminaries
Early work on slip boundary models was pioneered by Maxwell [7], who connected kinetic and continuum theory by deriving an expression for the shear stress at the wall, τ w , which may be expressed as where μ w represents the 'effective' viscosity at the wall, which Maxwell considered to be equal to that of the bulk flow. Maxwell considered that the wall brought some fraction, α, of fluid particles are brought into equilibrium with the wall per collision. In this model this fraction is directly equivalent to the TMAC. The accommodated particles are reflected from the wall with a Gaussian profile defined with the wall temperature, with the remaining fraction being perfectly reflected with no momentum or energy exchange. By equating the momentum flux and shear stress at the wall, Maxwell determined an expression for the velocity slip at the wall as where u sli p represents the non-zero tangential velocity value at the wall and the derivative term corresponds to the velocity gradient as shown in Fig. 1. Maxwell's solution to this problem leads to a linear profile near the wall, which does not account for the non-equilibrium Knudsen layer present near the wall. To account for this non-equilibrium region, the first order slip where U and Z are non-dimensional forms of the tangential velocity and wall normal distance respectively and C(α) is some function of the TMAC and is referred to as the first-order (or Maxwell) slip coefficient. Many theoretical and numerical analyses have been undertaken to determine the form of C over the last few decades. For a detailed comparison of forms of this function the reader is directed to the review paper by Zhang et al. [22]. Generalised models of boundary reflection in rarefied flows may be expressed as an integral transform of the particle distribution, f , with some reflection kernel, R, as in [2] where the particle distribution function, f , is defined over the domain of spatial and velocity coordinates x = (x, y, z) and u = (u, v, w) steady system in 3 physical dimensions. This integral is to be understood as the integral over all incident particle velocities at the boundary defined at x = x wall . The presence of the normal velocity w and w ensures zero mass flux through the wall. As outlined in [2] this kernel must obey the following conditions of: non-negativity for all u and u , and the law of detailed balance where g * has been introduced as the Gaussian distribution as where the thermal velocity of the particles is denoted as v T = √ k B T /m, with k B as Boltzmann's constant, T as the thermodynamic temperature, and m as the mass of each identical fluid particle. This final condition on the kernel enforces that the equilibrium distribution is the Gaussian profile. The Maxwell kernel may be expressed as where v w = √ k B T w /m, T w is the temperature of the wall and C w is a linear operator which reverses the normal component of the velocity as (u, v, w) → (u, v, −w). The Maxwell model is commonly used due to its simplicity, with only one accommodation coefficient required to define the whole kernel. Alongside the Maxwell kernel, the other most common forms of reflection kernel used are the Cercignani-Lampis (CL) model [23] and its derivatives. The CL kernel can be expressed as two independent kernels acting on the normal and tangential velocities separately, as functions on the normalised velocities these kernels read where α n and α t represent the normal and tangential energy accommodation coefficients (NEAC,TEAC) and I 0 is the modified Bessel function of the first kind of order zero. This boundary can model specular reflection for α t = α n = 0 and fully diffuse reflection for α t = α n = 1.

Governing Equations and Non-dimensionalization
The problem investigated is steady-state flow over a flat plate with small surface roughness, often referred to in the literature as Kramer's problem. The 2D plane is located at z = 0, and extends infinitely along the x-axis. The particles are defined in a 6D coordinate space of x * 3 = (x * , y * , z * ) and u * 3 = (u * , v * , w * ). For this work, the x-direction is to be defined as aligned to the bulk velocity, and as such the y * and v * are to be neglected, reducing the system to a quasi-2D geometry. The fluid particles may be defined by a probability distribution, f * , over the two vector coordinates x * = (x * , z * ) and u * = (u * , w * ). Here the fluid will be considered to be incompressible, and will be normalised such that its density is equal to 1, allowing for the definition of the bulk quantities where q is some microscopic property defined over the same spatial and velocity components as f * , and q is the mean continuum property defined only over the physical space. Conditional quantities in terms of the normal velocity w are defined as where f * w and q w (x * , w) represent the distribution function and mean of q conditioned on w respectively. The macroscopic fluid velocity is defined by the moment where u * = (u * , w * ). The evolution of this system is governed by the steady-state Boltzmann equation where J * represents the collision operator. This operator, in general, has a very complicated integral form, so without an approximate form analytical modelling generally becomes intractable [2]. Several collision models exist, with the most popular being the Bhatnagar-Gross-Krook (BGK) model [24] which was also presented independently about the same time by Welander [8]. Other models of this operator have been developed, such as the S-model [25] and the ellipsoidal-statistical [26] operators, which can account for systems with varying Prandtl number. Since this work investigates an isothermal system, the BGK relaxation operator can be used without concern of the non-physical Prandtl number. The BGK operator is written in dimensional form as where τ is the mean relaxation time and f * eq represents the equilibrium distribution, which for Maxwell particles in a 2D system is given by the Gaussian distribution f * eq (x * , u * ) = g * (u * (x * ) − u * ). The relaxation time, τ , which is treated as independent of u * in the BGK model [27], represents the characteristic time over which inter-particle collisions cause the system to develop towards equilibrium. For ideal gases, this parameter can be shown to be proportional to the ratio of the mean free path and the thermal velocity, τ ∝ λ/v T , where the proportionality constant is dependent on the model used for the fluid-fluid interactions [28]. Here the exact constant is not of great importance and it suffices that it is of order unity and remains constant over the domain. Since this relaxation time only depends on thermodynamic properties, a constant value for τ is a valid assumption for incompressible, isothermal flow. For Maxwell molecules the MFP may be expressed in terms of the Chapman-Enskog viscosity, μ, as λ = π 0. 5 To transform the governing equations into non-dimensional form, the following variables are introduced where acts as the scale variable for the following analysis. By considering the introduced variables, Eq. (17) becomes where f is the particle distribution in terms of the dimensionless variables. The dimensionless equilibrium distribution may be written in separable form in terms of two Gaussian distributions with unity variance, as f eq (x, u) = g(u(x) − u)g(w). In this non-dimensionalization process, a parameter, , has been introduced, this parameter characterises the relative rarefaction of a fluid flow (as such, it can be seen to be proportional to the Knudsen number). This will be the base equation for the analysis, which will performed by a series of asymptotic expansions. Asymptotic matching of the linearized Boltzmann equation in Couette flow has been performed before with standard textbooks on Boltzmann's equation (for an example, see [2]) outlining the general method. Despite this Aoki et al. [21] note that rigorous derivations of this behaviour are rare in peer-reviewed journals, with the exception being the work of Coron [20]. These two works focus on the deviation of the slip conditions from these equations for use with the Navier-Stokes equation. Here we are mainly interested in the form of the distribution for use as boundary conditions in a free-molecular roughness layer. Due to this, the full matching procedure will be outlined here.

Outer Expansion
The outer distribution, as in the far-field solution away from the wall, may be approximated by a Chapman-Enskog expansion, as the power series f = f 0 + f 1 + · · · . Introducing this expansion into Eq. (20) and collecting the zeroth and first-order terms leads to the distributions To the first order the moments of this distribution behaves as the Navier-Stokes equations [2]. The zeroth-order distribution is equivalent to the given equilibrium distribution. The outer distribution to the first-order may be written To determine the limiting behaviour of this distribution, a boundary layer expansion in the wall normal direction is performed. We consider the coordinate transform, Z = z/ , and expand this function in terms of . Taking the limit as → 0 determines the inner limit of the outer expansion as where-on the coordinate x will be absent from notation, as it does not play a role in this limit. In this layer, the fluid velocity is governed by the incompressible NS equations. Considering no external body forces and steady flow, the continuity and momentum equation in a dimensionless form reduce to where P = P * /(ρv 2 T ) is the dimensionless pressure, where it can be shown that that the dimensionless kinematic viscosity is proportional to [28]. Collecting the 0 terms in the momentum equation and enforcing zero mass flux through the solid boundary at w = 0 at Z = 0, the equations reduce as The tangential flow profile is linear in this limit, and there is zero normal pressure gradient as expected. The limiting velocity profile and distribution become where u 0 and du 0 d Z are constant values corresponding to the effective slip value and gradient, these values are determined through the matching procedure.

Knudsen Layer Expansion
Consider now the Knudsen layer, a fluid boundary layer with thickness of the order of the collision length. In this layer the fluid behaviour is equally affected by the inter-particle collisions and the particle-wall collisions. This layer is defined by introducing the strained wall distance, z = Z , as in the previous section. Substituting this strained coordinate into Eq. (20) leads to the expression Expanding the Knudsen layer distribution as a polynomial in , with . . , the defining ODE for the leading-order distribution is found as where-on the superscript 0 will be dropped, as only the leading order expression will be investigated. This PDE is linear, so integrating factors may be used to find a solution of form where Z 1 acts as a reference point from which the path of integration is defined. There are two cases to consider to track a single particles path: 1. Incoming particles, with w < 0 and contour of integration from Z 1 = ∞ to Z 2. Reflected particles, with w > 0 and contour of integration from Z 1 = 0 to Z This leads to a distribution with the following general form Using the Maxwell model at the boundary leads to a reflected distribution of the form for w > 0. Using this expression in Eq. (30) a full closed form solution for the distribution in the boundary layer may be found. A similar closed form solution may be found for the CL kernel in the case of perfect normal reflection, α n = 0. In this case the reflected distribution reads for w > 0. As before, this expression may be introduced to Eq. (30) to fully define the distribution in the Knudsen layer. This expression for the distribution differs from the Maxwell case, and alters the argument of the bulk velocity instead of taking a convex sum of the incoming distribution and the wall thermal distribution.

Tangential Velocity Profile
Considering the full distribution with reflected behaviour governed by either Eqs. (31) or (32), taking the conditional moment of u leads to Taking the remaining conditional moment of the above with g(w) leads to a Fredholm integral equation of the second kind governing the tangential bulk velocity as where the integral kernel has been defined as Welander [8] first derived a similar expression for considering a purely diffusive boundary, such that α = 1. This integral equation can be solved numerically to find equivalent continuum slip conditions from a boundary layer expansion. Despite the difference between the two kernels, if only the tangential component of the CL kernel is considered, the expressions for the conditional and bulk tangential velocities Eqs. (33) and (34) are identical for both the CL and Maxwell models. The difficulty of using arbitrary reflection kernels is in large part due to the flux conservation term, |w |/w, present in Eq. (6). If perfect reflection, α n = 0, can be safely assumed for the normal velocity, the analysis can still be performed simply as shown above. It is well known [2,8] that the limiting asymptotic behaviour of Eq. (34) is a linear profile, which will be denoted here as for some constant term B. Introducing the above into Eq. (33) leads to the asymptotic form for the conditional velocity which is linear in both Z and w as

Asymptotic Matching
Matching of the bulk velocity from the Knudsen layer to the free-stream solution is rather trivial. For boundary layer matching, the required equality is x z λ r Knudsen Layer Roughness Layer As Z → ∞, the velocity term in the Gaussian distribution becomes linear. This leads to the analytic solution of where erfc is the complementary error function. In the limit of Z → ∞ an asymptotic expression for erfc is used [29], leading to the limiting form as which can be seen to match asymptotically to Eq. (26) exactly.

Roughness Layer Model
In the following section, a model of surface roughness is developed such that the use of phenomenological accommodation coefficients may be avoided. This model assumes the solid boundary to be defined by some function of x. The deviation of this wall from the x-axis is to be of the order of the wall roughness, which has a roughness height of order λ r λ. To determine the form of the reflected distribution, a sub-layer within the Knudsen layer is developed, as shown in Fig. 2 . This boundary layer is of the order of the characteristic roughness and as such is referred to as the 'roughness layer'. This boundary layer is governed by the roughness parameter of r = λ r /L . To derive the governing transport equation, the wall normal distance is further strained, by introducing the strained variable Z = z/ r . Introducing this variable into Eq. (20), and collecting the leading order r terms leads to the equation In this layer, this distribution is independent of all spatial coordinates at the leading order. This independence is due to the physical scale length now being of lower order than the MFP, thus introducing a small region of free-molecular flow. Considering the free-molecular dynamics of particles in this layer, a reflection kernel for use as the boundary condition in the Knudsen layer may be derived. The function defining the wall height is considered to be Lipschitz continuous with respect to x. This requires the existence of a value M > 0 such that, for any x 1 , x 2 , the wall height, z, must satisfy the restriction |z( To investigate the behaviour in this layer the following polar coordinates are introduced where V is the velocity magnitude and γ is the incident angle, orientated to have positive angle as shown in Fig. 3. Two subspaces are now defined, with incident particles having angle in the space i = (−π, 0) and the reflected particles having angle in the space r = (0, π).

Wall Angle PDF
Consider a wall with angle defined probabilistically, such that for a uniformly selected sample along x, the angle may be defined by a PDF, P (θ ). For an incident particle with angle γ with respect to the wall, such that γ < θ lim , the likelihood of this particle striking a given 'effective' wall angle is dependent on the value of γ , and as such, the effective angle will be defined by a PDF denoted P(θ, γ ). The effective wall angle models the total change in angle of a particle due to further effects, such as multiple collisions, shadowing and screening, such as in [13]. For a section of wall x ∈ [a, b] it is assumed that the incident particles are uniformly distributed over x at some height y = H . As shown in Fig. 3, a small section of wall angled at θ , over a length dx 0 , experiences a flux of particles over a projected length, dx γ , as Neglecting the effects of screening and multiple collisions, this allows for the effective PDF to be approximated as which under the chosen coordinates is equivalent to the shadow effect model derived by Sommerfeld and Huber [12] if P is chosen to be normally distributed. A pivotal part of the following argument is the assumption that the limiting angle θ lim is small, thus reducing the effect of multiple collisions and screening. In the limit of θ lim → 0, it is clear that the distribution of the effective angle distribution approaches that of the physical angle.

PDF of Wall Angle for Given Wall Geometry
Now consider a physical wall defined by its height as a function, z(x). From this function, a PDF governing the corresponding wall angle may be determined. Consider a section of this wall over x ∈ [a, b]. If a point along x is chosen with some distribution, then the distribution of the wall height measured against x may be determined as P Z (z)|dz| = P X (x)|dx|. It is more convenient to use the PDF in terms of the local wall angle. This angle is defined as tan θ = d Z/dx with the limits due to the Lipschitz continuity as θ lim = arctan(M). The PDF of the physical wall angle may be found as A distinction is immediately made between the distribution defining the physical wall, which is independent of the incident molecules and the effective probability of the wall, which takes into account the incoming distribution. As before, the effective distribution is determined as While using a given PDF to define a wall clearly does not admit a unique solution, given a wall defined over some x, a unique PDF may be generated. If two walls can be generated by the same PDF, they are referred to as statistically similar.

Reflection Kernel
In this layer, all particles move along straight trajectories towards and away from the wall. The reflection of a given particle from the wall is given by the reflection from some 'effective angle' θ . In this work a Local Specular Reflection (LSR) model will be considered as in the works of Nicholson and Bhatia [9,10]. In this model the reflection of a particle in the local frame of the wall is perfectly specular, enforcing that no kinetic energy loss occurs per collision, with the Maxwell surface accommodation coefficient emerging as the net result when considering all incident angles and collisions over the surface. The local specular reflection model (LSR) is grounded on the fact that in a hard sphere system the interaction force between an incident atom and the solid surface atom with which it collides is directed along the line joining their centres and therefore has no tangential component. Consequently the LSR model is reasonable. Here we note that models considering energy transfer to the wall could be considered, but for this initial work, the LSR will be investigated to allow for simple expressions to be derived. In the local frame of a section of wall at angle, θ , the collision of a single particle is given by the operator C * (θ ) = cos 2θ sin 2θ sin 2θ − cos 2θ where the above is the operator acting on a single incident particle with velocity vector (u, w) as it collides with the wall. In the introduced polar coordinates the collision operator for a θ x z γ i γ r Fig. 4 Incident particle with angle γ ∈ i being processed by a virtual wall with effective wall angle θ to attain a reflected angle of γ ∈ r single particle acts as the mapping (V , γ ) → (V , 2θ − γ ). For collision against a wall with given angle, θ , the reflection kernel may be expressed where V = (V , γ ). An example collision is shown in Fig. 4. This operator only accounts for a distribution reflecting against a wall with a set angle θ . To model the effect of collisions of all possible wall angles, the weighted convex sum of the single-angle kernels can be introduced as where integration against θ will be taken over the domain θ ∈ [−θ lim , θ lim ] throughout this work. This expression cannot be simplified without first some discussion of the form of P(θ, γ ). This distribution must obey the constraints of non-negativity and normalisation. If the rough wall is statistically isotropic then a symmetry condition on γ and θ may be imposed. These constraints are written Some generic comments can be made of this, so far, mostly undetermined wall angle PDF. The dependence of this function on γ is to model shadowing effects as in the works by [12] and [13]. Thus, for small permissible effective wall angles, and for small incident angles, this dependence would be expected to vanish, as near-normally aligned particles reflecting from a smooth wall would not experience any of these prior effects. Neglecting higher order screening effects the reflection kernel may be expressed Using the restrictions enforced on P from Eq. (51), it can be shown that this kernel satisfies the restrictions outlined in Eqs. (6)- (9). The normal flux conservation term seen in Eq. (6) is a local effect, and since there is no total kinetic energy loss in an individual collision, this term vanishes, leaving the reflection kernel expressed as an integral over the wall angle for where the V term arises due to the change of basis. In the case where the limiting wall angle is small, the shadowing term can be neglected with small error. Using a norm with respect to P as where the reflected distribution is approximated by the mapping An important property of this introduced norm relies on the assumption of the symmetry of P . With this property, this norm maps all odd functions in θ to zero.

Molecular Beam Reflection
A motivating example for this analysis is the reflection of a molecular beam from a surface with some roughness. This is one of the primary methods of experiment that may be undertaken to determine the effective accommodation coefficients of a material with roughness at the micro/nanoscale [30,31]. Experimental data of molecular beam scattering reveals lobular patterns which the Maxwell kernel cannot model, with this concern leading to the development of the CL kernel [23]. For a molecular beam, the incoming distribution may be modelled as a delta function. Consider a stream of particles to be defined with incident velocity magnitude and angle (V 0 , γ 0 ), and by the distribution If the incident angle, γ 0 , is large enough such that shadowing effects are not present, then the effective wall PDF reduces to the leading order approximation. Thus, the reflected distribution becomes for γ ∈ r . This represents a 'smearing' of the reflected distribution around the specular case, and of course can be seen to reduce to the specular case when P is defined as a delta function. This is similar in form to the expression given in [12], which is shown in Eq. (2), but with the generalisation that the PDF term need not necessarily assumed to be normally distributed. This generalisation allows for walls with different statistics to be analysed, which will be explored in later sections.

Kramer's Problem
As the roughness layer is a sub-layer to the Knudsen layer, the incoming distribution towards the wall is determined by the asymptotic behaviour of the incoming particles in the Knudsen layer. If the bulk fluid velocity is much lower than the thermal velocity, implying a low Mach number, then the equilibrium distribution may be approximated with a small perturbation around a Gaussian profile with zero mean value as for w < 0. This may be expressed in polar coordinates as for γ ∈ i . It is clear from Eq. (33) that u w does not have a simple closed form. To keep the analysis tractable, the asymptotic form from Eq. (37) will be used, writing u w (w) = Aw + B for some constants A and B such that We use the relation sin(2γ ) = 2 sin γ cos γ and take the reflected distribution by introducing the above into Eq. (55). Expanding the brackets leads to the expression for γ ∈ r . Using simple trigonometric relations, the terms within the norm may be reduced, noting that the introduced norm maps any odd function of θ to zero sin γ = − cos 4θ cos θ sin 2γ − sin 4θ sin θ cos γ cos 2γ sin γ = − cos 4θ cos θ − 1 2 sin 4θ sin θ sin 2γ where the second step in Eq. (62) can be made when particles with angles nearly parallel to the wall are neglected. Introducing the above back into Eq. (60) and rearranging leads to the expression Regrouping the A and B terms back into the linearized approximation and returning to Cartesian coordinates leads to the final expression where the parameters α and β have been introduced above and read in full as α = 1 − cos 3θ β = 1 − cos 4θ cos θ − 1 2 sin 4θ sin θ cos 3θ (65) where α approximates the TMAC and β represents a new scattering coefficient, these terms model effective linear momentum loss due to fluid particle decorrelation caused by scattering from the rough wall collisions. At the leading order, the reflected distribution, as shown in Eq. (64), has a simple form in terms of the incident particle distribution. This allows for a closed form expression for the reflection kernel to be developed as This boundary condition is similar to the Maxwell kernel, but includes a straining coordinate β, where in the case of β = 0, the Maxwell kernel is recovered. Note, if the incoming distribution is Gaussian in w this simply reduces to the Maxwell kernel, which would be equivalent to matching the roughness sub-layer directly to the continuum layer. Such an approximation would enforce A = 0 which is increasingly accurate as the effect of the Knudsen layer diminishes, which occurs as α, K n → 0.

Knudsen Layer Profile with Rough Wall Boundary Condition
Consider the governing distribution equation as in the Knudsen layer and apply the derived boundary condition. After re-arranging, the reflected distribution may be written for w > 0. Where the equilibrium distribution has been written as f eq (u) = f eq,u (Z , u)g(w). This expression may be introduced into Eq. (30) to determine the total distribution for the new boundary condition. Taking the first conditional moment of the distribution with respect to u leads to the conditional velocity as Integrating this conditional velocity expression with respect to w leads again to a Fredholm integral equation of second kind, similar to Eq. (34) but now with the introduced straining parameter, β, present in the second integral as where the integral kernel has been defined as in Eq. (35). This leads to an equation to solve for the tangential velocity profile in the Knudsen layer only using coefficients defined from the physical statistics of the wall.

Numerical Results and Comparisons
To compare the derived model to the commonly used boundary conditions, numerical solutions are obtained. First the molecular beam case is considered and compared to the most common molecular beam scattering model, being the CL model. As mentioned throughout this report, the assumption that the wall roughness angle is defined by a Gaussian profile is quite common, especially when considering surfaces with microscale roughness [12,13].
Here we consider a wall PDF with a normal profile defined with some small roughness angle variance, θ 2 . Approximating the case as by Sommerfeld and Huber [12], then limiting angle is treated as θ lim = ∞. Under this approximation α and β have analytic solutions [29] Equating these expressions leads to the relation (1 − α) 16 9 (71) Considering the physical wall angle to be defined by a normal distribution, the reflected distribution as in Eq. (56) becomes which is equivalent to the case in [12] without shadowing under a different coordinate system. It is important to note that, for the case of large θ 2 , a truncated and normalised normal distribution must be used to satisfy the restrictions outlined in Eq. (51). These derived accommodation coefficients parameters may be related when the limiting θ values are small, in this case the integrands may be expanded as a Taylor series in θ , leading to the approximate leading order expressions This implies for a boundary with small permissible maximum collision angles, hence a very smooth boundary with low TMAC, that the value for β approaches 4/3 of the α value.

Effective Maxwell Model
Now attention is directed to the numerical solutions of the first-order slip coefficient. For the traditional Maxwell boundary problem, Eq. (34) may be solved numerically to determine the flow profile in the Knudsen layer for a given wall. There are numerous methods to perform this integration, but two will be considered here. First is the method of Neumann iteration, which may be performed since the integral kernel is L 2 integrable [32]. This method involves numerical iteration from an initial approximate, linear solution. This iteration method may be expressed as where U (i) represents the numerical form of the bulk velocity u after i iterations and M + and M − are (2N + 1) × (2N + 1) matrices corresponding to the integral operators introduced in Eq. (34). The introduced matrices approximate the continuous linear operators truncated to some wall distance Z = 2Z max . Since the values of M − corresponding to large values of Z approach zero, this matrix can approximate the integral form if Z max is large without further adjustment. The elements of this matrix are where * (Z i ; Z k ) is introduced as the kernel operator truncated to Z = Z k and normalised such that * (Z i ; Z k ) = 0 when i ≥ k and the singularity at zero is treated by truncating the kernel near Z = as * (0) = * ( Z ) for grid spacing Z . To account for the truncated singularity analytic integration to determine the elements in the matrix M −,sing . This integration is performed using the asymptotic form of the singularity around Z = 0 as where γ e is the Euler-Mascheroni constant. Considering the solution to be piecewise linear, analytic integration over the weakly singular point at Z = 0 can be performed, leading to the definition of the functions leading to elements of M −,sing as Since the M + values do not approach zero as Z becomes large, the operator must be adjusted accordingly. It is simple to see that as Z → ∞ the solution u(Z ) becomes linear, so this operator is modified such that it conserves linear behaviour for large Z . The elements of this matrix are where using the derived functions from before it is found that Alternatively this problem may be treated as an eigenvalue problem to find the eigenvector of the operator M + + (1 − α)M − with an eigenvalue of 1. Both methods were considered for the Maxwell and simplified CL case, but emphasis is placed on the iterative method, as in the next section it will be shown that the eigenvalue method is more difficult to apply to the derived model. Solving the numerical problem for a selection of α values and them finding the projected linear solution, as shown in Fig. 5, the form of C(α) may be obtained. Solving the iterative method from an approximate linear solution for the parameters Z max = 40 and Z = 0.01 with 1000 iterations leads to the values in Table 1. To ensure the validity of this method of solution, the calculated values are compared to the values found by Loyalka [33], which have been used frequently as a traditional benchmark in this field. As the errors between these calculations are within 0.3%, the method is considered valid, with the same fit as found in [33] of the form   [33] This form will be used in this report. These numerics are used only to compare to the later developed roughness boundary model, so this error is considered to be acceptable.

Numerical Knudsen Layer Profiles with Rough Wall Boundary
As in the previous section, Eq. (69) can be solved numerically via iteration to determine the flow profile in the Knudsen layer against wall with some roughness. The iteration is expressed here in a similar form to Eq. (74) as where U 1−β is the discretized form of u((1 − β)Z ). In the numerical solution U 1−β is found by interpolating U (i−1) , this interpolation is reversed after multiplication by M − so it may be added to the vector M + U (i) . Some example flow profiles for varying β are shown in Fig. 6. Solving Eq. (82) and finding the projected linear solution allows for a modified slip coefficient, C(α, β), to be obtained. A form similar to Eq. (5) is sought of the form Here an altered Knudsen number, K n * , is introduced. This term is introduced due to the ambiguity over the scale length this K n value should be defined. It is expected that this ratio should be approximately unity, yet, for confined flows, where the roughness layer does take up a non-negligible proportion of the total flow, then this value might deviate from this  approximation. In the case where K n * = K n, then clearly C(α, 0) = C(α) as expected.
Results for α and β values in the range [0, 0.5] are shown in Table 2. Since this model is expected to only be valid for relatively smooth surfaces the larger coefficients corresponding to higher roughness are omitted. A fit can be performed on the data in Table 2 to find the approximate relation where the values of the constants are determined numerically as A = 0.593 and B = 0.557. It appears that β acts as a higher order effect than α. The combined effect of α and β may be characterised by an 'effective' TMAC denoted as α(α, β) such that we enforce C(α)K n * = C(α, β)K n. This effective Knudsen number would lead to a tailoring for this method, allowing for the consideration of multiple collisions. Here, unless stated we will assume K n * = K n.
If it is assumed that C(α) follows the trend Loyalka derived, presented here in Eq. (81), then we may rearrange this expression to find, At very small α values, this simply reduces to no change, with α = α. For increasing α, the addition of the variable β acts to increase the effective TMAC, modelling the added friction due to the addition of particle-particle collisions in the Knudsen layer. For an example, here we consider a surface with normally distributed roughness, which is frequently investigated (1 − α) 16 9 . Case with with β = 0 added as reference [12,13,19] to model rarefied micro-flows, finding the slip coefficient as an adjustment to the case outlined in [33]. This profile is shown in Fig. 7.

Application to Model Nanoscale Wall: Close Packed Spheres
In this section, the behaviour of the derived model is investigated for use with a model nanoscale wall. To compare to nanoscale experiments and simulations, a given wall geometry is considered to represent a 2D analogue of an atomic lattice, such as an atomic Pt surface or a carbon lattice such as graphene. The model considers overlapping spheres of radius, R, with centre spacing, H , interacting with fluid particles of radius, r , as shown in Fig. 8. The centre spacing is allowed to be less than R such that that the wall is represented by the overlapping potential of each atom, with a cut-off radius which models the wall atoms as hard spheres. The expression of relative wall height over one period of x ∈ [−H , H ] of this locus may be written As in Sect. 5.2, the distribution of the physical wall angle may be expressed as Here a roughness parameter is chosen as s = H /(r + R) as in the works of Nicholson and Bhatia [9,10]. The limiting wall angle may be expressed as θ lim = arcsin(s). This expression readily leads to analytical expressions for the leading order α and β values in terms of this roughness parameter  It is seen that β → 4/3α as s → 0 as expected. The above expressions provide, to the authors' knowledge, the first presentation of a simple analytic approximation of the TMAC in nanoscale systems for arbitrary surface-gas combination. The radii used in the determination of s may be approximated by the Lennard-Jones (LJ) diameters of each species. A common test case for this style of surface interaction is the Ar-Pt system [6,30]. In Fig. 9 the effective value is compared to the theoretical analysis performed by [10] and some values extracted from molecular dynamics simulations of a Lennard-Jones system conducted by Arya et al. [34] and Cao et al. [30]. Since both numerical works consider the long range effects between the fluid and wall atoms, where applicable the case of low wall potential compared to particle energy is used. In [34] this is given as ε/k b T = 0.05. This limit is used as a comparison, as it is the closest case to the non-interacting surface investigated here. This analysis results in comparable values of α to the work of Nicholson and Bhatia [10] for very smooth surfaces where s < 0.5. We note that the comparisons available for use are strictly from numerical simulations and not directly from experiment. To the authors' knowledge, there are no appropriate experimental data available for the determination of the TMAC from the scattering of gas flows from atomically smooth surfaces. It is noted that, for very rough surfaces, this model and that of Nicholson and Bhatia no longer coincide. This is due to the increased possibility of multiple collisions as s approaches unity. Even with this note, it is seen that the differences in these models attains a maximum of the order of 10%, implying that multiple collisions do not contribute a large fraction of momentum accommodation, even at rough surfaces, which agrees with observations in other rough wall models [12,14,15]. The main advantages of this model are that: it may be used to approximate the TMAC for any fluid-surface system with the only data required being the LJ parameters which are well known and documented; that the resulting predictions have simple form, and the lower order expression is fully analytic.
In the right plot of Fig. 9 a comparison is made with an increased Knudsen ratio as outlined in Sect. 6.1. The reasoning behind this consideration is in the interest of modelling the effect of multiple collisions. When the roughness of a surface allows for multiple collisions, the effective mean free path is altered, changing the effective Knudsen number. No formal analysis is conducted here, but it can be shown that considering K n * /K n ≈ 1.2 leads to a good match between this work and the calculations performed by Nicholson and Bhatia [10] for all values of s.
Future work to extend this developed model could include: the effect of curving incident particle paths in the roughness layer due to a long range interatomic potential, the effect of thermal fluctuations and energy accommodation of the wall atoms and the behaviour when extended to a full 3D analysis. Investigations into these effects are planned.

Conclusions
This study has investigated the connection between the structure of a surface and its effect of the tangential momentum accommodation coefficient in rarefied flows of the nanoscale. This work has considered the asymptotic matching between three fluid layers: a roughness layer, the Knudsen layer and a continuum layer, effectively extending the classical works of Coron [20] and Aoki et al. [21] to no longer require the use of any accommodation coefficients a priori. This matching led to the derivation of a novel boundary condition and integral equation defining the Knudsen Layer velocity profile due to disturbances from the roughness layer. This work has provided a method for predicting the TMAC in real atomic surfaces, generalising the ideas presented by Tsuji [11] and Sommerfeld and Huber [12], allowing for scattering kernels to be described in scenarios where the assumption of normally distributed surface roughness is not valid. regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.