2.1 Thermophoresis on a sphere by a uniform temperature gradient in the far field

Consider a gas at rest with uniform, constant pressure $p_{0}^{\ast }$ and temperature $T_{0}^{\ast }$ surrounding a sphere of radius $a^{\ast }$ at the same temperature and motionless with respect to the laboratory frame. Under these conditions, the gas is in equilibrium with vanishing heat flux and deviatoric stress. In energy units, the temperature in this state is given by $\unicode[STIX]{x1D703}_{0}^{\ast }=R^{\ast }T_{0}^{\ast }$ , where $R^{\ast }$ is the gas specific constant. Suppose that this state of equilibrium is disturbed by imposing, far from the sphere, a temperature field, $T_{0}^{\ast }+z^{\ast }(\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast })_{\infty }$ , with $(\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast })_{\infty }$ constant, and where plane $z^{\ast }=0$ passes through the sphere’s centre (figure 1). The viscosity and thermal conductivity coefficient of the gas evaluated at the equilibrium state are denoted by $\unicode[STIX]{x1D707}_{0}^{\ast }$ and $k_{0}^{\ast }$ , respectively. The thermal conductivity coefficient for the sphere’s material is denoted by $k_{s}^{\ast }$ (constant). The gas is assumed to be ideal, so that in the equilibrium state the gas density is $\unicode[STIX]{x1D70C}_{0}^{\ast }=p_{0}^{\ast }/\unicode[STIX]{x1D703}_{0}^{\ast }$ . Assume also that the ratio of the gas molecules’ mean free path to the sphere radius is such that rarefaction effects cannot be ignored. On the basis of this ratio, we introduce a Knudsen number

(2.1) $$\begin{eqnarray}Kn=\frac{\unicode[STIX]{x1D707}_{0}^{\ast }\,\unicode[STIX]{x1D703}_{0}^{\ast \,1/2}}{p_{0}^{\ast }\,a^{\ast }}.\end{eqnarray}$$

Note that, from kinetic theory, a commonly used definition for the gas mean free path in the undisturbed state is $(\unicode[STIX]{x03C0}/2)^{1/2}\unicode[STIX]{x1D707}_{0}^{\ast }\,\unicode[STIX]{x1D703}_{0}^{\ast \,1/2}/p_{0}^{\ast }$ . In definition of (2.1), we have dropped the factor $(\unicode[STIX]{x03C0}/2)^{1/2}$ for simplicity.

In what follows, we consider the governing equations for a monatomic gas composed of Maxwell molecules. In this case, the Prandtl number $Pr=2/3$ . From the well-known definition of $Pr$ , a useful relationship between the gas thermal conductivity and dynamic viscosity can be obtained (Struchtrup Reference Struchtrup2005b )

(2.2) $$\begin{eqnarray}k_{0}^{\ast }=\frac{5}{2}\frac{\unicode[STIX]{x1D707}_{0}^{\ast }R^{\ast }}{Pr}=\frac{15}{4}\unicode[STIX]{x1D707}_{0}^{\ast }R^{\ast }.\end{eqnarray}$$

For a monatomic gas, the ratio of specific heats is $\unicode[STIX]{x1D6FE}=5/3$ (e.g. see Young Reference Young2011).

With the aim of modelling flow and heat transfer phenomena in a rarefied gas, we consider the conservation laws for mass, momentum and energy, supplemented by the constitutive equations for the deviatoric stress and heat flux from the R13 theory, and the associated augmented set of boundary conditions for a gas–solid interface. We are interested here in the steady-state gas flow and temperature fields, in the gas and solid, resulting from the far-field temperature gradient and gas rarefaction. Assuming that the dimensionless group $a^{\ast }(\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast })_{\infty }/T_{0}^{\ast }\ll 1$ , we can model the transport phenomena using the linearized version of the governing equations written in terms of deviations from the equilibrium state. The non-dimensional form of these deviations can be written as

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}p=p^{\ast }/p_{0}^{\ast },\quad \unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}^{\ast }/\unicode[STIX]{x1D703}_{0}^{\ast },\quad \unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}^{\ast }/\unicode[STIX]{x1D70C}_{0}^{\ast },\\ \text{u}=\text{u}^{\ast }/\unicode[STIX]{x1D703}_{0}^{\ast \,1/2},\quad \text{q}=\text{q}^{\ast }/(p_{0}^{\ast }\unicode[STIX]{x1D703}_{0}^{\ast \,1/2}),\quad \unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}^{\ast }/p_{0}^{\ast },\end{array}\right\}\end{eqnarray}$$

for the pressure, temperature, density, velocity, heat flux and deviatoric stress, respectively, in the gas. The deviatoric stress is symmetric and trace free. Length is non-dimensionalized with the sphere’s radius $a^{\ast }$ . The temperature deviation in the sphere $\unicode[STIX]{x1D703}_{s}^{\ast }$ is non-dimensionalized as that for the gas. The dimensionless form of the temperature gradient defines a new dimensionless group, the Epstein number – coined by Young (Reference Young2011) after P. S. Epstein, a pioneer in the study of thermophoresis of spherical particles, who presented the first theory on the subject (Epstein Reference Epstein1929). It is given by

(2.4) $$\begin{eqnarray}Ep=\frac{a^{\ast }(\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast })_{\infty }}{T_{0}^{\ast }}.\end{eqnarray}$$

Note that the non-dimensional pressure and temperature fields in the gas are given by $1+p$ and $1+\unicode[STIX]{x1D703}$ , respectively, whereas the temperature in the sphere is $1+\unicode[STIX]{x1D703}_{s}$ . On the other hand, $\text{ u}$ , $\text{q}$ and $\unicode[STIX]{x1D748}$ in (2.3) determine the actual velocity, heat flux and deviatoric stress in the gas.

The linearized, steady conservation equations of gas, momentum and energy are

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

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

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

(2.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \left(1-{\textstyle \frac{2}{3}}Kn^{2}\unicode[STIX]{x1D6E5}\right)\unicode[STIX]{x1D748}-{\textstyle \frac{4}{5}}Kn^{2}\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}}=-2Kn\overline{\unicode[STIX]{x1D735}\text{u}}-{\textstyle \frac{4}{5}}Kn\overline{\unicode[STIX]{x1D735}\text{q}}, & \displaystyle\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left(1-{\textstyle \frac{9}{5}}Kn^{2}\unicode[STIX]{x1D6E5}\right)\text{q}=-{\textstyle \frac{15}{4}}Kn\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-{\textstyle \frac{3}{2}}Kn\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6E5}$

(2.7) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{s}=0.\end{eqnarray}$$

Far from the sphere’s surface, the imposed disturbance is represented by the uniform temperature gradient $Ep\,\text{k}$ , where $\text{k}$ designates the unit vector pointing in the direction of the positive $z$ semi-axis. It is convenient to write the problem in such a way that deviations from the base (equilibrium) state vanish in the far field. For this purpose we introduce the transformations

(2.8a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\check{\unicode[STIX]{x1D703}}+Ep\,z,\quad \unicode[STIX]{x1D703}_{s}=\check{\unicode[STIX]{x1D703}}_{s}+Ep\,z,\quad \text{q}=\check{\text{q}}-{\textstyle \frac{15}{4}}Kn\,Ep\,\text{k},\end{eqnarray}$$

whereas the rest of the variables are left unchanged. With (2.8), expressions (2.5)–(2.7) remain invariant in form. On the other hand, the boundary conditions do change after using (2.8). Note that subscript ‘ $s$ ’ will be used to denote quantities in the solid sphere.

Equations (2.5)–(2.7) will be solved subjected to the following boundary conditions taken from Struchtrup et al. (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) in the absence of phase change (see also Struchtrup & Frezzotti Reference Struchtrup and Frezzotti2016) at the solid–gas interface $r=1$ , with $r=|\text{x}|$ , where $\text{x}$ is the position vector with origin at the sphere’s centre. The quantities at the interface corresponding to the liquid in Struchtrup et al. (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) are, in the present work, associated with the solid (sphere). Note that in this case, boundary conditions (37)–(41) in Struchtrup et al. (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) reduce to the wall boundary conditions (33)–(37) in Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2008) if the gas flow is two-dimensional. In all these equations for the boundary the accommodation coefficient has been set equal to one and $\unicode[STIX]{x1D703}$ , $\unicode[STIX]{x1D703}_{s}$ and $\text{q}$ have been substituted according to (2.8).

The generalized slip condition is

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{t_{\unicode[STIX]{x1D6FC}}\,n}=\,-\left(\frac{2}{\unicode[STIX]{x03C0}}\right)^{1/2}\left[(u_{t_{\unicode[STIX]{x1D6FC}}}-u_{s,t_{\unicode[STIX]{x1D6FC}}})+\frac{1}{5}\left(\check{q}_{t_{\unicode[STIX]{x1D6FC}}}-\frac{15}{4}Kn\,Ep\,\text{k}\boldsymbol{\cdot }\text{t}_{\unicode[STIX]{x1D6FC}}\right)+\frac{1}{2}\text{m}_{t_{\unicode[STIX]{x1D6FC}}\,nn}\right].\end{eqnarray}$$

The generalized temperature jump condition is

(2.10) $$\begin{eqnarray}\check{q}_{n}-\frac{15}{4}Kn\,Ep\,\text{k}\boldsymbol{\cdot }\text{n}=-\left(\frac{2}{\unicode[STIX]{x03C0}}\right)^{1/2}\left[2(\check{\unicode[STIX]{x1D703}}-\check{\unicode[STIX]{x1D703}}_{s})+\frac{1}{2}\unicode[STIX]{x1D70E}_{nn}+\frac{5}{28}\text{R}_{nn}\right].\end{eqnarray}$$

The generalized interface conditions for higher moments are

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \text{m}_{nnn}=\left(\frac{2}{\unicode[STIX]{x03C0}}\right)^{1/2}\left[\frac{2}{5}(\check{\unicode[STIX]{x1D703}}-\check{\unicode[STIX]{x1D703}}_{s})-\frac{7}{5}\unicode[STIX]{x1D70E}_{nn}-\frac{1}{14}\text{R}_{nn}\right], & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \text{m}_{t_{\unicode[STIX]{x1D6FC}}t_{\unicode[STIX]{x1D6FD}}n}=-\left(\frac{2}{\unicode[STIX]{x03C0}}\right)^{1/2}\left[\unicode[STIX]{x1D70E}_{t_{\unicode[STIX]{x1D6FC}}t_{\unicode[STIX]{x1D6FD}}}+\frac{1}{14}\text{R}_{t_{\unicode[STIX]{x1D6FC}}t_{\unicode[STIX]{x1D6FD}}}+\frac{1}{5}\{(\check{\unicode[STIX]{x1D703}}-\check{\unicode[STIX]{x1D703}}_{s})-\unicode[STIX]{x1D70E}_{nn}\}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\right], & \displaystyle\end{eqnarray}$$

and

(2.13) $$\begin{eqnarray}\text{R}_{t_{\unicode[STIX]{x1D6FC}}\,n}=\left(\frac{2}{\unicode[STIX]{x03C0}}\right)^{1/2}\left[(u_{t_{\unicode[STIX]{x1D6FC}}}-u_{s,t_{\unicode[STIX]{x1D6FC}}})-\frac{11}{5}\left(\check{q}_{t_{\unicode[STIX]{x1D6FC}}}-\frac{15}{4}Kn\,Ep\,\text{k}\boldsymbol{\cdot }\text{t}_{\unicode[STIX]{x1D6FC}}\right)-\frac{1}{2}\text{m}_{t_{\unicode[STIX]{x1D6FC}}\,nn}\right].\end{eqnarray}$$

Furthermore, the solution must also satisfy the non-penetration condition

(2.14) $$\begin{eqnarray}u_{n}=0,\end{eqnarray}$$

and the interfacial linearized energy balance (Young Reference Young2011)

(2.15) $$\begin{eqnarray}\check{q}_{n}+{\textstyle \frac{15}{4}}Kn\,(\unicode[STIX]{x1D6EC}-1)Ep\,\text{n}\boldsymbol{\cdot }\text{k}=-{\textstyle \frac{15}{4}}Kn\,\unicode[STIX]{x1D6EC}\,\text{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\check{\unicode[STIX]{x1D703}}_{s},\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}$ is the solid-to-gas thermal conductivity ratio $k_{s}^{\ast }/k_{0}^{\ast }$ . Here, indices $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ can take a value of 1 or 2; $\text{n}$ is the unit vector normal to the interface pointing into the gas; $\text{t}_{1}$ and $\text{t}_{2}$ represent two mutually orthogonal unit vectors tangential to the interface, and subscripts $n$ , $t_{1}$ and $t_{2}$ denote components in the corresponding normal or any of the two tangential directions, respectively. In these interfacial conditions one must set $u_{s,t_{1}}=u_{s,t_{2}}=0$ . Regarding the far field, all deviations from the basic equilibrium state must vanish as $r\rightarrow \infty$ .

Boundary conditions (2.14) and (2.15) hold for the interface between a fluid and an impenetrable solid. Their counterparts for a fluid–fluid interface are given by the conditions for mass and energy conservation presented, for instance, in appendix A of Struchtrup et al. (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) – note that in their expression (A5), the internal energy should be written instead of the enthalpy.

The interfacial conditions contain the components of the higher-order moments $\text{R}$ and $\text{m}$ , a rank-two and rank-three tensor, respectively. These are defined as

(2.16a ) $$\begin{eqnarray}\displaystyle & \text{R}=-{\textstyle \frac{24}{5}}Kn\,\overline{\unicode[STIX]{x1D735}\check{\text{q}}}, & \displaystyle\end{eqnarray}$$

(2.16b ) $$\begin{eqnarray}\displaystyle & \text{m}=-2Kn\,\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}}. & \displaystyle\end{eqnarray}$$

$\text{q}$

$\unicode[STIX]{x1D748}$

$\text{R}$

$\text{m}$

$\text{q}$

$\text{R}$

$\text{m}$

It is important to note that because $\unicode[STIX]{x1D748}$ , $\text{R}$ and $\text{m}$ are trace-free tensors, we have that $\unicode[STIX]{x1D70E}_{nn}=-\unicode[STIX]{x1D70E}_{t_{1}t_{1}}-\unicode[STIX]{x1D70E}_{t_{2}t_{2}}$ , $\text{R}_{nn}=-\text{R}_{t_{1}t_{1}}-\text{R}_{t_{2}t_{2}}$ and $\text{m}_{nnn}=-\text{m}_{t_{1}t_{1}n}-\text{m}_{t_{2}t_{2}n}$ . Using these constraints, we can easily show that boundary condition (2.11) can be obtained by writing (2.12) twice, first for $\text{m}_{t_{1}t_{1}n}$ and then for $\text{m}_{t_{2}t_{2}n}$ , and adding the resulting expressions. Therefore, it suffices for the solution to satisfy only one of these equations provided (2.11) is also satisfied.

Since these boundary conditions will be applied to a spherical interface, it is fitting to introduce a system of spherical-polar coordinates $(r,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})$ with its origin, $r=0$ , located at the centre of the solid sphere (figure 1). Here, $0\leqslant r<\infty$ , $0\leqslant \unicode[STIX]{x1D717}\leqslant \unicode[STIX]{x03C0}$ , and $0\leqslant \unicode[STIX]{x1D711}<2\unicode[STIX]{x03C0}$ . Semi-axis $\unicode[STIX]{x1D717}=0$ coincides with the positive $z$ semi-axis. To these coordinate directions correspond unit vectors $\hat{\text{r}}$ , $\hat{\unicode[STIX]{x1D751}}$ and $\hat{\unicode[STIX]{x1D753}}$ , respectively; these triplet forms an orthogonal set. Thus $\{\hat{\text{r}},\hat{\unicode[STIX]{x1D751}},\hat{\unicode[STIX]{x1D753}}\}$ take the place of the set $\{\text{n},\text{t}_{1},\text{t}_{2}\}$ when writing the boundary conditions. The problems considered here are axisymmetric, so quantities do not vary in the $\unicode[STIX]{x1D711}$ direction.

Finally, if needed, the gas density deviation from its value at equilibrium can be computed from the linearized form of the ideal gas equation of state, $p=\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D703}$ .

2.2 Uniform flow past a sphere

Suppose that instead of prescribing a far-field temperature gradient, the state of equilibrium of the gas described in the previous section is disturbed by imposing, far from the sphere, a uniform flow with constant velocity $U_{0}^{\ast }$ in the $z^{\ast }$ -direction. We can model this flow, including rarefaction effects, by using the linearized conservation laws, R13 constitutive relations and boundary conditions written for the disturbances, provided the dimensionless (Reynolds) number $\unicode[STIX]{x1D70C}_{0}^{\ast }U_{0}^{\ast }a^{\ast }/\unicode[STIX]{x1D707}_{0}^{\ast }\ll 1$ . Again, steady state will be assumed. Our work will extend Torrilhon’s (Reference Torrilhon2010) efforts on this problem by including heat conduction throughout the sphere. Furthermore, we will present an expression for the drag force over the sphere from our analytical solution – a closed-form expression for the drag was not given in Torrilhon’s work. Combining this expression for the drag with that for the thermophoretic force will result in the thermophoretic velocity when these two forces balanced each other.

For convenience, we transform the original problem to the equivalent one of the disturbance flow resulting from a sphere translating with dimensionless velocity $-Ma\,\text{k}$ in a fluid at rest far away from the sphere. Here, $Ma$ is a pseudo Mach number

(2.17) $$\begin{eqnarray}Ma=\frac{U_{0}^{\ast }}{\unicode[STIX]{x1D703}_{0}^{\ast \,1/2}}.\end{eqnarray}$$

The actual Mach number can be obtained by multiplying $Ma$ by $\unicode[STIX]{x1D6FE}^{-1/2}$ , with $\unicode[STIX]{x1D6FE}$ the specific heat ratio.

We thus set $\text{u}=\check{\text{u}}+Ma\text{k}$ , so that $\check{\text{u}}\rightarrow 0$ as $r\rightarrow \infty$ . The governing equations and boundary conditions for this problem are obtained from those in § 2.1 by writing $\check{\text{u}}$ instead of $\text{u}$ , dropping the ‘  $\check{~}\,$ ’ from $\check{\unicode[STIX]{x1D703}}$ , $\check{\unicode[STIX]{x1D703}}_{s}$ and $\check{\text{q}}$ and setting $Ep=0$ , ${\check{u}}_{s,t_{1}}={\check{u}}_{s,\unicode[STIX]{x1D717}}=-Ma\,\text{k}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D751}}$ and ${\check{u}}_{s,t_{2}}={\check{u}}_{s,\unicode[STIX]{x1D711}}=0$ . In addition, instead of (2.14) we must have

(2.18) $$\begin{eqnarray}{\check{u}}_{n}=-Ma\,\text{k}\boldsymbol{\cdot }\text{n}.\end{eqnarray}$$

Hereinafter, we drop symbol ‘  $\check{~}\,$ ’ unless otherwise noted.

2.3 Equivalent Stokes–Fourier system of equations

We shall proceed to rewrite the linearized R13-moment equations introduced previously as an equivalent Stokes–Fourier set of equations supplemented by non-homogeneous elliptic equations for the pressure, heat flux and deviatoric stress. To this alternative system of equations, we can apply, somewhat straightforwardly, analytical tools used successfully in scalar and, more notably, vector equations appearing in slow flow hydrodynamics. By introducing the auxiliary variables

(2.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D755} & = & \displaystyle \text{u}+{\textstyle \frac{2}{5}}\text{q},\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701} & = & \displaystyle \unicode[STIX]{x1D703}-{\textstyle \frac{2}{5}}p,\end{eqnarray}$$

it can be readily shown that expressions (2.5)–(2.6) can be written as the equivalent Stokes–Fourier system

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

(2.21b ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D72E}=0, & \displaystyle\end{eqnarray}$$

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

(2.22a ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D72E}=-2Kn\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D755}}, & \displaystyle\end{eqnarray}$$

(2.22b ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D734}=-{\textstyle \frac{15}{4}}Kn\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}, & \displaystyle\end{eqnarray}$$

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D706}_{1}^{2})p=-\unicode[STIX]{x1D706}_{1}^{2}\unicode[STIX]{x1D6F1},\quad \unicode[STIX]{x1D706}_{1}=(5/6)^{1/2}Kn^{-1}, & \displaystyle\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D706}_{2}^{2})\text{q}=-\unicode[STIX]{x1D706}_{2}^{2}\unicode[STIX]{x1D734},\quad \unicode[STIX]{x1D706}_{2}=(5/9)^{1/2}Kn^{-1}, & \displaystyle\end{eqnarray}$$

(2.25) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D706}_{3}^{2})\unicode[STIX]{x1D748}=-\unicode[STIX]{x1D706}_{3}^{2}\,\unicode[STIX]{x1D72E}+{\textstyle \frac{4}{5}}Kn^{2}\unicode[STIX]{x1D706}_{3}^{2}\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}p},\quad \unicode[STIX]{x1D706}_{3}=(3/2)^{1/2}Kn^{\,-1}. & \displaystyle\end{eqnarray}$$

Furthermore, because of (2.21a ) and (2.21c ), (2.21b ) and (2.22b ) lead to

(2.26) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D6F1}=0, & \displaystyle\end{eqnarray}$$

(2.27) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D701}=0, & \displaystyle\end{eqnarray}$$

and, with (2.22a ), expression (2.21b ) becomes

(2.28) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D755}=Kn^{-1}\unicode[STIX]{x1D735}\,\unicode[STIX]{x1D6F1}.\end{eqnarray}$$

Equations (2.21a ) and (2.21c ) result from (2.5a ), (2.5c ) and from the divergence of (2.6b ). Expression (2.21b ) is obtained by taking the divergence of (2.6a ), using (2.5b ) to eliminate $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}$ in favour of $\unicode[STIX]{x1D735}p$ , collecting the terms containing this vector and introducing $\unicode[STIX]{x1D6F1}$ defined in (2.23). Finally, expressions (2.24) and (2.25) are simply (2.6a ) and (2.6b ) written in terms of auxiliary variables $\unicode[STIX]{x1D734}$ and $\unicode[STIX]{x1D72E}$ , respectively. Rewriting the R13-moment equations in a way that includes the more familiar form of (2.21) or, more recognizably, (2.26)–(2.28), will facilitate the application to this set of equations of analytical methods used in the solution of the equations for Stokes flow in vector form. We shall extend these methods to (2.23)–(2.25), which include a rank-two tensor, second-order differential equation.

3.1 Particular solutions

Starting with (2.28), to find a particular solution we let

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D71B}_{i}^{(p)}=Kn^{-1}\unicode[STIX]{x1D6FC}_{i}\unicode[STIX]{x1D6F1},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{i}$ is a vector field to be determined and the subscript $i$ denotes Cartesian index notation. Computing

(3.2) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D71B}_{i}^{(p)}=Kn^{-1}(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{i})\unicode[STIX]{x1D6F1}+2Kn^{-1}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{i}}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}}{\unicode[STIX]{x2202}x_{j}},\end{eqnarray}$$

after using the fact that $\unicode[STIX]{x1D6F1}$ is harmonic, and comparing this result with the left-hand side of (2.28), namely, $Kn^{-1}\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}/\unicode[STIX]{x2202}x_{i}$ , we find $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{i}=0$ and $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{i}/\unicode[STIX]{x2202}x_{j}=\unicode[STIX]{x1D6FF}_{ij}/2$ , where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. Therefore, $\unicode[STIX]{x1D6FC}_{i}=x_{i}/2$ and (Leal Reference Leal2007)

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D755}^{(p)}={\textstyle \frac{1}{2}}Kn^{-1}\text{x}\,\unicode[STIX]{x1D6F1},\end{eqnarray}$$

with $\text{x}$ the position vector from the centre of the sphere. Next, for (2.23) and (2.24), solutions $p^{(p)}$ and $\text{q}^{(p)}$ can be readily found by inspection by using the fact that $\unicode[STIX]{x1D6F1}$ and $\unicode[STIX]{x1D701}$ are harmonic functions. This leads to

(3.4) $$\begin{eqnarray}p^{(p)}=\unicode[STIX]{x1D6F1},\end{eqnarray}$$

and

(3.5) $$\begin{eqnarray}\text{q}^{(p)}=-{\textstyle \frac{15}{4}}Kn\unicode[STIX]{x1D735}\unicode[STIX]{x1D701},\end{eqnarray}$$

respectively.

Finally, solution $\unicode[STIX]{x1D748}^{(p)}$ for (2.25), after eliminating $\unicode[STIX]{x1D72E}$ with (2.22a ), can be found by letting

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D748}^{(p)}=2Kn\,\unicode[STIX]{x1D706}_{3}^{2}\,\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D737}}+{\textstyle \frac{4}{5}}Kn^{2}\,\unicode[STIX]{x1D706}_{3}^{2}\,\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}},\end{eqnarray}$$

such that vector field $\unicode[STIX]{x1D737}$ and scalar field $\unicode[STIX]{x1D702}$ satisfy, respectively,

(3.7a ) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D706}_{3}^{2})\unicode[STIX]{x1D737}=\unicode[STIX]{x1D755}, & \displaystyle\end{eqnarray}$$

(3.7b ) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D706}_{3}^{2})\unicode[STIX]{x1D702}=p. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D706}_{1}^{2}\unicode[STIX]{x1D702}$

(3.8a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D737}=-\unicode[STIX]{x1D755}/\unicode[STIX]{x1D706}_{3}^{2}+m\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}, & \displaystyle\end{eqnarray}$$

(3.8b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=p/(\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2})+n\unicode[STIX]{x1D6F1}. & \displaystyle\end{eqnarray}$$

$m=-Kn^{-1}/\unicode[STIX]{x1D706}_{3}^{4}$

$n=-(\unicode[STIX]{x1D706}_{1}/\unicode[STIX]{x1D706}_{3})^{2}/(\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2})$

(3.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D737}=-{\textstyle \frac{1}{\unicode[STIX]{x1D706}_{3}^{4}}}(\unicode[STIX]{x1D706}_{3}^{2}\unicode[STIX]{x1D755}+Kn^{-1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}), & \displaystyle\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=\frac{1}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}\left(p-\frac{\unicode[STIX]{x1D706}_{1}^{2}}{\unicode[STIX]{x1D706}_{3}^{2}}\unicode[STIX]{x1D6F1}\right), & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D748}^{(p)}$

3.2 Solution of the homogeneous equations

The form of the solutions of the homogeneous equations associated with (2.28) and (2.23)–(2.25) is constructed by means of the method of multipole potentials. This method is presented in the context of classical hydrodynamics and applied to problems of low Reynolds number flows in Leal (Reference Leal2007, chap. 8) as well as thoroughly discussed and used in a wider range of problems of mathematical physics in Hess (Reference Hess2015, chap. 10). In the former, it is named as the method of superposition of vector harmonic functions. The most illustrative example discussed in those sources is perhaps the problem of Stokes flow past a sphere. As we shall see, a significant advantage of the method of multipole potentials concerning our problem is that we can obtain the solution of the deviatoric stress $\unicode[STIX]{x1D748}$ without going through the complexity of writing and solving differential equations for its scalar components in spherical coordinates. The general definitions of multipole potential tensors, the form of the first few of them and properties relevant to this work are presented in appendix A.

Because the governing equations and boundary conditions are linear in the perturbation fields, and the non-homogeneous terms in the boundary conditions must be linear in $G\text{k}$ , with $G=Ep$ or $G=Ma$ for the problems in §§ 2.1 or 2.2, respectively, the homogeneous solutions are constructed by adding products of multipole potentials with vector $G\text{k}$ . From the set of all multipole potential tensors, the method of multipole potential guides the selection of those potentials that should be used to construct the solution. Only those multipole potentials that conform to the rank of the unknown field (i.e. whether it is a scalar, a vector, a rank-two tensor, etc.), its symmetry, and its parity (i.e. whether it is a true scalar, true vector or true tensor or a pseudo-scalar, pseudo-vector or pseudo-tensor), after their product with $G\text{k}$ can be part of the solution. Regarding the parity attribute, in general terms, the components of a true scalar, vector or tensor change sign when subjected to an improper rotation, such as a reflection or inversion – e.g. changing the coordinate system from right-handed to left-handed – in order to continue properly describing an unchanged physical situation (Arfken et al. Reference Arfken, Weber and Harris2012, chap. 3). For a so-called pseudo-scalar, vector or tensor, their components do not change sign accordingly. For instance, the angular velocity is a pseudo-vector. A deeper look at this aspect and, in general, at the method of multipole potentials for constructing solutions of physical quantities is beyond the scope of the present article. For detailed discussions and a list of enlightening examples, the reader may referred to the cited books by Leal (Reference Leal2007) and Hess (Reference Hess2015).

Following these considerations, using Cartesian index notation, and letting $G_{i}$ be the components of vector $G\text{k}$ , we can write for the solutions

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}=A_{1}G_{j}X_{j}, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}=B_{1}G_{j}X_{j}, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D71B}_{i}=C_{1}G_{i}X_{0}+C_{2}G_{j}X_{ji}+\unicode[STIX]{x1D71B}_{i}^{(p)}, & \displaystyle\end{eqnarray}$$

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle p=a(r)\,G_{j}X_{j}+p^{(p)}, & \displaystyle\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle q_{i}=b(r)\,G_{i}X_{0}+c(r)\,G_{j}X_{ji}+q_{i}^{(p)}, & \displaystyle\end{eqnarray}$$

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}=f(r)\,\left[{\textstyle \frac{1}{2}}(G_{i}X_{j}+G_{j}X_{i})-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{ij}\,G_{k}X_{k}\right]+g(r)\,G_{k}X_{kij}+\unicode[STIX]{x1D70E}_{ij}^{(p)}, & \displaystyle\end{eqnarray}$$

where terms with superscript ‘ $(p)$ ’ are given by (3.3)–(3.6); $A_{1}$ , $B_{1}$ , $C_{1}$ and $C_{2}$ , are constants to be determined by the boundary conditions, and fields $X_{0}$ , $X_{i}$ , $X_{ij},\ldots \,$ , denote the descending multipole potentials listed in appendix A. Note that in (3.13)–(3.15) we use functions of the radial coordinate $a(r)$ , $b(r)$ , $c(r)$ , $f(r)$ and $g(r)$ , instead of constants as for (3.10)–(3.12), because the solutions of the homogeneous equations for $p$ , $\text{q}$ and $\unicode[STIX]{x1D748}$ , must satisfy modified Helmholtz equations rather than Laplace equations. The descending multipole potentials satisfy Laplace equation for $r\neq 0$ . In the present problem, $r=0$ is outside the gas domain.

For the interior of the sphere, the temperature profile is given by the ascending multipole potential

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{s}=D_{1}r^{3}G_{j}X_{j},\end{eqnarray}$$

which introduces another constant, $D_{1}$ . Both the equivalence in tensor rank and parity between the left- and right-hand sides of (3.10)–(3.16) determined the type of products employed between the external perturbation vector $G_{i}$ and the multipole potential tensors in these expressions.

Because $\unicode[STIX]{x1D701}$ is harmonic, equation (3.5) yields $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{q}^{(p)}=0$ . With the divergence-free vector $\text{q}=\text{q}^{(h)}+\text{q}^{(p)}$ , where superscript ‘( $h$ )’ denotes the homogeneous solution, we are left with $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{q}^{(h)}=0$ . Similarly, one can show that taking the divergence of $\unicode[STIX]{x1D748}^{(p)}$ using (3.6) results in $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(p)}=-\unicode[STIX]{x1D735}p$ . Therefore, since $\unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}^{(h)}+\unicode[STIX]{x1D748}^{(p)}$ , and because of expression (2.5b ), we must also have $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(h)}=0$ . Note in (3.15) that, by construction, $\unicode[STIX]{x1D748}^{(h)}$ is also symmetric and trace free.

Taking the divergence of $\unicode[STIX]{x1D755}$ using (3.12) and (3.3), setting it equal to zero and using the definitions of the multipole potentials from appendix A results in

(3.17) $$\begin{eqnarray}C_{1}={\textstyle \frac{1}{2}}Kn^{-1}A_{1}.\end{eqnarray}$$

With $\text{q}^{(h)}$ from (3.14), setting $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{ q}^{(h)}=0$ yields this constraint for functions $b(r)$ and $c(r)$

(3.18) $$\begin{eqnarray}r^{2}b^{\prime }+2c^{\prime }-rb=0.\end{eqnarray}$$

Furthermore, the vector equation $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(h)}=0$ , with $\unicode[STIX]{x1D748}^{(h)}$ from (3.15), leads to

(3.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle 2r^{2}f^{\prime }-rf+18g^{\prime }=0, & \displaystyle\end{eqnarray}$$

(3.19b ) $$\begin{eqnarray}\displaystyle & 32r^{2}f^{\prime }+rf-18g^{\prime }=0, & \displaystyle\end{eqnarray}$$

$f^{\prime }(r)\equiv 0$

$g^{\prime }(r)=rf(r)/18$

$p$

$\text{q}$

$\unicode[STIX]{x1D748}$

(3.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle a^{\prime \prime }-2r^{-1}a^{\prime }-\unicode[STIX]{x1D706}_{1}^{2}a=0, & \displaystyle\end{eqnarray}$$

(3.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle b^{\prime \prime }-\unicode[STIX]{x1D706}_{2}^{2}b=0, & \displaystyle\end{eqnarray}$$

(3.20c ) $$\begin{eqnarray}\displaystyle & \displaystyle c^{\prime \prime }-4r^{-1}c^{\prime }-\unicode[STIX]{x1D706}_{2}^{2}c=0, & \displaystyle\end{eqnarray}$$

(3.20d ) $$\begin{eqnarray}\displaystyle & \displaystyle f^{\prime \prime }-2r^{-1}f^{\prime }-\unicode[STIX]{x1D706}_{3}^{2}f=0, & \displaystyle\end{eqnarray}$$

(3.20e ) $$\begin{eqnarray}\displaystyle & \displaystyle g^{\prime \prime }-6r^{-1}g^{\prime }-\unicode[STIX]{x1D706}_{3}^{2}g=0. & \displaystyle\end{eqnarray}$$

$f^{\prime }(r)\equiv 0$

$f(r)\equiv 0$

$g^{\prime }(r)\equiv 0$

$g(r)\equiv 0$

$\unicode[STIX]{x1D748}^{(h)}(\text{x})\equiv 0$

$\unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}^{(p)}$

(3.21a ) $$\begin{eqnarray}\displaystyle & \displaystyle a(r)=E_{1}\exp [-\unicode[STIX]{x1D706}_{1}(r-1)](1+\unicode[STIX]{x1D706}_{1}r), & \displaystyle\end{eqnarray}$$

(3.21b ) $$\begin{eqnarray}\displaystyle & \displaystyle b(r)=F_{1}\exp [-\unicode[STIX]{x1D706}_{2}(r-1)], & \displaystyle\end{eqnarray}$$

(3.21c ) $$\begin{eqnarray}\displaystyle & \displaystyle c(r)=F_{2}\exp [-\unicode[STIX]{x1D706}_{2}(r-1)](1+\unicode[STIX]{x1D706}_{2}r+\unicode[STIX]{x1D706}_{2}^{2}r^{2}/3), & \displaystyle\end{eqnarray}$$

$r$

(3.22) $$\begin{eqnarray}F_{1}=-{\textstyle \frac{2}{3}}\unicode[STIX]{x1D706}_{2}^{2}F_{2}.\end{eqnarray}$$

Note that one could have also written explicit expressions for $f(r)$ and $g(r)$ from solving (3.20d ) and (3.20e ) in a manner similar to that followed to obtain (3.21), only to find out that, after enforcing $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(h)}=0$ , the associated integration constants would be equal to zero, thus yielding $f=g\equiv 0$ .

According to the R13 theory, the structure of the Knudsen layer near a boundary is partially determined by the factors $\exp (-\unicode[STIX]{x1D706}_{i}^{2}r)$ in the solution ( $i=1,2,3$ for the present case). The contribution with $\unicode[STIX]{x1D706}_{3}$ vanishes in this case resulting from the vanishing of $\unicode[STIX]{x1D748}^{(h)}$ . The factor containing the exponential with $\unicode[STIX]{x1D706}_{1}$ is present in the pressure, temperature and, through the pressure, in the deviatoric stress. On the other hand, the factor containing the exponential with $\unicode[STIX]{x1D706}_{2}$ occurs in the heat flux and velocity. All these features were already noted by Torrilhon (Reference Torrilhon2010) in the solution for the special case of slow flow past an isothermal sphere.

3.3 Final field expressions – system of equations for the integration constants

We can now write the final form of the exact solutions for the field variables. The pressure deviation in the gas $p$ can be computed from (3.13), using (3.4) and (3.10). The temperature deviation $\unicode[STIX]{x1D703}$ is found from (2.20), with (3.11) and (3.13). The expression for the gas velocity $\text{u}$ is obtained from (2.19), using (3.12) and (3.14), together with (3.3), (3.10), (3.5) and (3.11). The heat-flux vector $\text{q}$ is computed using (3.14), (3.5) and (3.12), whereas the stress tensor $\unicode[STIX]{x1D748}$ is determined from expressions (3.15) and (3.6), with (3.9), (3.10), (3.12) and (3.13). The temperature deviation in the solid $\unicode[STIX]{x1D703}_{s}$ is given by (3.16). After substitution of (3.17), (3.21) and (3.22), the final expressions for the deviations are

(3.23a ) $$\begin{eqnarray}\displaystyle p & = & \displaystyle A_{1}Gr^{-2}\cos \unicode[STIX]{x1D717}+E_{1}Gr^{-2}(\unicode[STIX]{x1D706}_{1}r+1)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.23b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703} & = & \displaystyle {\textstyle \frac{1}{5}}(2A_{1}+5B_{1})Gr^{-2}\cos \unicode[STIX]{x1D717}+{\textstyle \frac{2}{5}}E_{1}Gr^{-2}(\unicode[STIX]{x1D706}_{1}r+1)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\cos \unicode[STIX]{x1D717},\qquad\end{eqnarray}$$

(3.23c ) $$\begin{eqnarray}\displaystyle u_{r} & = & \displaystyle Kn^{-1}Gr^{-3}(A_{1}r^{2}-3Kn^{2}B_{1}+2KnC_{2})\cos \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle -\,{\textstyle \frac{4}{5}}F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.23d ) $$\begin{eqnarray}\displaystyle u_{\unicode[STIX]{x1D717}} & = & \displaystyle -{\textstyle \frac{1}{2}}Kn^{-1}Gr^{-3}(A_{1}r^{2}+3Kn^{2}B_{1}-2KnC_{2})\sin \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle -\,{\textstyle \frac{2}{5}}F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}^{2}r^{2}+\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.23e ) $$\begin{eqnarray}\displaystyle q_{r} & = & \displaystyle {\textstyle \frac{15}{2}}KnB_{1}Gr^{-3}\cos \unicode[STIX]{x1D717}+2F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.23f ) $$\begin{eqnarray}\displaystyle q_{\unicode[STIX]{x1D717}} & = & \displaystyle {\textstyle \frac{15}{4}}KnB_{1}Gr^{-3}\sin \unicode[STIX]{x1D717}+F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}^{2}r^{2}+\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.23g ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{rr} & = & \displaystyle \displaystyle \frac{2}{5}Gr^{-4}(5A_{1}r^{2}-12Kn^{2}A_{1}-30\unicode[STIX]{x1D706}_{3}^{-2}A_{1}+30KnC_{2})\cos \unicode[STIX]{x1D717}+\displaystyle \frac{8}{15}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}\nonumber\\ \displaystyle & & \displaystyle \times \,Kn^{2}E_{1}Gr^{-4}(\unicode[STIX]{x1D706}_{1}^{3}r^{3}+4\unicode[STIX]{x1D706}_{1}^{2}r^{2}+9\unicode[STIX]{x1D706}_{1}r+9)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.23h ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r} & = & \displaystyle -{\textstyle \frac{6}{5}}Gr^{-4}(2Kn^{2}A_{1}+5\unicode[STIX]{x1D706}_{3}^{-2}A_{1}-5KnC_{2})\sin \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{4}{5}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}Kn^{2}E_{1}Gr^{-4}(\unicode[STIX]{x1D706}_{1}^{2}r^{2}+3\unicode[STIX]{x1D706}_{1}r+3)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.23i ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}_{s} & = & \displaystyle D_{1}Gr\cos \unicode[STIX]{x1D717}.\end{eqnarray}$$

$\unicode[STIX]{x1D748}$

$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}\unicode[STIX]{x1D719}}=-\unicode[STIX]{x1D70E}_{rr}/2$

$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}r}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D719}}=0$

To compute the higher-order moment tensors $\text{R}$ and $\text{m}$ defined in (2.16a ) and (2.16b ), respectively, we first calculate rank-two and rank-three tensors $\unicode[STIX]{x1D735}\text{q}$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}$ , respectively. We compute these gradients by means of a popular symbolic algebra and calculus computer package (Wolfram Mathematica), leading to expressions that we do not reproduce here for the sake of space. An alternative path to $\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}$ is to extract components $\unicode[STIX]{x1D70E}_{rr}$ , $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$ and $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}}$ from $\unicode[STIX]{x1D748}$ , compute their gradients and combine them with the dyads and the gradients of the dyads formed by the spherical coordinate unit vectors. Another option is to apply the formulae for the components in spherical coordinates of the gradient of a rank-two tensor field given in Torrilhon (Reference Torrilhon2010). Once $\unicode[STIX]{x1D735}\text{q}$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}$ has been determined, their corresponding symmetric and trace-free tensors are computed from the formulae in appendix C. The expressions for the components of $\text{R}$ and $\text{m}$ needed in the boundary conditions are

(3.24a ) $$\begin{eqnarray}\displaystyle \text{R}_{rr} & = & \displaystyle 108Kn^{2}B_{1}Gr^{-4}\cos \unicode[STIX]{x1D717}+{\textstyle \frac{48}{5}}KnF_{2}Gr^{-4}\nonumber\\ \displaystyle & & \displaystyle \times \exp [-\unicode[STIX]{x1D706}_{2}(r-1)](\unicode[STIX]{x1D706}_{2}^{2}r^{2}+3\unicode[STIX]{x1D706}_{2}r+3)\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.24b ) $$\begin{eqnarray}\displaystyle \text{R}_{\unicode[STIX]{x1D717}r} & = & \displaystyle 54Kn^{2}B_{1}Gr^{-4}\sin \unicode[STIX]{x1D717}+\frac{12}{5}KnF_{2}Gr^{-4}\nonumber\\ \displaystyle & & \displaystyle \times \exp [-\unicode[STIX]{x1D706}_{2}(r-1)](\unicode[STIX]{x1D706}_{2}^{3}r^{3}+3\unicode[STIX]{x1D706}_{2}^{2}r^{2}+6\unicode[STIX]{x1D706}_{2}r+6)\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.24c ) $$\begin{eqnarray}\displaystyle \text{m}_{rrr} & = & \displaystyle -{\textstyle \frac{48}{5}}KnGr^{-5}(10\unicode[STIX]{x1D706}_{3}^{-2}A_{1}+4Kn^{2}A_{1}-r^{2}A_{1}-10KnC_{2})\cos \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{16}{25}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}Kn^{3}E_{1}Gr^{-5}\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\nonumber\\ \displaystyle & & \displaystyle \times \,(\unicode[STIX]{x1D706}_{1}^{4}r^{4}+7\unicode[STIX]{x1D706}_{1}^{3}r^{3}+27\unicode[STIX]{x1D706}_{1}^{2}r^{2}+60\unicode[STIX]{x1D706}_{1}r+60)\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

(3.24d ) $$\begin{eqnarray}\displaystyle \text{m}_{\unicode[STIX]{x1D717}rr} & = & \displaystyle -{\textstyle \frac{8}{5}}KnGr^{-5}(30\unicode[STIX]{x1D706}_{3}^{-2}A_{1}+12Kn^{2}A_{1}-r^{2}A_{1}-30KnC_{2})\nonumber\\ \displaystyle & & \displaystyle \times \sin \unicode[STIX]{x1D717}+\frac{32}{25}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}Kn^{3}E_{1}Gr^{-5}\nonumber\\ \displaystyle & & \displaystyle \times \exp [-\unicode[STIX]{x1D706}_{1}(r-1)](\unicode[STIX]{x1D706}_{1}^{3}r^{3}+6\unicode[STIX]{x1D706}_{1}^{2}r^{2}+15\unicode[STIX]{x1D706}_{1}r+15)\sin \unicode[STIX]{x1D717}.\end{eqnarray}$$

$\text{R}_{rr}$

$\text{R}_{\unicode[STIX]{x1D717}r}$

$\text{m}_{rrr}$

$\text{m}_{\unicode[STIX]{x1D717}rr}$

Regarding the boundary conditions, from the solutions we also find that $\text{R}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}}=\text{R}_{\unicode[STIX]{x1D719}\unicode[STIX]{x1D719}}$ , and $\text{m}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}r}=\text{m}_{\unicode[STIX]{x1D719}\unicode[STIX]{x1D719}r}$ . Because of this, and referring to the discussion at the end of § 2.1, we have that by enforcing (2.11), constraint (2.12) can be dropped entirely. Moreover, equations (2.9) and (2.13) provide constraints for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$ and $\text{R}_{\unicode[STIX]{x1D717}r}$ whereas they are trivially satisfied for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}r}$ and $\text{R}_{\unicode[STIX]{x1D719}r}$ . Therefore, including (2.10), (2.14) and (2.15), in total, we are left with six scalar boundary conditions.

Substitution of the solutions for the various fields into these boundary conditions leads to a system of six linear algebraic equations for the six unknowns $A_{1}$ , $B_{1}$ , $C_{2}$ , $D_{1}$ , $E_{1}$ and $F_{2}$ . For the problem of a sphere exposed to a uniform temperature gradient in the far field, this system may be written as

(3.25a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}-18Kn(256\sqrt{2}Kn^{4}+64\sqrt{\unicode[STIX]{x03C0}}Kn^{3}-8\sqrt{2}Kn^{2}+5\sqrt{2})A_{1}-135\sqrt{2}Kn^{3}B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,180Kn^{2}(24\sqrt{2}Kn^{2}+6\sqrt{\unicode[STIX]{x03C0}}Kn+\sqrt{2})C_{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,18Kn^{2} [\!216\sqrt{2}Kn^{3}+18(4\sqrt{15}+3\sqrt{\unicode[STIX]{x03C0}})Kn^{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,9\sqrt{2}(\sqrt{15\unicode[STIX]{x03C0}}+8)Kn+4\sqrt{15}+15\sqrt{\unicode[STIX]{x03C0}}\!]E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,4\sqrt{2}(9Kn^{2}+3\sqrt{5}Kn+5)F_{2}=-135\sqrt{2}Kn^{3},\end{eqnarray}$$

(3.25b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}-84\sqrt{2}Kn(32Kn^{2}-9)A_{1}+30Kn(270\sqrt{2}Kn^{2}+105\sqrt{\unicode[STIX]{x03C0}}Kn+28\sqrt{2})B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,2520\sqrt{2}Kn^{2}C_{2}-840\sqrt{2}KnD_{1}-42(54\sqrt{2}Kn^{3}+18\sqrt{15}Kn^{2}+12\sqrt{2}Kn-\sqrt{15})E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,40[54\sqrt{2}Kn^{2}+3(6\sqrt{10}+7\sqrt{\unicode[STIX]{x03C0}})Kn+10\sqrt{2}+7\sqrt{5\unicode[STIX]{x03C0}}]F_{2}=1575\sqrt{\unicode[STIX]{x03C0}}Kn^{2},\end{eqnarray}$$

(3.25c ) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}-42Kn(1280\sqrt{\unicode[STIX]{x03C0}}Kn^{3}+224\sqrt{2}Kn^{2}-120\sqrt{\unicode[STIX]{x03C0}}Kn-33\sqrt{2})A_{1}+30\sqrt{2}Kn(135Kn^{2}-7)B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +1260Kn^{2}(40\sqrt{\unicode[STIX]{x03C0}}Kn+7\sqrt{2})C_{2}+210\sqrt{2}KnD_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,21[\!2160\sqrt{\unicode[STIX]{x03C0}}Kn^{4}+18\sqrt{2}(20\sqrt{15\unicode[STIX]{x03C0}}+21)Kn^{3}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,18(7\sqrt{15}+45\sqrt{\unicode[STIX]{x03C0}})Kn^{2}+\sqrt{2}(35\sqrt{15\unicode[STIX]{x03C0}}+144)Kn+13\sqrt{15}+25\sqrt{\unicode[STIX]{x03C0}}\!]E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,40\sqrt{2}(27Kn^{2}+9\sqrt{5}Kn+5)F_{2}=0,\end{eqnarray}$$

(3.25d ) $$\begin{eqnarray}\displaystyle & & \displaystyle -18\sqrt{2}Kn(256Kn^{4}-8Kn^{2}-5)A_{1}+135Kn^{3}(72\sqrt{\unicode[STIX]{x03C0}}Kn+13\sqrt{2})B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,180\sqrt{2}Kn^{2}(24Kn^{2}-1)C_{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,72Kn^{2}(54\sqrt{2}Kn^{3}+18\sqrt{15}Kn^{2}+18\sqrt{2}Kn+\sqrt{15})E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,4[\!648\sqrt{\unicode[STIX]{x03C0}}Kn^{3}+9(13\sqrt{2}+24\sqrt{5\unicode[STIX]{x03C0}})Kn^{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,3(13\sqrt{10}+60\sqrt{\unicode[STIX]{x03C0}})Kn+65\sqrt{2}+20\sqrt{5\unicode[STIX]{x03C0}}\!]F_{2}=-1485\sqrt{2}Kn^{3},\end{eqnarray}$$

(3.25e ) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}A_{1}-3Kn^{2}B_{1}+2KnC_{2}-{\textstyle \frac{4}{15}}(3Kn+\sqrt{5})F_{2}=0,\end{eqnarray}$$

(3.25f ) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}{\textstyle \frac{15}{2}}Kn^{2}B_{1}+{\textstyle \frac{15}{4}}\unicode[STIX]{x1D6EC}Kn^{2}D_{1}+{\textstyle \frac{2}{3}}(3Kn+\sqrt{5})F_{2}=-{\textstyle \frac{15}{4}}(\unicode[STIX]{x1D6EC}-1)Kn^{2}.\end{eqnarray}$$

For the velocity problem, the left-hand sides of expressions (3.25) remain the same, whereas the right-hand sides for (3.25a ), (3.25b ), (3.25d ), (3.25e ) and (3.25f ) become, respectively, $180\sqrt{2}Kn^{2}$ , 0, $-180\sqrt{2}Kn^{2}$ , $-Kn$ and 0. The solution of this linear system of equations is found by means of a computer symbolic algebra package. The expressions for the coefficients are exceedingly large and are not reproduced here. They depend on the parameters $Kn$ and $\unicode[STIX]{x1D6EC}$ . Once these coefficients are determined, the field variables can be computed at any position $(r,\unicode[STIX]{x1D717})$ , for given values of parameters $Kn$ , $\unicode[STIX]{x1D6EC}$ and $Ep$ (or $Ma$ ), using (3.23).

In general terms, the main attributes of the solution method proposed here are twofold. First, it replaces an a priori knowledge of the exact dependency of the various scalar fields, including spherical-coordinate components of vectors and tensors, on the polar angle – e.g. from the classical solution of Stokes flow past a sphere – with the rather more general constraint of a dependency of the various tensor fields (scalar, vectors and rank-two tensors) on the external perturbation vector $\text{G}$ multiplied by a few multipole potentials appropriately chosen. Second, instead of dealing with the components of the Laplacian of the stress tensor in spherical coordinates and with the scalar differential equations carrying them, we solve a second-order differential equation for the stress in compact, tensor form.

4.1 Thermophoresis on a heat-conducting sphere

4.1.1 Flow and temperature fields

Examining the response of some of the variables modelled with the R13 equations to changes in the Knudsen number $Kn$ and solid-to-gas conductivity ratio $\unicode[STIX]{x1D6EC}$ near the surface of the sphere is of interest. In figure 2, we show the profiles of the velocity components, density and temperature deviations as functions of the radial coordinate starting at the surface of the sphere, $r=1$ , obtained with the R13 exact solution derived here, the numerical solutions of the linearized Boltzmann equation by Sone (Reference Sone2007) for a hard-sphere gas, and his asymptotic expression for small Knudsen number. Curves are presented for an isothermal sphere, i.e. $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , and $Kn=0$ , 0.090, 0.180, 0.269 and 0.539, corresponding to Sone’s $k=0$ , 0.1, 0.2, 0.3 and 0.6, respectively – the relationship between $Kn$ and $k$ is given in the figure’s caption. Predictions from the theories depict typical rarefaction effects, namely, velocity slip (figure 2 b), and temperature jump (figure 2 d) at the sphere’s surface. When $k\rightarrow 0$ , and $k=0.1$ , R13 agrees well with the results from kinetic theory; however, as $k$ increases the differences between the two models becomes noticeable, suggesting that R13’s quantitative description of the Knudsen layer near the surface of the sphere begins to become incomplete for $Kn\gtrsim 0.1$ .

Figure 2. Profiles of (a) radial velocity, (b) polar velocity, (c) density deviation and (d) temperature deviation in the gas as functions of the radial coordinate for the problem of thermophoresis of a sphere with a uniform temperature (i.e. $\unicode[STIX]{x1D6EC}\rightarrow \infty$ ). Results are obtained from the R13 exact solution, the numerical solution of Sone (Reference Sone2007) and his asymptotic expression for $k\rightarrow 0$ . Sone’s results are for a hard-sphere gas. Knudsen number $Kn$ is related to parameter $k$ as $Kn=\sqrt{2}\unicode[STIX]{x1D6FE}_{1}k/2$ , with $\unicode[STIX]{x1D6FE}_{1}=1.270042427$ from Sone (Reference Sone2007).

Figure 3. Speed contours and velocity streamlines in the case of thermophoresis of a sphere for $Kn=0.02$ and 0.2, and $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ computed with the exact solution from R13. Far-field temperature gradient points to the right.

Contour plots of gas speed (normalized by $Ep$ ) and velocity streamlines are presented in figure 3 for combinations of $Kn=0.02$ and 0.2 and $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . The temperature gradient points from left to right. Note that except for the case of $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , the gas flow near the sphere is in the direction of the temperature gradient. Following the reasoning by Sone (Reference Sone2007), this flow is driven by the force exerted by the solid surface on the gas in this direction. This force is the reaction to the momentum transferred onto the sphere’s surface by the gas in the opposite direction, that is, from the hot to the cold region. This corresponds to a scenario of normal or ‘positive’ thermophoresis. For the combination $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , on the other hand, the situation is inverted, and the gas flow is from hot to cold. This case is thus known as reversed or ‘negative’ thermophoresis. We shall discuss this phenomenon below with reference to the macroscopic transport equations of the previous section and, in particular, the boundary condition for slip. This is followed by a quantitative investigation of the net thermophoretic force acting on the sphere with special attention to its change in direction. Another feature shown in figure 3 is that for $Kn=0.02$ , the maximum speed in the case $\unicode[STIX]{x1D6EC}\rightarrow \infty$ is approximately one order of magnitude smaller than in the case $\unicode[STIX]{x1D6EC}=4$ , and nearly two orders of magnitude smaller than the other two cases corresponding to $Kn=0.2$ . One may thus expect that, among the four cases considered, the smallest net force on the sphere should be attained in the case $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ .

Contour plots of gas temperature deviation (normalized by $Ep$ ) with streamlines for the heat-flux vector in the gas and within the solid sphere are shown in figure 4 for four cases with $Kn=0.02$ and 0.2 and $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . The heat-flux vectors point to the left, opposite to the temperature gradient. In the cases where $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , figures show no temperature deviation inside the sphere. In contrast, for $\unicode[STIX]{x1D6EC}=4$ , a non-zero temperature gradient is observed in the spherical particle, with the temperature varying linearly with the axial coordinate. In the two cases with $\unicode[STIX]{x1D6EC}=4$ and also in the case where $\unicode[STIX]{x1D6EC}\rightarrow \infty$ and $Kn=0.2$ , temperature jumps across the spherical surface are evident. Enhancing the gas rarefaction by increasing $Kn$ for fixed $\unicode[STIX]{x1D6EC}$ increases the temperature jump across the spherical surface. For instance, with $\unicode[STIX]{x1D6EC}=4$ , there are more contour lines inside the sphere for $Kn=0.02$ than for 0.2, pointing toward a smoother temperature gradient in the latter, whereas the opposite takes place in the gas, as the isothermal lines are more bent for $Kn=0.02$ than for 0.2. This yields a greater temperature jump at a given point of the sphere’s surface ( $z\neq 0$ ) in the case with higher $Kn$ .

Figure 4. Temperature contours and heat-flux streamlines in the case of thermophoresis of a sphere for $Kn=0.02$ and 0.2, and $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ computed with the exact solution from R13. Between the contour levels 0 and 1.0 ( $-1.0$ ), the levels shown correspond to 0.1, 0.3 and 0.5 ( $-0.1$ , $-0.3$ and $-0.5$ ); then they increase by 0.5 (decrease by $-0.5$ ). Far-field temperature gradient points to the right.

Another notable feature in figure 4 is that lines of constant temperature intercept the surface of the sphere at various points in the two cases where $\unicode[STIX]{x1D6EC}=4$ as well as for $Kn=0.2$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , indicating the presence of a temperature gradient in the tangential direction. Under rarefied conditions, it is known that this temperature gradient induces gas motion along its direction (i.e. from cold to hot) near the solid surface, an effect known as thermal creep or thermal transpiration (Maxwell Reference Maxwell1879; Kennard Reference Kennard1938; Sone Reference Sone2007; Mohammadzadeh et al. Reference Mohammadzadeh, Rana and Struchtrup2015). This is the type of flow depicted by the streamlines in the corresponding plots of figure 3. On the other hand, for $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , the isothermal lines in the gas tend to wrap the sphere’s surface resulting in temperature gradients that essentially vanish in the tangential direction on the gas side of the sphere’s surface. This points toward a hindering of the thermal creep in comparison with the other cases, in accord with the small velocity magnitudes and, more importantly, the reversal in the direction of the streamlines shown for the same case in figure 3. In fact, this reversed flow, now from the hot to the cold region, occurs because another type of flow, the so-called thermal-stress slip flow (Sone Reference Sone2007; Young Reference Young2011), becomes dominant over the thermal creep. The thermal-stress slip flow represents a Knudsen number higher-order effect and, even though it is also induced by changes in the gas temperature distribution, is of a different nature from the thermal creep.

From his asymptotic analysis of the linearized Boltzmann equation for a gas that deviates slightly from a state of uniform equilibrium at rest, Sone (Reference Sone2007) identified that the slip flow in the slip boundary condition was determined by the term proportional to $-\text{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\text{n}$ multiplied by a positive constant when the boundary surface has a uniform temperature or in the absence of thermal creep (here $\text{t}$ and $\text{n}$ are unit vectors tangential and normal to the bounding surface, respectively). Because $\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$ multiplied by a constant is one of the terms in his expression for the stress tensor, he designated this flow as thermal-stress slip flow. As discussed by Sone (Reference Sone2007), if the boundary temperature is constant then $-\text{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\text{n}=-\text{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\text{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})=-\unicode[STIX]{x2202}(\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}n)/\unicode[STIX]{x2202}s$ , where $s$ is the arc length. Then, if the isothermal surfaces in the gas are not parallel to the solid boundary (also of constant temperature), i.e. the component of the temperature gradient normal to the boundary ( $\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}n$ ) changes along it, when the boundary temperature is higher (lower) than that in the gas, a flow is promoted in the direction in which the isothermal surfaces converge (diverge). This is exemplified by the contour plots in figures 3 and 4 for $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , and figure 5 sketches this behaviour for clarity. The thermal-stress slip flow is the primary cause of a negative thermophoretic force.

Figure 5. Sketch of the thermal-stress slip flow on the surface of a sphere (gas motion from hot to cold). The thin lines represent isothermal surfaces in the gas; the thick line represents the sphere’s surface (with uniform temperature $\unicode[STIX]{x1D703}_{s}$ ) and the dashed line is the axis of symmetry.

Recalling the generalized slip boundary condition (2.9), rewritten here for convenience in a slightly different form

(4.1) $$\begin{eqnarray}u_{\unicode[STIX]{x1D717}}-u_{s,\unicode[STIX]{x1D717}}=-\frac{1}{5}q_{\unicode[STIX]{x1D717}}-\left(\frac{\unicode[STIX]{x03C0}}{2}\right)^{1/2}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}-\frac{1}{2}\text{m}_{\unicode[STIX]{x1D717}rr},\end{eqnarray}$$

setting $u_{s,\unicode[STIX]{x1D717}}=0$ and using the expressions from the analytical solution obtained with the R13-moment equations, we investigate the contribution of each of the terms on the right-hand side of (4.1) to the gas slip velocity $u_{\unicode[STIX]{x1D717}}$ on the surface of the sphere, $r=1$ , and selected $\unicode[STIX]{x1D717}$ values. Considering for this exercise the conditions of figures 3 and 4, we have that for the cases of $\unicode[STIX]{x1D6EC}=4$ and also for $\unicode[STIX]{x1D6EC}\rightarrow \infty$ and $Kn=0.2$ , where thermal creep is predominant, the greater contribution in absolute value is from the first term on the right-hand side, giving a negative value (flow from cold to hot), whereas for the remaining case of $\unicode[STIX]{x1D6EC}\rightarrow \infty$ and $Kn=0.02$ , exhibiting thermal-stress slip flow, the larger value is positive and results from the shear-stress term $-(\unicode[STIX]{x03C0}/2)^{1/2}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$ . In all cases, the last term in (4.1), a Knudsen number higher-order term, results in a significantly smaller absolute value in comparison. Using now the balance equation for the heat flux (2.6b ) for $q_{\unicode[STIX]{x1D717}}$ when the flow has the direction of the temperature gradient or the balance equation for the deviatoric stress (2.6a ) for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$ when the flow has the opposite direction – after moving the terms with derivatives of $\text{q}$ and $\unicode[STIX]{x1D748}$ , of higher order in $Kn$ , to the right-hand side – we found that in the former, the term $(3/4)Kn\,\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}\unicode[STIX]{x1D717}$ carries the largest weight whereas in the latter, the term proportional to $Kn(\overline{\unicode[STIX]{x1D735}\text{ q}})_{\unicode[STIX]{x1D717}r}$ is predominant. Note that this term and not the shear stress from classical hydrodynamics, $-2Kn(\overline{\unicode[STIX]{x1D735}\text{ u}})_{\unicode[STIX]{x1D717}r}$ , is the leading contributor in this case. Examination of the term proportional to $Kn(\overline{\unicode[STIX]{x1D735}\text{ q}})_{\unicode[STIX]{x1D717}r}$ by substitution of (2.6b ) for $\text{q}$ when the slip flow is from hot to cold revealed that, as expected, the major role in this situation is played by the term $-3(\unicode[STIX]{x03C0}/2)^{1/2}Kn^{2}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})_{\unicode[STIX]{x1D717}r}$ ( ${\approx}-3(\unicode[STIX]{x03C0}/2)^{1/2}Kn^{2}\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}r\unicode[STIX]{x2202}\unicode[STIX]{x1D717}$ when $\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}\unicode[STIX]{x1D717}\approx 0$ ) in agreement with the argument of Sone (Reference Sone2007) for the thermal-stress slip flow and the illustration in figure 5. From a microscopic perspective, Sone (Reference Sone2007) emphasizes that for either thermal creep or thermal-stress slip flow, the gas motion is induced by the difference in the velocity distribution functions of the molecules colliding with the solid boundary and those leaving it. The different nature of the flow and force between thermal creep and thermal-stress slip flow is affected by the fact that, before impinging onto the solid boundary, particles starting within a mean free path essentially keep the attributes of their origins, and the temperature distribution in the gas surrounding the boundary is notably different in one case or the other.

4.1.2 Thermophoretic force

The thermophoretic force acting on the sphere is obtained by integrating over the surface of the sphere the projection of the total stress vector at $r^{\ast }=a^{\ast }$ onto direction $\text{k}$ , i.e. $(-p^{\ast }\hat{\text{r}}-\unicode[STIX]{x1D748}^{\ast }\boldsymbol{\cdot }\hat{\text{r}})\boldsymbol{\cdot }\text{k}$ . Young (Reference Young2011) introduced the dimensionless thermophoretic force

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=\frac{1}{Ep\,Kn}\frac{F_{T}^{\ast }}{\unicode[STIX]{x1D707}_{0}^{\ast }\,\unicode[STIX]{x1D703}_{0}^{\ast \,1/2}a^{\ast }},\end{eqnarray}$$

where $F_{T}^{\ast }$ denotes the dimensional thermophoretic force. From our solution of the R13 equations for the pressure and the deviatoric stress, we obtain the exact expression

(4.3) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=-12\unicode[STIX]{x03C0}\,\frac{\mathop{\sum }_{m=0}^{7}(\unicode[STIX]{x1D6FC}_{m}^{(0)}+\unicode[STIX]{x1D6FC}_{m}^{(1)}\unicode[STIX]{x1D6EC})Kn^{\,m}}{\mathop{\sum }_{m=0}^{9}(\unicode[STIX]{x1D6FD}_{m}^{(0)}+\unicode[STIX]{x1D6FD}_{m}^{(1)}\unicode[STIX]{x1D6EC})Kn^{\,m}},\end{eqnarray}$$

and coefficients $\unicode[STIX]{x1D6FC}_{m}^{(0)}$ , $\unicode[STIX]{x1D6FC}_{m}^{(1)}$ , $\unicode[STIX]{x1D6FD}_{m}^{(0)}$ and $\unicode[STIX]{x1D6FD}_{m}^{(1)}$ are given in table 1. For the particular case of an isothermal sphere, we simply compute the limit $\unicode[STIX]{x1D6EC}\rightarrow \infty$ in this expression, resulting in $\unicode[STIX]{x1D6F7}=-12\unicode[STIX]{x03C0}\sum \unicode[STIX]{x1D6FC}_{m}^{(1)}Kn^{m}/\sum \unicode[STIX]{x1D6FD}_{m}^{(1)}Kn^{m}$ , with index $m$ spanning the same ranges as in (4.3).

Table 1. Coefficients used in expression (4.3) for the thermophoretic force on a sphere modelled with R13.

Results from (4.3) are plotted as $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})$ versus $(\unicode[STIX]{x03C0}/2)^{1/2}Kn$ for $\unicode[STIX]{x1D6EC}=4$ , 10, and $22.4\times 10^{3}$ in figure 6. The highest $\unicode[STIX]{x1D6EC}$ value is chosen motivated by one of the experimental data sets depicted in the figure (see below). We also include predictions from various models based on the linearized Boltzmann equation, namely, by Sone & Aoki (Reference Sone and Aoki1983) for an isothermal sphere using the BGK equation – denoted in their work as the Boltzmann–Krook–Welander (BKW) equation – by Beresnev & Chernyak (Reference Beresnev and Chernyak1995) using the S model, and by Sone (Reference Sone2007) for a hard-sphere gas. Results from G13 and from its modification represented by Young’s (Reference Young2011) interpolation formula are also added. In addition, values from Waldmann (Reference Waldmann1959) formula, valid for large $Kn$ and insensitive to changes in $\unicode[STIX]{x1D6EC}$ , are presented. This formula as well as the expressions for the dimensionless thermophoretic force from G13 and Young (Reference Young2011) are presented, for completeness, in appendix D.

In figure 6 we included the experimental results on thermophoresis recently presented by Bosworth & Ketsdever (Reference Bosworth and Ketsdever2016) for spheres made of the polymer acrylonitrile butadiene styrene (ABS) and by Bosworth et al. (Reference Bosworth, Ventura, Ketsdever and Gimelshein2016) for spheres of copper. In the latter, data include observations of reversed or negative thermophoresis. In the experiments, they set the sphere in the mid-plane between two copper plates and placed the entire assembly in a vacuum chamber filled with argon. Rarefied conditions were attained by reducing the pressure to values significantly below the atmospheric pressure. The force was generated by keeping one plate at the ambient temperature while the other plate was heated uniformly to a higher temperature. To plot the data in figure 6 – originally presented as dimensional force versus pressure – we used an ambient temperature of 24.5 $^{\circ }$ C (Bosworth, R. & Ketsdever, A. 2017 Private communication.) and a temperature difference between plates of $35~~\text{K}~\text{m}^{-1}$ , a distance between plates of 0.40 m, a particle radius of 0.0254 m, a gas specific constant of $208~\text{J}~\text{kg}^{-1}\,\text{K}^{-1}$ and a gas viscosity of $2.295\times 10^{-5}~\text{Pa}~\text{s}$ at the mid-plane temperature (NIST 2017). The solid-to-gas heat conductivity ratios were estimated to be $\unicode[STIX]{x1D6EC}=10$ and $22.4\times 10^{3}$ for the ABS and copper spheres, respectively. In practice, in the copper case, such a large value would correspond to an isothermal sphere ( $\unicode[STIX]{x1D6EC}\rightarrow \infty$ ). To account for Knudsen number effects in the temperature gradient by the mid-plane between the plates, we adopted the model in chapter 4 of Sone (Reference Sone2007) based on the linearized Boltzmann equation for a hard-sphere gas. These gas rarefaction effects include a lower temperature gradient in the mid-plane with respect to the nominal value, gas–wall temperature jumps at the plates’ surfaces and noticeable deviations from linearity in the temperature profile in regions near the walls (Knudsen layers). These effects are evident, for instance, in figure 5 in Bosworth et al. (Reference Bosworth, Ventura, Ketsdever and Gimelshein2016), obtained using DSMC. This corrected value of the temperature gradient was used to compute the Epstein number $Ep$ for each experimental point. We obtain values of $Ep$ in the order of $3.0\times 10^{-3}$ for the experimental data sets.

Figure 6. Dimensionless thermophoretic force from various theories as a function of Knudsen number for a heat-conducting spherical particle in a gas with a constant temperature gradient in the far field. Plots correspond to particle-to-gas thermal conductivity ratios of (a) $\unicode[STIX]{x1D6EC}=4$ , (b) 10 and (c) $22.4\times 10^{3}$ . The experimental data of Bosworth & Ketsdever (Reference Bosworth and Ketsdever2016) for ABS spheres ( $\unicode[STIX]{x1D6EC}=10$ ) and of Bosworth et al. (Reference Bosworth, Ventura, Ketsdever and Gimelshein2016) for copper spheres ( $\unicode[STIX]{x1D6EC}=22.4\times 10^{3}$ ) are included in (b) and (c), respectively. The model by Sone & Aoki (Reference Sone and Aoki1983) assumes a uniform temperature in the entire sphere ( $\unicode[STIX]{x1D6EC}\rightarrow \infty$ ). Predictions from Waldmann’s (Reference Waldmann1959) model, valid for the free-molecule regime ( $Kn\gg 1$ ) and independent of $\unicode[STIX]{x1D6EC}$ , are also included. The thin-dashed line indicates the zero-force level.

For $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 0.1$ , figure 6 shows that the thermophoretic force from R13 tends very closely to the results from models of the linearized Boltzmann equation. In contrast, G13 significantly under-predicts these results. As demonstrated in the figure, this discrepancy is corrected in Young’s (Reference Young2011) interpolation formula, obtained from G13 after fitting the thermal creep, velocity slip and temperature jump coefficients from the so-called Maxwell–Smoluchowski values (e.g. see Kennard Reference Kennard1938; Nguyen & Wereley Reference Nguyen and Wereley2002; Gu & Emerson Reference Gu and Emerson2007; Young Reference Young2011) with those resulting from solutions of model Boltzmann equations compiled by Sharipov (Reference Sharipov2004). Young argues that changing these coefficients in the boundary conditions of G13 is needed to compensate for the inability of this model to reproduce the Knudsen layer. On the other hand, our analytical result from R13 gives the correct prediction without changing the original numerical factors appearing in the various terms of both the boundary conditions and bulk equations, i.e. with no fitting. Whereas the new set of coefficients adopted by Young (Reference Young2011) for G13 works well for the thermophoretic force on a sphere for very small $Kn$ , they might not be applicable to situations other than for which they were fitted (e.g. for different geometries).

When $Kn\rightarrow 0$ we have from (4.3) for R13 that

(4.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7} & \!=\! & \displaystyle -12\unicode[STIX]{x03C0}\!\left[\frac{1.13776419}{\unicode[STIX]{x1D6EC}+2}-1.03507976\frac{(\unicode[STIX]{x1D6EC}^{2}+8.34774566\unicode[STIX]{x1D6EC}+6.22407475\times 10^{-1})}{(\unicode[STIX]{x1D6EC}+2)^{2}}Kn\right]\nonumber\\ \displaystyle & & \displaystyle +\,O(Kn^{2}),\end{eqnarray}$$

whereas in the special case $\unicode[STIX]{x1D6EC}\rightarrow \infty$ (isothermal sphere), we find the following asymptotic result when $Kn\rightarrow 0$ ,

(4.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7} & = & \displaystyle 39.0215875\,Kn-1043.54657\,Kn^{2}+O(Kn^{3}).\end{eqnarray}$$

Note from these expressions that for very small $Kn$ and finite $\unicode[STIX]{x1D6EC}$ , $-\unicode[STIX]{x1D6F7}$ tends to a positive constant, indicating positive or normal thermophoresis. On the other hand, after passing the limit $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , the asymptotic result shows that $-\unicode[STIX]{x1D6F7}$ tends to zero linearly with $Kn$ through negative values corresponding to slightly reversed thermophoresis. A similar trend can be observed in the asymptotic formulae by Sone (Reference Sone2007).

Within the interval $0.02\lesssim \sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 0.2$ , figure 6(c) shows that the experimental data for the copper sphere exhibits negative or reversed thermophoresis, $-\unicode[STIX]{x1D6F7}<0$ . The reversed net force direction, such that it is now from the cold to the hot region for large values of $\unicode[STIX]{x1D6EC}$ , is mainly caused by the thermal-stress flow near the sphere’s surface. In this figure, R13, G13, Young’s formula and the models by Sone & Aoki (Reference Sone and Aoki1983) and Beresnev & Chernyak (Reference Beresnev and Chernyak1995) predict negative thermophoresis. In particular, the contour plots for $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ in figures 3 and 4 show the various features that characterized a scenario of negative thermophoresis, such as the reversed flow direction and the absence of temperature changes in the gas side along the solid surface, in comparison to the other cases considered. Nonetheless, whilst R13 agrees well with the benchmark Boltzmann solutions (for hard spheres and the S-model), both predict force magnitudes significantly smaller than the experimental values; a discrepancy which at present we cannot explain. Referring again to the values of $Kn$ and $\unicode[STIX]{x1D6EC}$ considered in figure 3, the force magnitude in the case $Kn=0.02$ and the largest $\unicode[STIX]{x1D6EC}$ , according to figure 6, is notably smaller than for the other cases, which seems to correspond with the fact that the velocity magnitude in this case is also one or more orders of magnitude smaller than in the others.

Expressions of thermophoretic force for R13 and G13, as well as Young’s interpolation formula exhibit a critical $\unicode[STIX]{x1D6EC}$ above which reversed thermophoresis occurs for some bounded $Kn$ interval. For R13, this threshold is $\unicode[STIX]{x1D6EC}=87.17$ , whereas for G13 it is a much lower ratio, perhaps less likely to actually occur, of 10.93. In the case of Young’s interpolation formula the threshold is $\unicode[STIX]{x1D6EC}=26.05$ . For larger, but finite values of $\unicode[STIX]{x1D6EC}$ , the corresponding curve of $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})$ will intercept the $\unicode[STIX]{x1D6F7}=0$ axis at two points, and in both $Kn>0$ . The maximum negative value of $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})$ occurs when $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . From the R13 model, this is $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})=-0.068$ at $Kn=0.024$ ; for G13 it is $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})=-0.321$ at $Kn=0.087$ , and for Young’s formula it is $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})=-0.180$ at $Kn=0.049$ . Negative thermophoresis cannot occur for $Kn\geqslant 0.054$ for R13, $Kn\geqslant 0.251$ for G13 and $Kn\geqslant 0.130$ for Young’s formula; these critical values are attained with $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . This upper limit of $Kn$ will decrease for $\unicode[STIX]{x1D6EC}<\infty$ . The critical magnitudes of Knudsen numbers listed in this paragraph were instrumental in choosing $Kn=0.02$ and 0.2 for the contour plots in figures 3 and 4.

For $0.2\lesssim \sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 1$ , R13 shows qualitatively the best performance among the models derived with the moments method for both the copper and polymer (ABS) data. In particular, for $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn~\lesssim ~0.4$ , its predictions agree with the experimental points for the copper’s sphere. It is likely coincidental that R13 performs better with respect to the experimental data than the benchmark models of Beresnev & Chernyak (Reference Beresnev and Chernyak1995) and Sone (Reference Sone2007) from kinetic theory. G13 depicts a rather unsatisfactory performance, as its maximum is significantly shifted to higher $Kn$ values. When $Kn>1$ , Waldmann’s formula matches the experiments, as expected. For larger $Kn$ values, the various model Boltzmann equations approximate well the experiments. G13 significantly over-predicts these values, although it tends to zero in the same asymptotic manner as Waldmann’s formula, that is, as $Kn^{-1}$ when $Kn\rightarrow \infty$ (see appendix D), but evidently with an incorrect coefficient. By construction, Young’s formula matches Waldmann’s for large $Kn$ . On the other hand, for R13, as $Kn\rightarrow \infty$ , we obtain from (4.3) that

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}=-3.80431806Kn^{-2}-\frac{5.46577139\,\unicode[STIX]{x1D6EC}-5.32067322\times 10^{-1}}{\unicode[STIX]{x1D6EC}\,Kn^{3}}+O(Kn^{-4}),\end{eqnarray}$$

and thus goes to zero faster than Waldmann’s.

For large $Kn$ ( $Kn>1$ , say), heat conductivity ratio $\unicode[STIX]{x1D6EC}$ plays a fairly inconsequential role. In this regime, according to figure 6, R13 under-predicts both the experimental force as well as Waldmann’s results. It is known that for such large values of $Kn$ , models such as G13 and R13 are not expected to give the correct quantitative result (Torrilhon Reference Torrilhon2010).

4.2 Uniform flow past a heat-conducting sphere

Profiles of the gas velocity components, density and temperature deviations as functions of the radial coordinate computed with the R13 exact solution of the previous section are presented in figure 7 for the case of a uniform flow past a stationary sphere with a thermal conductivity ratio $\unicode[STIX]{x1D6EC}\rightarrow \infty$ (isothermal sphere). For $Kn=0.05$ and 0.3, our solution is validated as its predictions match very well the exact solution of Torrilhon (Reference Torrilhon2010) obtained from his most recent Wolfram Mathematica code (see § 5.3 in Claydon et al. Reference Claydon, Shrestha, Rana, Sprittles and Lockerby2017). Good agreement is also observed, in general, between our solution and results from the asymptotic expressions of Sone (Reference Sone2007) when $Kn\rightarrow 0$ ( $k\rightarrow 0$ ) and his numerical solution of a kinetic theory model for $Kn=0.090$ ( $k=0.1$ ).

Figure 7. Profiles of (a) radial velocity, (b) polar velocity, (c) density deviation and (d) temperature deviation in the gas as functions of the radial coordinate for the problem of streaming flow past a sphere with a uniform temperature (i.e. $\unicode[STIX]{x1D6EC}\rightarrow \infty$ ). Results are from the R13 exact solutions obtained here and from Torrilhon (Reference Torrilhon2010), the numerical solution of Sone (Reference Sone2007) for $k=0.1$ , and his asymptotic expression for $k\rightarrow 0$ . Sone’s results are for a hard-sphere gas. Knudsen number $Kn$ is related to parameter $k$ as $Kn=\sqrt{2}\unicode[STIX]{x1D6FE}_{1}k/2$ , with $\unicode[STIX]{x1D6FE}_{1}=1.270042427$ from Sone (Reference Sone2007);  $k=0.055$ and 0.33 correspond to $Kn=0.05$ and 0.3, respectively.

The drag force acting on the sphere is computed from the same surface integral mentioned in the previous sub-section but using the pressure and deviatoric stress fields obtained from the solution of the problem of slow flow past a sphere. In dimensionless form, this drag force from R13 may be written as

(4.7) $$\begin{eqnarray}\frac{F_{D}^{\ast }}{F_{\text{Stokes}}^{\ast }}=\frac{2+\unicode[STIX]{x1D6EC}+\mathop{\sum }_{m=1}^{7}(\unicode[STIX]{x1D6FC}_{m}^{(0)}+\unicode[STIX]{x1D6FC}_{m}^{(1)}\unicode[STIX]{x1D6EC})Kn^{m}}{2+\unicode[STIX]{x1D6EC}+\mathop{\sum }_{m=1}^{9}(\unicode[STIX]{x1D6FD}_{m}^{(0)}+\unicode[STIX]{x1D6FD}_{m}^{(1)}\unicode[STIX]{x1D6EC})Kn^{m}},\end{eqnarray}$$

where $F_{\text{Stokes}}^{\ast }$ denotes Stokes formula for drag on a sphere, $F_{\text{Stokes}}^{\ast }=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{0}^{\ast }U_{0}^{\ast }a^{\ast }$ and coefficients $\unicode[STIX]{x1D6FC}_{m}^{(0)}$ , $\unicode[STIX]{x1D6FC}_{m}^{(1)}$ , $\unicode[STIX]{x1D6FD}_{m}^{(0)}$ and $\unicode[STIX]{x1D6FD}_{m}^{(1)}$ are given in Table 2. For the special case of the drag on an isothermal sphere due to a streaming flow, expression (4.7), in the limit of $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , reduces to $F_{D}^{\ast }/F_{\text{Stokes}}^{\ast }=(1+\sum \unicode[STIX]{x1D6FC}_{m}^{(1)}Kn^{\,m})/(1+\sum \unicode[STIX]{x1D6FD}_{m}^{(1)}Kn^{\,m})$ , where $m$ takes the same values as in (4.7). It should be noted that Torrilhon (Reference Torrilhon2010) showed the curve of drag versus Knudsen number resulting from his analysis, but did not present a closed-form expression for this force.

As $Kn\rightarrow 0$ , the drag can be computed from

(4.8) $$\begin{eqnarray}\frac{F_{D}^{\ast }}{F_{\text{Stokes}}^{\ast }}=1-1.30731306Kn+\frac{3.51286149+1.41122879\unicode[STIX]{x1D6EC}}{2+\unicode[STIX]{x1D6EC}}Kn^{2}+O(Kn^{3}).\end{eqnarray}$$

Furthermore, when $Kn\rightarrow \infty$ , we have from (4.7) the expression

(4.9) $$\begin{eqnarray}\displaystyle \frac{F_{D}^{\ast }}{F_{\text{S tokes }}^{\ast }} & = & \displaystyle 4.86486844\times 10^{-1}Kn^{-2}-2.50036697\times 10^{-1}Kn^{-3}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1.03874711\times 10^{-3}+6.79151242\times 10^{-2}\unicode[STIX]{x1D6EC}}{\unicode[STIX]{x1D6EC}Kn^{4}}+O(Kn^{-5}).\end{eqnarray}$$

Table 2. Coefficients used in expression (4.7) for the drag force on a sphere caused by a slow, streaming flow modelled with R13.

Figure 8 shows the dimensionless drag force as a function of $Kn$ for $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . We added results extracted from the work by Torrilhon (Reference Torrilhon2010) for an isothermal sphere. The curves match each other very well. This clearly indicates that the drag force is insensitive to changes in the solid-to-gas thermal conductivity ratio, a tendency already noted by Sone (Reference Sone2007) from his analysis based on the linearized Boltzmann equation in chapter 4 of his monograph. For small and large $Kn$ , this trend becomes evident in expressions (4.8) and (4.9), respectively. The matching with Torrilhon’s results serves as validation of the solution method adopted in our analysis. In addition, we included predictions from Young’s (Reference Young2011) solution for an isothermal sphere of the G13 equations using Maxwell–Smoluchowski values for the thermal creep, velocity slip and temperature jump coefficients, as well as with the coefficients extracted by Young from the work of Sharipov (Reference Sharipov2004) (see appendix D). For the interval shown, no significant difference is noted between these two sets. We included also predictions from kinetic theory by Sone (Reference Sone2007) for a hard-sphere gas and an isothermal sphere. Whereas G13 over-predicts the results by Sone, the R13 tendency is to under-predict them, although with a smaller difference. For $Kn\lesssim 0.4$ , R13 approximates Sone’s curve very well. For each of the models, as in the case of R13, changing $\unicode[STIX]{x1D6EC}$ from infinity to 4 produced no significant changes. Figure 8 also includes experimental data by Goldberg (Reference Goldberg1954) (extracted from Torrilhon, Reference Torrilhon2010) and by Allen & Raabe (Reference Allen and Raabe1982) (extracted from Claydon et al., Reference Claydon, Shrestha, Rana, Sprittles and Lockerby2017). Predictions from both R13 and Sone (Reference Sone2007) under-predict the drag, with the former resulting in the largest difference.

Figure 8. Drag force acting on a sphere due to a streaming flow non-dimensionalized with Stokes’ drag versus Knudsen number resulting from the R13 moment equations for particle-to-gas heat conductivity ratios $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . Results from Torrilhon (Reference Torrilhon2010) with R13; from Young (Reference Young2011) with G13 using (a) Maxwell–Smoluchowski’s set of coefficients and (b) an alternative set (see appendix D); from Sone (Reference Sone2007) using a model Boltzmann equation, and from experiments (Goldberg Reference Goldberg1954; Allen & Raabe Reference Allen and Raabe1982) are also included.

4.3 Thermophoretic velocity

The thermophoretic velocity is the velocity of a moving particle driven by the thermophoretic force, when this force is balanced by the drag; that is, when the net force on the particle vanishes. Setting $F_{T}^{\ast }+F_{D}^{\ast }=0$ , and using expression (4.2) for $F_{T}^{\ast }$ , we readily have for the dimensionless velocity

(4.10) $$\begin{eqnarray}Ma=-\frac{\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x1D6F9}}Ep\,Kn,\end{eqnarray}$$

expressed as the pseudo Mach number. Previously, $Ma$ was a prescribed parameter and now it is the unknown. Here, $\unicode[STIX]{x1D6F9}$ , a dimensionless quantity introduced by Young (Reference Young2011), is equal to $6\unicode[STIX]{x03C0}F_{D}^{\ast }/F_{\text{Stokes}}^{\ast }$ . In general, it is a function of $Kn$ and $\unicode[STIX]{x1D6EC}$ . The expression for $\unicode[STIX]{x1D6F9}$ from the G13 theory is given in appendix D.

With (4.3) and (4.7) for the non-dimensional thermophoretic force and drag, respectively, relation (4.10) provides an expression for the thermophoretic velocity according to R13. Figure 9 shows results from this expression normalized with factor $\sqrt{\unicode[STIX]{x03C0}/2}\,Ep\,Kn$ as a function of $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn$ for different values of $\unicode[STIX]{x1D6EC}$ . This figure includes predictions from G13 and from solutions of model Boltzmann equations (Beresnev & Chernyak Reference Beresnev and Chernyak1995; Sone Reference Sone2007), and from the combination in (4.10) of Waldmann’s formula for $\unicode[STIX]{x1D6F7}$ and Epstein’s formula for $\unicode[STIX]{x1D6F9}$ intended for the free-molecule regime of large $Kn$ (see appendix D). The velocity from this model is independent of $\unicode[STIX]{x1D6EC}$ .

Figure 9. Dimensionless thermophoretic velocity $Ma$ divided by product $\sqrt{\unicode[STIX]{x03C0}/2}KnEp$ as a function of Knudsen number from the R13 moment equations and from other models in the literature. The thin-dashed line indicates the zero-velocity level. Except for Waldmann–Epstein’s formula, the results are obtained for solid-to-gas heat conductivity ratios $\unicode[STIX]{x1D6EC}=4$ , 10 and $22.4\times 10^{3}$ .

From figure 9 we note that R13 predicts the correct values given by the model Boltzmann equations for $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 0.2$ . For greater values of $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn$ , the agreement deteriorates. For $Kn\rightarrow \infty$ , R13 predicts a constant value, as is the case with the Waldmann–Epstein formula and with the results from kinetic theory, but with a smaller magnitude. G13, on the other hand, although also tending to a constant value, shows significant discrepancies with the other theories considered. Even though R13 is outside its limits of applicability for $Kn>1$ , it provides a reasonable engineering tool in that it has the correct asymptotic behaviour and, quantitatively, its differences with the benchmark models from kinetic theory and Waldmann–Epstein’s formula are not significantly large.