2 Statement of the problem

Consider a semi-confined monoatomic gas layer of uniform density $\unicode[STIX]{x1D70C}_{0}^{\ast }$ and temperature $T_{0}^{\ast }$ surrounding an impermeable sphere of radius $r^{\ast }=r_{0}^{\ast }$ (hereafter asterisks denote dimensional quantities). The gas is initially at rest and in thermodynamic equilibrium with the sphere. At time $t^{\ast }\geqslant 0$ , radial surface pulsations are imposed with velocity

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D700}\boldsymbol{U}_{\boldsymbol{w}}^{\ast }(t^{\ast })=\unicode[STIX]{x1D700}U_{w}^{\ast }(t^{\ast })\hat{\boldsymbol{r}},\end{eqnarray}$$

together with normal heat-flux excitations

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D700}\boldsymbol{Q}_{\boldsymbol{w}}^{\ast }(t^{\ast })=\unicode[STIX]{x1D700}Q_{w}^{\ast }(t^{\ast })\hat{\boldsymbol{r}}.\end{eqnarray}$$

In (2.1) and (2.2), $\hat{\boldsymbol{r}}$ marks a unit vector in the radial direction, and $\unicode[STIX]{x1D700}\ll 1$ , so that the system description may be linearized about its initial equilibrium. In practice, the monitoring of time-varying heat-flux excitations requires the determination of the amount of energy that is transmitted to the gas. A similar issue, concerning the evaluation of energy transmission efficiency, has been considered in the context of the thermophone device (see Shinoda et al. (Reference Shinoda, Nakajima, Ueno and Koshida1999), Julius et al. (Reference Julius, Gold, Kleiman, Leizeronok and Cukurel2018), and works cited therein). Accordingly, a more comprehensive formulation would couple the thermoacoustic gas state to a heat-conduction problem within the sphere. This calculation would involve the replacement of the boundary condition (2.2) with appropriate matching conditions of heat fluxes and temperatures between the solid and fluid sides of the sphere. The solution of the coupled problem should yield a transfer function relating the amount of invested energy and the heat flux transmitted to the gas, affected by the thermal properties of the solid material. Such an analysis is not followed in the present work, which focuses on the effect of gas rarefaction on monopole sound radiation. In particular, a large portion of our discussion focuses on the generation of sound by an adiabatic (thermally insulating) sphere, which may be easily realized in experiments.

In what follows we analyse the gas response to sphere actuations in the entire range of sphere radii and excitation time scales relative to their counterpart molecular quantities. While the analysis is carried out in part for arbitrary small-amplitude actuations (see § 4), our results are presented for harmonic wall excitations,

(2.3a,b ) $$\begin{eqnarray}U_{w}^{\ast }(t^{\ast })=\overline{U}_{w}^{\ast }\cos (\unicode[STIX]{x1D714}^{\ast }t^{\ast })\quad \text{and}\quad Q_{w}^{\ast }(t^{\ast })=\overline{Q}_{w}^{\ast }\cos (\unicode[STIX]{x1D714}^{\ast }t^{\ast }),\end{eqnarray}$$

where $\unicode[STIX]{x1D714}^{\ast }$ denotes the common mechanical and heat-flux actuation frequency, and we focus on the final time-periodic state of the system. Here, $\unicode[STIX]{x1D700}$ has been omitted for convenience, yet the numerical values of $\overline{U}_{w}^{\ast }$ and $\overline{Q}_{w}^{\ast }$ will be taken small in their normalized form.

Further analysis requires the scaling of the problem. To this end we normalize the length and velocity by the nominal sphere radius $r_{0}^{\ast }$ and molecular mean thermal speed $U_{th}^{\ast }=\sqrt{2{\mathcal{R}}^{\ast }T_{0}^{\ast }}$ (where ${\mathcal{R}}^{\ast }$ marks the specific gas constant), respectively. The non-dimensional problem is then governed by the gas mean Knudsen number and scaled frequency,

(2.4a,b ) $$\begin{eqnarray}\mathit{Kn}=l_{0}^{\ast }/r_{0}^{\ast }\quad \text{and}\quad \unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}^{\ast }r_{0}^{\ast }/U_{th}^{\ast },\end{eqnarray}$$

respectively, where $l_{0}^{\ast }$ marks the molecular mean free path, specified after (3.8). In terms of the governing parameters, it is expected that free-molecular conditions should prevail when $\mathit{Kn}\gg 1$ or $\unicode[STIX]{x1D714}\mathit{Kn}\gg 1$ , for which either the length scale or time scale is short compared with the mean free path or mean free time, respectively. In contrast, as will be discussed in § 6, continuum-limit conditions require that both $\mathit{Kn}\ll 1$ and $\unicode[STIX]{x1D714}^{2}\mathit{Kn}\ll 1$ . For later convenience we denote

(2.5) $$\begin{eqnarray}\mathit{Kn}_{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D714}\mathit{Kn}=\unicode[STIX]{x1D714}^{\ast }l_{0}^{\ast }/U_{th}^{\ast },\end{eqnarray}$$

as the frequency-based Knudsen number, representing the ratio between the excitation frequency and the molecular collision frequency. The following analysis will be presented in terms of $\mathit{Kn}$ and $\unicode[STIX]{x1D714}$ , while $\mathit{Kn}_{\unicode[STIX]{x1D714}}$ will be used to facilitate the discussion in § 6.

To better assess the effect of gas rarefaction on the acoustic efficiency of the source, we recapitulate the known result for the acoustic field of a pulsating sphere at continuum inviscid ( $\mathit{Kn},\unicode[STIX]{x1D714}\mathit{Kn}\rightarrow 0$ ) conditions (Morse Reference Morse1948). Considering the case of harmonic mechanical actuations only ( $\overline{U}_{w}\neq 0,\overline{Q}_{w}=0$ ), the far-field acoustic pressure and radial velocity are given by

(2.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle p_{inviscid}(r,t)=\frac{\overline{U}_{w}\unicode[STIX]{x1D714}}{r(1+\unicode[STIX]{x1D714}^{2}/c_{0}^{2})}\left[\frac{\unicode[STIX]{x1D714}}{c_{0}}\cos (\unicode[STIX]{x1D714}t_{r})-\sin (\unicode[STIX]{x1D714}t_{r})\right]\quad \text{and}\\ \displaystyle u_{inviscid}(r,t)=\frac{\overline{U}_{w}\unicode[STIX]{x1D714}}{rc_{0}(1+\unicode[STIX]{x1D714}^{2}/c_{0}^{2})}\left[\frac{\unicode[STIX]{x1D714}}{c_{0}}\left(\frac{c_{0}^{2}}{\unicode[STIX]{x1D714}^{2}r}+1\right)\cos (\unicode[STIX]{x1D714}t_{r})+\left(\frac{1}{r}-1\right)\sin (\unicode[STIX]{x1D714}t_{r})\right]\!,\end{array}\right\}\end{eqnarray}$$

respectively, where $t_{r}=t-(r-1)/c_{0}$ is the acoustic retarded time (with $r$ denoting the scaled radial distance from the sphere centre), and $c_{0}$ marks the mean speed of sound scaled by the mean thermal speed, $c_{0}^{\ast }/U_{th}^{\ast }$ . For an ideal monatomic gas considered hereafter, $c_{0}^{\ast }=\sqrt{\unicode[STIX]{x1D6FE}{\mathcal{R}}^{\ast }T_{0}^{\ast }}$ (where $\unicode[STIX]{x1D6FE}=5/3$ marks the ratio of gas specific heats), and $c_{0}=\sqrt{5/6}$ . The overall source strength is evaluated through its acoustic power

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}(r,t)\equiv \oint _{\tilde{r}=r}p(\tilde{r},t)u(\tilde{r},t)\,\text{d}S,\end{eqnarray}$$

where the integration is carried out over a spherical surface $S$ of radius $\tilde{r}>1$ centred at $\tilde{r}=0$ . To obtain a time-independent result, we define the time-averaged acoustic power,

(2.8) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D6F1}}(r)=\frac{1}{T}\int _{0}^{T}\unicode[STIX]{x1D6F1}(r,t)\,\text{d}t,\end{eqnarray}$$

obtained by integrating $\unicode[STIX]{x1D6F1}(r,t)$ over a period $T=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$ . Substituting (2.6) into (2.7) and then into (2.8), we obtain for the inviscid far-field limit

(2.9) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D6F1}}_{inviscid}(r\gg 1)\approx 2\unicode[STIX]{x03C0}\overline{U}_{w}^{2}\sqrt{\frac{6}{5}}\left(\frac{\unicode[STIX]{x1D714}^{2}}{1+6\unicode[STIX]{x1D714}^{2}/5}\right).\end{eqnarray}$$

The result is independent of $r\gg 1$ due to the $r^{-1}$ leading-order decay of $p_{inviscid}$ and $u_{inviscid}$ , and thus remains unchanged at arbitrarily large distances from the source. As will be shown in §§ 3 and 4, non-zero gas rarefaction results in stronger decay rates which, in turn, lead to a decrease in (and far-field vanishing of) the acoustic power with increasing $r$ (see figure 7).

3 Near-continuum conditions

To analyse the problem at relatively small Knudsen numbers, we make use of the hydrodynamic R13 model mentioned in § 1 (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2003; Struchtrup Reference Struchtrup2005). The scheme has been derived for a Maxwell-type model of molecular interaction, where the three-dimensional formulation consists of balances for the scalar density and pressure, the vectors of velocity and heat flux and the stress tensor. Higher-order moments are coupled to the system through a Chapman–Enskog-type expansion of the probability density function in the kinetic Boltzmann equation. In this work, we apply the R13 model for an unsteady, one-dimensional ( $r$ -dependent) and linearized set-up. Focusing on the case of sinusoidal excitation specified in (2.3), we assume harmonic time dependence of all hydrodynamic perturbations,

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}(r,t)=Re\{\overline{\unicode[STIX]{x1D6F7}}(r)\exp [\text{i}\unicode[STIX]{x1D714}t]\}.\end{eqnarray}$$

Using spherical coordinates, the macroscopic balances of continuity, $r$ -momentum and energy are given by

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}\unicode[STIX]{x1D714}\overline{\unicode[STIX]{x1D70C}}+\frac{1}{r^{2}}\frac{\text{d}}{\text{d}r}(r^{2}\overline{u})=0, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}\unicode[STIX]{x1D714}\overline{u}=-\frac{\text{d}\overline{p}}{\text{d}r}-\frac{1}{r^{3}}\frac{\text{d}}{\text{d}r}(r^{3}\overline{\unicode[STIX]{x1D70E}}), & \displaystyle\end{eqnarray}$$

and

(3.4) $$\begin{eqnarray}\text{i}\unicode[STIX]{x1D714}\overline{T}=-\frac{2}{3r^{2}}\frac{\text{d}}{\text{d}r}[r^{2}(\overline{u}+2\overline{q})],\end{eqnarray}$$

respectively, and are supplemented by the linearized form of the equation of state,

(3.5) $$\begin{eqnarray}\overline{p}=(\overline{\unicode[STIX]{x1D70C}}+\overline{T})/2.\end{eqnarray}$$

Appearing in (3.2)–(3.5) are the amplitudes of the density perturbation $\overline{\unicode[STIX]{x1D70C}}(r)$ , radial velocity $\overline{u}(r)$ , pressure $\overline{p}(r)$ and temperature $\overline{T}(r)$ . Also appearing in (3.3)–(3.4) are the amplitudes of the radial normal stress $\overline{\unicode[STIX]{x1D70E}}(r)$ , and the radial heat flux $\overline{q}(r)$ , which satisfy the additional equations

(3.6) $$\begin{eqnarray}\text{i}\unicode[STIX]{x1D714}\overline{\unicode[STIX]{x1D70E}}+\frac{1}{r^{4}}\frac{\text{d}}{\text{d}r}(r^{4}\overline{m})=-\frac{\overline{\unicode[STIX]{x1D70E}}}{2\mathit{Kn}}-\frac{2r}{3}\frac{\text{d}}{\text{d}r}\left[r^{-1}\left(\overline{u}+\frac{4}{5}\overline{q}\right)\right]\end{eqnarray}$$

and

(3.7) $$\begin{eqnarray}\text{i}\unicode[STIX]{x1D714}\overline{q}+\frac{1}{2r^{3}}\frac{\text{d}}{\text{d}r}(r^{3}\overline{R})=-\frac{\overline{q}}{3\mathit{Kn}}-\frac{1}{2r^{3}}\frac{\text{d}}{\text{d}r}(r^{3}\overline{\unicode[STIX]{x1D70E}})-\frac{5}{8}\frac{\text{d}\overline{T}}{\text{d}r}-\frac{1}{6}\frac{\text{d}\overline{\unicode[STIX]{x1D6E5}}}{\text{d}r}.\end{eqnarray}$$

Equations (3.6) and (3.7) couple the higher-order moments’ amplitude functions $\overline{m}(r)$ , $\overline{R}(r)$ and $\overline{\unicode[STIX]{x1D6E5}}(r)$ , which are fixed through a linearized form of the closure relations derived by Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003),

(3.8a-c ) $$\begin{eqnarray}\overline{m}=-\frac{6\mathit{Kn}}{5}r^{2}\frac{\text{d}}{\text{d}r}(r^{-2}\overline{\unicode[STIX]{x1D70E}}),\quad \overline{R}=-\frac{16\mathit{Kn}}{5}r\frac{\text{d}}{\text{d}r}(r^{-1}\overline{q})\quad \text{and}\quad \overline{\unicode[STIX]{x1D6E5}}=-\frac{12\mathit{Kn}}{r^{2}}\frac{\text{d}}{\text{d}r}(r^{2}\overline{q}).\end{eqnarray}$$

The Knudsen number appearing in (3.6)–(3.8) is given by

(3.9) $$\begin{eqnarray}\mathit{Kn}=\unicode[STIX]{x1D708}_{0}^{\ast }/r_{0}^{\ast }U_{th}^{\ast },\end{eqnarray}$$

where $l_{0}^{\ast }=\unicode[STIX]{x1D708}_{0}^{\ast }/U_{th}^{\ast }$ in accordance with (2.4). At continuum-limit conditions, $\unicode[STIX]{x1D708}_{0}^{\ast }$ marks the coefficient of kinematic viscosity. Equations (3.2)–(3.8) are supplemented by the conditions of impermeability and normal heat flux at the sphere surface (see (2.1)–(2.2)),

(3.10a,b ) $$\begin{eqnarray}\overline{u}(1)=\overline{U}_{w}\quad \text{and}\quad \overline{q}(1)=\overline{Q}_{w},\end{eqnarray}$$

together with a high-order moment condition obtained from Gu & Emerson (Reference Gu and Emerson2007),

(3.11) $$\begin{eqnarray}-\frac{6}{\sqrt{\unicode[STIX]{x03C0}}}\overline{\unicode[STIX]{x1D70E}}(1)-\frac{6}{7\sqrt{\unicode[STIX]{x03C0}}}\overline{R}(1)-4\overline{m}(1)=\overline{U}_{w}+\frac{4}{5}\overline{Q}_{w},\end{eqnarray}$$

and far-field decay conditions. As expected, a temperature jump condition is missing from the present formulation, and is replaced by the heat-flux condition in (3.10). The condition (3.11) is derived by eliminating the unknown wall temperature from the linearized version of conditions (30) and (31) in Gu & Emerson (Reference Gu and Emerson2007) (or, alternatively, equations (35) and (36) in Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2008)). Note that in both Gu & Emerson (Reference Gu and Emerson2007) and Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2008), the boundary considered is not shifted in the normal direction, different from the present set-up. Hence, the wall normal velocity contribution should be added, similarly to Rana, Lockerby & Sprittles (Reference Rana, Lockerby and Sprittles2018) and Beckmann et al. (Reference Beckmann, Rana, Torrilhon and Struchtrup2018), in the context of evaporation-driven flows.

The system of equations may be recast to form a single triharmonic equation for  $\overline{T}(r)$ ,

(3.12) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathit{Kn}^{3}\left[1+{\textstyle \frac{36}{5}}\text{i}\unicode[STIX]{x1D714}\mathit{Kn}\right]\unicode[STIX]{x1D6FB}_{r}^{6}\overline{T}-{\textstyle \frac{5}{12}}\mathit{Kn}\left[1+{\textstyle \frac{62}{5}}\text{i}\unicode[STIX]{x1D714}\mathit{Kn}-{\textstyle \frac{888}{25}}(\unicode[STIX]{x1D714}\mathit{Kn})^{2}-{\textstyle \frac{2592}{125}}\text{i}(\unicode[STIX]{x1D714}\mathit{Kn})^{3}\right]\unicode[STIX]{x1D6FB}_{r}^{4}\overline{T}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\textstyle \frac{5}{18}}\text{i}\unicode[STIX]{x1D714}\left[1+{\textstyle \frac{48}{5}}\text{i}\unicode[STIX]{x1D714}\mathit{Kn}-{\textstyle \frac{864}{25}}(\unicode[STIX]{x1D714}\mathit{Kn})^{2}-{\textstyle \frac{864}{25}}\text{i}(\unicode[STIX]{x1D714}\mathit{Kn})^{3}\right]\unicode[STIX]{x1D6FB}_{r}^{2}\overline{T}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\textstyle \frac{1}{3}}\text{i}\unicode[STIX]{x1D714}^{3}[1+5\text{i}\unicode[STIX]{x1D714}\mathit{Kn}-6(\unicode[STIX]{x1D714}\mathit{Kn})^{2}]\overline{T}=0,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FB}_{r}^{2}=r^{-2}(\text{d}/\text{d}r)(r^{2}(\text{d}/\text{d}r))$ denotes the $r$ -dependent part of the Laplace operator. The equation is subject to the solution

(3.13) $$\begin{eqnarray}\overline{T}(r)=\mathop{\sum }_{n=1}^{6}\unicode[STIX]{x1D6FC}_{n}\unicode[STIX]{x1D6E9}_{n}(r),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{n}$ are constants and $\unicode[STIX]{x1D6E9}_{n}(r)$ satisfy the Helmholtz equation

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{r}^{2}\unicode[STIX]{x1D6E9}_{n}=\unicode[STIX]{x1D706}_{n}^{2}\unicode[STIX]{x1D6E9}_{n}.\end{eqnarray}$$

In (3.14), $\unicode[STIX]{x1D706}_{n}$ are the roots of the characteristic polynomial of (3.12). Expanding in $\mathit{Kn}\ll 1$ , the leading-order expressions for $\unicode[STIX]{x1D706}_{n}$ are

(3.15) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D706}_{1\pm }\approx \pm \sqrt{\frac{6}{5}}\left(\frac{7}{5}\unicode[STIX]{x1D714}^{2}\mathit{Kn}+\text{i}\unicode[STIX]{x1D714}\right),\quad \unicode[STIX]{x1D706}_{2\pm }\approx \pm \sqrt{\frac{1}{3}\frac{\unicode[STIX]{x1D714}}{\mathit{Kn}}}(1+\text{i})\quad \text{and}\\ \displaystyle \unicode[STIX]{x1D706}_{3\pm }\approx \pm \sqrt{15}\left(\frac{1}{6}\mathit{Kn}^{-1}+\frac{3}{10}\text{i}\unicode[STIX]{x1D714}\right),\end{array}\right\}\end{eqnarray}$$

and the general solution for $\overline{T}(r)$ may be put in the form

(3.16) $$\begin{eqnarray}\overline{T}(r)=\mathop{\sum }_{n=1}^{3}\left\{a_{n}\frac{\sin (\text{i}\unicode[STIX]{x1D706}_{n}r)}{\text{i}\unicode[STIX]{x1D706}_{n}r}+b_{n}\frac{\cos (\text{i}\unicode[STIX]{x1D706}_{n}r)}{\text{i}\unicode[STIX]{x1D706}_{n}r}\right\}.\end{eqnarray}$$

Expecting a casual outward radiating signal, we set $a_{n}=\text{i}b_{n}$ , which ensures that no incoming propagating wave is contained, to yield

(3.17) $$\begin{eqnarray}\overline{T}(r)=\mathop{\sum }_{n=1}^{3}C_{n}\frac{\exp [-\unicode[STIX]{x1D706}_{n}r]}{r}.\end{eqnarray}$$

For later reference, we denote the $\unicode[STIX]{x1D706}_{1},\unicode[STIX]{x1D706}_{2}$ and $\unicode[STIX]{x1D706}_{3}$ components of the solution the ‘compression’, ‘thermal’ and ‘Knudsen-layer’ modes, respectively. The compression mode is characterized by the slowest decay rate at low $\mathit{Kn}$ , and therefore dominates the far-field response, whereas the thermal and Knudsen-layer components decay much stronger. In particular, the Knudsen-layer mode decays at rate ${\sim}\exp [-\sqrt{15}\mathit{Kn}^{-1}r/6]$ , and is therefore significant only in an $O(\mathit{Kn})$ (Knudsen-layer) vicinity of the wall. Yet, apart from its direct effect on the gas behaviour inside the layer, it affects the values of the constants $C_{1}$ and $C_{2}$ multiplying the $\unicode[STIX]{x1D706}_{1,2}$ modes, which then modify the gas behaviour in the entire domain. For later discussion we note that the leading-order decay rate of $\unicode[STIX]{x1D706}_{3}$ is independent of the excitation frequency, and is therefore not affected by a change in the system time scale.

Having solved for $\overline{T}(r)$ , expressions for all other hydrodynamic fields follow by substitutions of (3.17) into the above set of balance equations, as specified in appendix A. Importantly, all fields combine the three decaying components contained in (3.17). The constants $C_{n}$ are fixed using the conditions in (3.10) and (3.11).

3.1 The Navier–Stokes–Fourier limit

To assess the differences between the R13 and continuum-limit solutions, it is instructive to derive the Navier–Stokes–Fourier (NSF) description of the system. Either by replacing (3.6) and (3.7) with the Newtonian and Fourier constitutive relations for $\overline{\unicode[STIX]{x1D70E}}(r)$ and $\overline{q}(r)$ , respectively, or by taking the consistent $Kn\ll 1$ leading-order limit of (3.12), the problem for $\overline{T}(r)$ reduces to the biharmonic equation

(3.18) $$\begin{eqnarray}\frac{5}{4}\mathit{Kn}\unicode[STIX]{x1D6FB}_{r}^{4}\overline{T}+\frac{\unicode[STIX]{x1D714}}{6}(-5\text{i}+23\unicode[STIX]{x1D714}\mathit{Kn})\unicode[STIX]{x1D6FB}_{r}^{2}\overline{T}-\text{i}\unicode[STIX]{x1D714}^{3}\overline{T}=0.\end{eqnarray}$$

Different from the R13 result, the reduction in the equation order leads to the solution

(3.19) $$\begin{eqnarray}\overline{T}(r)=\mathop{\sum }_{n=1}^{2}D_{n}\frac{\exp [-\unicode[STIX]{x1D706}_{n}^{NSF}r]}{r},\end{eqnarray}$$

where

(3.20a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{1\pm }^{NSF}\approx \pm \sqrt{\frac{6}{5}}\left(\frac{7}{5}\unicode[STIX]{x1D714}^{2}\mathit{Kn}+\text{i}\unicode[STIX]{x1D714}\right)\quad \text{and}\quad \unicode[STIX]{x1D706}_{2\pm }^{NSF}\approx \pm \sqrt{\frac{1}{3}\frac{\unicode[STIX]{x1D714}}{\mathit{Kn}}}(1+\text{i})\end{eqnarray}$$

are the roots of the characteristic polynomial corresponding to (3.18).

Inspecting the similarities between the NSF and R13 results, we observe that the leading orders of $\unicode[STIX]{x1D706}_{1,2}$ in (3.15) and $\unicode[STIX]{x1D706}_{1,2}^{NSF}$ in (3.20) are identical. The higher-order corrections, not specified here for brevity, differ in some quantitative details. Yet more importantly, the reduction in problem order in the NSF model results in the omission of the $\unicode[STIX]{x1D706}_{3}$ mode from the NSF description. At small enough $\mathit{Kn}$ , these differences have only a minor effect on the results, as the Knudsen layer (and its associated outer-layer impact) diminishes. At non-vanishing $\mathit{Kn}$ , however, the effect of the Knudsen layer becomes significant, contributing to the breakdown of the NSF description. As demonstrated by the results in § 6, the R13 model remains valid in part of these parameter combinations, yielding predictions that are in agreement with DSMC computations.

4 Free-molecular limit

Free-molecular conditions should prevail wherever $\mathit{Kn}\gg 1$ or $\mathit{Kn}_{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D714}\mathit{Kn}\gg 1$ , for which either the length scale or time scale is short compared with the mean free path or mean free time, respectively. For the spherically symmetric set-up considered, the gas state is governed by the probability density function $f=f(r,t,\unicode[STIX]{x1D743})$ of finding a gas molecule with velocity and radial position about $\unicode[STIX]{x1D743}$ and  $r$ , respectively, at time  $t$ . Expanding $f(r,t,\unicode[STIX]{x1D743})$ about its initial non-dimensional Maxwellian distribution $F=\unicode[STIX]{x03C0}^{-3/2}\exp [-\unicode[STIX]{x1D709}^{2}]$ , we write

(4.1) $$\begin{eqnarray}f(r,t,\unicode[STIX]{x1D743})=F[1+h(r,t,\unicode[STIX]{x1D743})],\end{eqnarray}$$

where $h(r,t,\unicode[STIX]{x1D743})$ denotes the unknown small ( $|h|\ll 1$ ) perturbation function due to sphere excitations. At free-molecular conditions, $h(r,t,\unicode[STIX]{x1D743})$ satisfies the collisionless one-dimensional ( $r$ -dependent) Boltzmann equation in spherical coordinates $(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ (Kogan Reference Kogan1969),

(4.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D709}_{r}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D703}}^{2}+\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D711}}^{2}}{r}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{r}}+\frac{1}{r}\left(\frac{\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D711}}^{2}}{\tan \unicode[STIX]{x1D703}}-\unicode[STIX]{x1D709}_{r}\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D703}}\right)\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D711}}}-\frac{\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D711}}}{r}\left(\frac{\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D703}}}{\tan \unicode[STIX]{x1D703}}+\unicode[STIX]{x1D709}_{r}\right)\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D711}}}=0,\end{eqnarray}$$

where the system parametrization is specified in figure 1. Performing the variable transformation

(4.3a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D709}=\sqrt{\unicode[STIX]{x1D709}_{r}^{2}+\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D703}}^{2}+\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D711}}^{2}},\quad \cos \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}=\frac{\unicode[STIX]{x1D709}_{r}}{\unicode[STIX]{x1D709}}\quad \text{and}\quad \tan \unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}=\frac{\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D711}}},\end{eqnarray}$$

equation (4.2) assumes the form

(4.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D709}\cos \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}-\frac{\unicode[STIX]{x1D709}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}{r}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}+\tan \unicode[STIX]{x1D703}\frac{\unicode[STIX]{x1D709}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\cos \unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}}{r}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}}=0.\end{eqnarray}$$

The equation is supplemented by the initial condition

(4.5) $$\begin{eqnarray}h(r,t=0^{-},\unicode[STIX]{x1D709},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}},\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}})=0,\end{eqnarray}$$

together with a linearized form of the diffuse boundary condition at the spherical surface,

(4.6) $$\begin{eqnarray}h(r=1,0\leqslant \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\leqslant \unicode[STIX]{x03C0}/2)=\unicode[STIX]{x1D70C}_{w}(t)+2\unicode[STIX]{x1D709}\cos \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}U_{w}(t)+(\unicode[STIX]{x1D709}^{2}-3/2)T_{w}(t),\end{eqnarray}$$

applied to the reflected $\unicode[STIX]{x1D709}_{r}>0$ molecules at $r=1$ . Here, $\unicode[STIX]{x1D70C}_{w}(t)$ and $T_{w}(t)$ are treated as unknown, with the latter denoting the time perturbation of the sphere temperature about $T_{0}=1$ .

Figure 1.

Schematic of the system parametrization in the free-molecular limit. The bold point marks the position of an arbitrary particle, and the outer-sphere volume confined by the blue lines and the spherical cap denotes the $\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}$ -interval acquired by the particle when reflected from the spherical surface.

Inspecting (4.6), it is noted that the boundary condition is independent of $\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}}$ . Since boundary excitation is the only cause for system deviation from equilibrium in the present problem, we assume that $h\neq h(\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}})$ , and (4.4) reduces to

(4.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D709}\cos \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}-\frac{\unicode[STIX]{x1D709}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}{r}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}=0.\end{eqnarray}$$

Equation (4.7) is amenable to the closed-form solution

(4.8) $$\begin{eqnarray}\displaystyle & & \displaystyle h(r,t,\unicode[STIX]{x1D709},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}},\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D709}})\nonumber\\ \displaystyle & & \displaystyle \quad =\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{w}(t_{w})+2\unicode[STIX]{x1D709}U_{w}(t_{w})\sqrt{1-r^{2}\sin ^{2}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}+T_{w}(t_{w})(\unicode[STIX]{x1D709}^{2}-3/2),\quad & \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\leqslant \unicode[STIX]{x1D702}\\ 0\quad & \unicode[STIX]{x1D702}<\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\leqslant \unicode[STIX]{x03C0},\end{array}\right.\qquad\end{eqnarray}$$

where

(4.9) $$\begin{eqnarray}t_{w}=\left.t-\left(r\cos \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}-\sqrt{1-r^{2}\sin ^{2}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}\right)\right/\!\unicode[STIX]{x1D709}\equiv t-x_{w}(r,\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})/\unicode[STIX]{x1D709}\end{eqnarray}$$

is the molecular-level retarded time in the ballistic limit. In addition,

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D702}=\sin ^{-1}(1/r),\end{eqnarray}$$

also indicated in figure 1.

To determine $\unicode[STIX]{x1D70C}_{w}(t)$ and $T_{w}(t)$ , the macroscopic conditions of impermeability,

(4.11) $$\begin{eqnarray}\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }h(r=1)\unicode[STIX]{x1D709}^{3}\sin 2\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709}=U_{w}(t),\end{eqnarray}$$

and heat flux,

(4.12) $$\begin{eqnarray}\frac{1}{2\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }h(r=1)\unicode[STIX]{x1D709}^{5}\sin 2\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709}-\frac{5}{4}U_{w}(t)=Q_{w}(t),\end{eqnarray}$$

are imposed at the walls, yielding

(4.13a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{w}(t)=\frac{7\sqrt{\unicode[STIX]{x03C0}}}{8}U_{w}(t)-\sqrt{\unicode[STIX]{x03C0}}Q_{w}(t)\quad \text{and}\quad T_{w}(t)=\frac{\sqrt{\unicode[STIX]{x03C0}}}{4}U_{w}(t)+2\sqrt{\unicode[STIX]{x03C0}}Q_{w}(t).\end{eqnarray}$$

Substituting (4.13) into (4.8) and then into (4.1), the density, radial velocity, pressure and radial heat-flux perturbations may be computed via velocity-space quadratures over the probability density function. Specifically,

(4.14) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70C}(r,t)=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }h(r,t,\unicode[STIX]{x1D709},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\unicode[STIX]{x1D709}^{2}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709},\\ \displaystyle u(r,t)=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }h(r,t,\unicode[STIX]{x1D709},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\unicode[STIX]{x1D709}^{3}\sin 2\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709},\\ \displaystyle p(r,t)=\frac{2}{3\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }h(r,t,\unicode[STIX]{x1D709},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\unicode[STIX]{x1D709}^{4}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709}\quad \text{and}\\ \displaystyle q(r,t)=\frac{1}{2\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }h(r,t,\unicode[STIX]{x1D709},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\unicode[STIX]{x1D709}^{5}\sin 2\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709}-\frac{5}{4}u(r,t),\end{array}\right\}\end{eqnarray}$$

whereas the temperature perturbation is given by the linearized form of the equation of state,

(4.15) $$\begin{eqnarray}T(r,t)=2p(r,t)-\unicode[STIX]{x1D70C}(r,t).\end{eqnarray}$$

Evaluation of the integrals in (4.14) is facilitated by carrying out a change in variables $\unicode[STIX]{x1D70F}=t-x_{w}(r,\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})/\unicode[STIX]{x1D709}$ , to yield

(4.16) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70C}(r,t)=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{t}{\hat{h}}(r,t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\frac{x_{w}^{3}}{(t-\unicode[STIX]{x1D70F})^{4}}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp \left[-\left(\frac{x_{w}}{t-\unicode[STIX]{x1D70F}}\right)^{2}\right]\text{d}\unicode[STIX]{x1D70F},\\ \displaystyle u(r,t)=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{t}{\hat{h}}(r,t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\frac{x_{w}^{4}}{(t-\unicode[STIX]{x1D70F})^{5}}\sin 2\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp \left[-\left(\frac{x_{w}}{t-\unicode[STIX]{x1D70F}}\right)^{2}\right]\text{d}\unicode[STIX]{x1D70F},\\ \displaystyle p(r,t)=\frac{2}{3\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{t}{\hat{h}}(r,t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\frac{x_{w}^{5}}{(t-\unicode[STIX]{x1D70F})^{6}}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp \left[-\left(\frac{x_{w}}{t-\unicode[STIX]{x1D70F}}\right)^{2}\right]\text{d}\unicode[STIX]{x1D70F}\quad \text{and}\\ \displaystyle q(r,t)=\frac{1}{2\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{t}{\hat{h}}(r,t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\frac{x_{w}^{6}}{(t-\unicode[STIX]{x1D70F})^{7}}\sin 2\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\exp \left[-\left(\frac{x_{w}}{t-\unicode[STIX]{x1D70F}}\right)^{2}\right]\text{d}\unicode[STIX]{x1D70F}-\frac{5}{4}u(r,t),\end{array}\right\}\end{eqnarray}$$

where

(4.17) $$\begin{eqnarray}\displaystyle & & \displaystyle \!\!\!\!{\hat{h}}(r,t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\nonumber\\ \displaystyle & & \displaystyle \!\!\!\!\!\quad =\left\{\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D70C}_{w}(\unicode[STIX]{x1D70F})-\frac{3}{2}T_{w}(\unicode[STIX]{x1D70F})+\frac{2x_{w}}{t-\unicode[STIX]{x1D70F}}U_{w}(\unicode[STIX]{x1D70F})\sqrt{1-r^{2}\sin ^{2}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}+\left(\frac{x_{w}}{t-\unicode[STIX]{x1D70F}}\right)^{\!2}T_{w}(\unicode[STIX]{x1D70F}),\quad & \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\leqslant \unicode[STIX]{x1D702}\\ 0\quad & \unicode[STIX]{x1D702}<\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\leqslant \unicode[STIX]{x03C0}.\end{array}\right.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The hydrodynamic fields can now be calculated via numerical integration. Towards this end, introduce

(4.18) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle I_{n}^{(m)}=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x1D702}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }U_{w}(t_{w})\sqrt{1-r^{2}\sin ^{2}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}}f^{(m)}(\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\unicode[STIX]{x1D709}^{n}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709},\\ \displaystyle J_{n}^{(m)}=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x1D702}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }\unicode[STIX]{x1D70C}_{w}(t_{w})f^{(m)}(\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\unicode[STIX]{x1D709}^{n}\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709}\quad \text{and}\\ \displaystyle K_{n}^{(m)}=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{\unicode[STIX]{x1D702}}\,\text{d}\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}\int _{0}^{\infty }T_{w}(t_{w})f^{(m)}(\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}})\exp [-\unicode[STIX]{x1D709}^{2}]\,\text{d}\unicode[STIX]{x1D709},\end{array}\right\}\end{eqnarray}$$

with $m=1,2$ and

(4.19) $$\begin{eqnarray}f^{(m)}=\left\{\begin{array}{@{}ll@{}}\sin \unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}},\quad & m=1\\ \sin 2\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}},\quad & m=2,\end{array}\right.\end{eqnarray}$$

to cast (4.16) in the form

(4.20) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70C}(r,t)=4I_{3}^{(1)}+2J_{2}^{(1)}-3K_{2}^{(1)}+2K_{4}^{(1)},\\ \displaystyle u(r,t)=2I_{4}^{(2)}+J_{3}^{(2)}-{\textstyle \frac{3}{2}}K_{3}^{(2)}+K_{5}^{(2)},\\ \displaystyle p(r,t)={\textstyle \frac{1}{3}}(4I_{5}^{(1)}+2J_{4}^{(1)}-3K_{4}^{(1)}+2K_{6}^{(1)}),\\ \displaystyle T(r,t)={\textstyle \frac{1}{3}}(-12I_{3}^{(1)}+8I_{5}^{(1)}-6J_{2}^{(1)}+4J_{4}^{(1)}+9K_{2}^{(1)}-9K_{4}^{(1)}+4K_{6}^{(1)})\quad \text{and}\\ \displaystyle q(r,t)={\textstyle \frac{1}{8}}(-20I_{4}^{(2)}+8I_{6}^{(2)}-10J_{3}^{(2)}+4J_{5}^{(2)}+15K_{3}^{(2)}-16K_{5}^{(2)}+4K_{7}^{(2)}).\end{array}\right\}\end{eqnarray}$$

Different from the analysis in § 3, the solution in the free-molecular regime, involving the numerical evaluation of $I_{n}^{(m)}$ , $J_{n}^{(m)}$ and $K_{n}^{(m)}$ , is valid for any small-amplitude $U_{w}(t)$ and $Q_{w}(t)$ input signals. Focusing on the case of harmonic excitation specified in (2.3), explicit approximations for the hydrodynamic fields may be obtained for $\unicode[STIX]{x1D714}(r-1)\gg 1$ . The approximation, presented in appendix B, applies the method of steepest descent to evaluate the $\unicode[STIX]{x1D709}$ -integrals, together with the method of stationary phase to estimate the $\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D709}}$ -integrals. These yield

(4.21) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle I_{n}^{(1)}\approx J_{n}^{(1)}\approx K_{n}^{(1)}\approx Re\left\{\frac{1}{2r}\sqrt{\frac{z_{0}^{n-2}}{3^{n-1}}}[1-\text{e}^{-(z_{0}r/3)\unicode[STIX]{x1D702}^{2}}]\left(1+\frac{a_{1}}{z_{0}}+\frac{a_{2}}{z_{0}^{2}}\right)\exp [\text{i}\unicode[STIX]{x1D714}t-z_{0}]\right\}\\ \displaystyle \text{and}\quad (I_{n}^{(2)},J_{n}^{(2)},K_{n}^{(2)})\approx 2(I_{n}^{(1)},J_{n}^{(1)},K_{n}^{(1)}),\end{array}\right\}\end{eqnarray}$$

where $z_{0}=3[\unicode[STIX]{x1D714}(r-1)/2]^{2/3}\exp [\text{i}\unicode[STIX]{x03C0}/3]\gg 1$ , and the expressions for the constants $a_{1}$ and $a_{2}$ are detailed in appendix B. Substituting (4.21) into (4.20), we find that the decay rate of all hydrodynamic perturbations is

(4.22) $$\begin{eqnarray}{\mathcal{D}}={\mathcal{D}}_{planar}{\mathcal{D}}_{curv},\end{eqnarray}$$

where ${\mathcal{D}}_{planar}$ marks the perturbations decay in the planar case (Manela & Pogorelyuk Reference Manela and Pogorelyuk2015),

(4.23) $$\begin{eqnarray}{\mathcal{D}}_{planar}\approx \sqrt{\frac{z_{0}^{n-2}}{3^{n-1}}}\left(1+\frac{a_{1}}{z_{0}}+\frac{a_{2}}{z_{0}^{2}}\right)\exp [-z_{0}],\end{eqnarray}$$

and ${\mathcal{D}}_{curv}$ denotes the additional decay caused by source sphericity,

(4.24) $$\begin{eqnarray}{\mathcal{D}}_{curv}\approx \frac{1}{2r}\left\{1-\exp \left[-\frac{z_{0}r}{3}(\sin ^{-1}(1/r))^{2}\right]\right\}.\end{eqnarray}$$

While ${\mathcal{D}}_{planar}$ reflects the dampening effect of gas rarefaction in the absence of boundary curvature, ${\mathcal{D}}_{curv}$ combines the impacts of gas rarefaction and geometric reduction on the signal decay. As specified by (4.24), ${\mathcal{D}}_{curv}$ is the product of an $r^{-1}$ decay rate (characteristic of continuum flow conditions; cf. (2.6)) with the $\{1-\exp [-(z_{0}r/3)(\sin ^{-1}(1/r))^{2}]\}$ component, caused by the narrowing in the volume occupied by the ballistically reflected molecules (see figure 1). Indeed, the number of emitted particles that reach a certain distance from the boundary (and therefore contribute to the far acoustic field) is reduced with  $r$ , resulting in a stronger decay rate compared with the planar free-molecular and sphere inviscid set-ups combined.

The above approximation should be valid at distances $(r-1)$ from the source for which $\unicode[STIX]{x1D714}(r-1)\gg 1$ . Since the signal at free-molecular conditions vanishes at short distances from $r=1$ , this requires that the actuation frequency is relatively large. Comparison between our asymptotic and full numerical solutions indicates that the present estimate holds for $\unicode[STIX]{x1D714}(r-1)\gtrsim 5$ .

5 Numerical scheme: DSMC method

The DSMC method proposed by Bird (Reference Bird1994) is a stochastic particle method commonly applied for the analysis of rarefied gas flows. In the present work, the method is used to compute the system acoustic field at arbitrary rarefied-flow conditions, and to validate our analytical ballistic- and continuum-limit solutions. We adopt Bird’s algorithm in one-dimensional spherical coordinates, together with the variable hard sphere (VHS) model of molecular interaction (Bird Reference Bird1994), to simulate the gas response. To enable comparison between the DSMC and R13 model for Maxwell molecules, the molecular collisional cross-section should to be accordingly defined. Towards this end, the Chapman–Enskog value for the dimensional dynamic viscosity of a VHS gas is given by (Bird Reference Bird1994)

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{0}^{\ast }=\frac{15(\unicode[STIX]{x03C0}m^{\ast }k_{B}^{\ast })^{1/2}(4k^{\ast }/m^{\ast })^{\unicode[STIX]{x1D71B}-1/2}T_{0}^{\ast \unicode[STIX]{x1D71B}}}{8\unicode[STIX]{x1D6E4}(9/2-\unicode[STIX]{x1D71B})\unicode[STIX]{x1D6F4}^{\ast }u_{rel}^{\ast 2\unicode[STIX]{x1D71B}-1}},\end{eqnarray}$$

where $k_{B}^{\ast }$ is the Boltzmann constant, $m^{\ast }$ and $\unicode[STIX]{x1D6F4}^{\ast }$ are the molecular mass and collisional cross-section, respectively, $u_{rel}^{\ast }$ denotes the relative velocity of the molecules, $\unicode[STIX]{x1D71B}$ marks the viscosity temperature index and $\unicode[STIX]{x1D6E4}(\cdot )$ is the gamma function. Substituting $\unicode[STIX]{x1D71B}=1$ for a Maxwell gas, and using the relation $l_{0}^{\ast }=\unicode[STIX]{x1D707}_{0}^{\ast }/\unicode[STIX]{x1D70C}_{0}^{\ast }U_{th}^{\ast }$ for the mean free path of a Maxwell molecule, we obtain

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D6F4}^{\ast }=\frac{m^{\ast }U_{th}^{\ast }}{\unicode[STIX]{x1D70C}_{0}^{\ast }u_{rel}^{\ast }l_{0}^{\ast }},\end{eqnarray}$$

exhibiting the required inverse proportion between the collisional cross-section and the molecular relative velocity. The number of molecular collisions within a DSMC cell occurring during a timestep $\unicode[STIX]{x0394}t^{\ast }$ is then given by the no-time-counter (NTC) scheme (Bird Reference Bird1994),

(5.3) $$\begin{eqnarray}\frac{N_{c}-1}{2}\frac{n_{c}^{\ast }U_{th}^{\ast }\unicode[STIX]{x0394}t^{\ast }}{n_{s}^{\ast }r_{0}^{\ast }\mathit{Kn}},\end{eqnarray}$$

where it is assumed that the collision probability of all particles in the cell is uniform. Here, $N_{c}$ marks the number of molecules in a cell, and $n_{s}^{\ast }$ and $n_{c}^{\ast }$ denote the mass densities of simulated and cell molecules, respectively.

In line with the problem statement in § 2, the spherical boundary in the simulation is modelled as fully diffuse, with prescribed radial velocity and heat flux. Thus, different from traditional realizations of gas–wall interactions, the wall temperature is treated as unknown, and application of the heat-flux condition requires modification of the conventional algorithm. In a recent contribution by the authors (Ben Ami & Manela Reference Ben Ami and Manela2017), a non-iterative procedure for the imposition of a heat-flux condition in a DSMC calculation has been presented. The algorithm has been assigned to analyse a one-dimensional cylindrical set-up, and is similarly applied in the present work.

To simulate the gas behaviour in the semi-confined gas expanse described while restricting a finite computational domain, we placed a stationary and adiabatic spherical boundary in the far field, so that it does not affect the system response. A known problem in DSMC simulations of spherical domains is in the increase in the number of particles with increasing $r$ , where, for a uniform cell width, $N_{c}\propto r_{c}^{\ast 2}\unicode[STIX]{x0394}r^{\ast }$ . In cases where the acoustic signal decays slowly with the distance from the source, this causes a vast non-uniformity in the number of particles in each cell, which leads to inaccuracies in DSMC predictions. To overcome this difficulty, the computational domain was divided into non-uniform cell sizes distributed in a geometric series, with a prescribed ratio between the farthest and closest cell size, denoted by ${\mathcal{R}}_{fc}$ (Bird Reference Bird1994). The value of ${\mathcal{R}}_{fc}$ was chosen such that the ratio between the largest and smallest number of molecules within a cell did not exceed 20. In each simulation, the one-dimensional computational grid was divided into $N_{cell}=100{-}200$ cells of radial dimension

(5.4) $$\begin{eqnarray}\unicode[STIX]{x0394}r_{m}^{\ast }=L^{\ast }a^{m-1}\left(\frac{1-a}{1-a^{N_{cell}}}\right),\end{eqnarray}$$

where $m=1,2,\ldots ,N_{cell}$ , $L^{\ast }$ is the distance between the spherical source and the ‘fictitious’ outer sphere and

(5.5) $$\begin{eqnarray}a={\mathcal{R}}_{fc}^{1/(N_{cell}-1)}.\end{eqnarray}$$

A timestep of $\unicode[STIX]{x0394}t^{\ast }=0.3\unicode[STIX]{x0394}r_{min}^{\ast }/U_{th}^{\ast }$ was assigned, where $\unicode[STIX]{x0394}r_{min}^{\ast }$ is the size of the smallest cell.

For the harmonic wall excitation studied hereafter, the computation was followed in time until the system had reached its final periodic state. A single run consisted of $5\times 10^{7}$ particles, and 16 realizations were carried out simultaneously to sufficiently reduce the numerical noise in the calculated fields. Each simulation lasted between 12 and 48 h (depending on the system Knudsen number and excitation frequency) using a ten core Intel i7-6950 machine.