2.1. Governing equations

The NSE for a perfect gas are written in conservative form as

(2.1) \begin{equation} \frac {\partial \boldsymbol{U}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot } ( \boldsymbol{F^c} - \boldsymbol{F^v} ) = 0, \end{equation}

where the solution vector is

(2.2) \begin{equation} U = \left [ \begin{matrix} \rho \\ \rho \boldsymbol{u} \\ E \end{matrix} \right ] \end{equation}

and the convective and viscous fluxes respectively are

(2.3) \begin{equation} F^c = \left [ \begin{matrix} \rho \boldsymbol{u} \\ \rho \boldsymbol{u} \otimes \boldsymbol{u} + \boldsymbol{I}\! P \\ (E+P) \boldsymbol{u} \end{matrix} \right ]\!, \end{equation}

(2.4) \begin{equation} F^v = \left [ \begin{matrix} 0 \\ \boldsymbol{\sigma } \\ \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\sigma } - \boldsymbol{q} \end{matrix} \right ]\!, \end{equation}

where $\rho$ is density, $ \boldsymbol{u}$ is the velocity vector, $P$ is pressure, $\boldsymbol{I}$ is the identity matrix, $E=\rho c_{v} T + \rho ( { \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u}}/{2})$ is the total energy, where $T$ is temperature and $c_v$ is the specific heat at constant volume, and $ \boldsymbol{\sigma }=2 \mu \boldsymbol{S} - ({2}/{3}) \mu \boldsymbol{I}\! \boldsymbol{S}$ is the viscous stress tensor, where $\mu$ is the dynamic viscosity of the fluid, $ \boldsymbol{S}$ is the strain tensor and $ \boldsymbol{q}=\kappa ({\partial T}/{\partial \boldsymbol{x}})$ is the heat flux vector with thermal conductivity $\kappa$ . Pressure is defined by the ideal gas equation of state $P=\rho \textit{RT}$ , where the specific gas constant for air is $R=287.051 \rm \: J \,kg^{-1} \,K^{-1}$ . Dynamic viscosity is calculated using Sutherland’s law (Sutherland Reference Sutherland1893) as $\mu =1.456 \times 10^{-6} ({T^{3/2}}/{T+110.3})$ . Thermal conductivity is calculated as $\kappa = {\mu c_{\textit{p}}}/\textit{Pr}$ , where $\textit{Pr}=0.72$ is the Prandtl number for air and the ratio of specific heats is $\gamma ={c_{\textit{p}}}/{c_v}=1.4$ .

2.2. Numerical scheme

For cases with steady flow, convective fluxes are calculated using the Modified-Steger–Warming (MSW) upwinding flux scheme (MacCormack & Candler Reference MacCormack and Candler1989). Viscous fluxes are calculated using weighted least-squares-error gradients with a deferred-correction approach (Kim, Makarov & Caraeni Reference Kim, Makarov and Caraeni2003). Time is advanced using the data parallel line relaxation (DPLR) time stepping scheme near walls (Wright, Candler & Bose Reference Wright, Candler and Bose1998) with the full matrix data parallel (FMDP) method (Wright, Candler & Prampolini Reference Wright, Candler and Prampolini1996) in regions of the flowfield not near a wall.

Cases with unsteady behaviour are simulated using a low-dissipation numerical method. The details of the method are beyond the scope of this work, so we provide only a brief summary here. The reader is encouraged to refer to the cited works for details. Time stepping is performed using an implicit second-order backward difference (BDF2) scheme. The convective flux $F^c$ is calculated from a fourth-order central kinetic-energy-consistent (KEC) flux $F^k$ (Subbareddy & Candler Reference Subbareddy and Candler2009) and a dissipative MSW upwind flux $F^d$ (MacCormack & Candler Reference MacCormack and Candler1989) as

(2.5) \begin{equation} F^c=F^k + \alpha\! F^d, \end{equation}

where the switch $\alpha$ determines the level of dissipation in the flux calculation.

For supersonic flows, $\alpha$ is set using a shock sensor. The sensor takes the form of a Larsson switch (Larsson & Lele Reference Larsson and Lele2009):

(2.6) \begin{equation} \alpha = \max \left [ \alpha _m, \min \left [2\frac {-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}}{|\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}|\, + \parallel\!\boldsymbol{\nabla }\times \boldsymbol{u}\!\parallel +\, \epsilon \frac {u_s}{\delta } }, 1 \right ] \right ]\!, \end{equation}

where $\alpha$ is the sensor value, $\epsilon$ is the tuning parameter, $u_s$ is the free stream velocity on the stagnation streamline, $\delta$ is a normalisation parameter defined as the approximate boundary layer thickness and $\alpha _m$ is the minimum value for the sensor to maintain numerical stability (set to 0.05 in this work). Here, $\epsilon$ is tuned so that $\alpha$ is highest inside shocks and low everywhere else.

For subsonic flows, a numerical sponge layer is created by setting $\alpha$ to be high near inflow/outflow boundaries and low everywhere else using the logistic function as

(2.7) \begin{equation} \alpha = \max \left [\alpha _m, \frac {1}{1 + \textrm {e}^{-2 k_h \left (L - \epsilon \right )}} \right ]\!, \end{equation}

where $L$ is the distance from face centroid to the nearest inflow/outflow, $k_h$ is the sharpening parameter which determines how quickly $\alpha$ increases, and $\epsilon$ is the distance from the boundary below which $\alpha$ increases to $1$ . As for supersonic cases, $\alpha _m$ is set to $0.05$ .

3.1. Computational grid

An overview of the subsonic computational grid is shown in figure 3. The grid extends 10 diameters upstream, 8 diameters in the free stream-normal direction and 15 diameters downstream from the particle centre. The supersonic grid takes the same form as the subsonic grid, but extends only 2 diameters upstream. The number of elements in each grid is listed in table 1.

Grid convergence is assessed for the subset of cases shown in table 2 using the lift and drag coefficients. Figure 4 plots the magnitude of the computed lift and drag coefficient for subsonic and supersonic flows. For subsonic flows, the medium and fine grids produce lift and drag coefficients that agree within 6 %, indicating that the medium grid is converged for subsonic flows. For supersonic flows, the coarse and medium grids produce lift and drag coefficients that agree within 2 %, indicating that the coarse grid is converged for supersonic flows. Furthermore, the measured drag coefficient agrees with at least one known drag model within 6 % in all cases except for a low supersonic case where $M_s=1.2$ , where the agreement is still within 20 %. Some deviation in this regime is not unexpected, as drag models in this regime are less accurate as compared with the subsonic and high supersonic flow regimes (Nagata et al. Reference Nagata, Nonomura, Takahashi and Fukuda2020).

The size of the computational domain can also influence the computed lift coefficient (Shi & Rzehak Reference Shi and Rzehak2019), especially for low-Mach-number flows. To ensure sufficient domain size, run A.23 from table 2 was re-run with a domain wherein the distance between the sphere and the domain boundaries was doubled. This changed the lift coefficient by only $3\,\%$ , indicating that the domain size is sufficient.

3.2. Lift coefficient scaling with non-dimensional dynamic pressure gradient

In this section, we study the direction of lift force and the magnitude of the lift coefficient relative to the dynamic pressure gradient. Drag force is always in the direction of the free stream flow, which by convention is in the $x$ -direction. Lift force is normal to the drag force. This means that the direction of lift force is not predetermined, since any vector with zero $x$ component is normal to the drag force. Therefore, we must determine the lift force direction ourselves.

The directional lift coefficient for all steady-flow cases in Appendix A.1 is shown in figure 5. Lift coefficient is plotted in the $y$ – $z$ plane; normal to the $x$ -direction free stream flow. The $y$ -axis is aligned with the direction of the dynamic pressure gradient. The $z$ -axis is in the cross-gradient direction. As shown in figure 5, the lift coefficient is only in the negative $y$ -direction. There is negligible lift in the $z$ -direction. This shows that shear-induced lift force is in the direction opposite the dynamic pressure gradient, with negligible cross-gradient lift. This means that particles tend to be lifted from regions of higher dynamic pressure to regions of lower dynamic pressure.

Figure 5. Lift polar for a spherical particle in a dynamic pressure gradient.

Next, we show that the lift coefficient is a function of the non-dimensional dynamic pressure gradient using the subset of cases in Appendix A.2. Figure 6 plots the lift coefficient for a non-dimensional dynamic pressure gradient $\widetilde {\boldsymbol{\nabla\! }Q} = 10^{-3}$ created from only a density gradient ( $\boldsymbol{\nabla}\kern-1pt \rho _{\infty }$ ), only a velocity gradient ( $\boldsymbol{\nabla }u_{\infty }$ ), and both density and velocity gradients ( $\boldsymbol{\nabla}\kern-1pt \rho _{\infty } + \boldsymbol{\nabla }u_{\infty }$ ). We have shown already that the lift coefficient is only in the negative $y$ -direction, so we plot the negative of the $y$ -direction lift coefficient. The total lift coefficient $C_{L,y}$ is plotted, as well as its pressure and viscous components $C_{L,y,P}$ and $C_{L,y,\tau }$ . We see that pressure lift dominates viscous lift at $\textit{Re}_{\textit{p}} \geqslant 200$ and pressure lift is similar in magnitude to viscous lift at $\textit{Re}_{\textit{p}}=100$ . This is consistent with the observations of Shi & Rzehak (Reference Shi and Rzehak2019). When the dynamic pressure gradient is created with only a velocity gradient, the measured lift coefficient at low Mach number ( $M_s=0.3,{} \textit{Re}_{\textit{p}}=100$ ) is within the variation of the DNS used to develop the shear lift correlation for incompressible flow proposed by Shi & Rzehak (Reference Shi and Rzehak2019).

As shown in figure 6, at $M_s \geqslant 0.8$ , $\textit{Re}_{\textit{p}} \geqslant 200$ , both the total lift coefficient and its pressure and viscous components do not change regardless of whether the dynamic pressure gradient was created from a density and/or velocity gradient. At $M_s \leqslant 0.6$ , $\textit{Re}_{\textit{p}} \leqslant 100$ , the lift coefficient starts to vary based on how the dynamic pressure gradient was created, though the variation is small. This shows that lift is dependent on the dynamic pressure gradient, not on the density or velocity gradient.

The physical mechanism behind the change in behaviour at lower Mach and Reynolds number is nuanced and requires careful examination. As Reynolds number decreases, viscous forces in the flowfield increase in strength relative to inertial (pressure) forces. This affects the behaviour of the entire flow, not just the boundary layer – so the changing Reynolds number affects both pressure and viscous lift at the particle surface. Recall that in this study, flows with a density gradient come with an opposing temperature gradient such that zero static pressure gradient is maintained in the free stream. A temperature gradient in the free stream indicates a viscosity gradient in the free stream (Sutherland Reference Sutherland1893). This means that for runs with a density gradient, there is an additional viscosity gradient affecting the behaviour of the flow. As Reynolds number decreases and viscous forces increase in strength relative to inertial forces, this ‘second’ gradient due to viscosity contributes further to asymmetries in the flowfield, as compared with runs with only a velocity gradient. This results in stronger lifting effects for particles in density gradients. In figure 6 (for which all runs have the same non-dimensional dynamic pressure gradient), this effect is reflected in stronger lift at $\textit{Re}_{\textit{p}}=100$ when the dynamic pressure gradient is created with a density gradient (red) than when the dynamic pressure gradient is created with a velocity gradient (blue); wherein runs with both density and velocity gradients (green) fall between runs with only a density or velocity gradient.

Now that we have shown lift coefficient depends on the dynamic pressure gradient, we must determine the scaling of lift coefficient magnitude with dynamic pressure gradient. We use the subset of cases in Appendix A.2. In figure 7, we plot lift coefficient versus non-dimensional dynamic pressure gradient; wherein the lift coefficient is normalised by its value at $\widetilde {\boldsymbol{\nabla\! }Q} = 10^{-3}$ . As shown in figure 7, lift coefficient scales linearly with non-dimensional dynamic pressure gradient. For all Mach and Reynolds numbers except at $M_s=0.8$ , $\textit{Re}_{\textit{p}}=200$ , if non-dimensional dynamic pressure gradient doubles, lift coefficient doubles. This relationship between $\widetilde {\boldsymbol{\nabla\! }Q}$ and $C\!_L$ indicates that as $\widetilde {\boldsymbol{\nabla\! }Q}$ approaches zero, $C\!_L$ also approaches zero. This makes sense, because if $\widetilde {\boldsymbol{\nabla\! }Q}=0$ , there is no shear in the flow, so we would expect there to be no shear-induced lift.

Figure 6. Lift coefficient in density and/or velocity gradients at $\widetilde {\boldsymbol{\nabla\! }Q} = 10^{-3}$ .

Figure 7. Scaling of lift coefficient with $\widetilde {\boldsymbol{\nabla\! }Q}$ .

Special behaviour is observed for the case at $M_s=0.8$ , $\textit{Re}_{\textit{p}}=200$ . This is because this combination of Mach and Reynolds number results in a planar-symmetric wake, whereas the other cases in figure 7 have axisymmetric wakes. For planar-symmetric wakes, some lift is induced along the symmetry plane (Nagata et al. Reference Nagata, Nonomura, Takahashi and Fukuda2020) due to counter-rotating vortices in the wake, even when $\widetilde {\boldsymbol{\nabla\! }Q}=0$ . The symmetry plane would tend to be randomly oriented when $\widetilde {\boldsymbol{\nabla\! }Q}=0$ . However, in the range of $\widetilde {\boldsymbol{\nabla\! }Q}$ considered in this work, the symmetry plane tends to align with the gradient, amplifying shear-induced lift (see § 3.4). This means that for planar-symmetric wakes, as $\widetilde {\boldsymbol{\nabla\! }Q} \rightarrow 0$ , nonlinear lifting effects are induced as the orientation of the symmetry plane changes from gradient-aligned to randomly oriented. However, figure 7 shows that even for planar-symmetric wakes, $C\!_L$  is linear with $\widetilde {\boldsymbol{\nabla\! }Q}$ in the range of $\widetilde {\boldsymbol{\nabla\! }Q}$ studied here.

This assessment of lift direction and magnitude is consistent with simple scaling arguments. Lift force is by definition normal to the free stream flow, so we ignore forces in the $\hat {i}$ direction as they would be included in drag force. We align the dynamic pressure gradient with the $y$ -axis, so in steady flow, there are no asymmetric forces with a component in the $z$ -direction. Therefore, lift force is one-dimensional along the $y$ -axis:

(3.10) \begin{equation} F_{L,y} = \int _{A}^{} ( P_w \sin {\theta } \cos {\phi } + \tau _{w,y} ) \,{\rm d}A, \end{equation}

where $P_w$ is the pressure force at the wall and $\tau _{w,y}$ is shear stress at the wall in the $y$ -direction. The angles $\theta$ and $\phi$ are defined in figure 2.

Consider first the lift force due to pressure. There exists no analytic model for surface pressure in a general compressible flow. However, Newtonian aerodynamics (Heybey Reference Heybey1964) provides a crude estimate of surface pressure at high Mach number. Newtonian aerodynamics neglects pressure forces on the downstream half of the sphere as small and approximates pressure coefficient on the upstream half of the sphere as

(3.11) \begin{equation} C_P \simeq 2 \cos ^2{\theta }. \end{equation}

Using the definition of pressure coefficient and the definition of dynamic pressure, we can approximate the wall pressure as

(3.12) \begin{equation} P_w \simeq P_{\infty } + Q_{\infty } \: 2 \cos ^2{\theta }. \end{equation}

Recall we have aligned our axes such that the free stream dynamic pressure gradient is parallel to and in the direction of the $y$ -axis, so $\boldsymbol{\nabla \!}Q_{\infty }$ only has one dimension and can be represented by the scalar $||\boldsymbol{\nabla \!}Q_{\infty }||$ . Therefore, the free stream dynamic pressure at a given $y$ -coordinate is

(3.13) \begin{equation} Q_{\infty } = Q_s - \frac {D}{2} ||\boldsymbol{\nabla \!}Q_{\infty }|| \sin {\theta } \cos {\phi }, \end{equation}

where $D$ is the diameter of the spherical particle. Recall also we assume no static pressure gradient in the free stream, so $P_{\infty } = P_s$ everywhere in the free stream. Thus, wall pressure scales as

(3.14) \begin{equation} P_w \sim P_{\infty } + 2 \left ( Q_s - \frac {D}{2} ||\boldsymbol{\nabla \!}Q_{\infty }|| \sin {\theta } \cos {\phi } \right ) \cos ^2{\theta }. \end{equation}

By substituting into the first term in (3.10) and evaluating the integral, we see that

(3.15) \begin{equation} F_{L,y,P} \sim -\frac {D^3}{4} \frac {2\pi }{15} ||\boldsymbol{\nabla \!}Q_{\infty }||. \end{equation}

The negative sign on the right-hand side in (3.15) means that lifting force is in the direction opposite the dynamic pressure gradient – i.e. particles will tend to be lifted down the gradient, from regions of higher dynamic pressure to regions of lower dynamic pressure.

Applying the definition of lift coefficient, using the definition of non-dimensional dynamic pressure gradient from (3.5), taking the reference area for a sphere to be $A_{\textit{ref}}=\pi R^2=\pi ( {D^2}/{4})$ and removing constant factors, we see that

(3.16) \begin{equation} C_{L,y,P} \sim -\widetilde {\boldsymbol{\nabla\! }Q}, \end{equation}

which shows that lift due to pressure forces scales with non-dimensional dynamic pressure gradient and is in the direction opposite the gradient. This is consistent with the DNS results.

While Newtonian aerodynamics can provide some limited insight into pressure forces on the upstream portion of the sphere, it does not provide any insight into the viscous forces on the upstream portion of the sphere, nor does it provide insight into the evolution of the vorticity field or forces on the downstream portion of the sphere. We explore in § 3.4 how wake structure influences lift coefficient, but no analytic map back to the vorticity field or viscous stresses exists for a general compressible flow.

3.3. Effects of Mach Number, Reynolds number and wall temperature ratio

In this section, we study effects of changing Reynolds number, Mach number and wall temperature ratio by varying each of these parameters independently. We consider a subsonic baseline case at $\textit{Re}_{\textit{p}}=100, \: M_s=0.3$ and a supersonic baseline case at $\textit{Re}_{\textit{p}}=1000, \: M_s=2.0$ . Both baseline cases use $ {T_w}/{T_{aw}}=0.5$ . Case details are listed in Appendix A.3. Figure 8 shows how increasing each parameter by a factor of 2 affects both pressure and viscous lift by plotting the lift coefficient against each parameter, normalised by their value for the baseline case. As shown in figure 8, Reynolds number has the largest effect on both pressure and viscous lift, especially at low Mach number. This is because changing the Reynolds number results in a change in wake structure (Nagata et al. Reference Nagata, Nonomura, Takahashi and Fukuda2020). Discussion on this matter can be found in § 3.4. Mach number has a smaller but non-negligible effect, primarily on the pressure lift coefficient at lower Reynolds number. Effects of wall temperature ratio are even smaller, especially for pressure lift, which drives the overall lift coefficient in an unbounded compressible flow.

These conclusions are consistent with parameters that drive drag forces. Drag coefficient on a spherical particle in continuum flow is a function of Mach number and Reynolds number, with wall temperature effects becoming important only for rarefied flows (Singh et al. Reference Singh, Kroells, Li, Ching, Ihme, Hogan and Schwartzentruber2020; Loth et al. Reference Loth, Despit, Jeong, Nagata and Nonomura2021). Thus, we would also expect Mach and Reynolds number to drive lift coefficient for the continuum flow studied in this work.

Figure 8. Effects of $\textit{Re}_{\textit{p}}$ , $M_s$ and $ {T_w}/{T_{aw}}$ on pressure lift (top row) and viscous lift (bottom row).

3.4. Effects of wake structure on lift coefficient

In this section, we examine the effect of wake structure on lift coefficient. Both steady and unsteady wake structures significantly affect lift.

Unsteady wakes can develop cross-gradient asymmetries. When the cross-gradient asymmetric behaviour is stronger than the dynamic pressure gradient, cross-gradient lifting forces in directions not readily predictable are induced. While this behaviour is not intuitive, it is well documented. Williams & Smits (Reference Williams and Smits2024) compiled dozens of observations regarding asymmetric behaviour in nominally symmetric flows. They concluded that asymmetric behaviours can and do exist, mainly because they evolve on time scales much longer than are currently possible to observe in most simulations or experiments.

Figure 9. Wakes at $M_s=0.3$ and $\widetilde {\boldsymbol{\nabla\! }Q}=10^{-3}$ : top, steady axisymmetric (SA); middle, steady planar-symmetric (SP); bottom, unsteady hairpin with azimuthal oscillations (HaWAO).

Figure 10. Effect of wake structure on lift coefficient for seven cases at $\widetilde {\boldsymbol{\nabla\! }Q}=10^{-3}$ , $\textit{Re}_{\textit{p}}$ from $100$ to $2000$ , $M_s$ from $0.3$ to $2.0$ .

Steady-wake structure also influences lift. We observe that the symmetry plane for planar-symmetric wakes tends to align with the direction of the dynamic pressure gradient. This creates gradient-aligned viscous forces on the downstream half of the sphere which would otherwise be axisymmetrically balanced. Therefore, at low Mach number when viscous forces have a significant influence on the lift coefficient, planar-symmetric wakes result in higher-magnitude lift coefficients than for axisymmetric wakes.

For example, consider the wake structures shown in figure 9. At $\textit{Re}_{\textit{p}}=100$ , the wake is steady axisymmetric (SA). (Wake classification nomenclature from Nagata et al. Reference Nagata, Nonomura, Takahashi and Fukuda2020.) At $\textit{Re}_{\textit{p}}=200$ , the wake is steady planar-symmetric (SP) across the $x$ – $y$ plane, which is aligned with the dynamic pressure gradient. At $\textit{Re}_{\textit{p}}=600$ , the wake becomes an unsteady hairpin wake with azimuthal oscillations (HaWAO). All of these structures match those predicted by Nagata et al. (Reference Nagata, Nonomura, Takahashi and Fukuda2020).

In figure 10, we plot the lift polar for three sets of cases. The first set consists of the same three cases shown in figure 9 at $M_s=0.3$ . The other sets are pairs of cases at $M_s=0.6$ and $M_s=2.0$ . One member of each pair is an SA wake and the other is an SP wake. Case details are listed in Appendix A.4. All six steady-flow cases show lift in the expected negative $y$ -direction with negligible cross-gradient ( $z$ -direction) lift. However, at $M_s=0.3$ and $M_s=0.6$ , the SP case has a lift coefficient nearly an order of magnitude higher than the SA case. This trend holds at $M_s=2$ , but is much less pronounced, since lift coefficient at high Mach number is dominated by pressure forces on the upstream half of the sphere and is therefore less affected by wake behaviour. The unsteady case demonstrates significant non-zero-average cross-gradient lift, which is an asymmetric behaviour across the $x$ – $y$ plane in a flow that is nominally symmetric across that plane. The unsteady case also demonstrates lift in the positive $y$ -direction – opposite the direction predicted in § 3.2. This shows that unsteady wakes drive lift force in directions not readily predicted. The changing lift coefficient with changing wake structure means that predicting lift requires predicting wake structure.

4.1. Computational grid modifications

For subsonic and low-supersonic flows, the computational grid take the same form as in § 3.1. However, for high-supersonic flows, the use of the low-dissipation numerical scheme with a shock sensor described in § 2.2 requires grid-shock alignment to minimise flickering of the shock between cells, which introduces numerical noise into the simulation (Knutson, Thome & Candler Reference Knutson, Thome and Candler2021). An overview of the grid is shown in figure 11. As shown in figure 12, the grid is approximately aligned with the shock using shock shape equations from Bedarev, Fedorov & Fomin (Reference Bedarev, Fedorov and Fomin2012):

(4.1) \begin{equation} \begin{array}{ll} R = \dfrac {D}{2} \left ( 1 + \dfrac {1}{M-1} \right )\!, \\[2pt]\\ a = \textit{RM} ( M-1 ),\\[2pt] \\ b = R \sqrt {M ( M-1)}, \\[2pt]\\ \varDelta = \dfrac {D}{2} \left ( 1 + \dfrac {M}{2 ( M^2 -1) } \right )\!. \end{array} \end{equation}

Figure 11. Unsteady supersonic grid overview, $M_s=2.4$ (heavily coarsened for visualisation purposes).

Figure 12. Sample shock-aligned grid with contours of the shock sensor value $\alpha$ , $M_s=2.4$ .

The number of elements in each grid is listed in table 3.

Table 3. Grid sizes for unsteady/steady simulations.

Grid convergence for the subset of cases shown in table 4 is assessed using the drag coefficient over a window of more than five unsteady periods. Mean drag coefficient is compared with existing models in figure 13. The coarse, medium and fine grids produce drag coefficients which match each other within 2 % and match at least one known drag model within 7 %. Frequency content of the drag coefficient time-history, generated using a fast Fourier transform (FFT), is shown in figure 14. The medium grid and fine grid produce frequency spectra with a maximum deviation of less than 5 %. The close agreement in both mean drag and unsteady drag frequency content indicates that the medium grid is converged.

Table 4. Convergence test set for delineating steady versus unsteady behaviour.

Figure 13. Mean drag force versus existing drag models.

Figure 14. Drag coefficient frequency spectrum in the Strouhal number domain, $St=f( {D}/{u_s})$ .

Furthermore, wake structures observed in this work at $\textit{Re}_{\textit{p}} \leqslant 1000$ match those predicted by Nagata et al. (Reference Nagata, Nonomura, Takahashi and Fukuda2020). Figure 15 shows wake structures over a range of Mach and Reynolds numbers. At $\textit{Re}_{\textit{p}}=500$ and $M_s=0.8$ , the wake demonstrates hairpin structures with azimuthal oscillation (HaWAO). At $\textit{Re}_{\textit{p}}=1000$ and $M_s=0.8$ , the wake demonstrates helical structures (HeW). At $\textit{Re}_{\textit{p}}=1000$ and $M_s=1.2$ , we observe one of the types of pure hairpin wakes (HaW2). The wake is steady and axisymmetric (SA) at $\textit{Re}_{\textit{p}}=1000$ and $M_s=2.0$ .

Figure 15. Sample wake structures in steady and unsteady flows.

4.2. Results

We use data from the prior work and the DNS to develop approximate delineations among steady axisymmetric, steady planar-symmetric and unsteady wakes. Delineations among wake types are represented using critical Mach numbers, each of which is a function of Reynolds number. A critical Mach number $M_{\textit{ts}}$ which delineates steady and unsteady behaviours is

(4.2)

where if $M_s \gt M_{\textit{ts}}$ , the wake will likely assume a steady behaviour and will likely be unsteady otherwise. Here, $\textit{Re}_{\textit{ts}}=1395$ is the Reynolds number at which the delineation switches from the linear Nagata et al. (Reference Nagata, Nonomura, Takahashi and Fukuda2020) curves to the proposed logarithmic curve; found by value- and slope-matching the logarithmic curve to the linear curve to eliminate discontinuity (within rounding error) between the two. The delineation is shown graphically in figure 16, along with the continuum flow limit defined by the Knudsen number $\textit{Kn} \leqslant 0.01$ and related to the Mach and Reynolds numbers by

(4.3) \begin{equation} M_s \leqslant 0.01 \textit{Re}_{\textit{p}} \sqrt {\frac {2}{\pi \gamma }}. \end{equation}

Note that for consistency with the prior work, we develop this delineation with $\widetilde {\boldsymbol{\nabla\! }Q} = 0$ . However, as shown by the cases at $M_s=2, \: \textit{Re}_{\textit{p}}=1000$ , $M_s=2, \: \textit{Re}_{\textit{p}}=2000$ and $M_s=4,{} \: \textit{Re}_{\textit{p}}=10\,000$ , the wake structure does not change when a gradient is introduced.

Figure 16. Approximate delineation between wake behaviours in continuum flow.

A critical Mach number $M_{\textit{tp}}$ which delineates steady axisymmetric (SA) and steady planar-symmetric (SP) wakes is

(4.4) \begin{align} M_{\textit{tp}} \simeq \left \{ \begin{array}{ll} 0, & \textit{Re}_{\textit{p}} \lt 175, \\ 0.0115\textit{Re}_{\textit{p}}-1.42, & 175 \leqslant \textit{Re}_{\textit{p}} \lt 210, \\ \big ( 0.72\textit{Re}_{\textit{p}}+846\big ) \times 10^{-3}, & 210 \leqslant \textit{Re}_{\textit{p}} \lt \textit{Re}_{\textit{tp}},\\ 0.78\ln {\textit{Re}_{\textit{p}}}-3.8243, & \textit{Re}_{\textit{p}} \geqslant \textit{Re}_{\textit{tp}}, \end{array} \right . \end{align}

where if $M_s \leqslant M_{\textit{tp}}$ and $M_s \gt M_{\textit{ts}}$ , the wake will likely be steady planar-symmetric, and steady axisymmetric if $M_s \gt M_{\textit{tp}}$ . Here, $\textit{Re}_{\textit{tp}}=1083$ is calculated by value- and slope-matching the linear curve to the logarithmic curve. This delineation is also shown in figure 16. Note that we observe a slight increase in $M_{\textit{tp}}$ at low Reynolds number as compared with the delineation drawn by Nagata et al. (Reference Nagata, Nonomura, Takahashi and Fukuda2020). This is not unexpected, as some variance around this delineation also occurs in the many works cited by Nagata et al. (Reference Nagata, Nonomura, Takahashi and Fukuda2020).

It is important to acknowledge that the proposed delineations should be treated as approximations only. Various complex factors including noise in the free stream (for experimental studies), properties of the numerical scheme (for computational studies), and other considerations can have non-negligible impacts on the Mach and Reynolds number at which wake behaviour changes. This is reflected in inconsistencies around delineations between wake types in the prior work. Note also that increasing wall temperature ratio tends to delay the transition from steady to unsteady wake behaviour as Reynolds number increases for a given Mach number (Nagata et al. Reference Nagata, Nonomura, Takahashi, Mizuno and Fukuda2018a ). Since this work uses a relatively low wall temperature ratio ( ${T_w}/{T_{aw}}=0.5$ ), the delineations developed herein would be considered conservative for the purposes of a steady-flow lift model, because wakes are more likely to be steady for a given Mach and Reynolds number for higher wall temperature ratios.

5.1. Pressure lift

First, we develop a high-Mach limit for pressure lift coefficient. The scaling arguments in § 3.2 show that in the high-Mach limit of Newtonian aerodynamics, pressure lift coefficient scales with non-dimensional dynamic pressure gradient. This would indicate that

(5.4) \begin{equation} \lim _{M_s\to \infty } C\!_{L,p}= \widetilde {\boldsymbol{\nabla \!}Q_{\perp }}. \end{equation}

However, in supersonic flows, a bow shock will form in front of the particle. The shock will change the effective dynamic pressure gradient seen by the particle. We must therefore account for the presence of a shock in our dynamic pressure gradient calculation. Dynamic pressure can be written as

(5.5) \begin{equation} Q = \frac {1}{2} \rho M^2 \gamma \textit{RT} .\end{equation}

For a perfect gas with constant $\gamma$ and $R$ , the change in dynamic pressure across a shock is then

(5.6) \begin{equation} \frac {Q_{2}}{Q_{1}} = \frac {\rho _{2}}{\rho _{1}} \frac {M_{2}^2}{M_{1}^2} \frac {T_{2}}{T_{1}}, \end{equation}

where subscript $1$ indicates pre-shock conditions and subscript $2$ indicates post-shock conditions. Approximating this change using the Rankine–Hugoniot normal shock relations for a perfect gas results in the dynamic pressure jump across the shock:

(5.7) \begin{equation} \frac {Q_{2}}{Q_{1}} = \frac {\left ( \gamma - 1 \right ) M_{1}^2 + 2}{\left ( \gamma + 1 \right ) M_{1}^2}. \end{equation}

Conveniently, this is the same relation as the inverse density ratio $\epsilon ={\rho _{s1}}/{\rho _{s2}}$ , which must be the case because momentum is conserved across a shock. Recall that we consider only dynamic pressure gradients normal to the free stream flow. Since a normal shock is also normal to the free stream flow, the gradient will change only in magnitude (not in direction) across the shock. So, we can take the free stream-normal gradient of $Q_1$ and $Q_2$ . Recall also that we aligned the $y$ -axis of our particle to the free stream-normal dynamic pressure gradient, so $\boldsymbol{\nabla \!}Q_{\perp }$ only has one dimension and can be represented by the scalar $||\boldsymbol{\nabla \!}Q_{\perp }||$ . Therefore,

(5.8) \begin{equation} ||\boldsymbol{\nabla \!}Q_{\perp 2}|| = \frac {\left ( \gamma - 1 \right ) M_{1}^2 + 2}{\left ( \gamma + 1 \right ) M_{1}^2} ||\boldsymbol{\nabla \!}Q_{\perp 1}|| .\end{equation}

Taking $M_1$ to be the free stream Mach number on the stagnation line $M_s$ , we can define a parameter $\sigma$ which characterises the change in dynamic pressure gradient across the shock:

(5.9) \begin{equation} \sigma = \left \{ \begin{array}{ll} \dfrac {\left ( \gamma - 1 \right ) M_{s}^2 + 2}{\left ( \gamma + 1 \right ) M_{s}^2}, & M \geqslant 1.0, \\ 1.0, & \mathrm{otherwise}. \end{array} \right . \end{equation}

A high-Mach limit for pressure lift coefficient is then

(5.10) \begin{equation} \lim _{M_s\to \infty } C\!_{L,p}=\sigma \widetilde {\boldsymbol{\nabla \!}Q_{\perp }}. \end{equation}

Next, a high-Mach correction parameter $C_{\textit{MH}}$ is fit to all steady axisymmetric cases where $M_s=6$ , so that the lift coefficient due to pressure forces is

(5.11) \begin{equation} \lim _{M_s\to \infty } C\!_{L,p}=C_{\textit{MH}} \sigma \widetilde {\boldsymbol{\nabla \!}Q_{\perp }}. \end{equation}

The high-Mach correction parameter is found to be

(5.12) \begin{equation} C_{\textit{MH}} = 0.8539, \end{equation}

which is consistent with the high-Mach drag coefficient for a sphere $C_D \sim 0.9$ (Hornung, Schramm & Hannemann Reference Hornung, Schramm and Hannemann2019), despite the theoretical Newtonian limit being $C_D=1.0$ .

Next, we find well-behaved functions that represent subsonic, low-supersonic and planar-symmetric wake lifting effects relative to the high-Mach axisymmetric limit. A subsonic representation should asymptote to a fixed limit as $M_s \rightarrow 0$ and approach 1 as $M_s \rightarrow \infty$ . A function that fits this form and the DNS is

(5.13) \begin{equation} \alpha _{\textit{SP}} = 1 - \frac {1}{A_1 M_s^2 + A_2}. \end{equation}

This leaves us with a need for a low supersonic representation that limits to 1 at both high and low Mach number, but provides a suitable adjustment for lift coefficient near $1 \lt M_s \lt 1.75$ . A function that fits this form and the DNS is

(5.14) \begin{equation} \alpha _{\textit{TP}} = 1 - A_3 \: e^{-(M_{s}-1.25)^2}, \end{equation}

which is effectively an inverted Gaussian centred at $M_s=1.25$ .

Finally, we need a representation for planar-symmetric wakes. We observe in the DNS data that this correction must be stronger at lower $\widetilde {\boldsymbol{\nabla \!}Q_\perp }$ and weaker at higher $M_s$ . Furthermore, as $\widetilde {\boldsymbol{\nabla \!}Q_\perp }$ decreases, this function must grow more slowly than $\widetilde {\boldsymbol{\nabla \!}Q_\perp }$ decreases, so lift coefficient limits to a fixed value when combined with the high-Mach scaling term. This function must also limit to fixed values as $M_s \rightarrow 0$ and $M_s \rightarrow \infty$ . A function that fits this form and the DNS is

(5.15) \begin{equation} C_{\textit{WP}}=\left ( 1 - \frac {1}{B_1M_s^4+B_1} \right ) \: \widetilde {\boldsymbol{\nabla \!}Q_\perp }^{\frac {-1}{B_2M_s^2+1}}. \end{equation}

We must also account for the step change in lift coefficient when the wake behaviour switches from steady axisymmetric to steady planar-symmetric. We can use a combination of hyperbolic tangent functions and the delineation between wake types $M_{\textit{tp}}$ from (4.4) to define

(5.16) \begin{equation} f^{\pm }= 1 \pm \tanh \left [ S \log \left ( \frac {M_s}{M_{\textit{tp}}}\right ) \right ]\!, \end{equation}

which can be used to smoothly switch between $C_{\textit{WP}}$ when $M_s \lt M_{\textit{tp}}$ and 1 when $M_s \gt M_{\textit{tp}}$ . This leaves us with the steady planar-symmetric wake term

(5.17) \begin{equation} C_{\textit{WPf}} = \frac {1}{2}\big ( C_{\textit{WP}}f^- + f^+ \big ). \end{equation}

The final lift coefficient due to pressure is then

(5.18) \begin{equation} C_{\textit{L,P}} = C_{\textit{MH}} \sigma \widetilde {\boldsymbol{\nabla \!}Q_{\perp }} \: \alpha _{\textit{SP}} \: \alpha _{\textit{TP}} \: C_{\textit{WPf}}. \end{equation}

The model is then fitted to the DNS data. Fit parameters are listed in table 5. The pressure lift coefficient model is shown in figure 17 and generally matches the DNS data. Small discrepancies exist for certain cases because lifting effects near changing wake behaviour are often erratic and therefore difficult to quantify with simple well-behaved analytic expressions.

Table 5. Coefficients for pressure lift model.

Figure 17. Fitted pressure lift coefficient ( $\widetilde {\textit{Re}_{\textit{p}}}$ defined in Appendix C).

5.2. Viscous lift

Similar to pressure lift, we first develop a high-Mach limit for viscous lift coefficient. Since viscous lift and pressure lift are both functions of non-dimensional dynamic pressure gradient as shown in § 3.2, we start with the high-Mach term from (5.11). However, we need a correction term to relate the high-Mach pressure lift limit to the high-Mach viscous lift limit. We observe in the DNS data that viscous lift coefficient scales inversely with $\widetilde {\textit{Re}_{\textit{p}}}$ and directly with $M_s$ at high Mach number. A function that meets these requirements and represents the DNS data is

(5.19) \begin{equation} C_{\textit{VH}} = C_1 M_s^{C_2log(M_s+0.1)+C_3} \: \widetilde {\textit{Re}_{\textit{p}}}^{-C_4}, \end{equation}

where $\widetilde {\textit{Re}_{\textit{p}}}$ is a transformed compressible Reynolds number derived by Singh et al. (Reference Singh, Kroells, Li, Ching, Ihme, Hogan and Schwartzentruber2020). The transformed Reynolds number is based on the Howarth (Reference Howarth1948) and Stewartson (Reference Stewartson1949) compressible boundary layer transformation and accounts for the changing relationship between compressibility and viscous effects as the Mach number increases. A summary of the transformed Reynolds number can be found in Appendix C. The addition of $0.1$ inside the $\log (M_s)$ term limits the function as $M_s \rightarrow 0$ . The high-Mach limit for viscous lift coefficient is then

(5.20) \begin{equation} \lim _{M_s\to \infty } C_{L,\tau } = C_{\textit{MH}} \sigma \widetilde {\boldsymbol{\nabla \!}Q_{\perp }} \: C_{\textit{VH}}. \end{equation}

Next, we incorporate a subsonic representation. The term should asymptote to a fixed value as $M_s \rightarrow 0$ and approach 1 as $M_s \rightarrow \infty$ . A function that fits this form and the DNS is

(5.21) \begin{equation} \alpha _{S\tau } = D_1e^{-D_2(M_s-0.5)^2}+1. \end{equation}

We observe in the DNS data that no additional low supersonic term is needed for viscous lift. Finally, we incorporate the planar-symmetric wake representation, which takes the same form as for pressure lift,

(5.22) \begin{equation} C_{W\tau }=\left ( 1 - \frac {1}{E_1M_s^4+E_1} \right ) \: \widetilde {\boldsymbol{\nabla \!}Q_\perp }^{\frac {-1}{E_2M_s^2+1}}, \end{equation}

and is then switched on or off using (5.16) to give

(5.23) \begin{equation} C_{W \tau f} = \frac {1}{2}\big ( C_{W\tau }f^- + f^+ \big ). \end{equation}

Combining these terms leaves us with the viscous lift coefficient

(5.24) \begin{equation} C_{L,\tau } = C_{\textit{MH}} \sigma \widetilde {\boldsymbol{\nabla \!}Q_{\perp }} \: C_{\textit{VH}} \: \alpha _{S\tau } \: C_{W \tau f}. \end{equation}

The model is then fitted to the DNS data. Fit parameters are listed in table 6. The viscous lift coefficient model is shown in figure 18 and generally matches the DNS data.

5.3. Total lift

Finally, we compute the total lift coefficient as

(5.25) \begin{equation} C\!_L = C_\textit{L,P} + C_{L,\tau }. \end{equation}

Table 6. Coefficients for viscous lift model.

Figure 18. Fitted viscous lift coefficient ( $\widetilde {\textit{Re}_{\textit{p}}}$ defined in Appendix C).

The fitted model was tested against five blind validation cases excluded from the fitting process. Case details are listed in Appendix A.5. As shown in figure 19, the fitted model matches the DNS within 11 % for all blind validation cases except at $M_s=0.6$ , where the fitted model still matches within 20 %. Note that the case at $M_s=0.8$ has a steady planar-symmetric wake, while other cases have a steady axisymmetric wake, indicating that the model does predict effects of changing wake behaviour.

Figure 19. Blind validation of total lift coefficient, $\widetilde {\boldsymbol{\nabla\! }Q}=15\times 10^{-3}$ .

5.4. Model limitations

The present work uses the Navier–Stokes equations to model fluid dynamics, which means that this study is limited to continuum flow. There exists little to no data in the published literature regarding particles in compressible shear flows, either in continuum or rarefied flows. Recall also that the present work only studies linear shear flow. Therefore, the proposed model comes with several limitations.

1. (i) Lifting effects at very low Reynolds or Mach number are not studied. As a result, there are no model terms that collapse to incompressible and/or Stokes flow limits.

2. (ii) The proposed model includes no term(s) for rarefaction effects.

3. (iii) The proposed model does not address lifting effects induced by unsteady wakes.

4. (iv) The proposed model applies only to linear or near-linear shear flows, i.e. flows where $\parallel\! ( {D^2}/{Q_{s}}){\nabla} ^2 Q_{\infty} \! \parallel \ll 1$ .

It is expected that as more data become available to address these limitations, the proposed lift model will and should evolve, just as spherical drag models have evolved with data availability.