## 1 Introduction

Simulating low-speed dilute gas flow around microscale particulates has a raft of potential applications: in health (e.g. modelling inhalation of fine particulate), in security (e.g. predicting the transport of bacteria-laden aerosolized droplets) and in atmospheric science (e.g. understanding the climatic impact of volcanic ash). The fluid dynamics of individual particles in these cases can often be characterised as having very low Reynolds number (due to both speed and scale) and a small, but non-negligible, Knudsen number. The Knudsen number, $Kn$ , the ratio of the mean free path in the gas to some characteristic length scale (see (2.4), below), describes the extent to which the flow is departed from local thermodynamic equilibrium. Only for Knudsen numbers ${\lesssim}0.01$ are the Navier–Stokes equations accurate (though their applicability can be extended to higher $Kn$ by a modification to the boundary conditions). For $Kn\,\gtrsim \,1.0$ , accurate simulation methods based on gas kinetic theory (e.g. Bird Reference Bird1994; Naris & Valougeorgis Reference Naris and Valougeorgis2005; Wu, Reese & Zhang Reference Wu, Reese and Zhang2014) begin to become computationally tractable. Flows falling in between these approximate limits ( $0.01\lesssim Kn\lesssim 1$ ) are described as being in the transition regime, and present a major simulation challenge: being beyond the physical model of classical fluid dynamics, but inaccessible (computationally speaking) to the accurate techniques developed from kinetic theory. Airborne particles in the size range 0.1– $5~\unicode[STIX]{x03BC}\text{m}$ , at sea-level pressures (where the mean free path is ${\sim}0.07~\unicode[STIX]{x03BC}\text{m}$ ), lie in this transition regime.

While kinetic theory is the gold standard for accuracy in modelling flows across the entire
$Kn$
range, numerical solutions of the Boltzmann equation (e.g. Wu *et al.*
Reference Wu, Reese and Zhang2014) and its models (e.g. Naris & Valougeorgis Reference Naris and Valougeorgis2005) are far more expensive than continuum-fluid calculations (the Boltzmann equation is seven-dimensional, for instance). In the transition regime, and for complex geometries (with no exploitable symmetries) or for modelling a collection of suspended particles, such computations are intractable. Stochastic solution techniques to the Boltzmann equation, such as the direct simulation Monte Carlo method (DSMC) developed by Bird (Reference Bird1994), are also prohibitively expensive for transitional Knudsen numbers at low speeds.

The development of continuum equations from moment approximations to the Boltzmann equation are a potential route to accessing flows in the transition regime. Benefitting from the computational efficiency of continuum methods, while incorporating a physical model valid at higher $Kn$ , moment methods offer a best-of-both-worlds solution. Grad developed the original moment method of approximation to the Boltzmann equation, expanding the distribution function using Hermite polynomials, and deriving the 13-moment equations (Grad Reference Grad1949; Chapman & Cowling Reference Chapman and Cowling1970). Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003) proposed a regularisation of these equations (the R13 equations), which overcome a number of issues, and have proven successful in predicting a raft of canonical flows in the transition regime (Torrilhon & Struchtrup Reference Torrilhon and Struchtrup2004; Lockerby & Reese Reference Lockerby and Reese2008; Taheri & Struchtrup Reference Taheri and Struchtrup2010). For reviews of recent progress in the development of moment methods the reader is referred to Struchtrup (Reference Struchtrup2005), Struchtrup & Taheri (Reference Struchtrup and Taheri2011), Torrilhon (Reference Torrilhon2016).

Tractable numerical methods for three-dimensional external flows do not currently exist for moment equations. The difficulty goes beyond the additional implementation and computational complexity of these higher-order equations in three dimensions. Finite-difference and finite-volume approaches (which have been implemented, for example, in two-dimensional lid-driven cavity flow for the R13 equations (Rana, Torrilhon & Struchtrup Reference Rana, Torrilhon and Struchtrup2013)) are, in general, very difficult to apply to low-speed external flows, because of the complexity of the meshing required (e.g. in suspensions of particles) and the far reaching influence of any internal boundary (e.g. orders of magnitude larger than the submerged body itself).

In Stokes flow analysis (
$Re\rightarrow 0$
,
$Kn=0$
), sometimes referred to as creeping flow, these problems are avoided by a class of technique that exploit the linearity of the Stokes equations, which the Navier–Stokes equations reduce to in such conditions. They are built on a fundamental solution to the Stokes equations, known as the Stokeslet (so-called by Hancock (Reference Hancock1953), but first derived by Lorentz (Reference Lorentz1897)). This fundamental solution (a Green’s function for the Stokes equations) is the flow response to a Dirac delta forcing term applied to the momentum equation. In the most straightforward use, a flow field is represented by a superposition of Stokeslets (positioned, say, inside a particle) that are given a combination of strengths chosen to satisfy the same number of conditions at nodes on the boundary. This approach is known as the method of fundamental solutions (MFS) or the superposition method, amongst other names (Kupradze & Aleksidze Reference Kupradze and Aleksidze1964; Karageorghis & Fairweather Reference Karageorghis and Fairweather1987; Fairweather & Karageorghis Reference Fairweather and Karageorghis1998; Young *et al.*
Reference Young, Jane, Fan, Murugesan and Tsai2006). Another use of the Stokeslet (and its variants, e.g. the Blakelet and Rotlet) is within the related, but more established, boundary integral method (Pozrikidis Reference Pozrikidis1992; Ying, Biros & Zorin Reference Ying, Biros and Zorin2006; Veerapaneni *et al.*
Reference Veerapaneni, Gueyffier, Zorin and Biros2009). Both approaches have the advantage of having a flow domain that is meshless, and a dimensionality that is reduced in order by one (a finite boundary is discretized rather than an infinite volume). The Stokeslet is also in itself of interest, as a fundamental solution to the Stokes equations, and can be used to conveniently derive certain analytical results. What prevents such methods being applied to higher
$Kn$
conditions is, of course, the validity of the Stokes equations from which the Stokeslet is derived.

This leads us to the main purpose of the current work: to derive the first fundamental solutions to moment equations for low-speed non-equilibrium gas flows. Specifically, fundamental solutions to the linearised moment equations originally proposed by Grad (Reference Grad1949), and their regularised form due to Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003), for monatomic dilute gases.

The paper is structured as follows. In § 2, the fundamental solutions are derived. All of the analytical solutions presented have been verified (using symbolic manipulation software) as satisfying the governing equations, in each case. In § 3 these solutions are compared, in particular to the Stokeslet. In § 4, to illustrate the utility of fundamental solutions, an implementation of the MFS is described (using the fundamental solutions to Grad’s linearised 13-moment equations). In § 5 this MFS scheme is applied to low-speed non-equilibrium gas flow around a sphere (for which an analytical solution and experimental data exist), and used to explore a new non-equilibrium gas configuration (the flow generated by stationary adjacent spheres held at a fixed uniform temperature difference) to demonstrate the approach’s geometrical flexibility. In § 6 conclusions are drawn and future opportunities and challenges identified.

## 2 Derivation of fundamental solutions

Our starting point for this creeping flow analysis (i.e. having low speed and negligible inertia) is the linearised steady-state conservation equations for mass, momentum and internal energy:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\hat{\boldsymbol{v}}=0, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D735}}\hat{p}+\hat{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D64E}}=0, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\hat{\boldsymbol{q}}=0, & \displaystyle\end{eqnarray}$$

where
$\hat{\unicode[STIX]{x1D703}}_{e}=\hat{R}\hat{T}_{e}$
, and
$\hat{R}$
is the specific gas constant (
$\hat{\unicode[STIX]{x1D703}}\ll \hat{\unicode[STIX]{x1D703}}_{e}$
and
$\hat{\unicode[STIX]{x1D70C}}\ll \hat{\unicode[STIX]{x1D70C}}_{e}$
). Constitutive closures to (2.1*b*
) and (2.1*c*
), introduced later, relate
$\hat{\boldsymbol{v}}$
,
$\hat{\unicode[STIX]{x1D703}}$
,
$\hat{\unicode[STIX]{x1D64E}}$
and
$\hat{\boldsymbol{q}}$
.

From here on, we non-dimensionalise variables using appropriate combinations of $\hat{\unicode[STIX]{x1D70C}}_{e}$ , $\hat{\unicode[STIX]{x1D703}}_{e}$ and $\hat{L}$ , e.g.:

*a*-

*d*) $$\begin{eqnarray}\boldsymbol{v}=\hat{\boldsymbol{v}}/\sqrt{\hat{\unicode[STIX]{x1D703}}_{e}},\quad \unicode[STIX]{x1D64E}=\hat{\unicode[STIX]{x1D64E}}/(\hat{\unicode[STIX]{x1D70C}}_{e}\hat{\unicode[STIX]{x1D703}}_{e}),\quad \boldsymbol{r}=\hat{\boldsymbol{r}}/\hat{L},\quad \boldsymbol{q}=\hat{\boldsymbol{q}}/[\hat{\unicode[STIX]{x1D70C}}_{e}{\hat{\unicode[STIX]{x1D703}}_{e}}^{3/2}],\end{eqnarray}$$

where $\hat{L}$ is some characteristic length scale and $\hat{\boldsymbol{r}}$ is a three-dimensional position vector. We take this opportunity to define a Knudsen number

where $\hat{\unicode[STIX]{x1D706}}_{e}$ and $\hat{\unicode[STIX]{x1D707}}_{e}$ are the mean free path and dynamic viscosity at the equilibrium state, respectively. As in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003), we define a scaled Knudsen number, adopted purely for notational convenience:

In §§ 2.1–2.6 we seek fundamental solutions to the conservation equations alongside three forms of closure: Navier–Stokes–Fourier (§§ 2.1 and 2.4); Grad’s 13-moment equations (§§ 2.2 and 2.5); and the regularised 13-moment equations (§§ 2.3 and 2.6).

In §§ 2.1–2.3 the fundamental solutions sought correspond to the steady-state flow response of the fluid to a point body force:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}=\boldsymbol{f}\unicode[STIX]{x1D6FF}(\boldsymbol{r}), & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}=0, & \displaystyle\end{eqnarray}$$

In §§ 2.4–2.6 the fundamental solutions sought correspond to the steady-state response of the fluid to a point heat source:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}=0, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}=g\,\unicode[STIX]{x1D6FF}(\boldsymbol{r}), & \displaystyle\end{eqnarray}$$

The derivation of the Stokeslet in § 2.1 is well documented (see, e.g. Pozrikidis Reference Pozrikidis1992; Lisicki Reference Lisicki2013), but which we have included for completeness and to aid the exposition of later derivations.

### 2.1 The Stokeslet

In this section we wish to find an isothermal solution to the conservation equations (2.6) closed by the Navier–Stokes–Fourier constitutive relations:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64E}=-2Kn^{\prime }\,\overline{\unicode[STIX]{x1D735}\boldsymbol{v}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}=-(15/4)Kn^{\prime }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$

We start by determining the pressure field. The constant temperature (
$\unicode[STIX]{x1D703}=0$
) directly provides solution to (2.6*c*
) and (2.8*b*
), i.e.
$\boldsymbol{q}^{\mathit{St}}=0$
. Now, combining the remaining equations, (2.6*a*
), (2.6*b*
) and (2.8*a*
), gives

where
$\unicode[STIX]{x0394}$
is the Laplacian. Taking the divergence of (2.9), and noting the commutative property of the Laplacian and the divergence-free velocity field (2.6*a*
), leads to:

The well-known fundamental solution to Laplace’s equation, provides the useful expression

which when substituted into (2.10) and commuting the Laplacian gives an explicit expression for the pressure field generated by a point forcing (the Stokeslet pressure):

*a*) $$\begin{eqnarray}\displaystyle p^{\mathit{St}} & = & \displaystyle -\boldsymbol{f}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{1}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert }\right),\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{\boldsymbol{f}\boldsymbol{\cdot }\boldsymbol{r}}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert ^{3}}.\end{eqnarray}$$

*a*) into (2.9):

We now choose an appropriate ansatz for $\boldsymbol{v}$ , guided by the right-hand side of (2.13):

where $\unicode[STIX]{x1D6FD}(\boldsymbol{r})$ is to be found. After substituting (2.14) into (2.13), and commuting the Laplacian, we find

Applying the Laplacian to both sides of (2.15), then substituting (2.11), leaves a biharmonic equation for $\unicode[STIX]{x1D6FD}$ with a Dirac delta function in the right-hand side:

We can see immediately, then, that $\unicode[STIX]{x1D6FD}$ is the fundamental solution to the biharmonic equation:

Substitution of (2.17) back into the ansatz (2.14) gives the Stokeslet velocity field:

Finally, the Stokeslet stress tensor is obtained directly from the constitutive relation (2.8*a*
):

*a*) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}^{\mathit{St}} & = & \displaystyle -2Kn^{\prime }\,\overline{\unicode[STIX]{x1D735}\boldsymbol{v}^{\mathit{St}}},\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & = & \displaystyle -\frac{\boldsymbol{f}\boldsymbol{\cdot }\boldsymbol{r}}{4\unicode[STIX]{x03C0}}\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert ^{3}}-\frac{3\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{5}}\right).\end{eqnarray}$$

### 2.2 The Gradlet

We now find a solution to the conservation equations (2.6), but with the linearised form of Grad’s 13-moment equations replacing the Navier–Stokes–Fourier closure:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64E}=-2Kn^{\prime }\overline{\unicode[STIX]{x1D735}\boldsymbol{v}}-{\textstyle \frac{4}{5}}Kn^{\prime }\overline{\unicode[STIX]{x1D735}\boldsymbol{q}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}=-{\textstyle \frac{15}{4}}Kn^{\prime }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-{\textstyle \frac{3}{2}}Kn^{\prime }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\,. & \displaystyle\end{eqnarray}$$

Adopting the naming convention of Hancock, we refer to the isothermal solution of (2.6) and (2.20) as the Gradlet, and use a superscript ‘ $Gr$ ’ to denote the various components of the solution.

Given the conservation equations are unchanged by the introduction of an alternative constitutive relation, a solution to the momentum equation (2.6*b*
) is as derived in § 2.1, i.e.

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64E}^{\mathit{Gr}}=\unicode[STIX]{x1D64E}^{\mathit{St}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle p^{\mathit{Gr}}=p^{\mathit{St}}. & \displaystyle\end{eqnarray}$$

What remains is to find velocity and heat-flux vector fields to satisfy equations (2.6*a*
), (2.6*c*
), (2.20*a*
) and (2.20*b*
). The heat flux comes directly from substitution of (2.19*b*
) and (2.21*a*
) into (2.20*b*
):

*a*) $$\begin{eqnarray}\displaystyle \boldsymbol{q}^{\mathit{Gr}} & = & \displaystyle -{\textstyle \frac{3}{2}}Kn^{\prime }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{\mathit{St}},\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{3Kn^{\prime }}{8\unicode[STIX]{x03C0}}\boldsymbol{f}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert ^{3}}-\frac{3\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{5}}\right).\end{eqnarray}$$

Now we seek the associated velocity field. Substitution of (2.19*a*
) and (2.21*a*
) into (2.20*a*
) produces

which given the velocity and heat-flux fields must be divergence free ((2.6*a*
) and (2.6*c*
)), gives us the Gradlet velocity:

*a*) $$\begin{eqnarray}\displaystyle \boldsymbol{v}^{\mathit{Gr}} & = & \displaystyle \boldsymbol{v}^{\mathit{St}}-{\textstyle \frac{2}{5}}\boldsymbol{q}^{\mathit{Gr}},\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & = & \displaystyle \boldsymbol{v}^{\mathit{St}}-\frac{3Kn^{\prime }}{20\unicode[STIX]{x03C0}}\boldsymbol{f}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert ^{3}}-\frac{3\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{5}}\right).\end{eqnarray}$$

### 2.3 The Reglet

In this section we consider the regularised form of Grad’s 13-moment closure, the R13 equations (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2003; Struchtrup Reference Struchtrup2005). In non-dimensional and linearised form these are

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64E}=-2Kn^{\prime }\overline{\unicode[STIX]{x1D735}\boldsymbol{v}}-{\textstyle \frac{4}{5}}Kn^{\prime }\overline{\unicode[STIX]{x1D735}\boldsymbol{q}}+{\textstyle \frac{2}{3}}K{n^{\prime }}^{2}(\unicode[STIX]{x0394}\unicode[STIX]{x1D64E}+{\textstyle \frac{6}{5}}\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}}), & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}=-{\textstyle \frac{15}{4}}Kn^{\prime }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-{\textstyle \frac{3}{2}}Kn^{\prime }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}+{\textstyle \frac{9}{5}}K{n^{\prime }}^{2}\unicode[STIX]{x0394}\boldsymbol{q}. & \displaystyle\end{eqnarray}$$

As in § 2.2, the constitutive model does not modify the conservation laws, and thus the Stokeslet solution for stress (2.21*a*
) and pressure (2.21*b*
) remain a solution to the momentum equation (2.6*b*
):

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64E}^{\mathit{Rg}}=\unicode[STIX]{x1D64E}^{\mathit{St}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle p^{\mathit{Rg}}=p^{\mathit{St}}. & \displaystyle\end{eqnarray}$$

Substituting (2.21*a*
), (2.22*a*
) and (2.27) into (2.25*b*
) gives:

Combining equations (2.6*b*
), (2.11) (2.12*a*
), (2.21), (2.22*a*
) gives

which we use to guide the following choice of ansatz for $\boldsymbol{q}^{+}$ :

where $\unicode[STIX]{x1D6F7}(\boldsymbol{r})$ is to be found. Substitution of (2.29) and (2.30) into (2.28), and combining with (2.11), gives

A quick inspection of (2.31) tells us that $\unicode[STIX]{x1D6F7}$ is the fundamental solution to the Helmholtz equation (which, again, is well known):

where $\unicode[STIX]{x1D6FE}=\sqrt{5}/(3Kn^{\prime })$ . Substituting this back into the ansatz of (2.30), and combining with (2.27), yields an explicit expression for the heat-flux vector field:

To find the velocity field, we substitute (2.19*a*
) and (2.21*a*
) into (2.25*a*
), to give

where we have made use of the fact that $\unicode[STIX]{x0394}\unicode[STIX]{x1D64E}^{\mathit{St}}=2\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{\mathit{St}}}$ . Because the velocity and heat-flux fields must be divergence free, and $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{\mathit{St}})=0$ , equation (2.34) allows us to write

Finally, to arrive at the velocity field of the Reglet, we combine (2.22), (2.24*a*
) and (2.35):

### 2.4 The Thermal Stokeslet

In this section we find solutions to the conservation equations subject to a point heat source, equations (2.7), rather than a point force. To be complimentary to the Stokeslet, Gradlet and Reglet, derived above, in this and the following two subsections we seek stationary isobaric fundamental solutions, i.e.

For the Navier–Stokes–Fourier constitutive relations (2.8) we refer to the point heat source solution as the Thermal Stokeslet and use the superscript ‘ $ThSt$ ’ to denote it (note, the solution has little to do with the Stokes equations, but this naming convention provides consistency in differentiating the fundamental solutions of §§ 2.1–2.3, corresponding to the steady-state response to a point force, from those in §§ 2.4–2.6, corresponding to steady-state response to a point heat source).

The derivation of the Thermal Stokeslet begins by seeing that equations (2.7*b*
) and (2.8*a*
) have the simple solution
$\unicode[STIX]{x1D64E}^{\mathit{ThSt}}=0$
, and (2.7*c*
) and (2.8*b*
) combine to give

From this it apparent that $\unicode[STIX]{x1D703}$ is the fundamental solution to Laplace’s equation, i.e.

On substitution of (2.39) back into (2.8*b*
) we get the Thermal–Stokeslet heat flux:

### 2.5 The Thermal Gradlet

Now we solve the conservation laws with point heat source, equations (2.7), alongside Grad’s linearised 13-moment closure (2.20), as introduced in § 2.2. We refer to the isobaric stationary fundamental solution we derive as the Thermal Gradlet, and use the superscript ‘ $ThGr$ ’ to denote it.

We start by combining (2.7*b*
), (2.20*b*
) and (2.37), which produces Fourier’s law (2.8*b*
); hence the temperature and heat flux of the Thermal Stokeslet can immediately be identified as solutions:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}^{\mathit{ThGr}}=\unicode[STIX]{x1D703}^{\mathit{ThSt}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}^{\mathit{ThGr}}=\boldsymbol{q}^{\mathit{ThSt}}. & \displaystyle\end{eqnarray}$$

*b*) into (2.20

*a*), provides the stress tensor of the Thermal Gradlet:

*a*) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}^{\mathit{ThGr}} & = & \displaystyle -{\textstyle \frac{4}{5}}Kn^{\prime }\,\overline{\unicode[STIX]{x1D735}\boldsymbol{q}^{\mathit{ThSt}}},\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & = & \displaystyle -\frac{Kn^{\prime }\,g}{5\unicode[STIX]{x03C0}}\,\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert ^{3}}-\frac{3\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{5}}\right).\end{eqnarray}$$

### 2.6 The Thermal Reglet

In this final subsection of derivation, we seek the isobaric and stationary solution to (2.7) alongside the linearised R13 closure. We refer to this solution as the Thermal Reglet and denote it by the superscript ‘ $ThRg$ ’.

Equations (2.7*b*
) and (2.37) ensure that the stress tensor is divergence free, thus the linearised R13 relations (2.25) reduce to

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64E}=-{\textstyle \frac{4}{5}}Kn^{\prime }\overline{\unicode[STIX]{x1D735}\boldsymbol{q}}+{\textstyle \frac{2}{3}}K{n^{\prime }}^{2}\unicode[STIX]{x0394}\unicode[STIX]{x1D64E}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}=-{\textstyle \frac{15}{4}}Kn^{\prime }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+{\textstyle \frac{9}{5}}K{n^{\prime }}^{2}\unicode[STIX]{x0394}\boldsymbol{q}. & \displaystyle\end{eqnarray}$$

*a*) tell us that $\unicode[STIX]{x0394}\boldsymbol{q}=0$ , which reduces (2.43

*b*) to Fourier’s law (2.8

*b*). As such, the temperature and heat flux of the Thermal Stokeslet remain solutions:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}^{\mathit{ThRg}}=\unicode[STIX]{x1D703}^{\mathit{ThSt}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}^{\mathit{ThRg}}=\boldsymbol{q}^{\mathit{ThSt}}. & \displaystyle\end{eqnarray}$$

We decompose the stress tensor of the Thermal Reglet ( $\unicode[STIX]{x1D64E}^{\mathit{ThRg}}$ ) into the Thermal Gradlet ( $\unicode[STIX]{x1D64E}^{\mathit{ThGr}}$ ) and a component due to regularisation ( $\unicode[STIX]{x1D64E}^{+}$ ):

which with (2.42*a*
) and (2.44), substitute into (2.43*a*
) to give

Combining (2.7*c*
), (2.11), (2.42*a*
) and (2.44*b*
), provides an alternate expression for the Thermal Gradlet stress:

which we use to choose an appropriate ansatz for $\unicode[STIX]{x1D64E}^{+}$ :

where $\unicode[STIX]{x1D6F9}(\boldsymbol{r})$ is an unknown field. Substituting (2.47) and (2.48) into (2.46) gives

Because $\unicode[STIX]{x1D64E}^{+}$ has to be divergence free, $\unicode[STIX]{x0394}\unicode[STIX]{x1D6F9}$ must equal zero, which leaves

and tells us that $\unicode[STIX]{x1D64E}^{+}=0$ for $\boldsymbol{r}\neq 0$ . As such, the Thermal–Reglet stress is equal to the Thermal–Gradlet stress:

## 3 Comparison of fundamental solutions

In this section we compare the fundamental solutions derived in § 2 for each of the three constitutive models; in § 3.1 we compare the isothermal, point force fundamental solutions (the Stokeslet, Gradlet and Reglet) and then in § 3.2 the isobaric, stationary, point heat source fundamental solutions (the Thermal Stokeslet, Thermal Gradlet and Thermal Reglet).

### 3.1 Isothermal point force response

For each of the three constitutive models, the pressure and stress fields generated by the point forcing are identical, and correspond to those associated with the Stokeslet. The velocity field, however, is markedly different for the three models (see (2.18), (2.24), and (2.36)); although evidently the Gradlet contains the Stokeslet, and the Reglet contains the Gradlet (and thus the Stokeslet, too). It is also easy to notice, that in the far field (i.e. as $\Vert \boldsymbol{r}\Vert \gg 1$ ) all solutions reduce to the Stokeslet. An equivalent observation being that as $Kn^{\prime }\rightarrow 0$ , the Stokeslet is reclaimed (as is expected). Another noteworthy point is the presence of the decaying exponential factor in the Reglet velocity (2.36), which is Knudsen-layer-like in form, though no boundary exists in the fundamental solution. The factor is identical (in terms of the spatial decay rate) to the exponential factor appearing in Knudsen-layer solutions of the linearised R13 equations (Lockerby, Reese & Gallis Reference Lockerby, Reese and Gallis2005).

Figure 1 compares streamline plots of velocity responses of the Stokeslet (*a*), the Gradlet (*b*) and the Reglet (*c*). (For the plots in this section the characteristic length scale
$L$
has been set to the mean free path
$\unicode[STIX]{x1D706}_{e}$
; thus
$Kn^{\prime }=1$
). What is most striking is that, unlike for the Stokeslet where the
$x$
-component of the flow is positive everywhere (i.e. has the same sign as the point force, indicated by the bold arrow at the origin), the Gradlet and Reglet show regions of recirculation: a three-dimensional vortex-ring structure is generated. Figure 2 provides a farther-field view. The Stokeslet velocity streamlines are unchanged between figures 1(*a*) and 2(*a*), because the Stokeslet has no intrinsic spatial scale. The Gradlet and Reglet, on the other hand, have a spatial scale manifesting from the mean free path. As the field of view becomes much larger than the mean free path, the velocity fields predicted by each constitutive relation tend towards the Stokeslet.

In the absence of a temperature gradient, Fourier’s law (2.8*b*
) does not predict a heat flux. However, isothermal heat flux is a well-known rarefied gas effect (sometimes called a mechanocaloric heat flux; Loyalka Reference Loyalka1971), seen for example in low-speed channel flows (Taheri & Struchtrup Reference Taheri and Struchtrup2010). In figure 3, heat-flux streamlines of the Gradlet (*a*;
$\boldsymbol{q}^{\mathit{Gr}}$
) and the Reglet (*b*;
$\boldsymbol{q}^{\mathit{Rg}}$
) are compared; both show a circulatory heat flux, which, at the point of force, is in opposition to the direction of the velocity fields and the force itself. A vortex-ring analogy can also be applied to describe these heat-flux fields. The Gradlet shows a heat-flux ‘vortex ring’ with no core region, whereas the Reglet predicts a heat-flux vortex ring with a clearly definable diameter; approximately
$4\unicode[STIX]{x1D706}_{e}$
.

### 3.2 Isobaric point heat source response

In this section, we briefly compare the three solutions of the isobaric stationary response to a point source of heat (the Thermal Stokeslet, Thermal Gradlet and Thermal Reglet). All of the solutions describe a temperature that varies with $\Vert \boldsymbol{r}\Vert ^{-1}$ , and a radial heat flux that varies with $\Vert \boldsymbol{r}\Vert ^{-2}$ . The major difference between the solutions is in the prediction of a thermally generated stress field; the Thermal Stokeslet predicts none.

Given the underlying symmetry, it is not surprising that the Thermal Gradlet and Thermal Reglet stress fields (2.42*b*
) and (2.51) can be described by a biaxial stress state: a principal normal stress in the radial direction (
$\unicode[STIX]{x1D70E}_{r}$
) and a principal normal stress in any perpendicular direction (
$\unicode[STIX]{x1D70E}_{\bot }$
). The orientation of the maximum shear-stress planes are found by a
$\pm 45^{\circ }$
rotation of any radial plane (a plane containing the radial axis) around that plane’s perpendicular axis. In the far field, both components of the biaxial stress decay with
$\Vert \boldsymbol{r}\Vert ^{-2}$
; for the Thermal Gradlet and Thermal Reglet. In these solutions, the ratio of the principal stress magnitudes is constant and equal throughout the field: the radial normal stress is compressive (
$\unicode[STIX]{x1D70E}_{r}>0$
), and three times the magnitude of the perpendicular normal stress, which is tensional (
$\unicode[STIX]{x1D70E}_{r}<0$
).

## 4 Superposition of the fundamental solutions

As a demonstration of the usefulness of the fundamental solutions derived in § 2, we implement a simple superposition scheme (known as the MFS, amongst other names; Fairweather & Karageorghis Reference Fairweather and Karageorghis1998; Young *et al.*
Reference Young, Jane, Fan, Murugesan and Tsai2006) for the prediction of creeping external flows over arbitrary geometries (though the general approach is valid for internal flows, too). The approach exploits the fact that, since the equations are linear, a valid flow field can be generated by a superposition of any number and distribution of the fundamental solution. We call the location of one fundamental solution (where its origin lies) a singularity site. Placing these sites within a solid object (as illustrated in figure 4) allows the degrees of freedom of each fundamental solution (e.g. the magnitude and direction of
$\boldsymbol{f}$
) to be determined by satisfying boundary conditions at nodes on the solid’s surface (boundary nodes); achieved by solving a set of linear algebraic equations. In so doing, the dimensions of the given problem are reduced from three to two, since the discretisation is performed over the bounding surface not the flow volume. Furthermore, MFS is particular convenient for external flows, since far-field equilibrium conditions are automatically satisfied by each fundamental solution, and therefore also by any superposition of them. The most notable drawback of the MFS is that singularity sites need to be placed in predetermined locations outside of the flow domain, introducing a degree of arbitrariness to the solution. We remind the reader that the adoption here of MFS is only for the purposes of providing a simple demonstration of the usefulness of the fundamental solutions we have derived; even so, we will show that our results are quite insensitive to the singularity site choice for the simple geometries we consider. The MFS can be thought of as a convenient and instructive, less mathematically rigorous, alternative to the boundary integral method (which also uses fundamental solutions); a description of MFS is given in the undergraduate text on viscous flow by Sherman, in the context of introducing the Stokeslet (Sherman Reference Sherman1990, pp. 282–283).

### 4.1 Implementation of the Gradlet–MFS

We use the subscript $i$ to denote the $i$ th of $N$ singularity sites, and the subscript $j$ the $j$ th of $N$ boundary nodes; $\boldsymbol{x}_{i}^{s}$ is the position of the $i$ th singularity site and $\boldsymbol{x}_{j}^{b}$ the position of the $j$ th boundary node. The vector from the $i$ th singularity site to a position $\boldsymbol{x}$ in the three-dimensional domain is then:

and the vector from the $i$ th singularity site to the $j$ th boundary node:

In this demonstration we choose the Gradlet and Thermal Gradlet fundamental solutions to construct global solutions. The R13 equations are known to have greater accuracy, but for the purposes of a simple and clear demonstration the fundamental solution to the original (unregularised) 13-moment equations has been used. Primarily, this is due to greater complexity and number of R13 boundary conditions (Gu & Emerson Reference Gu and Emerson2007; Rana & Struchtrup Reference Rana and Struchtrup2016). As such, this demonstration is the proof of concept, leaving the implementation with the Reglet, Thermal Reglet and appropriate boundary conditions, the focus of future work, as discussed in § 6.

Using the expression derived earlier for the Gradlet velocity (2.24), the velocity field generated by superposition of $N$ Gradlets (one at each singularity site) is

where

*a*,

*b*) $$\begin{eqnarray}\unicode[STIX]{x1D645}=\frac{1}{8\unicode[STIX]{x03C0}\,Kn^{\prime }}\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert }+\frac{\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{3}}\right)\quad \text{and}\quad \unicode[STIX]{x1D646}=-\frac{3Kn^{\prime }}{20\unicode[STIX]{x03C0}}\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert ^{3}}-\frac{3\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{5}}\right).\end{eqnarray}$$

The Oseen tensor, $\unicode[STIX]{x1D645}$ , can be identified as the Stokeslet component of the Gradlet velocity. Fields generated by the superposition of the Thermal Gradlets (one at each of the same $N$ sites) are obtained similarly, e.g. for temperature:

In summation, the $N$ Gradlets and $N$ Thermal Gradlets, describe a thermal flow field that is defined by 4 $N$ degrees of freedom: the magnitude of $g_{i}$ and three components of $\boldsymbol{f}_{i}$ , for each of the $N$ singularity sites. Once these unknowns are established, by appropriate boundary conditions, the entire flow field is described.

### 4.2 Boundary conditions

We adopt the standard boundary conditions for Grad’s 13-moment equations, taken from Gu & Emerson (Reference Gu and Emerson2007), for diffusely reflecting boundaries and with nonlinear and regularisation terms omitted. Evaluated at the $j$ th boundary node, they read:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}_{j}=\boldsymbol{v}_{\text{w},j}-\sqrt{\frac{\unicode[STIX]{x03C0}}{2}}\boldsymbol{n}_{j}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{j}\boldsymbol{\cdot }(\unicode[STIX]{x1D644}-\boldsymbol{n}_{j}\boldsymbol{n}_{j})-\frac{1}{5}\boldsymbol{q}_{j}\boldsymbol{\cdot }(\unicode[STIX]{x1D644}-\boldsymbol{n}_{j}\boldsymbol{n}_{j}), & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{j}=\unicode[STIX]{x1D703}_{\text{w},j}-\frac{1}{2}\sqrt{\frac{\unicode[STIX]{x03C0}}{2}}\boldsymbol{n}_{j}\boldsymbol{\cdot }\boldsymbol{q}_{j}-\frac{1}{4}\boldsymbol{n}_{j}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{j}\boldsymbol{\cdot }\boldsymbol{n}_{j}, & \displaystyle\end{eqnarray}$$

Using the relevant expressions derived in §§ 2.2 and 2.4 for the Gradlet and Thermal Gradlet, flow variables can be evaluated at the $j$ th boundary node:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}_{j}=\mathop{\sum }_{i=1}^{N}\boldsymbol{f}_{i}\,\boldsymbol{\cdot }(\unicode[STIX]{x1D645}(\boldsymbol{r}_{ij})+\unicode[STIX]{x1D646}(\boldsymbol{r}_{ij})), & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{j}=\mathop{\sum }_{i=1}^{N}\frac{g_{i}}{15\unicode[STIX]{x03C0}Kn^{\prime }\Vert \boldsymbol{r}_{ij}\Vert }, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64E}_{j}=\frac{5}{3Kn^{\prime }}\mathop{\sum }_{i=1}^{N}(\boldsymbol{f}_{i}\boldsymbol{\cdot }\boldsymbol{r}_{ij})\unicode[STIX]{x1D646}(\boldsymbol{r}_{ij})+\frac{4}{3}\mathop{\sum }_{i=1}^{N}g_{i}\unicode[STIX]{x1D646}(\boldsymbol{r}_{ij}), & \displaystyle\end{eqnarray}$$

*d*) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}_{j}=-\frac{5}{2}\mathop{\sum }_{i=1}^{N}\boldsymbol{f}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D646}(\boldsymbol{r}_{ij})+\frac{1}{4\unicode[STIX]{x03C0}}\mathop{\sum }_{i=1}^{N}\frac{g_{i}\boldsymbol{r}_{ij}}{\Vert \boldsymbol{r}_{ij}\Vert ^{3}}. & \displaystyle\end{eqnarray}$$

Substitution of (4.7) into (4.6), and evaluating at each of the
$N$
boundary nodes (i.e.
$j=1$
to
$N$
), produces a system of
$4N$
linear equations with
$4N$
unknowns; this can be solved, for example, by lower–upper (LU) decomposition. In the simple calculations presented in this paper we have adopted the default algorithms implemented in MATLAB^{®}.

For unbounded external flows, as considered in the following sections, the superposed solution automatically tends to an equilibrium stationary state in the far field ( $\boldsymbol{v}=\unicode[STIX]{x1D64E}=\boldsymbol{q}=\unicode[STIX]{x1D703}=0$ as $\Vert \boldsymbol{x}\Vert \rightarrow \infty$ ). This is a convenient feature of MFS for external flows.

## 5 Application to creeping non-equilibrium gas flows around spheres

In the following subsections we consider two configurations of spheres in creeping flow, to demonstrate both the accuracy of the Gradlet–MFS and the geometrical flexibility that the method affords. It is important to stress that the primary purpose of this section is not to investigate or comment on the physical accuracy of moment equations (a topic that is covered extensively, elsewhere: see Struchtrup & Taheri Reference Struchtrup and Taheri2011; Torrilhon Reference Torrilhon2016) but to illustrate the numerical accuracy and convenience of the MFS implementation of higher-order fundamental solutions. In short, it is about verification not validation. Even so, in § 5.1.3 we will present a comparison with experiment and results from kinetic theory, since they are readily available for this case, and it is useful to know how the basic equations (Grad’s) perform; while noting that the regularised (R13) equations are known to perform better.

### 5.1 Creeping flow around a single isothermal sphere

In this section we consider steady-state low-speed flow around a perfect sphere of radius $\hat{a}$ ; the sphere’s temperature is kept uniform and equal to that of the gas in the far field ( $\unicode[STIX]{x1D703}_{w}=0$ ); diffusely reflecting boundaries are assumed. The coordinate system is illustrated in figure 5, where $\unicode[STIX]{x1D6FC}$ is the polar angle, $\unicode[STIX]{x1D719}$ is the azimuthal angle and $r$ is the radial coordinate. We simulate a translating sphere in a fixed reference frame ( $\boldsymbol{v}_{w}=[00-U_{\infty }]$ ), with a stationary flow imposed at the far field (i.e. $\boldsymbol{v}=0$ as $r\rightarrow \infty$ ); however, as is conventional, we present the results in a reference frame moving with the sphere, as illustrated in figure 5 (i.e. $\boldsymbol{v}_{w}=0$ and $\boldsymbol{v}=[0\,0\,U_{\infty }]$ as $r\rightarrow \infty$ ).

For all results in the following sections, non-dimensionalisation has been performed as described above, in (2.3), using the far-field temperature and density, and using the sphere radius ( $\hat{L}=\hat{a}$ ).

The sphere’s boundary nodes are located at:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}=i\unicode[STIX]{x03C0}/(N_{\unicode[STIX]{x1D6FC}}+1), & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=2\unicode[STIX]{x03C0}(j-1)/N_{\unicode[STIX]{x1D719}}, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle r=1, & \displaystyle\end{eqnarray}$$

There are an equal number of singularity sites located inside the sphere. These points are defined by (5.1*a*
) and (5.1*b*
), but with
$r=0.2$
(i.e. they are radially inset from the surface). We will investigate the effect of this choice in § 5.1.3.

For the purposes of verification, we compare our numerical results for the isothermal sphere flow to an analytical solution derived by Young (Reference Young2011). As we have done, Young adopted the linearised form of Grad’s 13-moment equations (2.20) and used boundary conditions equivalent to those presented in (4.6).

#### 5.1.1 Young’s analytical solution

The analytical solution derived by Young (Reference Young2011) is for the general case of a conducting sphere within both a uniform flow and a temperature gradient. In this section we consider the specific solution pertaining to a monatomic gas, purely diffusive surfaces, an infinitely conducting sphere (producing a uniform temperature) and a uniform far-field temperature. These conditions correspond to the following values for the coefficients used in Young’s paper: $K_{tc}=3/4$ , $C_{m}=1$ , $C_{e}=15/8$ , $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . The solution is then:

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle u_{r}=U_{\infty }\left(1+\frac{B_{p}}{r}-\frac{1+B_{p}}{r^{3}}\right)\cos \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle u_{\unicode[STIX]{x1D6FC}}=-U_{\infty }\left(1+\frac{B_{p}}{2r}+\frac{1+B_{p}}{2r^{3}}\right)\sin \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}=\frac{U_{\infty }B_{t}}{r^{2}}\cos \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

and the standard Knudsen number $Kn$ (see (2.4)) is based on the sphere radius. The dimensional drag force $\hat{F_{d}}$ (taken from equation (32b) of Young’s paper) is as follows (two typographical errors have been corrected – a missing parenthesis in the denominator and a missing $Kn$ in the first parenthesis of the numerator):

Note, it was Goldberg (Reference Goldberg1954) who first attempted this solution, but his derivation contained errors.

#### 5.1.2 Velocity and temperature profiles

In figure 6 we compare velocity and temperature profiles from equation (5.2) with predictions from the Gradlet–MFS for a range of $Kn$ ; simulation parameters are given in the caption. For all cases, the agreement between the analytical and numerical solution is accurate to within the resolution of the image; the convergence characteristics of the numerical to the analytical solution are explored in § 5.1.3. Each of the calculations presented were computed in less than a second on a single 2.5 GHz Intel Core i7 processor.

Figure 6(*a*) is a plot of the radial component of velocity (
$u_{r}$
) against
$r$
at
$\unicode[STIX]{x1D6FC}=0$
(refer to figure 5 for the coordinate system). At the downstream pole of the sphere (
$r=1$
,
$\unicode[STIX]{x1D6FC}=0$
), the impermeability condition creates a stagnation point for all
$Kn$
. The extent of the downstream influence is significant for all cases, but, perhaps intuitively, it is reduced with increasing
$Kn$
. However, even at
$Kn=0.2$
, only at 130 radii downstream of the sphere does the velocity reach 1 % of the far-field velocity. This highlights a common difficulty in simulating external creeping flow within finite domains, i.e. the domain has to be very large; a problem that techniques based on fundamental solutions avoid. Figure 6(*b*) shows the polar component of velocity (
$u_{\unicode[STIX]{x1D6FC}}$
) against
$r$
at
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x03C0}/2$
. The non-zero velocity at the sphere’s surface (
$r=1$
) for
$Kn>0$
is velocity slip, which increases with
$Kn$
, as expected.

Figure 6(*c*) shows the radial temperature variation at
$\unicode[STIX]{x1D6FC}=0$
. (Note,
$\unicode[STIX]{x1D703}=0$
at
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x03C0}/2$
, because of the symmetry of the geometry and the linearity of the equations). The generation of a temperature variation around an isothermal translating sphere (held at the same temperature as the far field) is a counter-intuitive non-equilibrium gas effect; a manifestation of the stress and heat-flux coupling described by Grad’s family of equations. This effect has been observed in linearised R13-equation (Torrilhon Reference Torrilhon2010) and Boltzmann-equation (Takata, Sone & Aoki Reference Takata, Sone and Aoki1993) solutions to the same problem.

#### 5.1.3 Drag predictions and numerical convergence

Adopting the MFS here is for convenience and simplicity; fundamental solutions could be implemented in more accurate methods (e.g. the boundary integral method). Nonetheless, we briefly consider the convergence characteristics of the Gradlet–MFS in this section. As a measure of global error ( $\unicode[STIX]{x1D716}$ ), we use the closeness of the numerical prediction of total drag force to the analytical value (5.5):

where $F_{d,MFS}$ is the drag force predicted by the Gradlet–MFS, obtained simply by summing the force strengths of the individual Gradlets:

where $\boldsymbol{i}_{z}$ is a unit vector in the flow direction. Figure 7 shows the dependence of error on the number of boundary nodes, for four different singularity site locations. The result indicates very high accuracy, and rate of convergence, for all cases; even when the singularities sites are translated (in the $x$ -direction) relative to the boundary nodes, breaking the symmetry of the numerical solution. The results also illustrate that higher accuracy is obtained when the singularity sites are further away from the boundary nodes; this is at the expense of poorer conditioning of the equation system to solve.

As explained above, determining the physical accuracy of Grad’s family of equations is not the focus of this paper. Even so, having now established the numerical accuracy of the Gradlet–MFS, it is useful perspective to compare the drag force predictions to first-order slip analysis (Basset’s slip solution; Basset Reference Basset1888), more accurate results (from kinetic theory; Sone & Aoki Reference Sone and Aoki1977) and experimental data (Millikan Reference Millikan1923; Allen & Raabe Reference Allen and Raabe1982); figure 8 shows normalised drag force against $Kn$ . Grad’s linearised 13-moment equations show a major improvement on the Stokes drag prediction ( $F_{S}$ ), and a significant improvement to first-order slip flow analysis (Basset’s solution). However, the solution remains some way from both kinetic theory and experiment. There are two points to make, here: (i) implementation with Reglets (which would introduce a Knudsen-layer contribution) will improve the prediction further (as seen by Torrilhon Reference Torrilhon2010); (ii) an advantage of Grad’s family of equations lies in being able to explore qualitative flow behaviour in configurations that are more complex than an isothermal and translating sphere (particularly arising from heat flux and stress coupling). This second point is illustrated in the following section.

### 5.2 Flow generated between adjacent stationary spheres of fixed temperature

To illustrate the geometrical flexibility that an implementation of the fundamental solution affords, we consider the creeping flow between two identical stationary spheres (with centres fixed 2.5 sphere radii apart), held at different uniform temperatures (having a difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ ). Note, in Stokes flow, no mass flow is generated by the temperature difference, because heat flux and stress are uncoupled.

Figure 9 shows temperature contour plots on a plane running through the centre of both spheres ( $y=0$ ) at various $Kn$ . (For consistency with our previous definitions, the Knudsen number is based on sphere radius $a$ ; arguably, it could be based on the sphere separation distance, $a/2$ , doubling $Kn$ ). The contours show that, as to be expected, the temperature jump at the sphere surfaces (the difference between the gas temperature at the surface and the temperature of the surface itself) increases with $Kn$ , resulting in smaller gradients in the temperature field.

Figure 10 shows the velocity field generated between the two spheres (note, there is no velocity field for
$Kn=0$
). The flow field in figure 10(*a*) might, at first glance, be attributed to a thermal creep (also known as thermal transpiration: the tendency of mass to be transported in the direction of surface gradients of temperature.) However, thermal creep would tend to move the gas in a direction from cold to hot locations (positive in
$z$
); so thermal creep is not the explanation. Instead, this is likely the flow response to a thermally generated stress, known as a thermal stress flow (see Sone Reference Sone2002, for a range of interesting thermal stress flow configurations which run in opposition to thermal creep flows). As the Knudsen number increases further, something even more interesting happens: the flow around the spheres switches direction (at around
$Kn=0.2$
); see figure 10(*b*). This switch occurs at a Knudsen number that is low enough for us to be confident that Grad’s equations are at least qualitatively accurate. The reason for the switch in direction with increasing
$Kn$
can be explained by the weakening of temperature gradients (discussed above) leading to lower gradients in heat flux, and thus lower thermally generated stresses. Temperature gradients along the surface, on the other hand, appear to increase with
$Kn$
, so that thermal creep eventually dominates, and the flow is reversed. However, the surface temperature gradients tend to zero at the poles of each sphere, meaning thermal creep will always be negligible there. This explains the existence of a region near the poles of the spheres that appears to remain driven by thermally generated stress, moving mass from hot to cold. The combined effect of the counter-flowing thermal creep and thermal stress effects is a recirculating flow field. A similar competition of thermal stress and thermal creep is discussed in Mohammadzadeh, Rana & Struchtrup (Reference Mohammadzadeh, Rana and Struchtrup2015) in the context of lid-driven cavity flow.

Figure 11 shows the total drag force on the two spheres against $Kn$ (the force on each sphere is equal, and in the same direction; there is no attractive or repulsive component). The calculations have been performed with different numbers of boundary nodes, to demonstrate the solution’s convergence.

## 6 Conclusions and future work

Fundamental solutions to Grad’s steady-state linearised 13-moment equations, and their regularised equivalent, have been derived. Any linear partial differential operator, with constant coefficients, has a fundamental solution (according to the Malgrange–Ehrenpreis theorem). So, finding fundamental solutions, in the way we have illustrated, offers a new approach to solving higher-order continuum equations for low-speed non-equilibrium gas flows, in general. For example, fundamental solutions can be found to the steady linearised Burnett equations (Chapman & Cowling Reference Chapman and Cowling1970), and their many variants. Future work can also include finding fundamental solutions to R13-equation extensions for hard sphere molecules (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2013) and diatomic gas mixtures (Gupta, Struchtrup & Torrilhon Reference Gupta, Struchtrup and Torrilhon2016). Potentially, even, going beyond 13 moments, to the steady linearised regularised 26-moment equations (Gu & Emerson Reference Gu and Emerson2009), and other moment-equation families; e.g. the maximum entropy moment equations (McDonald & Groth Reference McDonald and Groth2013) and Eu’s equations (Myong Reference Myong2004).

As an illustration of the utility of such fundamental solutions, in the context of creeping rarefied flow modelling, we have implemented the Gradlet and Thermal Gradlet into a method of fundamental solutions (a linear superposition scheme). When applied to flow around a sphere, the scheme showed very fast convergence to a published analytical solution with increasing boundary nodes. A typical simulation presented in this paper (involving approximately 70 boundary nodes) is computed in less than a second on a modern laptop (in this case, using a single 2.5 GHz Intel Core i7 processor). Despite the numerical accuracy, a comparison with experimental data illustrates, as is known, that though Grad’s equations are able to capture qualitative phenomena beyond classical models (for instance mechanocaloric heat flux), they are not quantitatively accurate for transitional Knudsen numbers. Since the linearised R13 equations are formally accurate to higher $Kn$ , the natural next step, now the concept is proven, is to implement the Reglet and Thermal Reglet, with appropriate boundary conditions (Gu & Emerson Reference Gu and Emerson2007; Rana & Struchtrup Reference Rana and Struchtrup2016), into a method of fundamental solutions.

In § 5.2 the geometrical flexibility afforded by using a superposition of fundamental solutions has been demonstrated, by exploring the flow between two adjacent spheres (held at different uniform temperatures). In the longer term, implementation with more efficient methods than the MFS, such as the boundary integral method (see Pozrikidis Reference Pozrikidis1992), will allow multiple interacting particles to be simulated, each, potentially, with moving boundaries (e.g. deforming droplets); these are simulations that are becoming increasingly tractable with the boundary integral method (Ying *et al.*
Reference Ying, Biros and Zorin2006; Veerapaneni *et al.*
Reference Veerapaneni, Gueyffier, Zorin and Biros2009).

## Acknowledgements

The authors would like to thank J. Sprittles for useful discussions, and R. Claydon and A. Shrestha for corrections made to a draft of the manuscript. This work is financially supported in the UK by EPSRC Programme Grants EP/I011927/1, EP/N016602/1, and EPSRC grant EP/K038664/1.