2.1 Force on the sphere

The scaled force, $\tilde{F}$ , exerted on the sphere by the gas in the direction of motion, $\hat{\boldsymbol{z}}$ , is

(2.22) $$\begin{eqnarray}\displaystyle \tilde{F}=\int _{0}^{2{\rm\pi}}\int _{0}^{{\rm\pi}}\left(\tilde{T}_{rr}\cos {\it\theta}-\tilde{T}_{r{\it\theta}}\sin {\it\theta}\right)\sin {\it\theta}\,\text{d}{\it\theta}\,\text{d}{\it\phi}|_{r=1}, & & \displaystyle\end{eqnarray}$$

where the force scale, $2a^{2}U_{0}P/c_{m}$ , is used and $P$ is the ambient gas pressure. Importantly, (2.22) is valid for all Knudsen numbers, $Kn$ , and frequency ratios, ${\it\Theta}$ . The functions $\tilde{T}_{rr}$ and $\tilde{T}_{r{\it\theta}}$ are spatial components of the normal and shear stresses at the surface of the sphere, scaled by $2U_{0}P/c_{m}$ , and defined

(2.23a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{T}_{rr}=\frac{1}{{\rm\pi}^{3/2}}\int _{-\infty }^{\infty }c_{r}^{2}\exp (-|\boldsymbol{c}|^{2})\tilde{h}(r=1,\boldsymbol{c})\,\text{d}\boldsymbol{c}, & \displaystyle\end{eqnarray}$$

(2.23b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{T}_{r{\it\theta}}=\frac{1}{{\rm\pi}^{3/2}}\int _{-\infty }^{\infty }c_{r}c_{{\it\theta}}\exp (-|\boldsymbol{c}|^{2})\tilde{h}(r=1,\boldsymbol{c})\,\text{d}\boldsymbol{c}, & \displaystyle\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\tilde{h}(r=1,\boldsymbol{c})=\left\{\begin{array}{@{}ll@{}}\displaystyle 2\boldsymbol{c}\boldsymbol{\cdot }\tilde{\boldsymbol{U}}_{sp}+\sqrt{{\rm\pi}}\tilde{\boldsymbol{U}}_{sp}\boldsymbol{\cdot }\hat{\boldsymbol{r}}-\frac{2}{{\rm\pi}}\int _{c_{r}<0}\exp (-|\boldsymbol{c}|^{2})c_{r}\tilde{h}(r=1,\boldsymbol{c})\,\text{d}\boldsymbol{c},\quad & c_{r}>0,\\ \displaystyle \int _{0}^{\infty }\frac{1}{VKn}\left(\tilde{w}+2c_{r^{\prime }}\tilde{U} _{r^{\prime }}+2c_{{\it\theta}^{\prime }}\tilde{U} _{{\it\theta}^{\prime }}+\left[V^{2}-\frac{3}{2}\right]\tilde{{\it\tau}}\right)\quad \\ \displaystyle \quad \times \exp \left(-\frac{1-\text{i}{\it\Theta}}{VKn}|\boldsymbol{x}-\boldsymbol{x}^{\prime }|\right)\,\text{d}|\boldsymbol{x}-\boldsymbol{x}^{\prime }|,\quad & c_{r}<0.\end{array}\right.\end{eqnarray}$$

Equations (2.22)–(2.24) can be easily manipulated to yield a simplified expression for the scaled force on the sphere,

(2.25) $$\begin{eqnarray}\displaystyle \tilde{F} & = & \displaystyle \frac{2\sqrt{{\rm\pi}}}{3}\int _{1}^{\infty }\frac{1}{Kn}\sqrt{{\rm\pi}}r^{\prime 2}\left[W_{1}(r^{\prime })q_{1}(r^{\prime })-2W_{2}(r^{\prime })q_{2}(r^{\prime })\vphantom{\frac{3}{2}W_{4}(r^{\prime })}\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,4W_{3}(r^{\prime })q_{3}(r^{\prime })+\frac{3}{2}W_{4}(r^{\prime })\right]\,\text{d}r^{\prime }+\frac{{\rm\pi}}{3}\left(g+\frac{8}{\sqrt{{\rm\pi}}}\right),\end{eqnarray}$$

where the functions $W_{i}$ for $i=\{1,2,3,4\}$ are given in appendix B and the constant $g$ is given in (A 15).

3.1 Unsteady flow

Analytical solutions are now presented in the continuum, slip and free-molecular asymptotic limits for unsteady flow. These are used in § 5 to provide benchmark results for the numerical solution outlined in § 4.

3.1.1 Continuum flow $(Kn\rightarrow 0~\mathit{and}~{\it\Theta}\rightarrow 0)$

In the continuum limit, the linearised Navier–Stokes equations, together with the no-slip boundary conditions, yield the radial and angular components of the bulk velocity (Landau & Lifshitz Reference Landau and Lifshitz1987):

(3.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{U} _{r}(r,{\it\theta})=\left\{\frac{2\hat{A}}{r^{3}}+2\hat{B}\exp (\text{i}\sqrt{2\text{i}{\it\beta}}r)\left(\frac{1}{r^{2}}-\frac{1}{\text{i}\sqrt{2\text{i}{\it\beta}}r^{3}}\right)\right\}\cos {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{U} _{{\it\theta}}(r,{\it\theta})=\left\{\frac{\hat{A}}{r^{3}}-\hat{B}\exp (\text{i}\sqrt{2\text{i}{\it\beta}}r)\left(\frac{1}{\text{i}\sqrt{2\text{i}{\it\beta}}r^{3}}+\frac{\text{i}\sqrt{2\text{i}{\it\beta}}}{r}-\frac{1}{r^{2}}\right)\right\}\sin {\it\theta}, & \displaystyle\end{eqnarray}$$

where

(3.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{A}=\frac{1}{2}\left(1-\frac{3}{\text{i}\sqrt{2\text{i}{\it\beta}}}-\frac{3}{2\text{i}{\it\beta}}\right), & \displaystyle\end{eqnarray}$$

(3.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{B}=\frac{3\exp (-\text{i}\sqrt{2\text{i}{\it\beta}})}{2\text{i}\sqrt{2\text{i}{\it\beta}}} & \displaystyle\end{eqnarray}$$

and ${\it\Theta}=Kn^{2}{\it\beta}$ . The corresponding (continuum) force on the sphere is

(3.3) $$\begin{eqnarray}\tilde{F}_{con}={\rm\pi}Kn\left(3-3\text{i}\sqrt{2\text{i}{\it\beta}}-\frac{2\text{i}{\it\beta}}{3}\right),\end{eqnarray}$$

which has been scaled by $2a^{2}P/c_{m}$ . An identical force scale is used in all calculations.

3.1.2 Slip flow $(Kn\ll 1~\mathit{and}~{\it\Theta}\ll 1)$

It has been shown (Cercignani Reference Cercignani2000; Sone Reference Sone2002; Ivchenko, Loyalka & Tompson Reference Ivchenko, Loyalka and Tompson2007) that the leading-order effect of gas rarefaction is formally obtained by solving the linearised Navier–Stokes equation subject to a slip boundary condition. This yields the required results for the bulk velocity:

(3.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{U} _{r}(r,{\it\theta})=\left\{\frac{2\hat{A}}{r^{3}}+2\hat{B}\exp \left(\text{i}\sqrt{2\text{i}{\it\beta}}r\right)\left(\frac{1}{r^{2}}-\frac{1}{\text{i}\sqrt{2\text{i}{\it\beta}}r^{3}}\right)\right\}\cos {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{U} _{{\it\theta}}(r,{\it\theta})=\left\{\frac{\hat{A}}{r^{3}}-\hat{B}\exp \left(\text{i}\sqrt{2\text{i}{\it\beta}}r\right)\left(\frac{1}{\text{i}\sqrt{2\text{i}{\it\beta}}r^{3}}+\frac{\text{i}\sqrt{2\text{i}{\it\beta}}}{r}-\frac{1}{r^{2}}\right)\right\}\sin {\it\theta}, & \displaystyle\end{eqnarray}$$

where

(3.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{A}=\frac{1}{2}\left\{1-\frac{(3+6{\it\sigma}Kn)\left(1-\displaystyle \frac{1}{\text{i}\sqrt{2\text{i}{\it\beta}}}\right)}{\text{i}\sqrt{2\text{i}{\it\beta}}+{\it\sigma}Kn\left(3\text{i}\sqrt{2\text{i}{\it\beta}}+2\text{i}{\it\beta}\right)}\right\}, & \displaystyle\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{B}=\frac{3}{2}\left\{\frac{1+2{\it\sigma}Kn}{\text{i}\sqrt{2\text{i}{\it\beta}}+{\it\sigma}Kn\left(3\text{i}\sqrt{2\text{i}{\it\beta}}+2\text{i}{\it\beta}\right)}\right\}\exp \left(-\text{i}\sqrt{2\text{i}{\it\beta}}\right) & \displaystyle\end{eqnarray}$$

and ${\it\sigma}=1.016191\ldots$  . The corresponding scaled force is

(3.6) $$\begin{eqnarray}\displaystyle \displaystyle \tilde{F}_{slip} & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle -\frac{1}{3}{\rm\pi}Kn\nonumber\\ \displaystyle & & \displaystyle \displaystyle \times \,\left(\frac{4\text{i}{\it\beta}^{2}Kn\,{\it\sigma}-9\sqrt{2\text{i}{\it\beta}}(1+2{\it\sigma}Kn)+2\text{i}{\it\beta}(9\text{i}+\sqrt{2\text{i}{\it\beta}}+3[6\text{i}+\sqrt{2\text{i}{\it\beta}}]Kn{\it\sigma})}{\sqrt{2\text{i}{\it\beta}}+(3\sqrt{2\text{i}{\it\beta}}+2{\it\beta})Kn\,{\it\sigma}}\right).\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

There is no temperature jump present at $O(Kn)$ and hence no temperature variation in the gas for slip flow (Takata et al. Reference Takata, Sone and Aoki1993). For ${\it\beta}\gg 1$ , the flow is highly unsteady and vorticity is confined to a thin layer near the sphere surface; the thickness of this layer is $O({\it\delta})$ , where ${\it\delta}\ll a$ and ${\it\delta}=\sqrt{2{\it\mu}/({\it\rho}_{0}{\it\omega})}$ . In the opposite limit, ${\it\beta}\ll 1$ , we recover the steady solution, which is discussed in § 3.2.

3.1.3 Free-molecular flow $(Kn\gg 1~\mathit{and/or}~{\it\Theta}\gg 1)$

If either $Kn\gg 1$ or ${\it\Theta}\gg 1$ is satisfied (or both are), then the flow is free-molecular. In the free-molecular limit, the linearised Boltzmann equation in (2.7) becomes

(3.7) $$\begin{eqnarray}\frac{1-\text{i}{\it\Theta}}{Kn}\tilde{h}+\boldsymbol{c}\boldsymbol{\cdot }\frac{\partial \tilde{h}}{\partial \boldsymbol{x}}=0.\end{eqnarray}$$

Equation (3.7) is solved by integrating along its characteristics (discussed above) and using the boundary conditions in (2.12). This gives

(3.8) $$\begin{eqnarray}\tilde{h}(\boldsymbol{x},\boldsymbol{c})=\left\{\begin{array}{@{}ll@{}}\displaystyle \tilde{h}_{0}\exp \left(\text{i}\frac{1-\text{i}{\it\Theta}}{Kn}\left|\boldsymbol{x}-\boldsymbol{x}_{0}\right|\right),\quad & \boldsymbol{c}\in R1,\\ 0,\quad & \boldsymbol{c}\in R2,\end{array}\right.\end{eqnarray}$$

where

(3.9) $$\begin{eqnarray}\tilde{h}_{0}(\boldsymbol{x}_{0},\boldsymbol{c})=2\boldsymbol{c}\boldsymbol{\cdot }\tilde{\boldsymbol{U}}_{sp}+\sqrt{{\rm\pi}}\tilde{\boldsymbol{U}}_{sp}\boldsymbol{\cdot }\hat{\boldsymbol{r}},\end{eqnarray}$$

and the regions $R1$ and $R2$ are defined in (2.14). As before, $(V,{\it\chi},{\it\psi})$ are the spherical coordinates of the particle velocity $\boldsymbol{c}$ ; see (2.9) and (2.14). The density, bulk velocity and temperature are obtained by substituting (3.8) into (2.16), giving

(3.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{w}(r,{\it\theta})=\left[\bar{W}_{1}(r)+\sqrt{{\rm\pi}}\bar{R}_{1}(r)\right]\cos {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{U} _{r}(r,{\it\theta})=\left[\bar{W}_{2}(r)+\sqrt{{\rm\pi}}\bar{R}_{2}(r)\right]\cos {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{U} _{{\it\theta}}(r,{\it\theta})=\left[\bar{W}_{3}(r)+\sqrt{{\rm\pi}}\bar{R}_{3}(r)\right]\sin {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.10d ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{{\it\tau}}(r,{\it\theta})=\left[\bar{W}_{4}(r)+\sqrt{{\rm\pi}}\bar{R}_{4}(r)\right]\cos {\it\theta}, & \displaystyle\end{eqnarray}$$

where the functions $\bar{W}_{i}(r),\bar{R}_{i}(r)$ for $i\in \{1,2,3,4\}$ are defined in appendix C. The scaled force is obtained from (2.22) and (2.23):

(3.11) $$\begin{eqnarray}\tilde{F}_{fm}=\frac{\sqrt{{\rm\pi}}}{3}\left({\rm\pi}+8\right).\end{eqnarray}$$

For ${\it\Theta}\gg Kn$ , the flow is confined to a thin boundary layer of thickness $O(c_{m}/{\it\omega})$ , i.e. the average distance a gas particle travels in one oscillation cycle. When $Kn\ll {\it\Theta}$ , the unsteady term in (3.7) vanishes and steady flow is recovered; see § 3.2.

3.2 Steady flow $({\it\Theta}=0)$

The corresponding results are now summarised for steady flow. These are obtained from the results in the preceding section by taking the limit ${\it\Theta}\rightarrow 0$ , and are used to validate the numerical solution specified in § 4. The steady flow asymptotic solutions have been reported previously and are provided here for completeness.

3.2.1 Continuum flow $(Kn\rightarrow 0)$

The bulk velocity in the gas is

(3.12a,b ) $$\begin{eqnarray}\displaystyle U_{r}(r,{\it\theta})=\left(\frac{3}{2r}-\frac{1}{2r^{3}}\right)\cos {\it\theta},\quad \displaystyle U_{{\it\theta}}(r,{\it\theta})=-\left(\frac{3}{4r}+\frac{1}{4r^{3}}\right)\sin {\it\theta},\end{eqnarray}$$

and the corresponding scaled force is

(3.13) $$\begin{eqnarray}F_{con}=3{\rm\pi}Kn.\end{eqnarray}$$

3.2.2 Slip flow $(Kn\ll 1)$

The slip solution for steady flow is obtained by taking the same limit, ${\it\Theta}\rightarrow 0$ , in (3.4) while keeping $Kn$ finite. This gives the following result for the bulk velocity (Takata et al. Reference Takata, Sone and Aoki1993):

(3.14a ) $$\begin{eqnarray}\displaystyle & \displaystyle U_{r}(r,{\it\theta})=\left\{\frac{3}{2}\left(\frac{1+2{\it\sigma}Kn}{1+3{\it\sigma}Kn}\right)\frac{1}{r}-\frac{1}{2}\left(\frac{1}{1+3{\it\sigma}Kn}\right)\frac{1}{r^{3}}\right\}\cos {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.14b ) $$\begin{eqnarray}\displaystyle & \displaystyle U_{{\it\theta}}(r,{\it\theta})=-\!\left\{\frac{3}{4}\left(\frac{1+2{\it\sigma}Kn}{1+3{\it\sigma}Kn}\right)\frac{1}{r}+\frac{1}{4}\left(\frac{1}{1+3{\it\sigma}Kn}\right)\frac{1}{r^{3}}\right\}\sin {\it\theta}, & \displaystyle\end{eqnarray}$$

where the corresponding scaled force is

(3.15) $$\begin{eqnarray}F_{slip}=3{\rm\pi}Kn\left(\frac{1+2{\it\sigma}Kn}{1+3{\it\sigma}Kn}\right).\end{eqnarray}$$

Equation (3.15) coincides with the result of Ivchenko et al. (Reference Ivchenko, Loyalka and Tompson2007).

3.2.3 Free-molecular flow $(Kn\gg 1)$

The density, bulk velocity and temperature are (Takata et al. Reference Takata, Sone and Aoki1993)

(3.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle w=\displaystyle \frac{\sqrt{{\rm\pi}}}{6}\left\{\left(1+\frac{6}{{\rm\pi}}\right)\frac{1}{r^{2}}-(r^{2}-1)^{1/2}\left(2+\frac{1}{r^{2}}\right)+2r\right\}\cos {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle U_{r}=\displaystyle \left\{-\frac{1}{2}-\frac{1}{2}\left(1-\frac{1}{r^{2}}\right)^{3/2}+\frac{1}{4r^{3}}+\frac{1}{8}\left[1+\frac{1}{r^{2}}+\frac{(r^{2}-1)^{2}}{2r^{3}}\log \left(\frac{r-1}{r+1}\right)\right]\right\}\cos {\it\theta}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.16c ) $$\begin{eqnarray}\displaystyle \hspace{-15.00002pt} & \displaystyle U_{{\it\theta}}=\displaystyle \left\{\frac{1}{2}+\frac{1}{8r^{3}}+\frac{1}{4}\left(1-\frac{1}{r^{2}}\right)^{1/2}\left(2+\frac{1}{r^{2}}\right)\right. & \displaystyle \nonumber\\ \displaystyle \hspace{-15.00002pt} & \qquad \qquad \qquad \qquad \displaystyle \left.\vphantom{\left(2+\frac{1}{r^{2}}\right)^{1/2}}-\frac{1}{32}\left[6-\frac{2}{r^{2}}+\frac{3r^{4}-2r^{2}-1}{r^{3}}\log \left(\frac{r-1}{r+1}\right)\right]\right\}\sin {\it\theta}, & \displaystyle\end{eqnarray}$$

(3.16d ) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}=\displaystyle -\frac{\cos {\it\theta}}{3\sqrt{{\rm\pi}}r^{2}}. & \displaystyle\end{eqnarray}$$

At any instant in time in an unsteady flow, the distribution function for outgoing particle velocities is a Maxwellian with a mean velocity $U_{mean}\,\hat{\boldsymbol{z}}$ , where $U_{mean}\in [-U_{0},U_{0}]$ . Since particles do not collide in the free-molecular limit, the distribution function of outgoing particle velocities is identical to that generated by a sphere moving with constant velocity $U_{mean}\,\hat{\boldsymbol{z}}$ . The steady and unsteady free-molecular forces are therefore identical, and given by (3.11).

4.1 Steady flow

To solve (2.21) for steady flow, the region of integration is first truncated at a radius $r=R_{N}\gg 1$ , giving

(4.1) $$\begin{eqnarray}q_{k}(r)=\int _{1}^{R_{N}}\mathop{\sum }_{m=1}^{4}K_{km}(r,r^{\prime })q_{m}(r^{\prime })\,\text{d}r^{\prime }+T_{k}(r),\end{eqnarray}$$

where $k\in \{1,2,3,4\}$ is implied throughout. Asymptotic solutions for small $Kn$ in § 3.2 show that the bulk velocity asymptotically decays ${\sim}r^{-1}$ far from the sphere; $R_{N}$ must therefore be very large for a negligible truncation error (Takata et al. Reference Takata, Sone and Aoki1993). To reduce the computational domain, and hence computational time, the approach proposed by Takata et al. (Reference Takata, Sone and Aoki1993) is used. Because the flow must be near-continuum far from the sphere, asymptotic solutions for small $Kn$ are used in this region in place of the unknown functions $q_{k}$ in (4.1). We therefore choose some radius $R_{M}\in (1,R_{N})$ , beyond which the small-Kn asymptotic solution applies. Thus, (4.1) becomes

(4.2) $$\begin{eqnarray}q_{k}(r)=\int _{1}^{R_{M}}\mathop{\sum }_{m=1}^{4}K_{km}(r,r^{\prime })q_{m}(r^{\prime })\,\text{d}r^{\prime }+\int _{R_{M}}^{R_{N}}\mathop{\sum }_{m=1}^{4}K_{km}(r,r^{\prime })Q_{m}(r^{\prime })\,\text{d}r^{\prime }+T_{k}(r),\end{eqnarray}$$

where the near-continuum solutions $Q_{m}(r)$ have the general form (Takata et al. Reference Takata, Sone and Aoki1993)

(4.3a ) $$\begin{eqnarray}\displaystyle \displaystyle Q_{1}(r) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \frac{C_{1}Kn-C_{3}}{r^{2}},\end{eqnarray}$$

(4.3b ) $$\begin{eqnarray}\displaystyle \displaystyle Q_{2}(r) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \frac{C_{1}}{r}+\frac{C_{2}}{r^{3}},\end{eqnarray}$$

(4.3c ) $$\begin{eqnarray}\displaystyle \displaystyle Q_{3}(r) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle -\frac{C_{1}}{2r}+\frac{C_{2}}{2r^{3}},\end{eqnarray}$$

(4.3d ) $$\begin{eqnarray}\displaystyle \displaystyle Q_{4}(r) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \frac{C_{3}}{r^{2}}.\end{eqnarray}$$

$C_{1},C_{2}$

$C_{3}$

$Q_{m}(r)$

$r=R_{M}$

Next, we apply the singularity subtraction technique (Loyalka & Tompson Reference Loyalka and Tompson2009) to the first integral on the right-hand side of (4.2). This removes the singularity at $r=r^{\prime }$ in the kernel $K_{km}(r,r^{\prime })$ and improves convergence. Equation (4.2) yields

(4.4) $$\begin{eqnarray}\displaystyle \displaystyle \mathop{\sum }_{m=1}^{4}A_{km}(r)q_{k}(r) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \int _{1}^{R_{M}}\mathop{\sum }_{m=1}^{4}K_{km}(r,r^{\prime })\left(q_{m}(r^{\prime })-q_{m}(r)\right)\text{d}r^{\prime }\nonumber\\ \displaystyle & & \displaystyle \displaystyle +\int _{R_{M}}^{R_{N}}\mathop{\sum }_{m=1}^{4}K_{km}(r,r^{\prime })Q_{m}(r^{\prime })\,\text{d}r^{\prime }+T_{k}(r),\end{eqnarray}$$

where

(4.5) $$\begin{eqnarray}A_{km}(r)={\it\delta}_{km}-\int _{1}^{R_{M}}K_{km}(r,r^{\prime })\,\text{d}r^{\prime }\end{eqnarray}$$

and ${\it\delta}_{km}$ is the Kronecker delta function.

The coefficients $C_{1},C_{2},C_{3}$ in (4.3) are calculated by matching (4.3) to numerical solutions of the integral equations, (4.4), at $r=R_{M}$ , i.e.

(4.6a ) $$\begin{eqnarray}\displaystyle \displaystyle q_{1}(R_{M}) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \frac{C_{1}Kn-C_{3}}{R_{M}^{2}},\end{eqnarray}$$

(4.6b ) $$\begin{eqnarray}\displaystyle \displaystyle q_{2}(R_{M}) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \frac{C_{1}}{R_{M}}+\frac{C_{2}}{R_{M}^{3}},\end{eqnarray}$$

(4.6c ) $$\begin{eqnarray}\displaystyle \displaystyle q_{3}(R_{M}) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle -\frac{C_{1}}{2R_{M}}+\frac{C_{2}}{2R_{M}^{3}},\end{eqnarray}$$

$q_{1}(R_{M}),q_{2}(R_{M}),q_{3}(R_{M})$

Equation (4.4) is solved by first dividing the domain $r\in [1,R_{M}]$ into two subdomains, $[1,R_{L}]$ and $[R_{L},R_{M}]$ , where $R_{L}=Kn+1$ . The subdomain $[1,R_{L}]$ is further divided into $L_{1}$ intervals, and the resulting integral in each interval is evaluated using Gauss–Legendre (GL) quadrature of order $N_{1}$ . The second subdomain $[R_{L},R_{M}]$ is divided into $L_{2}$ intervals where the integral in each interval is evaluated with GL quadrature of order $N_{2}$ . This allows a dense mesh to be used close to the sphere where large velocity gradients are expected. Additionally, the integral with limits $[R_{M},R_{N}]$ is evaluated using GL quadrature of order $N_{c}$ . A schematic of the integration domain is shown in figure 4.

Figure 4. Schematic illustrating division of the integration domain into three subdomains. The subdomain $[1,R_{L}]$ is divided into $L_{1}$ intervals, each of which is subdivided into $N_{1}$ subintervals. The subdomain $[R_{L},R_{M}]$ is divided into $L_{2}$ intervals, each of which is subdivided into $N_{2}$ subintervals. The subdomain $[R_{M},R_{N}]$ is divided into $N_{c}$ intervals.

This produces the following system of linear equations:

(4.7) $$\begin{eqnarray}\displaystyle \displaystyle \mathop{\sum }_{m=1}^{4}A_{km}(r_{i})q_{k}(r_{i}) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \mathop{\sum }_{j=1}^{N_{L}}w_{j}\left\{\mathop{\sum }_{m=1}^{4}K_{km}(r_{i},r_{j})(q_{m}(r_{j})-q_{m}(r_{i}))\right\}\nonumber\\ \displaystyle & & \displaystyle +\,\displaystyle \mathop{\sum }_{p=1}^{N_{c}}w_{p}\frac{R_{N}-R_{M}}{2}\mathop{\sum }_{m=1}^{4}K_{km}(r_{i},r_{p})Q_{m}(r_{p})+T_{k}(r_{i})\end{eqnarray}$$

and

(4.8) $$\begin{eqnarray}\displaystyle \displaystyle \mathop{\sum }_{m=1}^{4}A_{km}(R_{M})q_{k}(R_{M}) & \hspace{-2.0pt}=\hspace{-2.0pt} & \displaystyle \displaystyle \mathop{\sum }_{j=1}^{N_{L}}w_{j}\left\{\mathop{\sum }_{m=1}^{4}K_{km}(R_{M},r_{j})\left(q_{m}(r_{j})-q_{m}(R_{M})\right)\right\}\nonumber\\ \displaystyle & & \displaystyle +\,\displaystyle \mathop{\sum }_{p=1}^{N_{c}}w_{p}\frac{R_{N}-R_{M}}{2}\mathop{\sum }_{m=1}^{4}K_{km}(R_{M},r_{p})Q_{m}(r_{p})+T_{k}(R_{M}),\end{eqnarray}$$

where

(4.9) $$\begin{eqnarray}N_{L}=N_{1}L_{1}+N_{2}L_{2}\end{eqnarray}$$

is the total number of meshpoints. The $j$ th discretised node to the right of $r=1$ in figure 4 is

(4.10) $$\begin{eqnarray}r_{j}=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{R_{L}-1}{2L_{1}}y_{1+(j-1)\,\text{mod}\,N_{1}}+\left(2\left\lfloor \frac{j-1}{N_{1}}\right\rfloor +1\right)\frac{R_{L}-1}{2L_{1}}+1,\quad & j=1,\ldots ,N_{1}L_{1},\\ \displaystyle \frac{R_{M}-R_{L}}{2L_{2}}z_{1+(j-N_{1}L_{1}-1)\,\text{mod}\,N_{2}}\quad \\ \quad \displaystyle +\left(2\left\lfloor \frac{j-N_{1}L_{1}-1}{N_{2}}\right\rfloor +1\right)\frac{R_{M}-R_{L}}{2L_{2}}+R_{L},\quad & j=N_{1}L_{1}+1,\ldots ,N_{L},\end{array}\right.\end{eqnarray}$$

where $r_{i},r_{j}\in (1,R_{M})$ and the corresponding weights are

(4.11) $$\begin{eqnarray}w_{j}=\left\{\begin{array}{@{}ll@{}}\displaystyle u_{1+(j-1)\,\text{mod}\,N_{1}}\frac{R_{L}-1}{2L_{1}},\quad & j=1,\ldots ,N_{1}L_{1},\\ \displaystyle v_{1+(j-N_{1}L_{1}-1)\,\text{mod}\,N_{2}}\frac{R_{M}-R_{L}}{2L_{2}},\quad & j=N_{1}L_{1}+1,\ldots ,N_{L}.\end{array}\right.\end{eqnarray}$$

Here, $y_{{\it\alpha}}$ and $u_{{\it\alpha}}$ are the usual GL abscissas and weights with $y_{{\it\alpha}}\in (-1,1)$ for ${\it\alpha}=1,2,\ldots ,N_{1}$ . Similarly, $z_{{\it\beta}}$ and $v_{{\it\beta}}$ are GL abscissas and weights of order $N_{2}$ for ${\it\beta}=1,2,\ldots ,N_{2}$ (note that ${\it\beta}$ here refers to an index, not the inverse Stokes number). The $p$ th discretised node to the right of $R_{M}$ in figure 4 is

(4.12) $$\begin{eqnarray}r_{p}=\frac{R_{N}-R_{M}}{2}b_{p}+\frac{R_{M}+R_{N}}{2},\quad p=1,\ldots ,N_{c},\end{eqnarray}$$

where $r_{p}\in (R_{M},R_{N})$ and $b_{p}$ are the usual GL abscissas with corresponding weights  $w_{p}$ .

Equations (4.6)–(4.8) specify a fully determined linear system of equations that can be solved numerically. This system is implemented and solved in Mathematica, the results of which are given in § 5.2; see § 4.3 for details of the numerical implementation.

4.2 Unsteady flow

The asymptotic solution for unsteady flow in the continuum limit (i.e. for infinitesimal $Kn$ ) shows that the bulk velocity decays far from the sphere as ${\sim}r^{-2}$ , which is faster than the rate for the steady case; cf. (3.1) and (3.12a,b ). This negates the utility of the asymptotic solution far from the sphere; we therefore choose the truncation radius $R_{N}=R_{M}$ ; see figure 4. Proceeding with the singularity subtraction method (Loyalka & Tompson Reference Loyalka and Tompson2009), (2.21) then becomes

(4.13) $$\begin{eqnarray}\mathop{\sum }_{m=1}^{4}A_{km}(r)q_{k}(r)=\int _{1}^{R_{M}}\mathop{\sum }_{m=1}^{4}K_{km}(r,r^{\prime })\left(q_{m}(r^{\prime })-q_{m}(r)\right)\,\text{d}r^{\prime }+T_{k}(r),\end{eqnarray}$$

where all variables are as defined in the preceding section.

Similarly to the steady case, (4.13) is solved by first dividing the domain $r\in [1,R_{M}]$ into two subdomains, $[1,R_{L}]$ and $[R_{L},R_{M}]$ , where $R_{L}=1+\min (Kn,Kn/{\it\Theta})$ . This choice of $R_{L}$ ensures that all steep gradients in the transport variables near the surface can be captured accurately. The subdomains $[1,R_{L}]$ and $[R_{L},R_{M}]$ are discretised as before, leading to the following system:

(4.14) $$\begin{eqnarray}\mathop{\sum }_{m=1}^{4}A_{km}(r_{i})q_{k}(r_{i})=\mathop{\sum }_{j=1}^{N_{L}}w_{j}\left\{\mathop{\sum }_{m=1}^{4}K_{km}(r_{i},r_{j})\left(q_{m}(r_{j})-q_{m}(r_{i})\right)\right\}+T_{k}(r_{i}),\end{eqnarray}$$

where the abscissas $r_{i},r_{j}$ and corresponding weights $w_{j}$ are given in (4.9)–(4.11). The integration domain for the unsteady problem is therefore a subset of that for the steady case.

Details of the numerical implementation of the above equations for the steady and unsteady problems are discussed in the next section.

4.3 Numerical implementation

The kernel and inhomogeneous terms of the integral equations defined in appendix B are written in terms of Meijer-G functions, which are evaluated in Mathematica. Integration of the kernel functions is required after applying the singularity subtraction method, and this is performed numerically using Mathematica’s NIntegrate function with a working precision of at least $45$ significant figures. Greater working precision is required as $R_{M}$ is increased. Additionally, the precision of all other calculations is set to $100$ significant figures in Mathematica because this setting ensures the Meijer-G functions are evaluated to at least $45$ significant figures.

Numerical results are calculated for $Kn\in [0.1,10]$ and for unsteady flow ${\it\Theta}\in [0.1,10]$ . For each Knudsen number and frequency ratio, different values for the computational parameters $R_{L}$ , $R_{M}$ , $R_{N}$ , $L_{1}$ and $L_{2}$ are used. These computational parameters are increased systematically to a point where they do not change the calculated force on the sphere. This numerical convergence procedure, used to obtain accurate solutions, is presented in appendix D.

5.1 Temperature jump in the free-molecular limit

Consider a stationary solid sphere in a uniform and steady gas flow in the free-molecular limit, i.e. $Kn\gg 1$ . Far from the sphere, the gas velocity is unidirectional in the Cartesian $z$ direction with constant magnitude $U_{0}$ , and the gas temperature $T_{\infty }$ is constant. Gas particles interact with the surface of the sphere via the diffuse boundary condition. The temperatures of the sphere and the free-stream gas are identical, as specified in § 2.

The mass distribution function of particles incoming to the sphere surface ( $c_{r}<0$ ) is a Maxwellian centred at a bulk velocity $\boldsymbol{U}_{\infty }=U_{0}\hat{\boldsymbol{z}}$ with temperature $T_{\infty }$ . In contrast, due to the diffuse boundary condition, the distribution function for particles leaving the sphere surface ( $c_{r}>0$ ) is a Maxwellian centred at a bulk velocity of zero, with the same temperature $T_{\infty }$ ; these are the conditions at the surface of the sphere. The complete mass distribution function for molecular velocities in the radial direction, $c_{r}$ , at the sphere surface, i.e. the one-dimensional distribution, is therefore

(5.1) $$\begin{eqnarray}\displaystyle F_{1D}=\sqrt{\frac{m}{2{\rm\pi}kT_{\infty }}}\left\{\begin{array}{@{}ll@{}}\displaystyle \exp \left(-\frac{m}{2kT_{\infty }}c_{r}^{2}\right),\quad & c_{r}>0,\\ \displaystyle \exp \left(-\frac{m}{2kT_{\infty }}\left[c_{r}-\boldsymbol{U}_{\infty }\boldsymbol{\cdot }\hat{\boldsymbol{r}}\right]^{2}\right),\quad & c_{r}<0,\end{array}\right.\end{eqnarray}$$

where $\boldsymbol{U}_{\infty }=U_{0}\hat{\boldsymbol{z}}=U_{0}(\cos {\it\theta}\,\hat{\boldsymbol{r}}-\sin {\it\theta}\,\hat{{\bf\theta}})$ ; see (2.2).

Consider the position $(r=a,{\it\theta}={\rm\pi})$ on the sphere surface; the gas impinges at normal incidence to the sphere at this position. The outward normal, $\hat{\boldsymbol{r}}$ , to the surface is therefore in the direction opposite to the uniform velocity $\boldsymbol{U}_{\infty }$ . In this case, the one-dimensional distribution function in (5.1) becomes

(5.2) $$\begin{eqnarray}\displaystyle F_{1D}=\sqrt{\frac{m}{2{\rm\pi}kT_{\infty }}}\left\{\begin{array}{@{}ll@{}}\displaystyle \exp \left(-\frac{m}{2kT_{\infty }}c_{r}^{2}\right),\quad & c_{r}>0,\\ \displaystyle \exp \left(-\frac{m}{2kT_{\infty }}\left[c_{r}+U_{0}\right]^{2}\right),\quad & c_{r}<0,\end{array}\right.\end{eqnarray}$$

which is illustrated in figure 5. It should be noted that this distribution function has a ‘width’ (or variance) that is larger than the original Maxwellian distribution of the gas far from the sphere. Because the gas temperature is directly proportional to the variance of the distribution of molecular velocities, this immediately establishes that the gas temperature at the surface of the sphere is higher than that of the free stream and solid surface. That is, the gas is hotter than the solid sphere at the position $(r=a,{\it\theta}={\rm\pi})$ .

Figure 5. The one-dimensional mass distribution function, $F_{1D}$ , at a position $(r=a,{\it\theta}={\rm\pi})$ on the surface of the stationary sphere in a uniform and steady free-molecular gas flow, $\boldsymbol{U}_{\infty }=U_{0}\hat{\boldsymbol{z}}$ . The distribution function is shown for incoming particles ( $c_{r}<0$ , red) and outgoing particles ( $c_{r}>0$ , blue). This distribution function has a larger variance than that of a Maxwellian at the equilibrium temperature, i.e. the gas temperature represented by this distribution is higher than equilibrium.

Identical arguments establish that the gas is cooler than $T_{\infty }$ at $(r=a,{\it\theta}=0)$ and its temperature is equal to that of the sphere at $(r=a,{\it\theta}={\rm\pi}/2)$ . This shows that the gas temperature around the sphere surface varies and is different from that of the sphere. Consequently, there is a temperature jump in the gas at the surface which results in a non-zero temperature variation around the sphere, i.e. the flow is non-isothermal. This phenomenon is frequently referred to as thermal polarisation, and always occurs if a section of the solid surface is oriented normal to the flow direction. It exists at all Knudsen numbers to varying degrees, vanishing in the continuum limit. Interestingly, the sign of the thermal polarisation for small Knudsen number is opposite to that of the free-molecular flow (Aoki & Sone Reference Aoki and Sone1987), i.e. the temperature of the gas is lower than that of the solid sphere at the position $(r=a,{\it\theta}={\rm\pi})$ and higher at ( $r=a$ , ${\it\theta}={\rm\pi}/2$ ).

Figure 6. Steady flow: numerical results for the radial dependence of the density $q_{1}(r)$ , velocities $q_{2}(r)$ , $q_{3}(r)$ and temperature $q_{4}(r)$ as a function of the radius $r$ , using the present method. Results are given for $Kn=0.1$ (red solid), $Kn=1$ (green dashed) and $Kn=10$ (black dotted). The solutions for $Kn=10$ are shown as solid circles, with the dotted line being the interpolated curve; the observed spacing of these points is due to the use of low-order quadrature. The dependences of the macroscopic variables on the polar angle ${\it\theta}$ are given in (2.17).

5.2 Validation for steady flow

We now return to the original problem of a sphere moving in a quiescent gas at finite Knudsen number. Before presenting the new results for unsteady flow, we validate the numerical method by comparing our calculations for steady flow with published literature results. This also provides benchmark data to investigate the effects of unsteadiness on the flow field and drag force. Due to axisymmetry, results are presented only for the radial dependences of the density, velocity and temperature field; the polar angle dependences are given in (2.17).

Figure 6 shows these radially dependent terms, i.e. $q_{1},q_{2},q_{3},q_{4}$ , for the near-continuum ( $Kn=0.1$ ), transition ( $Kn=1$ ) and near-collisionless ( $Kn=10$ ) flow regimes. As expected, the tangential component of the slip velocity, $q_{slip}=q_{3}(1)+1$ , increases with increasing Kn; note that $q_{3}(1)=-1$ coincides with the no-slip condition. The density and temperature fields are similar to their equilibrium solutions (i.e. no density or temperature variations) for near-continuum flow; see figures 6(a,d). Finite and increasing Knudsen number enhances the density and temperature variations, in line with the discussion of the previous section. Most importantly, flow at finite Knudsen number is non-isothermal and compressible.

Figure 7 presents a comparison of the numerical results (for finite Knudsen number, $Kn$ ) with the asymptotic solutions of § 3 for the bulk velocity; results are given in the appropriate near-continuum and near-collisionless regimes, with $Kn=0.1$ and $10$ . Strikingly, the numerical results for both $Kn=0.1$ and $10$ are in excellent agreement with the slip solution of (3.14) and the free-molecular solution of (3.16) respectively. While slight differences are clearly observable, this is expected because the analytical solutions in (3.14) and (3.16) are strictly valid in the asymptotic limits $Kn\rightarrow 0$ and $Kn\rightarrow \infty$ respectively. This provides a second validation of the numerical method.

Figure 7. Steady flow: comparison of current numerical results with asymptotic solutions for (a) $Kn=0.1$ and (b) $Kn=10$ . Solid lines: current numerical results. Dashed lines: asymptotic slip and free-molecular (FM) solutions, (3.14) and (3.16) respectively. The black (upper) and red (lower) lines are the radial dependences of the radial velocity, $q_{2}(r)$ , and angular velocity, $q_{3}(r)$ , respectively. The numerical solutions for (b) are solid circles with interpolated (solid line) curves. The dependence of the bulk velocity on the polar angle ${\it\theta}$ is given in (2.17).

Figure 8 compares the density, bulk velocity and temperature profiles obtained using the current numerical method with the results of Law & Loyalka (Reference Law and Loyalka1986) (for non-isothermal flows). Again, excellent agreement is found for the distinct cases $Kn=1$ and $10$ . It should be noted that the data in Law & Loyalka (Reference Law and Loyalka1986) appear to contain some typographical errors, which have been excluded from all comparisons in this article.

Next, we compare our results with both the non-isothermal and isothermal results of Law & Loyalka (Reference Law and Loyalka1986), which were reported for $Kn=0.1$ . Law & Loyalka (Reference Law and Loyalka1986) obtained isothermal results by ignoring the temperature equation; that is, by setting $q_{4}=0$ in the above formulation. As discussed above, the steady flow around a sphere is strictly non-isothermal, and this comparison is presented only to highlight the effect of imposing an (artificial) isothermal assumption. Figure 9 gives a comparison of these literature data with the results of the present numerical method.

Figure 8. Steady flow: comparison of current numerical results and the non-isothermal data of Law & Loyalka (Reference Law and Loyalka1986) for (a) $Kn=1$ and (b) $Kn=10$ . Solid lines: current numerical results. Points: Law & Loyalka (Reference Law and Loyalka1986). Radial dependences of the density $q_{1}(r)$ , radial velocity $q_{2}(r)$ , angular velocity $q_{3}(r)$ and temperature $q_{4}(r)$ are given. The current numerical results are obtained by interpolating the numerical solutions at GL abscissas. See (2.17) for the dependences of the macroscopic variables on the polar angle  ${\it\theta}$ .

Figure 9. Steady flow for $Kn=0.1$ : comparison of current numerical results for the radial dependences of the macroscopic variables and the data of Law & Loyalka (Reference Law and Loyalka1986). (a) Solid lines: current numerical results for the radial velocity $q_{2}(r)$ and angular velocity $q_{3}(r)$ . Points: Law & Loyalka (Reference Law and Loyalka1986) non-isothermal radial and angular velocity. (b) Solid lines: current numerical results for the density $q_{1}(r)$ and temperature $q_{4}(r)$ . Points: Law & Loyalka (Reference Law and Loyalka1986) non-isothermal density (squares) and temperature (diamonds). Dotted line: Law & Loyalka (Reference Law and Loyalka1986) isothermal density.

Excellent agreement between the non-isothermal results of Law & Loyalka (Reference Law and Loyalka1986) and the present method is observed for the bulk velocity; see figure 9(a). However, figure 9(b) shows that the density and temperature numerical results of the current study are significantly different from the non-isothermal data of Law & Loyalka (Reference Law and Loyalka1986). The non-isothermal density and temperature solutions from Law & Loyalka (Reference Law and Loyalka1986) appear to increase with radius from the sphere – this is unexpected because the temperature field must approach the equilibrium result far from the sphere, i.e. $q_{4}(r)\rightarrow 0$ as $r\rightarrow \infty$ . Interestingly, we find good agreement between our non-isothermal density numerical solutions and the isothermal density results of Law & Loyalka (Reference Law and Loyalka1986); see figure 9(b). Compared with larger Knudsen numbers (where excellent agreement is observed, see figure 8), the numerical simulations of flows at small Knudsen numbers converge slower and require more computational time to obtain solutions of similar accuracy; this may explain these deviations.

In table 1, we present tabulated numerical data for the normalised force, $\tilde{F}/\tilde{F}_{fm}$ , where $\tilde{F}_{fm}$ is the free-molecular force (3.11), obtained using the present method over a comprehensive range of Knudsen numbers. Available data from the literature are also presented in table 1. We find reasonable agreement with the data of Law & Loyalka (Reference Law and Loyalka1986) and excellent agreement with the data of Takata et al. (Reference Takata, Sone and Aoki1993), both of which utilise a finite-differencing scheme under non-isothermal conditions. Our results also approach the correct asymptotic limits for both near-continuum and free-molecular flow.

The new high-accuracy data set for steady flow generated by a sphere, listed under ‘current results’ in table 1, is accurate to four decimal places. The linearised Boltzmann equation for hard spheres has been solved over a wide range of Knudsen numbers for this steady flow problem (Takata et al. Reference Takata, Sone and Aoki1993). Table 1 presents an equivalent range for the Boltzmann–BGK equation, which is larger than reported previously for this collision model. This data set is expected to be of value for the benchmarking and validation of computational and analytical studies.

Table 1. The normalised force $\tilde{F}/\tilde{F}_{fm}$ (see (2.25) and (3.11)) obtained in the current study, and by Law & Loyalka (Reference Law and Loyalka1986) and Takata et al. (Reference Takata, Sone and Aoki1993), for a range of Knudsen numbers. The current results are accurate to four decimal places. Asymptotic results are included for the near-continuum (slip) (3.15) and free-molecular (FM) (3.11) regimes.

5.3 Unsteady flow

Next, we present results for unsteady flow over a broad range of Knudsen numbers, $Kn$ , and frequency ratios, ${\it\Theta}$ , encompassing the near-continuum to near-collisionless regimes. We first validate the numerical method for unsteady flow by comparing its solutions with the analytical asymptotic solutions presented in § 3.1.

Figure 10. Unsteady flow for $Kn=30,{\it\Theta}=100$ : comparison of current numerical solutions (solid line) and asymptotic free-molecular (FM) solutions in (3.10) (dotted line). The magnitude, $|q_{i}(r)|$ , and phase, ${\it\zeta}_{i}(r)$ , are given for (a,b) density, $q_{1}$ , and temperature, $q_{4}$ , and (c,d) radial, $q_{2}$ , and angular, $q_{3}$ , components of the bulk velocity. Here, $q_{i}(r)=|q_{i}(r)|\exp (\text{i}{\it\zeta}(r))$ , see (2.17).

5.3.1 Near-continuum and near-collisionless limits

In figures 10 and 11, the radial dependences of the density, bulk velocity and temperature fields are presented, together with asymptotic solutions for small and large Knudsen numbers. Figure 10 shows numerical solutions for $Kn=30,{\it\Theta}=100$ and the free-molecular asymptotic solution given by (3.11). The numerical data coincide with the asymptotic solutions in this limiting case, as required; the slight differences are due to finite Knudsen number and frequency ratio.

In figure 11, we compare numerical solutions for $Kn=0.3,{\it\Theta}=0.1$ with the slip asymptotic solution given in (3.4). We find reasonable agreement for the velocity amplitudes, see figure 11(a), with larger differences in the velocity phases. Better agreement between the numerical solutions and the slip asymptotic solutions is expected for smaller Knudsen numbers and frequency ratios – the value of $Kn=0.3$ used here is not very small and the flow is thus not strictly in the slip regime. The computational expense of the numerical solutions increases with decreasing Knudsen number, $Kn$ , which limits the range of $Kn$ that can be investigated.

Figure 11. Unsteady flow for $Kn=0.3,{\it\Theta}=0.1$ : comparison of current numerical solutions (solid line) and asymptotic slip solutions of (3.4) (dotted line). The amplitude and phase are given for the bulk velocity, where $q_{i}(r)=|q_{i}(r)|\exp (\text{i}{\it\zeta}(r))$ ,see (2.17).

Figure 12. Unsteady flow for $Kn=1,{\it\Theta}=1$ : real and imaginary parts of the (a) radial and (b) polar angular components of the velocity field. The flow is in the transition regime. Previous frequency-domain Monte Carlo data for the radial velocity are shown by solid dots (Ladiges & Sader Reference Ladiges and Sader2015).

5.3.2 Transition flows

In figures 12–14, we present velocity profiles for three representative cases: (i) $Kn=1,{\it\Theta}=1$ ; (ii) $Kn=3,{\it\Theta}=0.1$ ; (iii) $Kn=0.1,{\it\Theta}=1$ . These flows are selected because they specify strongly non-equilibrium flows that are far from the limits considered in the previous section, i.e. they are in the transition regime where analytical solution is difficult. Here, the gas mean free path is comparable to the sphere radius, and/or the oscillation frequency is similar to the interparticle collision frequency. The first case was recently reported by Ladiges & Sader (Reference Ladiges and Sader2015) using a frequency-domain Monte Carlo method, for which we provide a direct comparison, whereas the other two cases are yet to be investigated; they all represent important benchmark solutions for unsteady non-continuum three-dimensional gas flows.

Figure 13. Unsteady flow for $Kn=3,{\it\Theta}=0.1$ : real and imaginary parts of the (a) radial and (b) polar angular velocity fields. The flow is quasi-steady.

We first consider the case where both the mean free path and the interparticle collision frequency are equal to the sphere radius and oscillation frequency respectively, i.e. $Kn=1,{\it\Theta}=1$ . Figure 12 shows the real and imaginary parts of the radial and polar angular components of the velocity field for this case. Here, the viscous penetration depth is comparable to the sphere radius (and gas mean free path), and thus the flow is highly unsteady. Spatial oscillations in the velocity field, as a function of the radius $r$ , are clearly visible over approximately 10 radii from the surface of the sphere; see figure 12. A large slip velocity is also evident in figure 12(b); this is expected because the flow is far from the continuum limit. The radial component of the bulk gas velocity at the surface coincides with the sphere velocity, as required by the no-penetration condition; see figure 12(a). The current numerical results also display excellent agreement with the independent Monte Carlo simulation data of Ladiges & Sader (Reference Ladiges and Sader2015), thus providing further validation of its robustness and accuracy. Ladiges & Sader (Reference Ladiges and Sader2015) presented results only for the radial component of the velocity as a function of $r$ . Similar agreement is found for the density and temperature fields (not shown). It should be noted that the current solution exhibits far greater accuracy (results correct to four significant figures in the drag force, see above) than that achieved by the Monte Carlo method, given finite computational resources.

Next, we examine the case where the interparticle collision frequency greatly exceeds the oscillation frequency, but the sphere radius is significantly smaller than the gas mean free path, i.e. $Kn=3,{\it\Theta}=0.1$ . Figure 13 gives the real and imaginary parts of the radial and polar angular components of the velocity field for this case. In contrast to figure 12, this flow profile presents minimal spatial oscillations in the radial direction, and the imaginary part of the velocity field is very close to zero. This behaviour is expected because the inverse Stokes number is small, i.e. ${\it\beta}\equiv {\it\Theta}/Kn^{2}\approx 0.01$ , indicating that the flow is quasi-steady. However, the slip velocity is large because $Kn=3$ , giving a strongly non-continuum flow; see figure 13(b).

Finally, we examine the opposite limit where the gas mean free path is small relative to the sphere radius, but the interparticle collision frequency and oscillation frequency of the sphere are identical, i.e. $Kn=0.1,{\it\Theta}=1$ . The density and temperature fields for this flow are given in figure 14. Here, the viscous penetration depth is a small fraction of the sphere radius, i.e. ${\it\delta}/a=0.1\ll 1$ , or equivalently ${\it\beta}\gg 1$ , establishing that unsteadiness can significantly affect the flow. As such, the velocity field exhibits strong spatial oscillations in the vicinity of the sphere; these decay to zero at approximately 10 times the viscous penetration depth; see figure 14.

Figure 14. Unsteady flow for $Kn=0.1,{\it\Theta}=1$ : real and imaginary parts of the (a) density and (b) temperature fields. The flow is highly unsteady.

5.3.3 Effect of systematically increasing the oscillation frequency

We now explore the effect of increasing unsteadiness on the flow field and the resulting drag force in the transition regime. Figure 15 shows magnitude and phase numerical data using the present method for the density, radial velocity, angular velocity and temperature, for $Kn=1$ and ${\it\Theta}=0,0.3,1,3,10$ , thus systematically tuning the effects of unsteadiness. As the frequency ${\it\Theta}$ increases, the flow becomes increasingly confined to the surface of the sphere. This is evident from figure 15(a–d), where the listed amplitudes for all transport variables decrease more rapidly with radius as ${\it\Theta}$ increases. Since ${\it\Theta}/Kn=a\,{\it\omega}/c_{m}$ , the flow approaches the collisionless limit for large ${\it\Theta}$ , where the dominant length scale is the average distance, $c_{m}/{\it\omega}$ , travelled by gas particles in one oscillation cycle. As such, the dominant length scale decreases relative to the sphere radius with increasing frequency, resulting in a thin boundary layer near the sphere surface. As expected, figure 15(e–h) shows that the phase of all four hydrodynamic quantities is zero for steady flow ( ${\it\Theta}=0$ ), and increases with ${\it\Theta}$ due to the effects of inertia.

Figure 15. Unsteady flow for $Kn=1$ : amplitude (a–d) and phase (e–h) of the density, temperature, radial velocity and angular velocity fields respectively, for ${\it\Theta}=0$ , dotted; ${\it\Theta}=0.3$ , short dashed; ${\it\Theta}=1$ , medium dased; ${\it\Theta}=3$ , long dashed; ${\it\Theta}=10$ , solid.

Figure 16. Magnitude of (a) the velocity slip and (b) the temperature jump at the surface of the sphere for steady ( ${\it\Theta}=0$ ) and unsteady ( ${\it\Theta}=0.1$ ) flow as a function of $Kn$ . The dots are calculated results obtained by extrapolating GL points near the surface and the lines are interpolation curves.

Figure 15(b,d) shows that the temperature jump and slip velocity (at $r=1$ ) are non-zero and increase with frequency. This is expected, because increasing the oscillation frequency of the sphere drives the flow further from equilibrium, accentuating finite-Knudsen-number effects; both the temperature jump and the velocity slip are strictly rarefied gas dynamics phenomena. We find that at small ${\it\Theta}$ , the temperature profile decays monotonically with increasing radius from the sphere, consistent with a steady flow. On the other hand, as ${\it\Theta}$ is increased, the temperature profile becomes non-monotonic and exhibits a maximum near the sphere surface before decaying with increasing $r$ . This maximum in temperature is due to high gas density near the surface, where a thin boundary layer forms. While the Knudsen number is not very large, flow in this high-frequency regime ( ${\it\Theta}\gg 1$ ) is quasi-collisionless because gas particles do not have time to collide within a single oscillation period. The non-monotonic dependence of temperature on radius observed in figure 15(b) is thus a free-molecular effect; indeed, it resembles the temperature profile in figure 10(a) obtained for very large Knudsen number (and frequency), as expected.

The effect of unsteadiness on the slip velocity and temperature jump at the surface of the sphere ( $r=1$ ) is shown in figure 16, where results as a function of both $Kn$ and ${\it\Theta}$ are given. For low frequency ( ${\it\Theta}\ll 1$ ), asymptotic analysis of the Boltzmann–BGK equation shows that the temperature jump is an $O(Kn^{2})$ effect, while the slip velocity occurs at $O(Kn)$ (Nassios & Sader Reference Nassios and Sader2012). The results in figure 16 highlight this asymptotic behaviour: the slip velocity varies linearly with $Kn$ near the continuum limit ( $Kn\approx 0$ ), whereas the temperature dependence on $Kn$ vanishes in the vicinity of $Kn\approx 0$ . Unsteadiness (non-zero ${\it\Theta}$ ) increases the magnitude of these effects, which are driven by departures from gas equilibrium. This dominates the behaviour of the slip velocity and temperature jump as $Kn\rightarrow 0$ , with these parameters approaching finite values, i.e. no-slip and isothermal conditions are not valid for finite ${\it\Theta}$ , regardless of the value of $Kn$ studied. This observed behaviour cannot be captured using the asymptotic theory of Nassios & Sader (Reference Nassios and Sader2012), which implicitly assumes ${\it\Theta}\ll 1$ . The finite value of ${\it\Theta}=0.1$ chosen thus produces a strong departure from equilibrium, even for small Knudsen number.

Figure 17. Magnitude of the drag force on the sphere. Comparison of numerical (points) and asymptotic slip (dotted lines) solutions for steady and unsteady flow. Blue: steady flow, ${\it\Theta}=0$ . Red: unsteady flow, ${\it\Theta}=0.1$ . Yellow: unsteady flow, ${\it\Theta}=0.3$ . Equation (3.15) is used for the steady slip solution, ${\it\Theta}=0$ , and (3.6) is used for the unsteady asymptotic solutions, ${\it\Theta}=0.1$ and 0.3.

5.3.4 Normalised force

Figure 17 presents a comparison of numerical data and analytical slip asymptotic solutions for the magnitude of the normalised force on the sphere, $|\tilde{F}/\tilde{F}_{fm}|$ , as a function of the Knudsen number. Since the asymptotic slip solutions are valid for small Knudsen numbers and small frequency ratios only, i.e. weak departures from equilibrium, we consider $0.1\leqslant Kn\leqslant 1$ and ${\it\Theta}=0$ , 0.1 and 0.3; see (3.6) and (3.15).

For steady flow ( ${\it\Theta}=0$ ), excellent agreement is observed between the slip and numerical solutions for $0.1\leqslant Kn<0.3$ , with the solutions diverging with increasing Knudsen number, as expected. For an unsteady flow at moderately low frequency ( ${\it\Theta}=0.1$ ), the asymptotic slip solution agrees well with the numerical solution for $0.1\leqslant Kn\leqslant 0.2$ , with the deviation occurring at lower $Kn$ than the steady case. The asymptotic solution still predicts qualitatively correct behaviour, exhibiting a drag minimum with increasing Knudsen number, albeit at lower $Kn$ than observed in the accurate numerical solution. It is interesting to note that while the asymptotic theory for small ${\it\Theta}$ does not correctly predict the behaviour of the slip velocity for $Kn=0.1$ , ${\it\Theta}=0.1$ (see figure 16), the overall force on the sphere is captured accurately. This is because the force is the integral of the stress at the surface, which averages deviations from equilibrium. As such, non-equilibrium effects are manifested less strongly in the force in comparison to local effects, such as the slip velocity. For the largest frequency studied in figure 17, i.e. ${\it\Theta}=0.3$ , representing the strongest deviation from equilibrium, the slip solution for the force disagrees both quantitatively and qualitatively with the numerical solution, and even predicts an (incorrect) force minimum with increasing Knudsen number – the real force displays a force maximum. This latter case generates a strongly non-equilibrium flow that cannot be properly described by the (near-equilibrium) slip theory despite the small Knudsen number of $Kn=0.1$ studied.

Table 2. Unsteady force: real and imaginary parts of the normalised force on the sphere $\tilde{F}/\tilde{F}_{fm}$ , see (2.25) and (3.11), for a range of Knudsen numbers and frequency ratios. The steady solutions obtained in Table 1 are also included. The accuracy of the data is given in decimal places (dp) or relative percentage error (%). See appendix D for convergence procedures.

Comprehensive tabulated numerical data for the normalised force, $\tilde{F}/\tilde{F}_{fm}$ , are given in table 2 over a broad range of $Kn$ and ${\it\Theta}$ . This data set can be used in a host of applications, including those involving the unsteady motion of small spheres (such as nanoparticles) in gas at ambient temperature and pressure, and general unsteady motion of spheres in a rarefied atmosphere (in conjunction with a Fourier analysis), and for benchmarking future analytical and numerical work. The provided solutions span the range $Kn\in [0.1,10]$ and ${\it\Theta}\in [0.1,10]$ , encompassing the near-continuum to near-collisionless limits. Because the numerical method presented here exhibits slow convergence at small Knudsen numbers and frequency ratios, i.e. near-equilibrium flows where strong gradients in the transport variables occur near the sphere surface, the accuracy in these cases is restricted. In all cases, the accuracies of individual results are specified in table 2. This table provides the first complete data set of the unsteady force generated by a sphere executing rectilinear oscillations in an unbounded gas, as a function of $Kn$ and ${\it\Theta}$ .