3.1 Importance of volume expansion in gases

Let us first quantify the role played by volume expansion of the fluid when heat addition is present. Consider the general form of the equation of state $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}(p,T)$ . This leads to

(3.1) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\text{D}\unicode[STIX]{x1D70C}}{\text{D}t}}=\left({\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}p}\bigg|_{T}\right){\displaystyle \frac{\text{D}p}{\text{D}t}}-\left(-{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}T}}\bigg|_{p}\right){\displaystyle \frac{\text{D}T}{\text{D}t}}=\unicode[STIX]{x1D6FC}{\displaystyle \frac{\text{D}p}{\text{D}t}}-\unicode[STIX]{x1D6FD}{\displaystyle \frac{\text{D}T}{\text{D}t}}. & & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FC}$ is the isothermal compressibility coefficient of the fluid while $\unicode[STIX]{x1D6FD}$ is the bulk expansion coefficient (or thermal expansion coefficient) of the fluid. Transforming this equation using viscous scaling of the variables and noting that the scaling of time $t$ does not matter in our analysis, we introduce $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{\infty }$ , $\unicode[STIX]{x1D6FD}^{\ast }=\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FD}_{\infty }$ , $\unicode[STIX]{x1D6FC}^{\ast }=\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FC}_{\infty }$ , $p^{\ast }=p/(\unicode[STIX]{x1D707}_{\infty }U_{\infty }/D)$ , $t^{\ast }=t/t_{ref}$ , and $T^{\ast }=(T-T_{\infty })/(T_{p}-T_{\infty })=(T-T_{\infty })/\unicode[STIX]{x1D706}T_{\infty }$ , to get

(3.2) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{1}{\unicode[STIX]{x1D70C}^{\ast }}}{\displaystyle \frac{\text{D}\unicode[STIX]{x1D70C}^{\ast }}{\text{D}t^{\ast }}}=\left({\displaystyle \frac{\unicode[STIX]{x1D707}_{\infty }U_{\infty }\unicode[STIX]{x1D6FC}_{\infty }}{D}}\right)\unicode[STIX]{x1D6FC}^{\ast }{\displaystyle \frac{\text{D}p^{\ast }}{\text{D}t^{\ast }}}-(\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty })\unicode[STIX]{x1D6FD}^{\ast }{\displaystyle \frac{\text{D}T^{\ast }}{\text{D}t^{\ast }}}. & & \displaystyle\end{eqnarray}$$

Noting that the speed of sound in an ideal gas is $c_{\infty }=\sqrt{\unicode[STIX]{x1D6FE}_{\infty }p_{\infty }/\unicode[STIX]{x1D70C}_{\infty }}$ , the isothermal compressibility can be used along with the ideal gas law, $p_{\infty }=\unicode[STIX]{x1D70C}_{\infty }RT_{\infty }$ ( $R$ being the specific gas constant), to derive $\unicode[STIX]{x1D6FC}_{\infty }=(1/\unicode[STIX]{x1D70C}_{\infty })(\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}p)_{T_{\infty }}=1/(\unicode[STIX]{x1D70C}_{\infty }RT_{\infty })$ implying $c_{\infty }^{2}=\unicode[STIX]{x1D6FE}_{\infty }RT_{\infty }=\unicode[STIX]{x1D6FE}_{\infty }/\unicode[STIX]{x1D70C}_{\infty }\unicode[STIX]{x1D6FC}_{\infty }$ . Then, $\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }\gg \unicode[STIX]{x1D707}_{\infty }U_{\infty }\unicode[STIX]{x1D6FC}_{\infty }/D$ if

(3.3) $$\begin{eqnarray}\displaystyle \left({\displaystyle \frac{\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }}{\unicode[STIX]{x1D6FE}_{\infty }}}\right)\left({\displaystyle \frac{\unicode[STIX]{x1D70C}_{\infty }U_{\infty }D}{\unicode[STIX]{x1D707}_{\infty }}}\right)\gg {\displaystyle \frac{\unicode[STIX]{x1D6FC}_{\infty }U_{\infty }^{2}\unicode[STIX]{x1D70C}_{\infty }}{\unicode[STIX]{x1D6FE}_{\infty }}}={\displaystyle \frac{U_{\infty }^{2}}{c_{\infty }^{2}}}\quad \;\Longrightarrow \;\quad Re\left({\displaystyle \frac{\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }}{\unicode[STIX]{x1D6FE}_{\infty }}}\right)\gg Ma^{2}.\qquad & & \displaystyle\end{eqnarray}$$

Note that this scaling holds when $\unicode[STIX]{x1D706}\neq 0$ . Since $\unicode[STIX]{x1D6FE}_{\infty }=O(1)$ and if $\unicode[STIX]{x1D706}ReT_{\infty }$ is fixed and $T_{\infty }=283$  K, then $\unicode[STIX]{x1D6FD}_{\infty }\approx 3.55\times 10^{-3}$ (K $^{-1}$ ) for air, while it is approximately $8.88\times 10^{-5}$ (K $^{-1}$ ) for water, which is approximately 44 times smaller. Even though this is a specific example, the trend of $\unicode[STIX]{x1D6FD}_{\infty }$ being larger for gases than liquids holds. In a recent remarkable and significant paper, Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016) used the Lorentz reciprocal theorem and perturbation analysis to analytically study the motion of hot particles in a viscous liquid at very low $Re$ (e.g. Brownian motion of hot nanoscopic particles) and low $Pe$ by accounting for viscosity variations while assuming incompressibility of the fluid. Their study is the first to derive analytical expressions for the force and the torque on an arbitrarily shaped, translating and rotating particle in this regime. The error incurred due to the assumption of incompressibility in such a scenario is small, as is shown by the scaling analysis above. However, for larger particles in gases like air, the correction due to volume expansion of the fluid is significant and non-negligible, as is evidenced by our simulations. The results are presented in terms of a small parameter $\unicode[STIX]{x1D716}=((T_{p}-T_{\infty })/\unicode[STIX]{x1D707}_{0})(\text{d}\unicode[STIX]{x1D707}/\text{d}T)|_{T_{\infty }}$ with a sign that depends on the behaviour of the viscosity distribution with temperature. For example, in ideal gases, $\text{d}\unicode[STIX]{x1D707}/\text{d}T>0$ , while in common liquids $\text{d}\unicode[STIX]{x1D707}/\text{d}T<0$ . The analysis that follows this section will inherit this characteristic due to the temperature-dependent viscosity in the momentum equation. In other words, if the prescribed viscosity variation in our solver had the characteristic that $\text{d}\unicode[STIX]{x1D707}/\text{d}T<0$ , we would observe a corresponding decrease of drag with increasing temperature. Our simulations corroborate Oppenheimer et al.’s (Reference Oppenheimer, Navardi and Stone2016) assumption that distortions induced by the flow on the temperature field can be neglected in the regime when $Pe$ is very small, which corresponds to $Re\ll 1$ and $\unicode[STIX]{x1D706}\ll 1$ as is shown in figure 5. Owing to this assumption, in addition to assuming constant-density incompressibility of the fluid, it is sufficient for Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016) to solve $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ and $-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(2\unicode[STIX]{x1D707}\unicode[STIX]{x1D64E})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+(\unicode[STIX]{x1D735}\unicode[STIX]{x1D707})\boldsymbol{\cdot }[\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}]=0$ , where $\unicode[STIX]{x1D64E}$ is the rate-of-strain tensor given by $\unicode[STIX]{x1D64E}=\frac{1}{2}(\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}})$ , thereby assuming that temperature-induced density variations, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , do not significantly affect the flow. However, in the case of interest to us, where (3.3) is true, $\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }=O(1)$ and $Re(\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }/\unicode[STIX]{x1D6FE}_{\infty })\gg Ma^{2}$ , we must include variable-density effects. Figure 11 shows a comparison of the rate of dilatation for different values of $\unicode[STIX]{x1D706}$ at two values of $Re=0.01$ and $0.1$ . It can be seen that the effects of variable density are more pronounced as $Re$ and $\unicode[STIX]{x1D706}$ increase. In figure 7, we show a comparison between the numerically obtained drag and the drag that is predicted by the Oppenheimer–Navardi–Stone (Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016) first-order drag correction $F_{drag}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{\infty }RU_{\infty }(1+\frac{5}{12}\unicode[STIX]{x1D716})+O(\unicode[STIX]{x1D716})$ . This correlation evaluates to $F_{drag}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{\infty }RU_{\infty }(1+\frac{10}{36}\unicode[STIX]{x1D706})$ following the notation in this paper.

Figure 7.

Comparison of the analytically predicted drag modification from Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016) and the observed drag modification from numerical simulations presented in table 2.

3.2 An explanation of drag modification due to heating as $Re\rightarrow 0$

Clearly, when the quantity $Re(\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }/\unicode[STIX]{x1D6FE}_{\infty })$ is not small, variable-density effects are important. To account for this and inspired by Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016), we derive a set of equations simpler than the full Navier–Stokes equations that help explain the large modification in drag by leveraging the results from our numerical simulations. This formulation should not be thought of as a complete and independent problem in itself but rather as a physical explanation for the sources of the observed drag enhancement in the presence of heat addition in ideal gases. The near-sphericity of the temperature field suggests that the convective terms in the energy equation take the form of linearized Oseen-like correction terms. The proposed simplification couples the energy equation with the continuity and momentum equations only via the convective Oseen correction terms. Since the fully resolved simulations solve the full Navier–Stokes equations (2.2)–(2.5) in three dimensions, we have the computed (true) values of the scalar and vector fields at all points in the domain at steady state. The simplified set of equations at steady state are

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

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}U_{\infty }\hat{\boldsymbol{\imath }}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}, & \displaystyle\end{eqnarray}$$

(3.4c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}C_{p}U_{\infty }\hat{\boldsymbol{\imath }}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}T) & \displaystyle\end{eqnarray}$$

$p=\unicode[STIX]{x1D70C}RT$

$\unicode[STIX]{x1D749}=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}-\frac{2}{3}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})\unicode[STIX]{x1D644})$

$\unicode[STIX]{x1D644}$

$\unicode[STIX]{x1D707}$

$Re\ll 1$

$Re$

$\unicode[STIX]{x1D706}$

$\unicode[STIX]{x1D706}$

$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}\neq 0$

$Re$

$Re$

$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}\neq 0$

$\unicode[STIX]{x1D706}$

$Re$

$-(1/\unicode[STIX]{x1D70C})(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C})\approx 0$

$\unicode[STIX]{x1D70C}U_{\infty }\hat{\boldsymbol{\imath }}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\approx 0$

In what follows, we show that the change in drag of a hot particle from the isothermal value of drag is explained well by adding a suitable irrotational and dilatational velocity field to the Stokes–Oseen flow. We seek a perturbation solution of (3.4) in the limit of small temperature (implies small density) variations and write (denoting Stokes–Oseen by the suffix SO) $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{\infty }(1+\unicode[STIX]{x1D706}\unicode[STIX]{x1D70C}^{(1)}+\unicode[STIX]{x1D706}^{2}\unicode[STIX]{x1D70C}^{(2)}+\cdots \,)$ , $T=T_{\infty }(1+\unicode[STIX]{x1D706}T^{(1)}+\unicode[STIX]{x1D706}^{2}T^{(2)}+\cdots \,)$ and $\boldsymbol{u}=\boldsymbol{u}_{SO}+\unicode[STIX]{x1D706}\boldsymbol{u}^{(1)}+\unicode[STIX]{x1D706}^{2}\boldsymbol{u}^{(2)}+\cdots \,$ . At first order ( $O(\unicode[STIX]{x1D706})$ ), the fluid continuity equation becomes $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{(1)}=-(\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{(1)}\boldsymbol{\cdot }\boldsymbol{u}_{SO})$ . Mass rebalancing is solely due to density differentials that are produced when there are temperature differentials. We can represent the modified velocity field as a correction to the Stokes–Oseen velocity field, which itself can be rewritten as the Helmholtz sum of a rotational and irrotational component that correspond to a scalar potential $\unicode[STIX]{x1D719}$ and a vector potential $\boldsymbol{A}$ as $\boldsymbol{u}=\boldsymbol{u}_{SO}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}+\unicode[STIX]{x1D735}\times \boldsymbol{A}=\boldsymbol{u}_{SO}+\unicode[STIX]{x1D706}\boldsymbol{u}^{(1)}$ . The rotational correction has been described in detail in Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016). As will be seen from our analysis, we are able to recover a significant portion (more than 90 %) of the drag modification by accounting for only the irrotational correction. Thus, $\boldsymbol{u}\approx \boldsymbol{u}_{mod}=\boldsymbol{u}_{SO}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ . The canonical approach to solving a problem using a perturbation expansion is to obtain the partial differential equations and the boundary conditions at each order of $\unicode[STIX]{x1D706}$ and solve them successively. The above approximation allows us to partially bypass this process. The irrotational part captures the correction to the velocity field to conserve mass when the flow changes due to heat addition. Taking the divergence of this approximation and recognizing that the Stokes–Oseen solution is divergence-free yields

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}\approx \unicode[STIX]{x1D706}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{(1)})=-\unicode[STIX]{x1D706}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{(1)}\boldsymbol{\cdot }\boldsymbol{u}_{SO}). & & \displaystyle\end{eqnarray}$$

In the low- $Re$ or the Stokes–Oseen limit, the contribution from the higher-order expansion terms in the perturbation expansion is small. The true solutions to $\unicode[STIX]{x1D70C}$ , $T$ and $\boldsymbol{u}$ are already known from the fully resolved simulations. Thus, as $\unicode[STIX]{x1D706}\rightarrow 0$ , we have that $\unicode[STIX]{x1D70C}^{(1)}\approx (\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{\infty })/\unicode[STIX]{x1D706}$ , where $\unicode[STIX]{x1D70C}$ is the density that is known from the simulations.

An analytical expression for $\boldsymbol{u}_{SO}$ is known. We solve the Poisson equation (3.5) to obtain $\unicode[STIX]{x1D719}$ . On the surface of the sphere, we impose a no-penetration boundary condition which is equivalent to imposing a Neumann boundary condition. Far away from the sphere, the velocity correction goes to zero as the velocity field becomes the imposed uniform velocity field. Thus, equation (3.5) is subject to $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\hat{\boldsymbol{r}}=0$ if $r=R$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0$ as $r\rightarrow \infty$ , where $r$ is the radial coordinate centred at the centre of the sphere and $\hat{\boldsymbol{r}}$ is the unit vector normal to the sphere. Since we are solving a Poisson equation with Neumann boundary conditions, the obtained solution $\unicode[STIX]{x1D719}$ is unique up to a constant; its value does not affect the rest of the solution. Calculating $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ yields the approximation to the modified velocity field denoted by $\boldsymbol{u}_{mod}=\boldsymbol{u}_{SO}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ . On substituting this modified velocity in the momentum equation, we obtain an equation for the approximation to the modified hydrodynamic pressure, denoted as $p_{mod}$ . Taking the divergence of both sides of this equation yields a Poisson equation for $p^{(1)}$ (where $p_{mod}=p_{SO}(1+\unicode[STIX]{x1D706}p^{(1)}+\cdots \,)$ ), which can be solved with appropriate boundary conditions. The equation at $O(\unicode[STIX]{x1D706})$ for $p^{(1)}$ is

(3.6) $$\begin{eqnarray}\displaystyle -p_{\infty }\unicode[STIX]{x1D6FB}^{2}p^{(1)}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{\infty }U_{\infty }\hat{\boldsymbol{\imath }}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{(1)}+\unicode[STIX]{x1D70C}_{\infty }\unicode[STIX]{x1D70C}^{(1)}U_{\infty }\hat{\boldsymbol{\imath }}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{SO}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{(1)})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D733}(\boldsymbol{x}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D749}=\unicode[STIX]{x1D749}_{SO}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D749}^{(1)}$ and is subject to the boundary conditions $-p_{\infty }\unicode[STIX]{x1D735}p^{(1)}\boldsymbol{\cdot }\hat{\boldsymbol{r}}=\unicode[STIX]{x1D733}(\boldsymbol{x})\boldsymbol{\cdot }\hat{\boldsymbol{r}}$ if $r=R$ and $p^{(1)}=0$ if $r\rightarrow \infty$ . It must also be noted that, since $\unicode[STIX]{x1D70C}$ is obtained from the numerical solution and fully known, the temperature field is also known due to the equation of state. Now $\boldsymbol{u}_{mod}$ and $p_{mod}$ can be used to calculate the viscous, pressure and total drag on the sphere, which can be compared to the values obtained from the full 3-D solver.

Table 5.

Comparison of the total drag and its components (pressure (P Drag) and viscous (V Drag) components) computed using the 3-D solver with those computed from the perturbation methodology as a function of $\unicode[STIX]{x1D706}$ at $Re=0.001$ and $0.01$ . The last column shows the error incurred in the perturbation approach compared to the 3-D solver.

The results of this low-Mach-number perturbation analysis are tabulated for the range of Boussinesq parameters $\unicode[STIX]{x1D706}$ under study at $Re=0.001$ and $Re=0.01$ in table 5, which compares the total coefficient of drag and its components obtained from solving the full 3-D Navier–Stokes equations to that obtained from the perturbation approach. We see that the perturbation approach recovers the true drag to within 10 % error even when $\unicode[STIX]{x1D706}$ is as large as 2 and strong non-Boussinesq effects are present. To develop a physical sense of the magnitude of heat supplied, consider the case when the ambient temperature of the gas is 300 K. The value $\unicode[STIX]{x1D706}=2$ corresponds to a particle temperature of 900 K. The Stokes–Oseen solution is obtained as a linearized approximation to the Navier–Stokes equations. As $\unicode[STIX]{x1D706}$ increases, the heat addition into the fluid increases causing strong density gradients, which increases the nonlinear coupling of the momentum and energy equations. As these nonlinearities rise, the explanation power of the linear model decreases. With the perturbation methodology, the drag decomposition maintains the approximate ratio of the pressure and viscous drag of $1:2$ consistent with the numerical solution. This further vindicates the physical validity of the perturbation methodology in the low- $Re$ limit. Comparison of the velocity contours resulting from the full 3-D solver and the contours resulting from the perturbation methodology is shown for two cases of $Re$ and $\unicode[STIX]{x1D706}$ in figure 8. Just as with drag, better agreement is seen at lower $Re$ and lower $\unicode[STIX]{x1D706}$ when the reduced equations are a good approximation to the nonlinear coupled equations. At higher Reynolds numbers, this approach fails to recover the drag even at very small $\unicode[STIX]{x1D706}$ due to the presence of nonlinear terms which cannot be neglected.

Figure 8.

Comparison of normalized $x$ -velocity,  $U/U_{\infty }$ , contours resulting from the 3-D solver (a–c) and the contours resulting from the perturbation methodology (d–f), each for three cases of the Boussinesq parameter: one with low heat addition ( $\unicode[STIX]{x1D706}=0.1$ ), one with moderate heat addition ( $\unicode[STIX]{x1D706}=0.5$ ), and one with large heat addition ( $\unicode[STIX]{x1D706}=1$ ), for $Re=0.01$ .

3.3 Suggested correlations for $C_{D}$ and $Nu_{D}$ via least-squares regression

The explanation of the observed drag modification as the effect of an irrotational correction to the Stokes–Oseen solution due to mass rebalancing and the effect of the fluid’s variable properties is encouraging. This approach is also consistent with the form of the correlations for the drag coefficient of a sphere at low Reynolds numbers in the literature as seen in Oseen (Reference Oseen1910), Goldstein (Reference Goldstein1929), Chester, Breach & Proudman (Reference Chester, Breach and Proudman1969), Chester et al. (Reference Chester, Breach and Proudman1969) and others stated in tables 3.1 and 5.1 in Clift et al. (Reference Clift, Grace and Weber1978), where the first term of the correlation is Stokes’ drag (independent of inertial effects and hence independent of the density of the fluid) and the subsequent terms account for inertial corrections to the Stokes’ drag and is a function of the density of the fluid. As $Re\rightarrow O(1)$ , inertial effects begin to become significant and the correction terms account for this effect. Inspired by the observations from our perturbation analysis and in order to collapse the data obtained from the parametric numerical study over the $\unicode[STIX]{x1D706}$ and $Re$ ranges studied in this paper ( $10^{-3}\leqslant Re\leqslant 10$ , $0\leqslant \unicode[STIX]{x1D706}\leqslant 3$ ), we seek a drag correlation composed of two terms: a first term that resembles the Stokes–Oseen drag (and thus uses non-dimensional parameters defined at the ambient fluid temperature) and a correction term that captures the mass rebalance due to the variable density of the fluid caused by the heating of the sphere (which uses non-dimensional parameters defined at an appropriate temperature, which is a hyper-parameter to be found by tuning). Taking inspiration from the correlations in table 5.1 in Clift et al. (Reference Clift, Grace and Weber1978), we model the second term to be of power-law form. We use all the data from the fully resolved simulations of the full 3-D Navier–Stokes equations over the two parameters $Re$ and $\unicode[STIX]{x1D706}$ . Statistical methods of feature selection and hypothesis testing are used to choose the optimal form of the power-law terms in the nonlinear regression procedure, while hold-out cross-validation is used to choose between multiple hypotheses (Vapnik Reference Vapnik1998). The regression has multiple hyper-parameters that must be chosen; for example, the density and viscosity used to define the Reynolds number used in the correction terms (denoted as the film temperature Reynolds number ( $Re_{T_{film}}$ )) can be evaluated at film temperatures which are a linear combination of $T_{p}$ and $T_{\infty }$ . A full description of the process of building these accurate correlations along with the algorithm used to tune the hyper-parameters is described in Ganguli (Reference Ganguli2018).

The Clift–Grace–Weber correlation (Clift et al. Reference Clift, Grace and Weber1978), valid for $Re\leqslant 20$ and stated in table 1, is extended to take the form $C_{D}=C_{D}^{CGW}+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D706}^{\unicode[STIX]{x1D6FD}}Re_{T_{film}}^{\unicode[STIX]{x1D6FE}}$ . The adjusted- $R^{2}$ is used as the metric to evaluate these regression fits. The best fit to the data is obtained when $Re_{T_{film}}=\unicode[STIX]{x1D70C}_{\infty }U_{\infty }D/\unicode[STIX]{x1D707}(T_{film})$ and $T_{film}=T_{\infty }$ . It has an adjusted- $R^{2}$ value of 99.56 % and the values of the hyper-parameters along with their 95 % confidence intervals are $\unicode[STIX]{x1D6FC}=10.7672\pm 0.0416$ , $\unicode[STIX]{x1D6FD}=0.9673\pm 0.0130$ and $\unicode[STIX]{x1D6FE}=-0.9529\pm 0.0095$ . Thus the correlation for the drag of a heated sphere when $10^{-3}\leqslant Re\leqslant 10$ with the definitions of the hyper-parameters defined above is

(3.7) $$\begin{eqnarray}\displaystyle C_{D}={\displaystyle \frac{24}{Re}}(1+0.1315Re^{0.82-0.05\log _{10}Re})+10.7672\unicode[STIX]{x1D706}^{0.9673}Re^{-0.9529}. & & \displaystyle\end{eqnarray}$$

A feature ablation study reveals that if the feature $\unicode[STIX]{x1D706}^{\unicode[STIX]{x1D6FD}}$ is left out, the best adjusted- $R^{2}$ drops to 57.70 %. If the feature $Re_{T_{film}}^{\unicode[STIX]{x1D6FE}}$ is left out, the adjusted- $R^{2}$ drops to a paltry 13.58 %. Rewritten in terms of $C_{DV}$ , the correlation becomes

(3.8) $$\begin{eqnarray}\displaystyle C_{DV}=6(1+0.1315Re^{0.82-0.05\log _{10}Re})+2.6918\unicode[STIX]{x1D706}^{0.9673}Re^{0.0471}. & & \displaystyle\end{eqnarray}$$

(As observed in (3.8), the exponent of $\unicode[STIX]{x1D706}$ is ${\approx}1$ and the exponent of $Re$ is ${\approx}0$ . This, and the insight from § 3.2 that the drag modification can be captured quite well with a linear perturbation expansion in $\unicode[STIX]{x1D706}$ , suggests an alternative correlation for $C_{DV}$ of a heated particle with an additive correction $\unicode[STIX]{x1D6FC}_{M}\unicode[STIX]{x1D706}$ . Linear regression on the data from tables 3 and 4 reveals the representation of the best fit as $C_{DV}=6(1+0.1315Re^{0.82-0.05\log _{10}Re})+2.3801\unicode[STIX]{x1D706}$ with an adjusted- $R^{2}$ value of 96.08 %. The 95 % confidence interval on this coefficient is $\unicode[STIX]{x1D6FC}_{M}=2.3801\pm 0.0423$ .)

It is worth mentioning that when $\unicode[STIX]{x1D706}=0$ , the correlations for drag in (3.7) and (3.8) reduce to the Clift–Grace–Weber correlation. A similar approach is used to derive a correlation for the Nusselt number from the data in table 4. When the velocity is everywhere zero (which is the case when $Re\rightarrow 0$ ) and the temperature difference between the fluid and sphere is small ( $\unicode[STIX]{x1D706}<O(10^{-2})$ ), diffusion of heat is in the radial direction only and thus $Nu_{D}=2$ . From the work of Brenner (Reference Brenner1963) and Masliyah & Epstein (Reference Masliyah and Epstein1972), we know that for $\unicode[STIX]{x1D706}\rightarrow 0$ and very small $Re$ , $Nu_{D}$ approaches 2, while for large $Re$ , $Nu_{D}$ becomes proportional to $Re^{1/3}$ . The neglected convection terms are not negligible far away from the sphere since the ratio of convective to diffusive terms is $O(Re(r/a))$ . Asymptotic solutions are available for $Nu_{D}$ using the matching procedure of Proudman & Pearson (Reference Proudman and Pearson1957). Analogous to the Oseen correction for drag, Brenner (Reference Brenner1963) suggests that for small $\unicode[STIX]{x1D706}$ , we have $\frac{1}{2}Nu-1=\frac{1}{4}Re$ as $Re\rightarrow 0$ and is accurate to $O(Re)$ . Acrivos & Taylor (Reference Acrivos and Taylor1962) extended this to higher-order terms with only slight improvement. Clift et al. (Reference Clift, Grace and Weber1978) state that for small $\unicode[STIX]{x1D706}$ , $Nu_{D}=1+(1+Re)^{1/3}$ , which is valid for all $Re$ in creeping flow within 2 % error. This may be more generally written (Clift et al. Reference Clift, Grace and Weber1978) as $Nu_{D}=1+(1+Pe)^{1/3}$ , where $Pe=RePr$ . This correlation is extended to take the form (Ganguli Reference Ganguli2018), $Nu_{D}=[1+(1+Re)^{1/3}]\text{e}^{-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D706}-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D706}^{2}-\unicode[STIX]{x1D6FE}(Re-Re_{T_{film}})}$ , where $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ and $T_{film}$ are the hyper-parameters. The statistical procedure outlined in Ganguli (Reference Ganguli2018) yields the optimal definition of $Re_{T_{film}}$ for the drag correlation to be $Re_{T_{film}}=\unicode[STIX]{x1D70C}_{\infty }U_{\infty }D/\unicode[STIX]{x1D707}(T_{film})$ such that $T_{film}$ is $T_{film}=0.080T_{p}+0.920T_{\infty }$ . The best fit has an adjusted- $R^{2}$ value of 97.59 %. If, instead, the conventional definition of the film temperature is used, i.e. $T_{film}=0.5(T_{p}+T_{\infty })$ , the fit has an adjusted- $R^{2}$ of ${\sim}97.05\,\%$ . It must be noted that the correlations proposed here depend on the specific form of the power law used, and using other power laws may induce discrepancies that may surmount this small difference in the values of adjusted- $R^{2}$ . We will therefore use the conventional definition of $T_{film}$ . The values of the hyper-parameters along with their 95 % confidence intervals are $\unicode[STIX]{x1D6FC}=0.2201\pm 0.0082$ , $\unicode[STIX]{x1D6FD}=-0.0331\pm 0.0005$ and $\unicode[STIX]{x1D6FE}=-0.0046\pm 0.0001$ . A feature ablation study reveals that if the feature $\text{e}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D706}}$ is suppressed by setting $\unicode[STIX]{x1D6FC}=0$ , the adjusted- $R^{2}$ drops to 56.59 %. Suppressing $\text{e}^{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D706}^{2}}$ by setting $\unicode[STIX]{x1D6FD}=0$ yields an adjusted- $R^{2}$ of 91.24 %. Suppressing $\text{e}^{\unicode[STIX]{x1D6FE}(Re-Re_{T_{film}})}$ by setting $\unicode[STIX]{x1D6FE}=0$ yields an adjusted- $R^{2}$ of 95.94 %. Thus the correlation for the Nusselt number of a heated sphere becomes

(3.9) $$\begin{eqnarray}\displaystyle Nu_{D}=[1+(1+Re)^{1/3}]\text{e}^{-0.2201\unicode[STIX]{x1D706}+0.0331\unicode[STIX]{x1D706}^{2}+0.0046(Re-Re_{T_{film}})} & & \displaystyle\end{eqnarray}$$

when $10^{-3}\leqslant Re\leqslant 10$ .

In both of the above correlations, the value of $\unicode[STIX]{x1D707}$ is obtained from the power law for viscosity stated in § 2. It is worth mentioning that, when $\unicode[STIX]{x1D706}=0$ , the correlation for $Nu_{D}$ in (3.9) reduces to the Clift–Grace–Weber correlation for $Nu_{D}$ . Figure 9 shows the performance of the inferred correlations for drag and the Nusselt number. Figure 9(a) shows the plot of the drag correlation while figure 9(b) shows the plot of the Nusselt number correlation as a function of $\unicode[STIX]{x1D706}$ for the five decades of $Re$ . The correlations are plotted as lines and the values obtained from the particle-resolved simulations are plotted as filled diamonds. Since these correlations have been inferred from simulations of fixed particles, they do not include the effect of mobility of particles. However, we expect our correlations to hold as long as a particle moves at nearly constant velocity. One way to measure the impact of the drag correlation is to measure its effect on the terminal velocity of a hot particle falling due to gravity in a fluid, starting at rest. If the fluid density is $\unicode[STIX]{x1D70C}_{f}$ , particle density is $\unicode[STIX]{x1D70C}_{p}$ , particle volume is $V$ , acceleration due to gravity is $g$ , particle diameter is $D$ , time is $t$ and particle velocity is $u$ , the equation governing the falling particle is

(3.10) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{p}V\frac{\text{d}u}{\text{d}t}=(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})Vg-\frac{3V\unicode[STIX]{x1D70C}_{f}}{4D}u^{2}C_{D}.\end{eqnarray}$$

(For metal particles (e.g. nickel) falling in air, $\unicode[STIX]{x1D70C}_{p}\gg \unicode[STIX]{x1D70C}_{f}$ , so that the influence of temperature on the density ratio $\unicode[STIX]{x1D70C}_{f}/\unicode[STIX]{x1D70C}_{p}$ can be safely ignored for the analysis shown here. Justifiably, we can also assume zero or negligibly small contributions from fluid inertia, added mass and Basset–Boussinesq history terms in the Maxey–Riley–Gatignol equation. When $Re\ll 1$ , the added-mass factor is $\unicode[STIX]{x1D6FD}=3\unicode[STIX]{x1D70C}_{f}/(\unicode[STIX]{x1D70C}_{f}+2\unicode[STIX]{x1D70C}_{p})\ll 1$ .)

At steady state, we have

(3.11) $$\begin{eqnarray}u^{2}C_{D}=\left(1-\frac{\unicode[STIX]{x1D70C}_{f}}{\unicode[STIX]{x1D70C}_{p}}\right)g\left(\frac{4\unicode[STIX]{x1D70C}_{p}D}{3\unicode[STIX]{x1D70C}_{f}}\right).\end{eqnarray}$$

If $\unicode[STIX]{x1D70C}_{f}$ , $\unicode[STIX]{x1D70C}_{p}$ , $g$ and $D$ are constant, then $u_{terminal}^{2}C_{D}$ is a constant, allowing us to obtain the ratio of the terminal velocity of a hot particle ( $u_{terminal}^{heated}$ ) and the terminal velocity of an unheated particle ( $u_{terminal}^{CGW}$ ) as $u_{terminal}^{heated}/u_{terminal}^{CGW}=\sqrt{C_{D}^{CGW}/C_{D}^{heated}}$ . This ratio is plotted as a function of $\unicode[STIX]{x1D706}$ for $Re=0.001$ , 0.01 and 0.1 in figure 9(c).

Figure 9.

Performance of the suggested correlations for (a) $C_{DV}$ and (b) $Nu_{D}$ stated in (3.8) and (3.9), respectively, as a function of $Re$ and $\unicode[STIX]{x1D706}$ . Lines denote correlations and diamonds denote simulation results. (c) The ratio of the terminal velocity of a heated particle to an unheated particle falling under gravity starting from rest.

3.4 A scaling for the rate of dilatation for $\unicode[STIX]{x1D706}\rightarrow 0$ and $Re\rightarrow 0$

We saw in previous sections that when $Re(\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }/\unicode[STIX]{x1D6FE}_{\infty })\gg Ma^{2}$ , $\unicode[STIX]{x1D6FD}_{\infty }\unicode[STIX]{x1D706}T_{\infty }=O(1)$ and $(\unicode[STIX]{x1D707}_{\infty }U_{\infty }\unicode[STIX]{x1D6FC}_{\infty })/D\ll 1$ , we must include variable-density effects. These effects result in a finite rate of dilatation. In this section, we analytically derive a scaling for the rate of dilatation ( $\unicode[STIX]{x1D701}$ ) in the limit of $Re\rightarrow 0$ and $\unicode[STIX]{x1D706}\rightarrow 0$ . Following the analysis in Leal (Reference Leal2007), assuming $\unicode[STIX]{x1D705}$ and $\unicode[STIX]{x1D707}$ as constants, and assuming viscous dissipation is negligible, the steady-state thermal energy equation for the fluid around a heated sphere is $Pe(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E9})=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6E9}$ subject to $\unicode[STIX]{x1D6E9}(|\boldsymbol{r}|=R)=1$ and $\unicode[STIX]{x1D6E9}(|\boldsymbol{r}|\rightarrow \infty )\rightarrow 0$ . Here, $Pe=RePr$ and $\unicode[STIX]{x1D6E9}=(T_{g}-T_{\infty })/(T_{p}-T_{\infty })$ . Seeking the solution to this equation in the limit $Pe\rightarrow 0$ , a regular perturbation expansion of $\unicode[STIX]{x1D6E9}$ as $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D6E9}_{0}(\boldsymbol{r};Re)+Pe\unicode[STIX]{x1D6E9}_{1}(\boldsymbol{r};Re)+O(Pe^{2})$ yields at first order $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6E9}_{0}=0$ subject to $\unicode[STIX]{x1D6E9}_{0}(|\boldsymbol{r}|=R)=1$ and $\unicode[STIX]{x1D6E9}_{0}(|\boldsymbol{r}|\rightarrow \infty )\rightarrow 0$ . The solution to this is $\unicode[STIX]{x1D6E9}_{0}=R/r\;\Longrightarrow \;T_{g}=T_{\infty }(\unicode[STIX]{x1D706}R/r+1)$ so that, under the low- $Ma$ approximation with constant thermodynamic pressure, the equation of state implies a density field prescription of

(3.12) $$\begin{eqnarray}p_{\infty }=\unicode[STIX]{x1D70C}R_{g}T_{g}\quad \;\Longrightarrow \;\quad \unicode[STIX]{x1D70C}=\frac{p_{\infty }}{R_{g}T_{\infty }}\,\frac{1}{\displaystyle \frac{\unicode[STIX]{x1D706}R}{r}+1}.\end{eqnarray}$$

With reference to (3.5), the quantity $\unicode[STIX]{x1D70C}^{(1)}$ may be approximated as $\unicode[STIX]{x1D70C}^{(1)}\approx -R/r$ . Assuming the base flow is Stokes flow, we also have the analytical solutions (Clift & Gauvin Reference Clift and Gauvin1970) for the velocity components $u_{r}$ and $u_{\unicode[STIX]{x1D703}}$ . Then, we obtain up to first order

(3.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}=-{\displaystyle \frac{\unicode[STIX]{x1D706}U_{\infty }R\cos \unicode[STIX]{x1D703}}{2r^{4}}}\left({\displaystyle \frac{R^{3}}{r}}-3rR+2r^{2}\right) & & \displaystyle\end{eqnarray}$$

so that $\unicode[STIX]{x1D701}$ scales as $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}\sim O(\unicode[STIX]{x1D706}U_{\infty }D^{3}/D^{4})=O(\unicode[STIX]{x1D706}U_{\infty }/D)=O(\unicode[STIX]{x1D706}Re\unicode[STIX]{x1D708}_{\infty }/D^{2})$ . This is consistent with the observation that $\unicode[STIX]{x1D701}$ decreases as $\unicode[STIX]{x1D706}$ and $Re$ decreases. Figure 10 shows the collapse of $\unicode[STIX]{x1D701}$ under this scaling at $Re=0.001$ and $0.01$ for small values of $\unicode[STIX]{x1D706}$ . This scaling does not hold as $Re$ or $\unicode[STIX]{x1D706}$ values increase. Contour plots of scaled $\unicode[STIX]{x1D701}$ at $Re=0.01$ and $0.1$ as $\unicode[STIX]{x1D706}$ varies from $0.1$ to $1$ are shown in figure 11. Another approach to obtain this scaling for $\unicode[STIX]{x1D701}$ motivated by a balance of terms in the governing equations at steady state is provided in Ganguli (Reference Ganguli2018). Using the methods of separation of variables and variation of parameters, the particular solution to $\unicode[STIX]{x1D719}$ satisfying (3.13) and its boundary conditions is $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D706}U_{\infty }R\cos \unicode[STIX]{x1D703}(-(7R^{3}/24r^{3})+(13R^{2}/16r^{2})-(3R/4r)+1/6)$ . Using the analytical solution for $\unicode[STIX]{x1D719}$ in the numerical evaluation of the excess pressure drag resulting from $p^{(1)}$ (governed by (3.6)) also agrees well with the numerically obtained values when $Re=0.001$ and $\unicode[STIX]{x1D706}=O(10^{-3})$ and lower. The derived expression also explains the observation in figure 7 where the drag increases almost linearly with $\unicode[STIX]{x1D706}$ as $Re\rightarrow 0$ and $\unicode[STIX]{x1D706}\rightarrow 0$ .

Figure 10.

Collapse of the scaled rate of dilatation, $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})/(\unicode[STIX]{x1D706}Re(\unicode[STIX]{x1D708}_{\infty }/D^{2}))$ , at small $Re$ ( $Re=0.001$ and $0.01$ ) and small $\unicode[STIX]{x1D706}$ ( $\unicode[STIX]{x1D706}=0.01$ , 0.02, 0.03 and 0.04).

3.5 Comparison with compressible low- $Re$ flow

The governing equations solved in this paper are the low- $Ma$ , variable-density Navier–Stokes equations. When the flow is compressible, $Ma$ effects would also need to be considered. The governing parameters of the problem then become $Ma$ , $Re$ and $\unicode[STIX]{x1D706}$ . As mentioned in the introduction, Kassoy et al. (Reference Kassoy, Adamson and Messiter1966) derived analytical solutions for compressible low- $Re$ flow around a sphere in the regimes of $\unicode[STIX]{x1D706}=O(Re)$ and $\unicode[STIX]{x1D706}=O(1)$ . It is an interesting exercise to qualitatively and quantitatively compare the results from the low- $Ma$ , variable-density case presented in this paper with those from the compressible case. It is instructive to note that in Kassoy et al. (Reference Kassoy, Adamson and Messiter1966), when $\unicode[STIX]{x1D706}=O(Re)$ , the expansions are carried out far enough to include the first Mach number effect. On the other hand, when $\unicode[STIX]{x1D706}=O(1)$ , only the zeroth-order term is calculated owing to the complexity of the equations. Note that the resulting relationship for drag in this strong heating regime is independent of $Re$ in the compressible case while, in the low- $Ma$ , variable-density case, there is a weak $Re$ dependence. Table 6 presents this comparison. Quantitatively, we observe that the deviation of the drag computed in the compressible case from the drag computed in the low- $Ma$ , variable-density case is not very large. This means that compressibility effects are present but are not significant in this regime. When $\unicode[STIX]{x1D706}=O(1)$ , the deviation is slightly larger owing to strong compressibility effects. However, the qualitative trends observed in the low- $Ma$ , variable-density case are consistent with those in the compressible case. Specifically, both the drag and Nusselt number are greatly affected by the impressed temperature difference $\unicode[STIX]{x1D706}$ . The pressure and viscous stress distributions are also altered considerably from the corresponding incompressible values. The drag coefficient increases almost linearly with $\unicode[STIX]{x1D706}$ ; its value when $\unicode[STIX]{x1D706}$ is near unity is substantially greater than the classical Stokes value.

Figure 11.

Comparison of the contours of the scaled rate of dilatation, $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})/$ $(\unicode[STIX]{x1D706}Re(\unicode[STIX]{x1D708}_{\infty }/D^{2}))$ , at $Re=0.01$ (a–c) and $0.1$ (d–f) as the value of $\unicode[STIX]{x1D706}$ varies from 0.1 to 1.

Table 6.

Comparison of the drag on a heated sphere in a low-Mach-number variable-density flow with that of the heated sphere in compressible flow as computed in Kassoy et al. (Reference Kassoy, Adamson and Messiter1966).

3.6 An analogy to incompressible flow over a spheroid

The drag experienced by a non-spherical body in a low-Reynolds-number incompressible flow has been successfully calculated for a variety of shapes via analytical, semi-analytical and numerical methods. Some examples include an ellipsoid (Oberbeck Reference Oberbeck1876), lens (Payne & Pell Reference Payne and Pell1960), hemisphere (Payne & Pell Reference Payne and Pell1960) and hemispherical cap (Collins Reference Collins1963). When bodies are volumes of revolution of an axisymmetric locus about the axis of symmetry, it makes solving the Stokes flow more amenable to mathematical analysis. Solutions using analytical methods for slightly distorted spheroids have also been worked out (Brenner Reference Brenner1964; Taylor & Acrivos Reference Taylor and Acrivos1964). However, for moderately or largely distorted spheroids created by deforming a sphere of radius $R$ , exact solutions are known only for ellipsoids (Oberbeck Reference Oberbeck1876; Payne & Pell Reference Payne and Pell1960). Recall that the Stokes drag on an undeformed sphere of radius $R$ is $6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{\infty }U_{\infty }R$ and is valid when the fluid has constant density, flow is incompressible and $Re\rightarrow 0$ . Semi-analytical solutions for the drag in Stokes flow for a wide variety of mathematical shapes have been obtained by O’Brien (Reference O’Brien1965a ,Reference O’Brien b , Reference O’Brien1966, Reference O’Brien1968). Results are expressed in terms of the drag of the deformed body ( $F_{DB}$ ), drag ratio ( $r_{D}$ ) and form factor ( $FF$ ) defined as $r_{D}=F_{DB}/(6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{\infty }U_{\infty }R)$ and $FF=r_{D}-1$ . If the volume of the deformed sphere (spheroid) is denoted as $V_{DS}$ and the volume of the sphere is denoted as $V$ , O’Brien (Reference O’Brien1968) demonstrates the collapse of the form factors obtained numerically for a wide variety of body shapes over a wide range of deformations in the Stokes limit into the empirical correlation $FF=\frac{1}{5}(V_{DS}-V)/V$ . It is seen that $FF$ for most of the shapes falls within $\pm 0.01$ (1 % of drag ratio) of this linear relationship. Figure 4 motivates an analogy between the increased drag on a spheroid when compared to a sphere and the increased drag on the spherical particle when heated. As $\unicode[STIX]{x1D706}$ increases, the flow is slower and has a spatially wider scale of variation in front of the sphere whereas the flow is faster and has a spatially smaller scale of variation behind the sphere. In other words, the effective shape of the sphere that the fluid ‘sees’ changes to a pear-shaped or egg-shaped spheroid where the base of the pear/egg is orthogonal to the direction of the free-stream velocity, the base serves as the forward stagnation region, and the axis of symmetry of the pear is aligned parallel to the direction of the free-stream flow. The dimensions of this spheroid can be calculated by equating the form factor obtained from (3.8) and the correlation suggested by O’Brien (Reference O’Brien1968). Let $R_{eq}$ be the radius of the sphere whose volume is equal to the volume of the pear-shaped spheroid ( $V_{eq}$ ) such that the base of the pear is on the forefront. Let $V$ be the volume of the undeformed spherical particle with radius $R$ . Equating the form factor suggested by O’Brien with the form factor for the heated sphere obtained from (3.8) leads to

(3.14) $$\begin{eqnarray}\displaystyle FF={\displaystyle \frac{C_{DV}}{6}}-1=0.1315Re^{0.82-0.05\log _{10}Re}+0.4486\unicode[STIX]{x1D706}^{0.9673}Re^{0.0471}={\displaystyle \frac{1}{5}}{\displaystyle \frac{V_{eq}-V}{V}}. & & \displaystyle\end{eqnarray}$$

The equivalent radius (radius of a sphere such that its volume is equal to the volume of the pear-shaped body) of the spheroid inferred from this equality can be written as

(3.15) $$\begin{eqnarray}\displaystyle R_{eq}=(1+0.6575Re^{0.82-0.05\log _{10}Re}+2.2430\unicode[STIX]{x1D706}^{0.9673}Re^{0.0471})^{1/3}R. & & \displaystyle\end{eqnarray}$$

In conclusion, the volume expansion of the fluid associated with heat transfer from the particle to the fluid makes the particle appear aspherical and to have a larger volume; the fore–aft asymmetry causes the effective shape of the body to have a wider head than its tail. The volume of this effective shape is characterized in (3.15). This qualitative change in the particle shape is consistent with the increased drag.