Fundamental solutions to moment equations for the simulation of microscale gas flows

Abstract

Fundamental solutions (Green’s functions) to Grad’s steady-state linearised 13-moment equations for non-equilibrium gas flows are derived. The creeping microscale gas flows, to which they pertain, are important to understanding the behaviour of atmospheric particulate and the performance of many potential micro/nano technologies. Fundamental solutions are also derived for the regularised form of the steady-state linearised 13-moment equations, due to Struchtrup & Torrilhon (Phys. Fluids, vol. 15 (9), 2003, pp. 2668–2680). The solutions are compared to their classical and ubiquitous counterpart: the Stokeslet. For an illustration of their utility, the fundamental solutions to Grad’s equations are implemented in a linear superposition approach to modelling external flows. Such schemes are mesh free, and benefit from not having to truncate and discretise an infinite three-dimensional domain. The high accuracy of the technique is demonstrated for creeping non-equilibrium gas flow around a sphere, for which an analytical solution exists for comparison. Finally, to demonstrate the method’s geometrical flexibility, the flow generated between adjacent spheres held at a fixed uniform temperature difference is explored.

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 1994; Naris & Valougeorgis 2005; Wu, Reese & Zhang 2014) 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. 2014) and its models (e.g. Naris & Valougeorgis 2005) 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 (1994), 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 1949; Chapman & Cowling 1970). Struchtrup & Torrilhon (2003) 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 2004; Lockerby & Reese 2008; Taheri & Struchtrup 2010). For reviews of recent progress in the development of moment methods the reader is referred to Struchtrup (2005), Struchtrup & Taheri (2011), Torrilhon (2016).

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 2013)) 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 (1953), but first derived by Lorentz (1897)). 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 1964; Karageorghis & Fairweather 1987; Fairweather & Karageorghis 1998; Young et al. 2006). 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 1992; Ying, Biros & Zorin 2006; Veerapaneni et al. 2009). 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 (1949), and their regularised form due to Struchtrup & Torrilhon (2003), 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:

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

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

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

where $\boldsymbol{v}$ is the velocity field, $p$ is the pressure, $\unicode[STIX]{x1D64E}$ is the stress tensor and $\boldsymbol{q}$ is the heat-flux vector; hats denote dimensional quantities. The conservation equations have been linearised around an equilibrium state, described by a constant density $\hat{\unicode[STIX]{x1D70C}}_{e}$ and temperature $\hat{T}_{e}$ and vanishing velocity, stress and heat-flux fields. The ideal gas law linearised about these conditions is

(2.2) $$\begin{eqnarray}\hat{p}=\hat{\unicode[STIX]{x1D70C}}_{e}\hat{\unicode[STIX]{x1D703}}+\hat{\unicode[STIX]{x1D70C}}\hat{\unicode[STIX]{x1D703}}_{e},\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.1b ) and (2.1c ), 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.:

(2.3a-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

(2.4) $$\begin{eqnarray}Kn=\frac{\hat{\unicode[STIX]{x1D706}}_{e}}{\hat{L}}=\sqrt{\frac{\unicode[STIX]{x03C0}}{2\unicode[STIX]{x1D703}_{e}}}\frac{\hat{\unicode[STIX]{x1D707}}_{e}}{\hat{\unicode[STIX]{x1D70C}}_{e}\hat{L}},\end{eqnarray}$$

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 (2003), we define a scaled Knudsen number, adopted purely for notational convenience:

(2.5) $$\begin{eqnarray}Kn^{\prime }=\sqrt{2/\unicode[STIX]{x03C0}}Kn.\end{eqnarray}$$

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:

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

(2.6b ) $$\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}$$

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

where $\unicode[STIX]{x1D6FF}$ is the Dirac delta function and $\boldsymbol{f}\unicode[STIX]{x1D6FF}(\boldsymbol{r})$ represents a point force ( $\hat{\boldsymbol{f}}=\hat{\unicode[STIX]{x1D70C}}_{e}\hat{\unicode[STIX]{x1D703}}_{e}\hat{L}^{2}\boldsymbol{f}$ ) at the origin ( $\boldsymbol{r}=0$ ). If all flow variables tend to their equilibrium values as $\Vert \boldsymbol{r}\Vert \rightarrow \infty$ an analytical solution can be obtained – if the Navier–Stokes closure is adopted, this solution is the Stokeslet (see § 2.1). To arrive at the simplest non-trivial solution, a constant temperature is prescribed, such that (2.2) becomes $p=\unicode[STIX]{x1D70C}$ .

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

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

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

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

where $g\,\unicode[STIX]{x1D6FF}(\boldsymbol{r})$ represents a point heat source at the origin ( ${\hat{g}}=\hat{\unicode[STIX]{x1D70C}}_{e}\hat{\unicode[STIX]{x1D703}}_{e}^{3/2}\hat{L}^{2}g$ ). Again, it is assumed that all flow variables tend to their equilibrium values as $\Vert \boldsymbol{r}\Vert \rightarrow \infty$ . In these cases, to be complementary to the solutions obtained in §§ 2.1–2.3, constant pressure and stationary solutions are sought (which from (2.2) implies $\unicode[STIX]{x1D703}=-\unicode[STIX]{x1D70C}$ ).

The derivation of the Stokeslet in § 2.1 is well documented (see, e.g. Pozrikidis 1992; Lisicki 2013), 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:

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

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

The overline denotes the traceless and symmetrical component of the tensor, i.e. $\overline{\unicode[STIX]{x1D735}\boldsymbol{v}}=(1/2)(\unicode[STIX]{x1D735}\boldsymbol{v}+\unicode[STIX]{x1D735}\boldsymbol{v}^{T})-(1/3)\unicode[STIX]{x1D644}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}$ , where $\unicode[STIX]{x1D644}$ is the identity matrix. The solution we derive in this section is the famous Stokeslet; a name coined by Hancock (1953). This term is commonly used to refer to the velocity field, only; but for generality and clarity, we use the phrase to encompass all variables in the solution (e.g. pressure and stress) and denote these using the superscript ‘ $St$ ’.

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

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D735}p-Kn^{\prime }\,\unicode[STIX]{x0394}\boldsymbol{v}=\boldsymbol{f}\unicode[STIX]{x1D6FF}(\boldsymbol{r}),\end{eqnarray}$$

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.6a ), leads to:

(2.10) $$\begin{eqnarray}\unicode[STIX]{x0394}p=\boldsymbol{f}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}(\boldsymbol{r}).\end{eqnarray}$$

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

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}(\boldsymbol{r})=-\unicode[STIX]{x0394}\left(\frac{1}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert }\right),\end{eqnarray}$$

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):

(2.12a ) $$\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}$$

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

Turning our attention to the velocity field, we substitute (2.11) and (2.12a ) into (2.9):

(2.13) $$\begin{eqnarray}\unicode[STIX]{x0394}\boldsymbol{v}=\frac{1}{Kn^{\prime }}\,\boldsymbol{f}\boldsymbol{\cdot }(\unicode[STIX]{x1D644}\unicode[STIX]{x0394}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D735})\left(\frac{1}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert }\right).\end{eqnarray}$$

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

(2.14) $$\begin{eqnarray}\boldsymbol{v}=\frac{1}{Kn^{\prime }}\,\boldsymbol{f}\boldsymbol{\cdot }(\unicode[STIX]{x1D644}\unicode[STIX]{x0394}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D735})\unicode[STIX]{x1D6FD}(\boldsymbol{r}),\end{eqnarray}$$

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

(2.15) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}=\frac{1}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert }.\end{eqnarray}$$

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:

(2.16) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}=-\unicode[STIX]{x1D6FF}(\boldsymbol{r}).\end{eqnarray}$$

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

(2.17) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}=\frac{\Vert \boldsymbol{r}\Vert }{8\unicode[STIX]{x03C0}}.\end{eqnarray}$$

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

(2.18) $$\begin{eqnarray}\boldsymbol{v}^{\mathit{St}}=\frac{\boldsymbol{f}}{8\unicode[STIX]{x03C0}\,Kn^{\prime }}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert }+\frac{\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{3}}\right).\end{eqnarray}$$

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

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

(2.19b ) $$\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:

(2.20a ) $$\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}$$

(2.20b ) $$\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.6b ) is as derived in § 2.1, i.e.

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

(2.21b ) $$\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.6a ), (2.6c ), (2.20a ) and (2.20b ). The heat flux comes directly from substitution of (2.19b ) and (2.21a ) into (2.20b ):

(2.22a ) $$\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}$$

(2.22b ) $$\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}$$

Note, this is an isothermal heat flux (since $\unicode[STIX]{x1D703}=0$ ): a heat flow generated by the divergence of stress rather than a temperature variation.

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

(2.23) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D735}\boldsymbol{v}^{\mathit{Gr}}}=\overline{\unicode[STIX]{x1D735}\boldsymbol{v}^{\mathit{St}}}-{\textstyle \frac{2}{5}}\overline{\unicode[STIX]{x1D735}\boldsymbol{q}^{\mathit{Gr}}},\end{eqnarray}$$

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

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

(2.24b ) $$\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 2003; Struchtrup 2005). In non-dimensional and linearised form these are

(2.25a ) $$\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}$$

(2.25b ) $$\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}$$

Continuing Hancock’s naming convention, we refer to the isothermal solution of (2.6) and (2.25) as the Reglet, denoting components of the solution using the superscript ‘ $Rg$ ’.

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

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

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

We start by finding the Reglet heat flux ( $\boldsymbol{q}^{\mathit{Rg}}$ ), which for convenience we decompose into the Gradlet ( $\boldsymbol{q}^{\mathit{Gr}}$ ) and an additional component due to regularisation ( $\boldsymbol{q}^{+}$ ):

(2.27) $$\begin{eqnarray}\boldsymbol{q}^{\mathit{Rg}}=\boldsymbol{q}^{\mathit{Gr}}+\boldsymbol{q}^{+}.\end{eqnarray}$$

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

(2.28) $$\begin{eqnarray}\unicode[STIX]{x0394}\boldsymbol{q}^{+}-\frac{5}{9K{n^{\prime }}^{2}}\boldsymbol{q}^{+}=-\unicode[STIX]{x0394}\boldsymbol{q}^{\mathit{Gr}}.\end{eqnarray}$$

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

(2.29) $$\begin{eqnarray}\boldsymbol{q}^{\mathit{Gr}}=-\frac{3}{2}Kn^{\prime }\boldsymbol{f}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}-\unicode[STIX]{x1D644}\unicode[STIX]{x0394})\left(\frac{1}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert }\right),\end{eqnarray}$$

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

(2.30) $$\begin{eqnarray}\boldsymbol{q}^{+}=-{\textstyle \frac{3}{2}}Kn^{\prime }\boldsymbol{f}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}-\unicode[STIX]{x1D644}\unicode[STIX]{x0394})(\unicode[STIX]{x1D6F7}(\boldsymbol{r})),\end{eqnarray}$$

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

(2.31) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}-\frac{5}{9K{n^{\prime }}^{2}}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6FF}(\boldsymbol{r}).\end{eqnarray}$$

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):

(2.32) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=-\frac{\text{e}^{-\unicode[STIX]{x1D6FE}\Vert \boldsymbol{r}\Vert }}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert },\end{eqnarray}$$

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:

(2.33) $$\begin{eqnarray}\boldsymbol{q}^{\mathit{Rg}}=\boldsymbol{q}^{\mathit{Gr}}-\frac{3Kn^{\prime }}{8\unicode[STIX]{x03C0}}\text{e}^{-\unicode[STIX]{x1D6FE}\Vert \boldsymbol{r}\Vert }\boldsymbol{f}\boldsymbol{\cdot }\left[\left(1+\unicode[STIX]{x1D6FE}\Vert \boldsymbol{r}\Vert +\unicode[STIX]{x1D6FE}^{2}\Vert \boldsymbol{r}\Vert ^{2}\right)\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert ^{3}}-\frac{3\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{5}}\right)+\frac{2\unicode[STIX]{x1D6FE}^{2}\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{3}}\right].\end{eqnarray}$$

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

(2.34) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D735}\boldsymbol{v}^{\mathit{Rg}}}=\overline{\unicode[STIX]{x1D735}\boldsymbol{v}^{\mathit{St}}}-{\textstyle \frac{2}{5}}\overline{\unicode[STIX]{x1D735}\boldsymbol{q}^{\mathit{Rg}}}+{\textstyle \frac{16}{15}}Kn^{\prime }\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{\mathit{St}}},\end{eqnarray}$$

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

(2.35) $$\begin{eqnarray}\boldsymbol{v}^{\mathit{Rg}}=\boldsymbol{v}^{\mathit{St}}-{\textstyle \frac{2}{5}}\boldsymbol{q}^{\mathit{Rg}}+{\textstyle \frac{16}{15}}Kn^{\prime }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{\mathit{St}}.\end{eqnarray}$$

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

(2.36) $$\begin{eqnarray}\displaystyle \boldsymbol{v}^{\mathit{Rg}} & = & \displaystyle \boldsymbol{v}^{\mathit{Gr}}-\frac{4Kn^{\prime }}{15\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)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{3Kn^{\prime }}{20\unicode[STIX]{x03C0}}\text{e}^{-\unicode[STIX]{x1D6FE}\Vert \boldsymbol{r}\Vert }\boldsymbol{f}\boldsymbol{\cdot }\left[(1+\unicode[STIX]{x1D6FE}\Vert \boldsymbol{r}\Vert +\unicode[STIX]{x1D6FE}^{2}\Vert \boldsymbol{r}\Vert ^{2})\left(\frac{\unicode[STIX]{x1D644}}{\Vert \boldsymbol{r}\Vert ^{3}}-\frac{3\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{5}}\right)+\frac{2\unicode[STIX]{x1D6FE}^{2}\boldsymbol{r}\boldsymbol{r}}{\Vert \boldsymbol{r}\Vert ^{3}}\right].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

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.

(2.37) $$\begin{eqnarray}\boldsymbol{v}=p=0.\end{eqnarray}$$

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.7b ) and (2.8a ) have the simple solution $\unicode[STIX]{x1D64E}^{\mathit{ThSt}}=0$ , and (2.7c ) and (2.8b ) combine to give

(2.38) $$\begin{eqnarray}\frac{15Kn^{\prime }}{4g}\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=-\unicode[STIX]{x1D6FF}(\boldsymbol{r}).\end{eqnarray}$$

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

(2.39) $$\begin{eqnarray}\unicode[STIX]{x1D703}^{\mathit{ThSt}}=\frac{g}{15Kn^{\prime }\,\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert }.\end{eqnarray}$$

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

(2.40) $$\begin{eqnarray}\boldsymbol{q}^{\mathit{ThSt}}=\frac{g\,\boldsymbol{r}}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert ^{3}}.\end{eqnarray}$$

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.7b ), (2.20b ) and (2.37), which produces Fourier’s law (2.8b ); hence the temperature and heat flux of the Thermal Stokeslet can immediately be identified as solutions:

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

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

What remains to be established is the thermal stress generated (i.e. the stress field that arises, in this stationary flow, due to variations in heat flux). Substitution of (2.37), (2.40) and (2.44b ) into (2.20a ), provides the stress tensor of the Thermal Gradlet:

(2.42a ) $$\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}$$

(2.42b ) $$\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.7b ) and (2.37) ensure that the stress tensor is divergence free, thus the linearised R13 relations (2.25) reduce to

(2.43a ) $$\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}$$

(2.43b ) $$\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}$$

Taking the divergence of (2.43a ) tell us that $\unicode[STIX]{x0394}\boldsymbol{q}=0$ , which reduces (2.43b ) to Fourier’s law (2.8b ). As such, the temperature and heat flux of the Thermal Stokeslet remain solutions:

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

(2.44b ) $$\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}^{+}$ ):

(2.45) $$\begin{eqnarray}\unicode[STIX]{x1D64E}^{\mathit{ThRg}}=\unicode[STIX]{x1D64E}^{\mathit{ThGr}}+\unicode[STIX]{x1D64E}^{+},\end{eqnarray}$$

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

(2.46) $$\begin{eqnarray}\frac{3}{2K{n^{\prime }}^{2}}\unicode[STIX]{x1D64E}^{+}-\unicode[STIX]{x0394}\unicode[STIX]{x1D64E}^{+}=\unicode[STIX]{x0394}\unicode[STIX]{x1D64E}^{\mathit{ThGr}}.\end{eqnarray}$$

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

(2.47) $$\begin{eqnarray}\unicode[STIX]{x1D64E}^{\mathit{ThGr}}=\frac{4g}{5}Kn^{\prime }\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}-\frac{1}{3}\unicode[STIX]{x1D644}\unicode[STIX]{x0394}\right)\left(\frac{1}{4\unicode[STIX]{x03C0}\Vert \boldsymbol{r}\Vert }\right),\end{eqnarray}$$

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

(2.48) $$\begin{eqnarray}\unicode[STIX]{x1D64E}^{+}=\frac{4g}{5}Kn^{\prime }\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}-\frac{1}{3}\unicode[STIX]{x1D644}\unicode[STIX]{x0394}\right)\unicode[STIX]{x1D6F9}(\boldsymbol{r}),\end{eqnarray}$$

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

(2.49) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D6F9}-\frac{3}{2K{n^{\prime }}^{2}}\unicode[STIX]{x1D6F9}=\unicode[STIX]{x1D6FF}(\boldsymbol{r}).\end{eqnarray}$$

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

(2.50) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}=-\frac{2\mathit{Kn}^{\prime 2}}{3}\unicode[STIX]{x1D6FF}(\boldsymbol{r}),\end{eqnarray}$$

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:

(2.51) $$\begin{eqnarray}\unicode[STIX]{x1D64E}^{\mathit{ThRg}}=\unicode[STIX]{x1D64E}^{\mathit{ThGr}}.\end{eqnarray}$$

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 2005).

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.

Figure 1. Streamline plots of the velocity field generated by a point force at the origin: a comparison of the Stokeslet (a, equation (2.18)), the Gradlet (b, equation (2.24)) and the Reglet (c, equation (2.36)). Plots are in a plane containing the point force ( $z=0$ ), with a range $\pm 1.8\unicode[STIX]{x1D706}_{e}$ in the $x$ and $y$ directions. The bold arrow indicates the direction of the point force.

Figure 2. Streamline plots, as in figure 1, but with a farther-field view: range $\pm 4.25\unicode[STIX]{x1D706}_{e}$ in the $x$ and $y$ directions.

In the absence of a temperature gradient, Fourier’s law (2.8b ) does not predict a heat flux. However, isothermal heat flux is a well-known rarefied gas effect (sometimes called a mechanocaloric heat flux; Loyalka 1971), seen for example in low-speed channel flows (Taheri & Struchtrup 2010). 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}$ .

Figure 3. Streamlines of the isothermal heat flux generated by a point force at the origin (note, $\unicode[STIX]{x1D703}=0$ ): a comparison of the Gradlet (a, equation (2.22)) and the Reglet (b, equation (2.33)). Plots are in a plane containing the point force ( $z=0$ ), with a range $\pm 2.7\unicode[STIX]{x1D706}_{e}$ in the $x$ and $y$ directions. The bold arrow indicates the direction of the point force.

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.42b ) 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 1998; Young et al. 2006) 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 1990, pp. 282–283).

Figure 4. Illustration of the method of fundamental solutions (MFS).

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:

(4.1) $$\begin{eqnarray}\boldsymbol{r}_{i}=\boldsymbol{x}-\boldsymbol{x}_{i}^{s},\end{eqnarray}$$

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

(4.2) $$\begin{eqnarray}\boldsymbol{r}_{ij}=\boldsymbol{x}_{j}^{b}-\boldsymbol{x}_{i}^{s}.\end{eqnarray}$$

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 2007; Rana & Struchtrup 2016). 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

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

where

(4.4a,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:

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\mathop{\sum }_{i=1}^{N}\frac{g_{i}}{15\unicode[STIX]{x03C0}Kn^{\prime }\Vert \boldsymbol{r}_{i}\Vert }.\end{eqnarray}$$

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 (2007), for diffusely reflecting boundaries and with nonlinear and regularisation terms omitted. Evaluated at the $j$ th boundary node, they read:

(4.6a ) $$\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}$$

(4.6b ) $$\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}$$

where the subscript ‘ $w$ ’ denotes a wall temperature/velocity boundary condition (distinct from the gas temperature/velocity) and $\boldsymbol{n}_{j}$ is the unit vector normal to the solid surface at the $j$ th boundary node.

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:

(4.7a ) $$\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}$$

(4.7b ) $$\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}$$

(4.7c ) $$\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}$$

(4.7d ) $$\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^®.

Figure 5. Coordinate system for flow around a sphere: $\unicode[STIX]{x1D6FC}$ is the polar angle, $\unicode[STIX]{x1D719}$ is the azimuthal angle and $r$ ( $=\Vert \boldsymbol{x}\Vert$ ) is the radial coordinate. Flow direction is in $z$ , and the sphere is centred at the origin, with its boundary at $r=1$ (non-dimensionalisation is performed with the sphere radius $\hat{a}$ ).

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 2011; Torrilhon 2016) 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:

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

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

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

and at the two poles, where $i$ is the polar node index ( $i=1$ to $N_{\unicode[STIX]{x1D6FC}}$ ) and $j$ is the azimuthal node index ( $j=1$ to $N_{\unicode[STIX]{x1D719}}$ ). In total there are $N=(N_{\unicode[STIX]{x1D6FC}}N_{\unicode[STIX]{x1D719}}+2)$ boundary nodes.

There are an equal number of singularity sites located inside the sphere. These points are defined by (5.1a ) and (5.1b ), 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 (2011). 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 (2011) 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:

(5.2a ) $$\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}$$

(5.2b ) $$\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}$$

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

where

(5.3) $$\begin{eqnarray}B_{p}=-\frac{60\unicode[STIX]{x03C0}+345\unicode[STIX]{x03C0}Kn+450(1+\unicode[STIX]{x03C0})Kn^{2}}{40\unicode[STIX]{x03C0}+270\unicode[STIX]{x03C0}Kn+18(25\unicode[STIX]{x03C0}+18)Kn^{2}+648Kn^{3}},\end{eqnarray}$$

(5.4) $$\begin{eqnarray}B_{t}=-\sqrt{\unicode[STIX]{x03C0}/2}\frac{30Kn^{2}+180(1+1/\unicode[STIX]{x03C0})Kn^{3}}{20\unicode[STIX]{x03C0}+135\unicode[STIX]{x03C0}Kn+9(25\unicode[STIX]{x03C0}+18)Kn^{2}+324Kn^{3}},\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):

(5.5) $$\begin{eqnarray}\hat{F_{d}}=6\unicode[STIX]{x03C0}\hat{\unicode[STIX]{x1D707}}\hat{a}\left[\frac{\left(1+2Kn\right)(1+\frac{15}{4}Kn)+\frac{15}{2\unicode[STIX]{x03C0}}Kn^{2}}{(1+3Kn)\left(1+\frac{15}{4}Kn+\frac{9}{\unicode[STIX]{x03C0}}Kn^{2}\right)-(\frac{3}{4}+9Kn)\frac{6}{5\unicode[STIX]{x03C0}}Kn^{2}}\right].\end{eqnarray}$$

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

Figure 6. Radial variation of: (a) radial velocity at $\unicode[STIX]{x1D6FC}=0$ ; (b) polar velocity at $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x03C0}/2$ ; and (c) temperature at $\unicode[STIX]{x1D6FC}=0$ . Comparison of Young’s analytical solution ( $\boldsymbol{\cdot }\,\boldsymbol{\cdot }\,\boldsymbol{\cdot }$ ) to the Gradlet–MFS solution with $N=66$ boundary nodes (—), for three Knudsen numbers ( $Kn=0$ , $0.2$ , and $1$ ). Note, the surface of the sphere is at $r=1$ and all plots are non-dimensional and normalised with the far-field velocity $U_{\infty }$ .

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 2010) and Boltzmann-equation (Takata, Sone & Aoki 1993) 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):

(5.6) $$\begin{eqnarray}\unicode[STIX]{x1D716}=\left|\frac{F_{d,MFS}-F_{d}}{F_{d}}\right|,\end{eqnarray}$$

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

(5.7) $$\begin{eqnarray}F_{d,MFS}=\mathop{\sum }_{i=1}^{N}\boldsymbol{f}_{i}\boldsymbol{\cdot }\boldsymbol{i}_{z},\end{eqnarray}$$

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.

Figure 7. Numerical error of Gradlet–MFS drag prediction, $\unicode[STIX]{x1D716}$ , against number of boundary nodes, $N$ . Four singularity site distributions are used: sites located on the surface of an origin-centred sphere with radius $a=0.2$ (–*–), $a=0.3$ (– $\circ$ –) and $a=0.4$ (–▵–); and sites located on the surface of a sphere centred at $\boldsymbol{x}=[0.1,0,0]$ with radius $a=0.2$ (– $+$ –).

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 1888), more accurate results (from kinetic theory; Sone & Aoki 1977) and experimental data (Millikan 1923; Allen & Raabe 1982); 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 2010); (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.

Figure 8. Normalised drag force versus Knudsen number. Comparison of Basset’s slip solution (Basset 1888) (— ⋅ —); Young’s solution of Grad’s linearised 13-moment equations, equation (5.5) (——); the Gradlet–MFS with $N=66$ boundary nodes ( $\boldsymbol{\cdot }$ ); the kinetic theory solution due to Sone & Aoki (1977) (–  –); and the experimental data of Millikan (1923) fitted by Allen & Raabe (1982)  (*).

Figure 9. Temperature contours (at $y=0$ ) around adjacent and identical stationary spheres having a fixed and uniform temperature difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ : (a)  $Kn=0$ ; (b)  $Kn=0.1$ ; (c)  $Kn=0.5$ . The spheres are centred at $\boldsymbol{x}=[0,0,\pm 1.25]$ ; length scale non-dimensionalisation and $Kn$ definition uses the sphere radii (i.e. $\hat{L}=\hat{a}$ ). Temperature contour values are normalised with $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}/2$ . The number of boundary nodes used in each case is $N=396$ (198 for each sphere).

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 2002, 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 (2015) 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.

Figure 10. Velocity field (at $y=0$ ) generated between stationary spheres having a fixed and uniform temperature difference: (a)  $Kn=0.1$ ; (b)  $Kn=0.5$ . See figure 9 caption for other details. Glyph lengths are normalised to the maximum velocity magnitude for each figure.

Figure 11. Total non-dimensional force on identical, adjacent and stationary spheres generated by a difference in sphere temperature ( $\unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D703}}$ ), against Knudsen number (based on sphere radii). Gradlet–MFS solutions with $N=132$ (– –), $N=204$ (—) and $N=292$ ( $\boldsymbol{\cdot }\boldsymbol{\cdot }\boldsymbol{\cdot }$ ) total boundary nodes.

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 1970), and their many variants. Future work can also include finding fundamental solutions to R13-equation extensions for hard sphere molecules (Struchtrup & Torrilhon 2013) and diatomic gas mixtures (Gupta, Struchtrup & Torrilhon 2016). Potentially, even, going beyond 13 moments, to the steady linearised regularised 26-moment equations (Gu & Emerson 2009), and other moment-equation families; e.g. the maximum entropy moment equations (McDonald & Groth 2013) and Eu’s equations (Myong 2004).

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 2007; Rana & Struchtrup 2016), 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 1992), 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. 2006; Veerapaneni et al. 2009).