Thermophoresis of a spherical particle: modelling through moment-based, macroscopic transport equations

Abstract

We consider the linearized form of the regularized 13-moment equations (R13) to model the slow, steady gas dynamics surrounding a rigid, heat-conducting sphere when a uniform temperature gradient is imposed far from the sphere and the gas is in a state of rarefaction. Under these conditions, the phenomenon of thermophoresis, characterized by forces on the solid surfaces, occurs. The R13 equations, derived from the Boltzmann equation using the moment method, provide closure to the mass, momentum and energy conservation laws in the form of constitutive, transport equations for the stress and heat flux that extend the Navier–Stokes–Fourier model to include non-equilibrium effects. We obtain analytical solutions for the field variables that characterize the gas dynamics and a closed-form expression for the thermophoretic force on the sphere. We also consider the slow, streaming flow of gas past a sphere using the same model resulting in a drag force on the body. The thermophoretic velocity of the sphere is then determined from the balance between thermophoretic force and drag. The thermophoretic force is compared with predictions from other theories, including Grad’s 13-moment equations (G13), variants of the Boltzmann equation commonly used in kinetic theory, and with recently published experimental data. The new results from R13 agree well with results from kinetic theory up to a Knudsen number (based on the sphere’s radius) of approximately 0.1 for the values of solid-to-gas heat conductivity ratios considered. However, in this range of Knudsen numbers, where for a very high thermal conductivity of the solid the experiments show reversed thermophoretic forces, the R13 solution, which does result in a reversal of the force, as well as the other theories predict significantly smaller forces than the experimental values. For Knudsen numbers between 0.1 and 1 approximately, the R13 model of thermophoretic force qualitatively shows the trend exhibited by the measurements and, among the various models considered, results in the least discrepancy.

1 Introduction

Thermophoresis refers to the force and, potentially, motion experienced by solid particles or surfaces exposed to gas under rarefied conditions and in the presence of a temperature gradient. This phenomenon seems to have been first noted by Tyndall (1870) while observing the spatial redistribution of ambient dust in the proximity of a heated surface. Commonly, the thermophoretic force has been assumed to point from the hot to the cold region, that is, opposite to the temperature gradient. This condition has been labelled as ‘positive’ thermophoresis. On the other hand, a reversal of the force direction, known as ‘negative’ thermophoresis, has also been predicted and – although rarely – observed, as discussed below. Thermophoresis belongs to a class of phenomena promoted in a gas away from thermodynamic equilibrium, which occurs when the collisions between gas molecules are insufficient. Among the various effects in this class are the velocity slip and temperature jump at solid walls and liquid–gas interfaces, transpiration flow, thermal stress, Knudsen layers and heat flux without temperature gradients (Sone 2007; Struchtrup & Torrilhon 2008). Review papers (Talbot et al. 1980; Bakanov 1991; Zheng 2002) and book sections (Sone 2007) have been written on the subject of thermophoresis. In particular, Young (2011) presented a comprehensive examination of the various theories on particle thermophoresis at arbitrary Knudsen numbers under the light of experimental data, and concluded that the accuracy of the measurements and the interval of Knudsen numbers explored are such that confirmation of the validity of the theories could not be achieved.

The main measure of departure from local thermodynamic equilibrium in a gas is the Knudsen number, defined as the ratio of the mean free path between collisions and a characteristic macroscopic length scale. Typically, it is accepted that for Knudsen numbers below 0.01 the Navier–Stokes equations and the Fourier law for conductive heat flux are reliable. Beyond this threshold, because departure from equilibrium may be significant, predictions from these models become doubtful and models from kinetic theory, based on the Boltzmann equation for the molecular velocity distribution function, take precedence. For Knudsen numbers of order unity or higher, solution of the Boltzmann equation either directly or by means of stochastic techniques such as the direct simulation Monte Carlo method (DSMC) of Bird (1994) are computationally tractable, in general. On the other hand, in the interval of Knudsen numbers between, roughly, 0.01 and 1 such computations become increasingly expensive (Torrilhon 2016).

Since for gas flows confined in micro- or nano-devices the Knudsen number may lie in the transition regime or beyond, it is important to resort to models capable of describing, at least qualitatively, rarefaction effects when designing or analysing such systems. Other examples include the transport of particles or droplets by a gaseous stream, or of bubbles by a liquid, when the size of these solid or fluid objects lies in the micrometre or nanometre range. Although the Navier–Stokes–Fourier equations of classical hydrodynamics can be extended into the transition regime for Knudsen numbers beyond 0.01 by including slip and jump effects in the boundary conditions, not all rarefaction effects occurring in the bulk can be captured by this approximation (Mohammadzadeh et al. 2015; Torrilhon 2016). Another aspect relevant to modelling micro- or nano-flows is that, for such small scales, the flow Reynolds number is typically smaller than one, so that neglecting inertia is an admissible assumption.

Efforts to model thermophoresis or related phenomena can be traced back for more than a century. The first theory on gas motion induced by temperature gradients was due to Maxwell (1879). His analysis stems from kinetic theory. Later, Epstein (1929) presented a theory for the thermophoretic force on a spherical particle for small Knudsen numbers based on the continuum approach that takes into account thermal slip as well as the particle’s thermal conductivity. Waldmann (1959) derived a model valid in the free-molecule regime (large Knudsen numbers) that has been widely applied. Brock (1962) improved Epstein’s theory to develop an expression for the thermophoretic force for Knudsen numbers ${\lesssim}$ 0.1. Talbot et al. (1980) proposed a correlation for the transition regime by modifying the coefficients in Brock’s expression such that Waldmann’s free-molecule regime expression was approached at large Knudsen numbers.

A commonly used route to study gas flow in the transition regime from first principles, since numerical solutions of the Boltzmann equation can be computationally very costly, has been the analysis of linearized versions of the Boltzmann equation with simpler models for the collision term. Some of the simpler models applied to the problem of thermophoresis of a spherical particle include the Bhatnager, Gross, and Krook model (BGK) (Sone & Aoki 1983; Yamamoto & Ishihara 1988; Takata et al. 1993; Sone 2007) the S model (Beresnev & Chernyak 1995) and the hard-sphere model (Sone 2007).

An alternative approach to model rarefied gas flow is the use of macroscopic transport equations derived using moment methods. In the original moment method, introduced by Grad (1949), the distribution function in the Boltzmann equation was expanded in Hermite polynomials and the macroscopic variables describing the flow were represented as moments (integrals) of this distribution. For the first approximation beyond the Navier–Stokes–Fourier equations, in total, thirteen moments are needed for the same number of fields, namely, mass density, macroscopic velocity vector, temperature, heat-flux vector and deviatoric stress tensor (symmetric and trace free), yielding Grad’s 13-moment equations (G13). The pressure appears through an equation of state, typically, the ideal gas law. Dwyer (1967) presented a theory for the thermophoretic force on a spherical particle based on Grad’s moment method predicting, for the first time, reversed thermophoresis. Recently, Young (2011) noted that Dwyer did not account for the totality of the stress and heat-flux coupling terms in the temperature jump boundary condition. Young corrected this and rederived the expression for thermophoretic force from the G13 equations. He then modified this result to present an interpolation formula that fits Waldmann’s expression for large Knudsen numbers and, by introducing values for the thermal creep, velocity slip, and temperature jump coefficients cited by Sharipov (2004) based on solutions of model Boltzmann equations, also matches results from kinetic theory for small Knudsen numbers.

A notorious deficiency of the G13 equations is their inability to describe Knudsen layers, that is, regions adjacent to solid surfaces where rarefaction effects are conspicuous. Struchtrup & Torrilhon (2003) regularized Grad’s 13-moment equations (R13) by adding second-order derivatives to the closures. Starting with the Boltzmann equation, R13 equations are best derived using the order of magnitude method (Struchtrup 2005a ; Struchtrup et al. 2017). The R13 equations are equipped with a set of boundary conditions which are valid at walls or gas–liquid interfaces without mass transfer (Torrilhon & Struchtrup 2008) or at evaporating and condensing interfaces (Struchtrup & Frezzotti 2016; Struchtrup et al. 2017). In contrast to G13, these equations can partially capture the structure and effects of the Knudsen layer. Because they involve the typical variables describing fluid flow and heat transfer at the macroscopic level, interpretation of results may be facilitated by inspecting specific terms in the differential equations, contributing to a good physical understanding. In addition, the extension of tested, well-established numerical techniques in computational fluid dynamics to this system of equations may be achieved straightforwardly. Also, because the boundary conditions have been derived for sharp surfaces of discontinuity, R13 equations can be applied to two-phase systems involving a rarefied gas and a liquid sharing an interface whose instantaneous position is another unknown in the problem. With R13, the limit of continuum models to give meaningful results when rarefaction effects are important has been pushed to a Knudsen number of approximately 0.5 in the transition regime (Torrilhon 2016; Struchtrup et al. 2017). Early advances in the development of moment methods, with emphasis on R13, can be found in the textbook by Struchtrup (2005b ); more recent developments and applications of R13 have been compiled in the reviews by Struchtrup & Taheri (2011) and Torrilhon (2016).

To the best of our knowledge, application of R13 to investigate transport phenomena involving spherical or near spherical particles in rarefied environments is limited to the analytical work of Torrilhon (2010) on the slow flow of gas past a sphere, and to the numerical treatment of the same problem by Claydon et al. (2017) using a mesh-free method. The application of R13 to model thermophoresis on a spherical particle has not yet been pursued. After recognizing the complexity of the R13 equations in comparison to G13, Young (2011) carried out the modelling of this problem with the latter. In the conclusions of his article, he recommends ‘solving the R13-moment equations in order to study reversed thermophoresis in greater detail’.

The aim of this work is to obtain the thermophoretic force acting on a sphere surrounded by a rarefied gas exposed to a uniform temperature gradient far from the sphere using the R13 equations and taking into account heat conductivity inside the sphere. In addition, we compute the thermophoretic velocity of the sphere when this is free to move under the thermophoretic force. This velocity corresponds to the balance between the thermophoretic force and the drag caused on the sphere in its motion by the surrounding gas. To model the drag, we expand the scope in Torrilhon’s (2010) work with R13 to include the thermal conductivity of the solid (i.e. a non-isothermal sphere), even though the drag has been shown to be fairly insensitive to changes in this parameter (Sone 2007). Instead of using the form of the solution for classical Stokes flow past a sphere to write the ansatz for the equations (cf. Torrilhon 2010), we apply a somewhat more general approach, namely, the method of multipole potentials. We present closed-form expressions for the thermophoretic force and the drag resulting from R13. We compare predictions for the thermophoretic force from this model with results from simplified models of the Boltzmann equation used in kinetic theory and from other systems of moment equations. We include in the comparison the new experimental data by Bosworth et al. (2016) for both positive and negative thermophoresis.

The content of this paper is organized as follows. In the next section, we separately formulate the problems of thermophoresis of a spherical particle and of slow flow of a rarefied gas past a sphere and introduce the main tool of analysis, the regularized 13-moment equations, or R13, and their boundary conditions for the gas–solid interface in linearized form. We then rewrite the system of equations in a different form, more convenient for the solution method of the following section. In § 3, we describe the solution of the system of equations; the procedure involves the method of multipole potentials. Next, § 4 begins by discussing results for the problem of thermophoresis on a sphere, including spatial profiles for the macroscopic field variables, contour plots and streamline patterns and the thermophoretic force, exploring the effect of the solid-to-gas thermal conductivity ratio. For the force, we compare results from R13 with recent experimental data showing reversed thermophoresis and with other models from the literature. Then, we present results for the drag force arising from the gas flow past a sphere from R13 considering the sphere’s heat conductivity and compare with theoretical predictions from the literature, including Torrilhon’s (2010) R13 results in the case of an isothermal particle, which serves as validation of the solution method implemented here. The section closes presenting results for the thermophoretic velocity from R13 and other models. Finally, § 5 contains some concluding remarks.

Figure 1. Sketch for the problem of thermophoresis on a sphere ( $G^{\ast }=\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast }|_{\infty }$ , the far-field temperature gradient) or uniform flow past a sphere ( $G^{\ast }=U_{0}^{\ast }$ , the far-field gas velocity). The sphere’s radius is denoted by $a^{\ast }$ . In both cases, the flow field is axisymmetric with respect to the $z^{\ast }$ axis. Unit vector $\text{k}$ points in the direction of the positive $z^{\ast }$ semi axis. The spherical coordinate system $\{r^{\ast },\unicode[STIX]{x1D717},\unicode[STIX]{x1D719}\}$ , with origin at the sphere’s centre, is depicted.

2 Problem formulation

In this section we formulate mathematically two problems involving a sphere in a rarefied gas. First, we consider the problem of thermophoresis of a spherical particle with the gas far from the sphere at rest with a uniform temperature gradient. Second, we detail the problem of a uniform flow past a sphere with no temperature gradient imposed in the far field. The results of these two problems will be combined later to obtain the thermophoretic velocity of a sphere.

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

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

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 (2011) after P. S. Epstein, a pioneer in the study of thermophoresis of spherical particles, who presented the first theory on the subject (Epstein 1929). 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}$$

whilst the R13 constitutive equations for the deviatoric stress and heat flux are given by (Struchtrup 2005b ; Lockerby & Collyer 2016; Torrilhon 2016)

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

The overbar in expressions (2.6) denotes a symmetric and trace-free tensor; $\unicode[STIX]{x1D6E5}$ denotes the Laplacian operator. The temperature in the spherical particle satisfies the classical steady heat equation

(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. (2017) in the absence of phase change (see also Struchtrup & Frezzotti 2016) 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. (2017) are, in the present work, associated with the solid (sphere). Note that in this case, boundary conditions (37)–(41) in Struchtrup et al. (2017) reduce to the wall boundary conditions (33)–(37) in Torrilhon & Struchtrup (2008) 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 2011)

(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. (2017) – 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}$$

The full nonlinear expressions for the heat flux $\text{q}$ , deviatoric stress $\unicode[STIX]{x1D748}$ and higher-order moments $\text{R}$ and $\text{m}$ that lead to (2.6) and (2.16) can be found in Struchtrup (2005b , see chapters 7 and 9) – there, a third (scalar) higher-order moment is also present that contributes nothing to the linear equations when $\text{q}$ is divergence free. Once these equations are linearized, relations (2.16) are then used to eliminate $\text{R}$ and $\text{m}$ resulting in (2.6).

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

where

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

coupled with the set of equations

(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 Solution applying the method of multipole potentials

In this section, we seek solutions for the flow disturbances in the gas that vanish in the far field alongside temperature profiles inside the sphere which are coupled through the set of conditions prescribed at the sphere’s surface.

Consider partial differential equations (2.26)–(2.28) and (2.23)–(2.25) for the gas, and (2.7) for the temperature in the solid particle. Based on this, scalar fields $\unicode[STIX]{x1D6F1}$ , $\unicode[STIX]{x1D701}$ and $\unicode[STIX]{x1D703}_{s}$ are harmonic functions. For fields $\unicode[STIX]{x1D755}$ , $p$ , $\text{q}$ , $\unicode[STIX]{x1D748}$ , their solutions will be written as the sum of the solution of the homogeneous equation plus a particular solution of (2.28) and (2.23)–(2.25), respectively (e.g. see Lamb 1932; Leal 2007). The solutions of the homogeneous equations will be found with the method of multipole potentials. After these fields have been determined, the velocity and temperature fields $\text{u}$ and $\unicode[STIX]{x1D703}$ can be obtained from expressions (2.19) and (2.20), respectively.

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

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

After adding and subtracting $\unicode[STIX]{x1D706}_{1}^{2}\unicode[STIX]{x1D702}$ to the left-hand side of (3.7b ), we seek solutions of (3.7) of the form

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

Substitution of (3.8) into (3.7) results in $m=-Kn^{-1}/\unicode[STIX]{x1D706}_{3}^{4}$ and $n=-(\unicode[STIX]{x1D706}_{1}/\unicode[STIX]{x1D706}_{3})^{2}/(\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2})$ , so that

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

and $\unicode[STIX]{x1D748}^{(p)}$ is given by (3.6) with (3.9); it is symmetric and trace free.

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 (2007, chap. 8) as well as thoroughly discussed and used in a wider range of problems of mathematical physics in Hess (2015, 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. 2012, 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 (2007) and Hess (2015).

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

so that $f^{\prime }(r)\equiv 0$ and $g^{\prime }(r)=rf(r)/18$ . Substitution of expressions (3.13)–(3.15) for $p$ , $\text{q}$ and $\unicode[STIX]{x1D748}$ into (2.23)–(2.25), respectively, and invocation of property (A 7) from appendix A, lead to the following ordinary differential equations

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

Using $f^{\prime }(r)\equiv 0$ in (3.20d ) yields $f(r)\equiv 0$ ; hence, from above, we have that $g^{\prime }(r)\equiv 0$ and (3.20e ) results in $g(r)\equiv 0$ . We thus find that $\unicode[STIX]{x1D748}^{(h)}(\text{x})\equiv 0$ . Therefore, $\unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}^{(p)}$ , given in (3.6) which, in turn, makes use of expressions (3.10), (3.12) and (3.13). The solution of (3.20b ) is well known. Obtaining exact solutions for (3.20a ) and (3.20c ) is discussed in appendix B. These solutions can be written as

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

where only the solution that decays for large $r$ has been retained. Substitution of (3.21b ) and (3.21c ) into the divergence-free condition for the heat flux (3.18) results in

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

From the solution for $\unicode[STIX]{x1D748}$ , we find that $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}\unicode[STIX]{x1D719}}=-\unicode[STIX]{x1D70E}_{rr}/2$ , a result already obtained by Torrilhon (2010) from an entirely different route. We also find that $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}r}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D719}}=0$ , because of the problem symmetry.

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 (2010). 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}$$

Torrilhon (2010) also gives formulae for components $\text{R}_{rr}$ , $\text{R}_{\unicode[STIX]{x1D717}r}$ , $\text{m}_{rrr}$ and $\text{m}_{\unicode[STIX]{x1D717}rr}$ directly in terms of derivatives of the components of the heat flux and stress, which can be used to derive expressions (3.24).

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 Results and discussion

This section contains three parts. First, we discuss predictions from our solution of the R13 moment equations for the problem of thermophoresis of a spherical particle presented in § 2.1. In particular, we show flow and temperature profiles for the sphere’s surroundings contrasted with results from kinetic theory and then compare predictions for the thermophoretic force from our solution to those from other theories as well as recently published experimental measurements. Second, we present results from R13 and other theories for the problem described in § 2.2, that is, the drag force due to a uniform flow past a sphere taking into account the particle-to-gas heat conductivity ratio. In both sub-sections, we present profiles of velocity components, density and temperature in a neighbourhood of the spherical surface. Third, the solutions from these two problems are coupled in the balance of the thermophoretic force by the drag acting on the sphere when this is free to move giving rise to the sphere’s thermophoretic terminal velocity. We show predictions from R13 and other theories for this quantity.

Note that the results presented here for the temperature deviations in the gas and solid, $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D703}_{s}$ and the gas heat flux, $\text{q}$ , correspond to the variables on the left-hand sides of expressions in (2.8). These are obtained after transforming back from the fields given in (3.23), which represent the variables with ‘  $\check{\,}$  ’.

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

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 (2007), 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 1879; Kennard 1938; Sone 2007; Mohammadzadeh et al. 2015). 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 2007; Young 2011), 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 (2007) 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 (2007), 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 (2007) for the thermal-stress slip flow and the illustration in figure 5. From a microscopic perspective, Sone (2007) 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 (2011) 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 (1983) for an isothermal sphere using the BGK equation – denoted in their work as the Boltzmann–Krook–Welander (BKW) equation – by Beresnev & Chernyak (1995) using the S model, and by Sone (2007) for a hard-sphere gas. Results from G13 and from its modification represented by Young’s (2011) interpolation formula are also added. In addition, values from Waldmann (1959) 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 (2011) are presented, for completeness, in appendix D.

In figure 6 we included the experimental results on thermophoresis recently presented by Bosworth & Ketsdever (2016) for spheres made of the polymer acrylonitrile butadiene styrene (ABS) and by Bosworth et al. (2016) 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 (2007) 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. (2016), 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 (2016) for ABS spheres ( $\unicode[STIX]{x1D6EC}=10$ ) and of Bosworth et al. (2016) for copper spheres ( $\unicode[STIX]{x1D6EC}=22.4\times 10^{3}$ ) are included in (b) and (c), respectively. The model by Sone & Aoki (1983) assumes a uniform temperature in the entire sphere ( $\unicode[STIX]{x1D6EC}\rightarrow \infty$ ). Predictions from Waldmann’s (1959) 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 (2011) 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 1938; Nguyen & Wereley 2002; Gu & Emerson 2007; Young 2011) with those resulting from solutions of model Boltzmann equations compiled by Sharipov (2004). 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 (2011) 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 (2007).

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 (1983) and Beresnev & Chernyak (1995) 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 (1995) and Sone (2007) 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 2010).

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 (2010) obtained from his most recent Wolfram Mathematica code (see § 5.3 in Claydon et al. 2017). Good agreement is also observed, in general, between our solution and results from the asymptotic expressions of Sone (2007) 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 (2010), the numerical solution of Sone (2007) 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 (2007);  $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 (2010) 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 (2010) 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 (2007) 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 (2011) 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 (2004) (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 (2007) 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 (1954) (extracted from Torrilhon, 2010) and by Allen & Raabe (1982) (extracted from Claydon et al., 2017). Predictions from both R13 and Sone (2007) 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 (2010) with R13; from Young (2011) with G13 using (a) Maxwell–Smoluchowski’s set of coefficients and (b) an alternative set (see appendix D); from Sone (2007) using a model Boltzmann equation, and from experiments (Goldberg 1954; Allen & Raabe 1982) 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 (2011), 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 1995; Sone 2007), 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.

5 Concluding remarks

In this work, we investigated theoretically an instance of the phenomenon of thermophoresis, a term that refers to the forces on and motions of objects caused by temperature gradients when these objects are exposed to rarefied gases. In particular, we considered the problem of thermophoresis of a spherical particle, obtaining an analytical solution by solving the R13 moment equations, a model that provides a macroscopic description of rarefied gas flows up to the transition regime for Knudsen numbers smaller than one. Besides writing expressions for the field variables describing the fluid flow and heat transfer problem, by integration of the total stress on the surface of the sphere, we obtained a closed-form expression for the thermophoretic force as a function of the Knudsen number, dimensionless temperature gradient (Epstein number), and solid-to-gas heat conductivity ratio.

Employing the closed-form expressions for the field variables in the gas obtained here, we plotted their profiles in a neighbourhood of the sphere alongside predictions from a model Boltzmann equation from the literature. This comparison revealed that up to a Knudsen number of about $0.1$ , based on the sphere’s radius, the agreement between these solutions is very good thereby showing that the solution based on R13 is capable of fully describing the Knudsen layer for this particular physical situation in this range of Knudsen numbers. We also showed contour plots of speed and temperature including velocity and heat flux streamlines for various conditions. These figures exhibit the interplay between the various mechanisms involved in the complex process of thermophoresis, namely, the more common gas thermal creep toward the hot region caused by a temperature gradient in the gas parallel to the solid surface and, under very specific conditions in the slip regime that include a highly thermally conductive solid, the reversal of the flow direction resulting from the now dominant mechanism of slip flow driven by thermal stresses.

We extended the same modelling approach with the R13 equations to obtain an analytical solution for the drag by considering the problem of a heat-conducting spherical particle in a uniform, slow gas stream, thereby extending the analysis by Torrilhon (2010) for an isothermal sphere. Predictions from the new solution are insensitive to changes in the solid-to-gas thermal conductivity ratio and agree very well with the results of Torrilhon, who used a different method to obtain his solution. For Knudsen numbers smaller than $0.4$ , R13 approximates reasonably well kinetic theory results for hard spheres. We then computed the thermophoretic velocity of the spherical particle when this can move driven by the thermophoretic force. The expression for the velocity results from balancing the thermophoretic force with the drag resistance exerted by the surrounding gas on the sphere. For Knudsen numbers smaller than $0.2$ , approximately, values for the thermophoretic velocity from R13 show good agreement with those from two models of the Boltzmann equation.

Results from the new thermophoretic force model derived with the R13 equations were compared with results from G13 and from various models based on the linearized Boltzmann equation, such as BGK and the S model, and with data from very recent experiments, for a wide range of Knudsen numbers. For Knudsen numbers below approximately $0.1$ , R13 results match predictions from the model Boltzmann equations. In this interval, the various theories considered, including R13, under-predict significantly the experimental measurements of reversed thermophoretic force for a metallic sphere (highly thermally conductive). For Knudsen numbers lying between $0.1$ and $1$ , approximately, the graphs from R13 follow qualitatively the experimental curves for both low and high solid-to-gas thermal conductivity ratios, and in general its predictions show, although important, the least differences with the experimental data. Surprisingly, solutions of models of the Boltzmann equation taken as benchmark demonstrate larger discrepancies with the experiments.

There are more involved macroscopic models based on the moments method, such as the R26 equations of Gu & Emerson (2009), that exhibit higher accuracy in the transition regime. Nevertheless, modelling phenomena in rarefied gas dynamics with the R13 equations may provide perhaps the best compromise between a sufficiently complex mathematical description capable of capturing the most significant features not only in the bulk but also near the boundaries (Knudsen layer) and an analytically tractable set of equations. In any case, the exact solution presented here can be useful in the step of validation of numerical tools developed for the R13 equations simulating flow and heat transfer phenomena occurring either in the interior of or external to complicated geometries.