Fundamental solutions to the regularised 13-moment equations: efficient computation of three-dimensional kinetic effects

Fundamental solutions (Green’s functions) are derived for the regularised 13-moment system (R13) of rarefied gas dynamics, for small departures from equilibrium; these solutions show the presence of Knudsen layers, associated with exponential decay terms, that do not feature in the solution of lower-order systems (e.g. the Navier–Stokes–Fourier equations). Incorporation of these new fundamental solutions into a numerical framework based on the method of fundamental solutions (MFS) allows for efficient computation of three-dimensional gas microflows at remarkably low computational cost. The R13-MFS approach accurately recovers analytic solutions for low-speed flow around a stationary sphere and heat transfer from a hot sphere (for which a new analytic solution has been derived), capturing non-equilibrium flow phenomena missing from lower-order solutions. To demonstrate the potential of the new approach, the influence of kinetic effects on the hydrodynamic interaction between approaching solid microparticles is calculated. Finally, a programme of future work based on the initial steps taken in this article is outlined.


Introduction
Accurate modelling of the flow of gases in micro-and nanoscale environments is known to be key to understanding and optimising numerous industrial and where the velocity field decay is slow so that the computational mesh must extend far from the object in order to achieve domain-independent results.
In classical Stokes flow problems, there is a huge body of research into 'singularity methods' based on fundamental solutions that represent the response of a flow to a singular point force. The primary singularity is the Stokeslet, and from its derivatives a host of other singularities can be derived (rotlets, stresslets, etc.). These solutions are routinely exploited in a number of ways; for example, they can be used to construct flow profiles over various body shapes (Chwang & Wu 1975) or as the basis of the powerful reduced-dimensionality boundary integral method (Pozrikidis 1992). These methods have been applied to numerous liquid and gas microflows, such as in the locomotion of micro-organisms (Lauga & Powers 2009), but have only recently been mentioned in the context of gas microflows.
The first steps towards harnessing the power of the aforementioned methods for rarefied gas flows were taken by Lockerby & Collyer (2016), who derived fundamental solutions for the G13 system and utilised them in the method of fundamental solutions (MFS) to solve a number of canonical problems. This paper also suggested a starting point for a derivation of the R13 system, but a full solution was not presented. The G13 'Gradlets' can actually be derived as a natural extension of the Stokeslets, with the same number of boundary conditions required at a solid surface (three on velocity/stress and one on temperature/heat flux). In contrast, the R13 system represents a step change in complexity due to the appearance of higher-order derivatives in the equations alongside the requirement of satisfying an additional five boundary conditions at a solid surface (nine in total). These are the challenges that have been overcome in this article, culminating in a new efficient simulation tool for 3D rarefied gas microflows.

Problem formulation
The steady linearised equations consider small deviations from an equilibrium state of densityρ e , temperatureT e , vanishing velocityv, deviatoric stressŜ and heat flux q, with dimensional quantities denoted by a hat. In this setting, introducingθ =RT withR the specific gas constant, the linearised ideal gas law (valid forρ ρ e and θ θ e ) for pressure isp =ρ eθ +ρθ e . Taking characteristic scales for lengthL, velocity θ e , stressρ eθe and heat flux ρ eθ 3/2 e , the conservation laws of mass, momentum and energy are ∇ · v = 0, ∇p + ∇ · S = 0 and ∇ · q = 0. (2.1a−c) The differences between the frameworks lie in the formulation of the constitutive relations for the heat flux and deviatoric stress. In order to highlight the additional terms that G13 and R13 add to the usual NSF closure, these can be written conveniently as (Torrilhon 2016) where the overbar operator, (·), means that the trace-free symmetric component of a tensor should be taken (for example (∇q ij is the Kronecker delta), and a scaled Knudsen number has been introduced as whereλ e andμ e are the mean free path and dynamic viscosity of the equilibrium state respectively. Notably, we have eliminated the higher-order moments R (scalar), R (rank-2 tensor) and m (rank-3 tensor) of the R13 system, which appear on the right-hand side of (2.2a) as −Kn ∇ · m and the right-hand side of (2.2b) as −(1/4)Kn(3∇ · R + ∇R), using the R13 closure relations, R = −12Kn ∇ · q, R = − 24 5 Kn ∇q and m = −2Kn ∇S. (2.4a−c) These terms, which were excluded by Grad, introduce higher derivatives than found in the G13 system and hence require the formulation of additional boundary conditions.
2.1. The R13 boundary conditions Consider an impermeable solid surface with outward unit normal n and tangent vectors t 1 , t 2 , so that subscripts n, t 1 , t 2 correspond to the projection of a particular moment onto these directions. Only the fluxes with odd normal components need to be prescribed at the wall (Torrilhon 2016), so that the NSF/G13 system requires four conditions on the variables v n , S nt i , q n . These are derived from Maxwell's accommodation model with diffuse reflection. For a surface with wall velocity v w and temperature θ w , the wall conditions are given by where β = √ 2/π. Higher-order contributions (from G13) are often retained in the NSF formulation to create velocity slip and temperature-jump conditions. For R13, Maxwell's model gives an additional five boundary conditions associated with the odd normal components R nt 1 , R nt 2 , m nt 1 t 2 , m nnn and m nt 1 t 1 (as m is trace-free, a condition on the odd normal moment m nt 2 t 2 is redundant), given by (2.5i) Notably, the boundary conditions have all been written in a form where the inputs (v w and θ w ) are isolated on the left-hand side, which is convenient for the MFS implementation.
2.2. Summary The R13 system requires the solution of 13 bulk equations, the usual five from conservation laws (2.1) alongside eight constitutive evolution equations (2.2), with nine boundary conditions applied at impermeable solid walls (2.5). Our approach to solving these equations requires us to now derive the fundamental solutions of this system.

R13 fundamental solutions
For the linearised NSF system, the fundamental solutions for the momentum and energy equation are well known, namely the Stokeslet and the Green's function for the Laplacian. Furthermore, higher-order singularities based on the Stokeslet's derivatives (doublet, quadrupole, etc.) or point sources of mass (source, source doublet, etc.) have also been derived. In Lockerby & Collyer (2016), this analysis was extended to derive the first fundamental solutions to the G13 system; however, attempts to derive R13 fundamental solutions were not completed due to the increased complexity of the R13 system.
Fundamental solutions to the NSF/G13 systems are obtained by calculating the response of the conservation equations of momentum and energy (2.1) to point sources of force f δ(r) and heat gδ(r) respectively, resulting in four degrees of freedom (three from f and one from g) that can be used in the MFS to satisfy the required four boundary conditions. For the R13 system, not only will these basic fundamental solutions differ from the NSF/G13 systems, but an additional fundamental solution will be obtained from finding the response of the constitutive relation for the deviatoric stress (2.2a) to a point source Gδ(r), yielding the required five more degrees of freedom from the trace-free symmetric rank-2 tensor G. A deeper understanding of how to structure point sources in higher-order moment equations remains an open problem, but building on the NSF/G13 system in this manner will be shown to be a sensible starting point. The full driven system to be investigated is therefore 3.1. Deriving fundamental solutions The fundamental solutions were derived using two independent methods, each yielding the same result, one utilising Fourier transforms and the other directly manipulating the equations by exploiting the properties of known fundamental solutions (see Lisicki (2013) for a range of methods). It was assumed that all quantities vanish as r → ∞, where r = r is the distance from the origin.
Interestingly, the R13 fundamental solutions can be decomposed into contributions from classical solutions and their derivatives. For the G13/NSF system, one requires the fundamental solution to Laplace's equation, and the biharmonic equation, Notably, the R13 system also requires solutions to the Helmholtz equation, which brings into the system exponential terms that will form the Knudsen layer, with γ the parameter characterising the decay rate. The full set of fundamental solutions will be given in § 3.2, but before doing so, we outline the method based on Fourier transforms in order to calculate the response of the system to a point source of heat only (i.e. f = G = 0 in (3.1)). Due to linearity, the responses to each source can be calculated separately and then superimposed.

Derivation of response to a point heat source
The starting point is to take the divergence of (3.1d) twice, to give where = ∇ k ∇ k is the Laplace operator, which can be rearranged into the form of a Helmholtz equation, with γ 2 2 = 5/6Kn 2 . Defining the Fourier transform and corresponding inverse as where k = k . Noting that the Fourier transform of the Helmholtz fundamental solution (3.4) gives (3.10) and that F[ f (x)] = −k 2f (k), the inverse transform of (3.9) is found to be Taking the divergence of the momentum equation in (3.1a-c) and using (3.11) then immediately gives the required solution for pressure, The other variables are found using the same method; that is, by taking the Fourier transform of the required equations to reduce them to an algebraic equation for the corresponding variable and making appropriate substitutions using the known Fourier transforms of (3.4).

The fundamental solutions
The fundamental solutions are where the coefficients for the Helmholtz solutions are γ 2 1 = (5/9)Kn −2 , γ 2 2 = (5/6)Kn −2 and γ 2 3 = (3/2)Kn −2 . These are the exponential functions that superpose to form the 833 R4-7 Knudsen layer in the R13 system. It is therefore critical that all decay rates are obtained, and, notably, γ 3 can only be obtained as a response to the point source in (3.1d). At the NSF level, this source creates Stokes doublet terms in (3.13a) and (3.13b), i.e. directional derivatives of the Stokeslet; however, these are not required in an NSF-MFS scheme (or G13), so G = 0 can be used in this case. The solutions both extend the R13 solutions derived in Lockerby & Collyer (2016) in number and also generalise the approach by making no assumptions on the flow class in order to obtain the complete solution set. For example, the R13 solutions in Lockerby & Collyer (2016) for isobaric stationary flow miss out the thermally driven pressure and stress fields. Notably, however, all computed solutions in their article were for the G13 system, and these agree with our results.
The fundamental solutions derived above have been obtained using two independent methods and have been verified using Wolfram Mathematica.

Method of fundamental solutions
The MFS (Fairweather & Karageorghis 1998;Young et al. 2006) will be used as a simple-to-implement numerical scheme to demonstrate how the fundamental solutions obtained in § 3 allow us to simulate 3D micro-gas-flows. Future work will develop these ideas towards a boundary integral method; however, the MFS is simpler to implement as fundamental solutions are placed outside the flow domain (at singularity sites), so that inside the fluid their contributions are non-singular and satisfy the equations of motion (by design). The solution is then represented as a superposition of fundamental solutions, and the weights of each solution (the f , g and G) are obtained by insisting that the relevant boundary conditions are satisfied (at boundary nodes).
As our focus is on the development of the fundamental solutions, the simplest variant of the MFS is implemented, in which the N singularity sites are (a) distributed over a surface and (b) equal in number to the boundary nodes (creating a square matrix). Notably, the MFS is particularly convenient for external flows (Lockerby & Collyer 2016), as the fundamental solutions automatically satisfy far-field boundary conditions, so that the fluid domain does not have to be artificially truncated, as in many classical numerical methods.

Implementation
The N singularity sites and boundary nodes are at {x s i } N i=1 and {x b j } N j=1 respectively; the choice of their placement will be discussed in § 5.1. The displacements r i and r ij from the ith singularity site to a position x and jth boundary node x b j respectively are then The 13 flow variables will be represented in an array Y = {p, v 1 , v 2 , v 3 , θ , S 11 , S 12 , S 13 , S 22 , S 23 , q 1 , q 2 , q 3 } T of dimension 13 × 1, while the weights associated with the fundamental solution at each singularity site are given as an array i } T of size 9 × 1. The flow variables are then represented as a superposition of contributions from all singularity sites as where A is a coefficient matrix of size 13 × 9. The entry A mn is obtained from the equations for the flow variables (3.13a)-(3.13e) by extracting the coefficient of the nth weight from the mth equation. For example, A 16 would be the coefficient of weight G 12 in (3.13a), i.e. ∇ 1 ∇ 2 (φ + µ γ 2 ). It should be noted that, up to this point, we have not specified a coordinate system, but in our implementation we use Cartesian coordinates, so that ∇ 1 ∇ 2 = ∂ x ∂ y .
The flow variables at each boundary node are obtained from (4.2), which is inserted into the right-hand side of (2.5) to construct a linear system of equations of size 9N × 9N, i.e.
where B contains the 9N boundary values (i.e. the left-hand side of (2.5)), X contains the 9N singularity weights, and the 9N × 9N matrix C connects the contribution of each singularity weight to each boundary value. The singularity weights are the unknowns, which can be found using standard inversion methods.

R13-MFS results
Given that analytic solutions to the NSF system are rare, it is of little surprise that there are few R13 analytic results that can be used as benchmarks for our new method. However, for two canonical problems, heat transfer from a sphere and creeping flow around a sphere (figure 1a), analytic solutions are available to establish the accuracy of the R13-MFS scheme. In the former case, in § 5.2, a new analytic solution is derived for R13, which exhibits a Knudsen layer not present for NSF/G13. In the latter case, in § 5.3, an analytic solution has been obtained by Torrilhon (2010), although its complexity is such that, as a by-product, our computations will verify its accuracy also. Having established the accuracy of the new approach, the flow and force generated by the hydrodynamic interaction between two approaching spheres is calculated in § 5.4.

MFS parameters
To simulate a single sphere, with relevant length scaleL equal to the radius of the sphere, both boundary nodes and singularity sites were evenly distributed over spheres of dimensionless radius 1 (the surface) and r s , respectively. It was found, as in Ramachandran (2002), Young et al. (2008) and Lockerby & Collyer (2016), that faster convergence to the analytic solutions was obtained when the singularities were furthest from the computational domain r s 1, but that this can lead to ill-conditioning. A compromise of r s = 0.1 was found to work well and was used in all that follows. For the R13-MFS solver, the number of nodes was chosen to be 102, with further increases giving results that were graphically indistinguishable. Notably, at low Kn, R13 requires more nodes than G13 in order to resolve sharp exponentially decaying Knudsen layers.

Heat transfer from a stationary sphere
We consider a stationary sphere (v w = 0) whose temperature (θ w ) is maintained above that of the far field (θ → 0 as r → ∞); see figure 1(a) with θ w = 1 and U ∞ = 0. In a spherical coordinate system (r, α, φ), because of the spherical symmetry, the flow properties are independent of the azimuthal and polar directions. Accordingly, all non-zero components can be expressed in terms of the temperature θ (r), radial heat flux q r (r) and radial-radial deviatoric stress S rr (r). The analytic expressions for these quantities are obtained from (2.1)-(2.2) as where κ * (effective heat conductivity) and k 1 are constants of integration. The macroscopic quantities in (5.1)-(5.2) are superpositions of NSF contributions (first order in Kn), G13 contributions (second order in Kn) and from R13 exponential functions describing Knudsen layers. It should be noted that for the NSF and G13 solutions, terms with k 1 are absent; therefore, only one boundary condition is required (to fix κ * ). For example, if the no temperature-jump condition is specified at the boundary (i.e. θ = θ w ), then κ * = 1, but this is known to only be valid for Kn → 0. In the free-molecular case of Kn → ∞, the effective heat conductivity κ * = 0.

Flow around a stationary sphere
We consider a sphere moving relative to the fluid with dimensionless speed 1 at a temperature equal to the far-field value θ w = 0, and sit in the frame moving with the sphere (see figure 1a with θ w = 0 and U ∞ = 1). The most natural coordinate system is cylindrical polars with the z-axis in the direction of the uniform far-field velocity, α the polar angle and φ the azimuthal angle about which axisymmetry is assumed. For this problem, an analytic solution was derived for R13 by Torrilhon (2010).
As can be seen in figure 2(a), at higher Kn, the disturbance to the uniform flow caused by the sphere is significantly reduced. This is confirmed in figure 3(a), where the polar velocity at the midpoint, α = π/2, is shown to barely deviate from the free stream speed of −1 at higher Kn, in contrast to the behaviour at lower Kn, where no slip, v α (r = 1) = 0, is approached and the bulk flow is significantly disturbed. Notably, our R13-MFS predictions in figure 3 are shown to agree perfectly with the analytic results. (It should be noted that the analytic results for temperature were obtained by a code provided by Torrilhon, updated since 2010.) Comparisons of the drag force on the sphere, in figure 2(b), show that the R13-MFS prediction (i) agrees with the analytic solution and (ii) provides a better agreement with experimental results from Millikan (1923) (fitted by Allen & Raabe 1982), aligning with kinetic theory solutions (Sone & Aoki 1977;Beresnev, Chernyak & Fomyagin 1990), than NSF is able to. The R13 is within 10 % of the experimental data up to Kn = 0.39, compared with the NSF, which fails at Kn = 0.11. As discussed in detail in Torrilhon (2010), in contrast to the NSF, the R13 system also predicts temperature variations, as also observed in kinetic solutions (Takata, Sone & Aoki 1993)

Hydrodynamic interactions between approaching solid spheres
To demonstrate the potential of our new method, we consider how the quasistatic hydrodynamic interaction between two equal-sized spheres each moving with unit speed towards the other, along their lines of centre, depends on the distance between their centres l (which could be used to define a second Knudsen number of interest) and Kn (see figure 4a). Calculation of the induced force on neighbouring spheres is key to the dynamics of particulate flows, and the Knudsen number associated with such flows is often in the transition regime, meaning that traditional approaches, e.g. that of Happel & Brenner (1965), become inaccurate. Figure 4(b) shows that as Kn → 0, R13 and G13 approach a common limit, and that limit, reassuringly, at larger separations (l = 10), agrees with the analytic result in Happel & Brenner (1965) obtained for the NSF. However, even at moderate Kn, G13 dramatically overpredicts the force on a sphere compared with R13. Taking the l = 10 case as an example, at Kn = 0.5, G13 predicts F/Kn = 15.01 while the R13 prediction is F/Kn = 10.16, giving an overprediction of 47.7 %. The NSF result is far worse, being insensitive to Kn, and figure 4(a) indicates that this occurs because at higher Kn the spheres disturb the flow less, so that their interaction forces are reduced.
Although this particular set-up happens to be axisymmetric, the method can, at no additional computational cost, obtain the interaction forces for all possible three-dimensional configurations in order to, for example, provide the microscopic information required for the construction of macroscopic particulate flow models.

Discussion
The derivation of the R13 fundamental solutions represents just the first step in a programme of research whose possibilities are outlined below. The drag force generated on one of the spheres at two different separations. The force predicted by R13-MFS (--) is compared with G13-MFS predictions (· · · · · ·). As Kn → 0, the curves for each separation approach a value corresponding to the Kn-independent NSF solution. For large separations (l = 10 here), this value is close to the analytic solution (horizontal dashed line) derived in Happel & Brenner (1965).
6.1. Exploitation of the R13-MFS Within this article, the full potential for the new computational tool has only been touched upon, as our focus here was on methodology and benchmarking. The R13-MFS has been shown to describe fully three-dimensional flows that include non-equilibrium gas effects in an incredibly efficient manner, giving it a significant advantage for flows in the range Kn ≈ 0.01-0.5 over traditional methods such as DSMC, which can become computationally intensive for small-Mach-number flows. This range includes numerous physical systems characterised by scalesL = 0.1-7 µm at atmospheric pressure, and larger at reduced pressure, e.g. for recent applications in free surface flows (Sprittles 2017). Of initial interest are canonical flow problems, studied intensively for Stokes flow, including solid objects of various shapes in classical flows (shear, extensional, etc.) which will provide insight into possible effects in more complex phenomena (e.g. particulate flows).

Extension to higher orders
The R13 system is the lowest-order set of regularised moments that captures some basic features of the Knudsen layer. For Kramer's problem of shear flow past a plate, it is known that R13 predicts the existence of one exponential decay rate to capture the behaviour of the layer. In reality, the Knudsen layer is formed from the superposition of many such exponential functions, so that approximation of this behaviour with one function can lead to poor accuracy as Kn is increased. This was part of the motivation for developing the R26 system in Gu & Emerson (2009, 2014, where three exponential functions are superimposed to form the Knudsen layer and significantly better accuracy is achieved in comparison with benchmark linearised Boltzmann solutions. These advantages of R26 appear to be sufficient to render an attempt to derive fundamental solutions worthwhile. It should be noted that while the complexity of the R13 fundamental solutions is a step up from those obtained for NSF and G13, once these solutions have been obtained and implemented into a numerical approach, the reduction in computational cost compared with other approaches more than compensates the initial investment. It is likely that the same will be true for the R26 system, where the approach developed here should also apply once a deeper understanding of the structure of the fundamental solutions for higher-order moment equations emerges. Looking further into the future, one may look to automate the process of deriving higher-order fundamental solutions so that one can pick the correct number of moments for a given problem; i.e. use the m-refinement suggested by Torrilhon & Sarna (2017).