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}$$

$\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}$

$\unicode[STIX]{x1D644}$

$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}$$

(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}$$

$\unicode[STIX]{x1D703}=0$

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 Reference Struchtrup and Torrilhon2003; Struchtrup Reference Struchtrup2005). 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}$$

$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}$$

$\boldsymbol{q}^{\mathit{Rg}}$

$\boldsymbol{q}^{\mathit{Gr}}$

$\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}$$

(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}$$

$\unicode[STIX]{x0394}\boldsymbol{q}=0$

(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.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.

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 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}$ .

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

(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 (Reference Gu and Emerson2007), 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}$$

$w$

$\boldsymbol{n}_{j}$

$j$

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.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}$$

$i$

$i=1$

$N_{\unicode[STIX]{x1D6FC}}$

$j$

$j=1$

$N_{\unicode[STIX]{x1D719}}$

$N=(N_{\unicode[STIX]{x1D6FC}}N_{\unicode[STIX]{x1D719}}+2)$

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

(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 (Reference Goldberg1954) 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.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 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.

Figure 8. Normalised drag force versus Knudsen number. Comparison of Basset’s slip solution (Basset Reference Basset1888) (— ⋅ —); 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 (Reference Sone and Aoki1977) (–  –); and the experimental data of Millikan (Reference Millikan1923) fitted by Allen & Raabe (Reference Allen and Raabe1982)  (*).

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

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.