3.1. Effect of density ratio

In this section, we investigate the dynamics of freely ascending or descending particles with thermal buoyancy $({\textit{Ga}}=60, {\textit{Ri}}=0.5)$ across a wide range of density ratios $(0.0 \leqslant \rho _r \leqslant 4)$ focusing on the effect of particle inertia. Our analysis begins from examining the trajectory patterns of particles under same buoyant conditions, as shown in figure 2. A significant dependence of path stability on particle inertia is observed. For inertialess particles $(\rho _r=0.0)$ , stable zigzag oscillations dominate the ascent dynamics. Moderately inertial particles $(\rho _r = 0.5, 0.9)$ transition to 3-D spiral rising patterns. Conversely, particles slightly heavier than the fluid $(\rho _r= 1.1,2)$ exhibit spiral settling trajectories with progressively diminishing oscillation amplitudes. Notably, oscillatory motions become nearly imperceptible at $\rho _r=4$ , indicating inertial stabilization of the settling path. In conclusion, as the density ratio increases, path instabilities are suppressed, resulting in reduced oscillation or spiralling amplitudes and delayed the onset of instability.

To quantify the deviation of particle motion from a vertical path, we perform an elliptical fitting (Fitzgibbon, Pilu & Fisher Reference Fitzgibbon, Pilu and Fisher1999) of the non-vertical portions of the particle trajectories projected onto the horizontal $XZ$ plane. From the fitted ellipses, we extract the semimajor axis length $a$ , the aspect ratio (AR), defined as the ratio of the major to minor axis, and the root-mean-square error (RMSE) as a measure of fit quality, defined as

(3.1) \begin{equation} {\textit{RMSE}}=\sqrt {\frac {1}{N} \sum _{i=1}^N d_i^2}, \end{equation}

where $d_i$ is the minimum Euclidean distance between the $i$ th data point and the fitted ellipse, and $N$ is the total number of points. As shown in figure 3, increasing the density ratio leads to reductions in both the lateral oscillation amplitude (quantified by $a$ ) and the aspect ratios, consistent with the qualitative observations discussed earlier.

Figure 3.

Particle trajectories projected onto the horizontal $XZ$ plane (orange lines) and corresponding ellipse fits (red lines) for various density ratios: (a) $\rho _r = 0.0$ , (b) $\rho _r = 0.5$ , (c) $\rho _r = 1.1$ , (d) $\rho _r = 4.0$ .

Table 2.

Ellipse fitting results for particle trajectories in the $XZ$ plane. The semimajor axis length $a$ , $\textit{AR}$ , $\textit{RMSE}$ and the normalized error $\textit{RMSE}$ / $a$ are reported for particles with different density ratios $\rho _r$ . Note that $\textit{RMSE}$ and $\textit{RMSE}$ / $a$ are presented in units of $10^{-3}$ .

The fitted ellipses effectively capture the dominant geometric features of the motion. For the $\rho _r=0.0$ particle, the projected trajectory is nearly linear, exhibiting an extremely large aspect ratio (AR=214.70), characteristic of zigzagging motion. For particles with higher density ratios $(\rho _r=0.5, 1.1)$ , the projected trajectories form ellipses with moderate AR and low fitting residuals (RMSE), indicative of spiralling motion. In contrast, the $\rho _r=4.0$ particle displays only very weak lateral instability. The oscillation amplitude is minimal ( $a = 0.02$ ), and the motion remains close to vertical. Since the RMSE represents the average deviation of the particle positions from the fitted ellipse, we also computed the dimensionless quantity RMSE/ $a$ to facilitate comparison across particles with different oscillation amplitudes. A lower value of RMSE/ $a$ indicates better agreement between the particle trajectory and the fitted ellipse. Overall, the fitting results are satisfactory, and the relevant parameters for various density ratios are summarized in the table 2. For clarity, the values of RMSE and RMSE/ $a$ reported in the table are scaled by $10^3$ .

All particles attain identical dimensionless terminal velocities prior to the onset of instability as shown in figure 4, indicating that steady-state motion under fixed ${\textit{Ga}}$ and ${\textit{Ri}}$ conditions is independent of particle inertia. However, the effects of particle inertia manifest in the duration required to reach terminal velocity and the time it begins to oscillate. Figure 4(b) provides a zoomed-in view to reveal subtle effects of particle inertia. Particles with low inertia exhibit amplified velocity oscillations and slightly reduced vertical velocities. In contrast, high-inertia particles demonstrate delayed vortex-induced oscillations with attenuated amplitudes, reflecting the enhanced stabilization through the dominance of particle inertia.

Figure 4.

Time evolution of particle vertical velocity with varying density ratio. Panel (b) provides a zoomed-in view of panel (a).

Figure 5.

Spectrum of the fluctuation of horizontal velocity of particles with different $\rho _r$ . Panel (b) provides a zoomed-in view of panel (a) near the frequency peak. At $\rho _r=0.0$ , the peak frequency is $\bar {f}=0.12$ , with the dominant Strouhal number ${\textit{St}}\approx 0.1207$ . The Fourier transform (FFT) was evaluated as $X[k] = \sum_{n=0}^{N-1} x[n],w[n],e^{-j2\pi kn/N}$ ($j^2=-1$) , with a Hanning window $w[n]=\tfrac{1}{2}[1-\cos(2\pi n/(N-1))]$ applied.

Complementary to vertical velocity analysis, we systematically investigated the horizontal velocity components in the XZ-plane during particle oscillations. Spectral analysis in figure 5, employing detrending and a Hanning window to avoid spectral leakage, reveals a fundamental frequency peak at $\bar {f}=0.12$ across cases except $\rho _r=4$ , particularly pronounced for light particles. The dimensionless frequency is defined as $\bar {f}=\textit{fL}_C/|U_C|$ , where $f$ is the oscillation frequency. Owing to the clearly periodic nature of the motion observed in this set of cases, the Strouhal number is employed to characterize the dynamics of particles. It is defined as ${\textit{St}}=\textit{fL}_C/|\widetilde {U_v}|=\bar {f} / (|\widetilde {U_v}|/U_C)$ . Under the present conditions, the Strouhal number exhibits a slight decreasing trend with increasing density ratio, ranging from 0.1207 at $\rho _r=0\,{\rm to}\,0.1142$ at $\rho _r=2$ . These values are of the same order of magnitude as those reported by Kotouč et al. (Reference Kotouč, Bouchet and Dušek2009), though slightly higher, likely due to differences in Prandtl number and the fact that the particle is freely moving in the present study. Furthermore, as illustrated in figure 5(b), the amplitude of the primary frequency progressively attenuates with increasing density ratios. Notably, a secondary low-frequency component ( $\bar {f}\lt 0.05$ ) develops at mediate density ratios, as indicated by the inset in figure 5(a). This phenomenon is reminiscent of the tertiary instability observed in the wake instability of particles. This subharmonic signature is associated with spiralling motion superposed on the primary oscillatory trajectory, suggesting inertial modification of 3-D instability modes.

Figure 6.

(a) Time evolution of the horizontal component of the hydrodynamic force $\boldsymbol{F}^H$ acting on particles with different density ratios. (b) Spectrum of the horizontal hydrodynamic force and (c) its zoomed-in view.

We also examine the hydrodynamic force $\boldsymbol{F}^H$ , defined as

(3.2) \begin{equation} \boldsymbol{F}^H = \oint _{\partial{\kern-1pt}{V_p}} {{\boldsymbol{\tau }} \boldsymbol{\cdot }{\boldsymbol{n}}{\textrm{d}} A}. \end{equation}

Although the vertical component of the hydrodynamic force is approximately 20 times greater than the horizontal component (not shown), it primarily acts as resistance to settling or rising. The horizontal component, in contrast, plays a direct role in influencing the lateral motion of the particle. Therefore, here we focus on the horizontal hydrodynamic force $\boldsymbol{F}^H_h$ , as illustrated in figure 6(a). Once the particle becomes unstable, this force initially increases and then exhibits periodic oscillations, closely corresponding to the fluctuations in the horizontal particle velocity. A spectral analysis of the horizontal hydrodynamic force is presented in figures 6(b) and 6(c), where figure 6(c) provides an enlarge view of the spectrum in the range $0\leqslant \bar {f} \leqslant 0.15$ . A dominant spectral peak at $\bar {f} \approx 0.12$ is observed, consistent with the primary oscillation frequency of the velocity signal. Particles with density ratios closer to unity experience stronger horizontal hydrodynamic forces, while those with density ratios further from one exhibit weaker horizontal forces. For the lightest case ( $\rho _r = 0$ ), the hydrodynamic force is minimal, yet the horizontal oscillation is most pronounced (see figure 5 c) due to the very low particle inertia. Conversely, for the heaviest case ( $\rho _r = 4$ ), its settling motion is dominated by particle inertia in the vertical direction, resulting in both weak horizontal forces and minimal lateral oscillations. In addition, a secondary low-frequency peak appears within $0\leqslant \bar {f} \leqslant 0.05$ , particularly for heavy particles $(\rho _r \gt 1)$ . This component is associated with the spiralling motion observed during particle settling.

Figure 7.

Wake visualized using the Q criterion (Iso-Q surface, $Q=0.12$ , coloured by the local value of vertical velocity): (a,b) $\rho _r = 0.5$ ; (c,d) $\rho _r = 2.0$ . The black spheres denote particle positions, while (a) and (c) present cross-sectional views in the Z–Y plane, and (b) and (d) present cross-sectional views in the X–Y plane.

To further investigate particle–fluid interactions, we visualize the wake dynamics of a light particle $(\rho _r=0.5)$ and a heavy particle $(\rho _r=2.0)$ using a $Q=0.12$ isosurface (vortex identification criterion), coloured by the vertical velocity component, as shown in figure 7. Note that these are given the same buoyancy conditions, i.e. ${\textit{Ga}}=60$ , and ${\textit{Ri}}=0.5$ . Under thermal buoyancy $({\textit{Ri}}=0.5)$ , both particles exhibit periodic vortex shedding at ${\textit{Ga}}=60$ , directly driving oscillatory motion. Spectral analysis of fixed-probe velocity signals within the particle wakes reveals a fundamental frequency $\bar {f}=0.06$ , precisely half the horizontal oscillation frequency of particles. This oscillation behaviour is similar to the horizontal oscillation during the settlement of a single cold particle in a two-dimensional channel, as described by Gan et al. (Reference Gan, Chang, Feng and Hu2003). The natural convection stemming from thermal buoyancy induces high pressure at the back stagnation point of the particles. The high-pressure counterbalance between the front and back stagnation points causes a ‘buckling instability’, which in turn prompts vortex shedding and subsequent wake oscillations (Gan et al. Reference Gan, Chang, Feng and Hu2003). Initially, the particles commence horizontal rotation, as illustrated in figure 8(a). When the rotational degrees of freedom of the particles are restricted, the onset of oscillation is postponed and the oscillation amplitude is amplified, as depicted in figure 8(b).

Figure 8.

(a) Particle rotation velocity $(\rho _r=0.5,{\textit{Ga}}=60,{\textit{Ri}}=0.5)$ . (b) Comparison of the time evolution of particles vertical velocity in the presence or absence of particle rotation.

Under thermal buoyancy with positive Richardson number, particle inertia significantly influences path stability, oscillation characteristics and wake dynamics even for particles with ${\textit{Ga}}=60$ . This may be attributed to the fact that thermal buoyancy enhances the particle velocity and promotes the onset of wake instability. As discussed in Dan & Wachs (Reference Dan and Wachs2010), the acceleration or deceleration of particle velocities is fundamentally influenced by the alignment of natural convection flows with the gravity-dominated motion. In our simulations with ${\textit{Ri}}=0.5$ , hot particles generate ascending natural convection plumes when rising, whereas cold particles exhibit descending convective flows when falling. Velocity amplification occurs when the thermal-buoyancy-driven flow direction aligns synergistically with the rising or settling velocity, leading to wake instability and complex kinematic patterns.

3.2. Effect of thermal buoyancy and gravity

Given the sensitivity of light particles to thermal-buoyancy effects, we examine a series of cases of light particles $(\rho _r =0.5)$ under various ${\textit{Ga}}$ and ${\textit{Ri}}$ to elucidate the effects of thermal buoyancy and gravity. Figure 9 presents the particle trajectories for the cases with ${\textit{Ga}} = 40$ and ${\textit{Ga}} = 100$ as examples. In general, increasing the Richardson number $({\textit{Ri}})$ induces progressive transitions of the ascent patterns of particles: from stable vertical rising to 3-D spiralling or planar zigzagging, and ultimately chaotic rising. For particles with ${\textit{Ga}}=40$ , stronger buoyancy forcing $({\textit{Ri}} \gtrsim 0.5)$ is required to initiate path instabilities, marked by the emergence of horizontal velocity fluctuations. Whereas ${\textit{Ga}}=100$ particles exhibit lower destabilization thresholds, entering fully developed chaotic ascent regimes at ${\textit{Ri}}=0.5$ . This contrast highlights the dual control of thermal buoyancy and gravity: higher ${\textit{Ga}}$ amplifies instability susceptibility by enhancing momentum advection, while increased ${\textit{Ri}}$ accelerates shear layer destabilization through the generation of thermal convection.

Figure 9.

The trajectories of light particles $(\rho _r=0.5)$ with varying buoyant forces: (a–c) ${\textit{Ga}}=40$ and ${\textit{Ri}}=0.5/ 0.75/ 1.0$ ; (d–f) ${\textit{Ga}}=100$ and ${\textit{Ri}}=0.1/ 0.25/ 0.5$ .

Similarly, we fit ellipses to the horizontal projection of trajectories ( $XZ$ plane) of light particles $(\rho _r =0.5)$ under varying gravitational and thermal-buoyancy effects. The corresponding parameters, including the semimajor axis length $a$ , AR, RMSE and the normalized error RMSE/ $a$ are summarized in table 3. For the classification purposes, the stable vertical regime is defined as the cases with $a\lt 0.02$ , the zigzagging regime is characterized by AR ${\gt }\, 30$ and RMSE $\times\, 10^3$ / $a\lt 0.1$ , the spiralling regime corresponds to AR between 1 and 30 with RMSE $\times\, 10^3$ / $a\lt 100$ , and the chaotic regime is defined by AR ${\lt }\, 5$ and RMSE $\times\, 10^3$ / $a\gt 100$ . These four regimes are mapped in figure 10, where the black regions denote stable vertical rising and green regions correspond to 3-D chaotic rising. It is worth noting that the transition between the zigzagging and spiralling regimes is not sharply resolved within the current parameter space. Furthermore, flow regimes such as the stable oblique regime, commonly observed under isothermal conditions, are not observed in the present study, in agreement with the numerical results reported by Majlesara et al. (Reference Majlesara, Abouali, Kamali, Ardekani and Brandt2020).

Table 3.

Ellipse fitting results for light particles ( $\rho _r=0.5$ ) rising under different effects of gravity ( ${\textit{Ga}}$ ) and thermal buoyancy ( ${\textit{Ri}}$ ). The semimajor axis length $a$ , $\textit{AR}$ , $\textit{RMSE}$ and normalized error $\textit{RMSE}$ / $a$ are reported. Note that $\textit{RMSE}$ and $\textit{RMSE}$ / $a$ are presented in units of $10^{-3}$ .

To further investigate the relationship between particle motion and horizontal hydrodynamic forces, we examine the temporal evolution of the angle between the horizontal hydrodynamic force vector and the horizontal particle velocity vector during the rising motion of light particles ( $\rho _r = 0.5$ ). This angle, denoted as $\theta _{\textit{FV}}$ , is defined as

(3.3) \begin{equation} \theta _{\textit{FV}}(t)=\cos ^{-1}\left (\frac {\boldsymbol{F}^H_h \boldsymbol{\cdot }\boldsymbol{U}_h}{\left |\boldsymbol{F}^H_h\right |\left |\boldsymbol{U}_h\right |}\right )\!, \end{equation}

Figure 10.

Regime map of the rising behaviour of a light sphere $(\rho _r=0.5)$ in an infinite domain with thermal buoyancy in the ( ${\textit{Ga}},{\textit{Ri}}$ ) parameter plane. The annotations are as follows: SV, stable vertical; ZZ, zigzagging; SP, spiralling; 3DC, 3-D chaotic.

where $\boldsymbol{F}^H_h$ and $\boldsymbol{U}_h$ represent the horizontal components of the hydrodynamic force and particle velocity, respectively. In the stable vertical regime, where the horizontal hydrodynamic force is negligible, the angle $\theta _{\textit{FV}}$ appears randomly distributed between 0 $^\circ$ to 180 $^\circ$ , reflecting the absence of lateral forcing. In contrast, in the zigzagging and spiralling regimes, $\theta _{\textit{FV}}$ exhibits periodic or quasiperiodic variations, as illustrated in figure 11. For the zigzagging regime, the angle oscillates periodically between 0 $^\circ$ and 180 $^\circ$ , indicating that the horizontal force primarily acts as drag – resisting motion without inducing deviation from the main oscillation direction. Interestingly, even within this regime, subtle differences are evident. For instance, the case $({\textit{Ga}}, {\textit{Ri}}) = (60, 0.4)$ shows a slight angular offset compared with $({\textit{Ga}}, {\textit{Ri}}) = (40, 0.75)$ , indicating a tendency towards the spiralling behaviour. This transition is further evidenced by the corresponding horizontal trajectory aspect ratios: AR=32.30 for $(60, 0.4)$ is nearly half that of $(40, 0.75)$ , for which AR=54.47. In the spiralling regime, the angle $\theta _{\textit{FV}}$ oscillates between approximately 20 $^\circ$ and 160 $^\circ$ , indicating that the hydrodynamic force includes both resistive and deflective components in the horizontal plane. This periodic deflection may be driven by the combined effects of vortex shedding, path instability and asymmetries in the incoming flow. Notably, although the case $({\textit{Ga}}, {\textit{Ri}}) = (60, 0.5)$ initially exhibits spiralling motion, it gradually transitions to a zigzag-like rising pattern at later times. We interpret this evolution as the result of a competition between thermal-buoyancy-induced asymmetric instabilities and the planar symmetry imposed by vortex shedding. At early stages, thermal buoyancy dominates, promoting vortex shedding and lateral deviations. However, as the vortex field intensifies, it tends to restore planar symmetry, thereby altering the long-term particle trajectory. Furthermore, as ${\textit{Ga}}$ increases, the periodicity of $\theta _{\textit{FV}}$ becomes progressively disrupted, suggesting a transition towards the chaotic regime. For the chaotic regime (not shown), the angle $\theta _{\textit{FV}}$ no longer displays clear periodicity but instead shows a broad, near-random distribution, revealing a complex and irregular interaction between horizontal force and motion.

Figure 11.

Temporal evolution of the angle $\theta _{\textit{FV}}$ for light particles $(\rho _r=0.5)$ under different combinations of $({\textit{Ga}}, {\textit{Ri}})$ , categorized as (a) zigzagging and (b) spiralling regimes.

Figure 12.

Instantaneous contours of vertical fluid velocity in the X–Y plane of ${\textit{Ga}}=40$ light particles $(\rho _r=0.5)$ with different ${\textit{Ri}}$ : (a) ${\textit{Ri}}=0.5$ , (b) ${\textit{Ri}}=0.75$ , (c) ${\textit{Ri}}=1.0$ . White solid lines indicate zero relative vertical velocity ( $u_{v,r}=|U_v-v|=0$ ), while black solid lines represent temperature contours ranging from 0.1 to 0.9.

Table 4.

The time-averaged recirculation length $\widetilde {L_v}$ , the normalized length $\widetilde {L_v}/\widetilde {L_{v0}}$ , and its standard deviation Std( $\widetilde {L_v}$ ) for light particles $(\rho _r=0.5)$ at various ${\textit{Ga}}$ and ${\textit{Ri}}$ . Here $\widetilde {L_{v0}}$ denotes the value at ${\textit{Ri}} = 0$ .

Figure 13.

Instantaneous iso-Q surface ( $Q=0.12$ ) of ${\textit{Ga}}=40$ light particles $(\rho _r=0.5)$ with different ${\textit{Ri}}$ : (a) ${\textit{Ri}}=0.5$ , (b) ${\textit{Ri}}=0.75$ , (c) ${\textit{Ri}}=1.0$ . Panel (d) provides a vertical side view of case (c). The black spheres denote particle positions, and the surfaces coloured by the local value of vertical velocity.

We analyse instantaneous flow fields during the ascent of ${\textit{Ga}}=40$ particles to elucidate thermal-buoyancy effects. Figure 12 presents XY-plane slices at increasing three buoyancy intensities $({\textit{Ri}}=0.5/0.75/1.0)$ with vertical velocity colormaps. The corresponding 3-D vortex structures in figure 13 are visualized using $Q = 0.12$ isosurfaces coloured by vertical velocity. Thermal buoyancy drives convection around the hot particle: heated fluid is advected from the rear of the sphere to the shear region, while colder ambient fluid is entrained into the wake. This mechanism intensifies the recirculation with increasing ${\textit{Ri}}$ , as evidenced by the stronger vertical velocities within the high-velocity regions (outlined by white contours) with higher ${\textit{Ri}}$ in figure 12. To further quantify this effect, we computed the time-averaged recirculation length $\widetilde {L_v}$ , defined as the maximum distance from the sphere surface to the point where the relative vertical velocity ( $u_{v,r}=|U_v-v|=0$ ) is zero. The results are summarized in table 4. For the isothermal cases ( ${\textit{Ri}} = 0$ ), the length $\widetilde {L_v}$ increases monotonically with ${\textit{Ga}}$ without vortex shedding. For a given ${\textit{Ga}}$ , $\widetilde {L_v}$ also increases with ${\textit{Ri}}$ , indicating that thermal buoyancy strengthens and extends the recirculation zone. An exception is observed at $({\textit{Ga}}, {\textit{Ri}}) = (40, 0.5)$ , where the length $\widetilde {L_v}$ is longer than those for $(40, 0.75)$ and even $(40, 1.0)$ . This irregularity results from the onset of unsteady vortex shedding, which introduces fluctuations in the wake and affects the recirculation length under the present definition. A similar trend was reported by Chouippe et al. (Reference Chouippe, Krayer, Uhlmann, Dušek, Kiselev and Leisner2019). The normalized length $\widetilde {L_v}/\widetilde {L_{v0}}$ , where $\widetilde {L_{v0}}$ corresponds to the isothermal case $({\textit{Ri}}=0)$ , also increases significantly with increasing ${\textit{Ri}}$ . However, the relative enhancement becomes less pronounced at higher ${\textit{Ga}}$ . Specifically, for ${\textit{Ri}} = 0.25$ and $0.5$ , the normalized length $\widetilde {L_v}/\widetilde {L_{v0}}$ decreases as ${\textit{Ga}}$ increases. Interestingly, the length $\widetilde {L_v}$ shows ${\textit{Ri}}$ -dependent trends: it increases with ${\textit{Ga}}$ at ${\textit{Ri}} = 0.25$ but decreases at ${\textit{Ri}} = 0.5$ , again reflecting the influence of wake instability. Moreover, the standard deviation of $\widetilde {L_v}$ also increases with both ${\textit{Ga}}$ and ${\textit{Ri}}$ , reflecting the development of stronger fluctuations and enhanced wake unsteadiness. These trends emphasize the destabilizing role of thermal buoyancy. Vorticity fields in figure 13 reveal that vortex shedding is modulated with ${\textit{Ri}}$ , deviating from the classical dynamics observed in isothermal cases at ${\textit{Ga}}=40$ . As thermal buoyancy increases $({\textit{Ri}}\gt 0)$ , vortex shedding becomes more intense and increasingly non-axisymmetric. The cases with ${\textit{Ri}} = 0.5$ and ${\textit{Ri}} = 1.0$ exhibit more pronounced asymmetries, corresponding to spiralling and chaotic regimes, respectively. In particular, at ${\textit{Ri}} = 0.5$ , the symmetry of the vortex ring is disrupted, although the shedding remains relatively organized. In contrast, at ${\textit{Ri}} = 1.0$ , even the planar symmetry of the shedding is disrupted, indicating a fully chaotic regime. By comparison, the case with ${\textit{Ri}} = 0.75$ maintains an approximately planar-symmetric structure, indicating that thermal buoyancy has a weaker influence on the vortex configuration in this intermediate regime. Enhanced vertical flow also leads to elevated temperatures within vortex cores, as indicated by black isotherms in figure 12. Collectively, these observations highlight the critical role of thermal buoyancy-induced convection in modifying wake dynamics and triggering path instabilities.

In addition, we present instantaneous iso-Q surface for the ${\textit{Ga}}=60$ cases to further investigate the transitions between different regimes. As illustrated in the regime map (figure 10), an inversion in the sequence of regimes is observed as ${\textit{Ga}}$ increases. At ${\textit{Ga}}=40$ , the particle transitions from stable vertical motion to spiralling, then to zigzagging and finally to chaotic motion (3-D chaotic) with increasing ${\textit{Ri}}$ . In contrast, for ${\textit{Ga}}=60$ , spiralling appears before zigzagging motion. This inversion can be explained by the interplay between thermal buoyancy-induced asymmetric instabilities and the planar symmetries enforced by vortex shedding. At low ${\textit{Ga}}$ , thermal buoyancy triggers a buckling instability that breaks the axial symmetry of the vortex ring, as illustrated in figure 13(a). Therefore, the particle exhibits spiralling motion which becomes more pronounced with increasing ${\textit{Ri}}$ , as evidenced by the increasing semimajor axis length $a$ and AR in the ${\textit{Ga}}=40$ cases (table 3). For the case with ${\textit{Ri}}=0.75$ , vortex shedding becomes more pronounced, but thermal buoyancy has weak impact on the planar symmetry of the shedding vortices. As a result, the motion is predominantly zigzagging, as shown in figure 13(b). At higher ${\textit{Ga}}$ , vortex shedding develops earlier and enforcing planar symmetry at moderate ${\textit{Ri}}$ , as illustrated in figures 14(a) and 14(b). These two cases thus exhibit transitional behaviour between spiralling and zigzagging, as discussed above. As ${\textit{Ri}}$ further increases, thermal buoyancy eventually breaks the symmetry, inducing more chaotic motion with increasing $a$ but decreasing AR. This is further supported by instantaneous vorticity structure in figure 14(c).

Figure 14.

Instantaneous iso-Q surface ( $Q=0.12$ ) of ${\textit{Ga}}=60$ light particles $(\rho _r=0.5)$ with different ${\textit{Ri}}$ : (a) ${\textit{Ri}}=0.4$ , (b) ${\textit{Ri}}=0.5$ , (c) ${\textit{Ri}}=0.75$ .

The vertical velocities of particles at ${\textit{Ga}}=40$ and $100$ modulated by thermal buoyancy with varying ${\textit{Ri}}$ are compared in figure 15. Increasing ${\textit{Ri}}$ , which means stronger positive thermal buoyancy, enhances ascent velocities while inducing progressive transitions from stable vertical to zigzagging or spiralling, and ultimately chaotic rising. For the ${\textit{Ga}}=100$ cases, we also present the results for negative thermal-buoyancy effect $({\textit{Ri}}=-0.5)$ , corresponding to a light particle that is colder than the surrounding fluid. In this configuration, thermal buoyant suppression manifests through inhibited cold fluid entrainment into shear layers, attenuating recirculation intensity and reducing the terminal velocity to 0.72 $U_{g,0}$ (versus $U_{g,0}$ at ${\textit{Ri}}=0$ ). This counterflow mechanism contrasts with positive ${\textit{Ri}}$ behaviour, demonstrating bidirectional control of dynamics of particles via the orientation of thermal buoyancy-induced convection.

Figure 15.

Vertical velocity of particles for different ${\textit{Ri}}$ : (a) ${\textit{Ga}}=40$ ; (b) ${\textit{Ga}}=100$ .

Figure 16.

Time-averaged vertical velocity $|\widetilde {U_v}|$ of light particles ( $\rho _r=0.5$ ) under negative thermal buoyancy ( ${\textit{Ri}}\lt 0$ ), plotted against ${\textit{Ga}}$ for various ${\textit{Ri}}$ . Green diamonds represent the cases without thermal buoyancy ( ${\textit{Ri}} = 0$ ), blue circles correspond to ${\textit{Ri}} = -0.5$ , black markers show corresponding results for heavy particles ( $\rho _r = 1.5$ ) at ${\textit{Ga}} = 60$ and the purple star represents ${\textit{Ri}} = -0.75$ at ${\textit{Ga}} = 100$ .

Under negative ${\textit{Ri}}$ , particles consistently reach a steady regime within the explored parameter space. Figure 16 summarizes the variation of the time-averaged vertical velocity $|\widetilde {U_v}|$ of light particles ( $\rho _r=0.5$ ) with respect to both ${\textit{Ga}}$ and ${\textit{Ri}}$ . It can be seen that as ${\textit{Ga}}$ decreases, the suppressive effect of negative thermal buoyancy on particle vertical velocity becomes more pronounced. We also compute the settling velocity of heavy particles ( $\rho _r=1.5$ ) at ${\textit{Ga}} = 60$ . The results closely mirror those observed for light particles ( $\rho _r=0.5$ ), as indicated by the black markers in figure 16. Additionally, simulations with stronger negative thermal buoyancy ( ${\textit{Ri}} = -0.75$ ) for the light particle ( $\rho _r=0.5$ ) at ${\textit{Ga}}=100$ confirm further suppression of particle velocity, though the motion remains steady (see the purple star in figure 16). To further illustrate the flow behaviour, figure 17 presents instantaneous contours of vertical fluid velocity during the rising process. With increasing ${\textit{Ga}}$ , convective heat transfer intensifies, resulting in a thinner and more elongated thermal wake. The zero-relative vertical velocity isosurface tightly concentrated around the particle surface, indicating that the downward motion of cold fluid nearly eliminates the recirculation zone. Furthermore, a larger $|{\textit{Ri}}|$ leads to faster descent of the cold fluid. Particles attain a steady and axisymmetric regime in all the cases studied. We did not observe the 3-D regimes ( $m \geqslant 1$ modes) reported by Kotouč et al. (Reference Kotouč, Bouchet and Dušek2009) in their §7, likely due to the relatively low particle Reynolds number in the present simulations.

Figure 17.

Instantaneous contours of vertical fluid velocity in the X–Y plane of light particles $(\rho _r=0.5)$ under negative thermal buoyancy at various ${\textit{Ga}}$ and ${\textit{Ri}}$ : (a) ${\textit{Ga}}=40$ and ${\textit{Ri}}=-0.5$ ; (b) ${\textit{Ga}}=60$ and ${\textit{Ri}}=-0.5$ ; (c) ${\textit{Ga}}=100$ and ${\textit{Ri}}=-0.5$ ; (d) ${\textit{Ga}}=100$ and ${\textit{Ri}}=-0.75$ . White solid lines indicate zero relative vertical velocity ( $u_{v,r}=|U_v-v|=0$ ), which are closely aligned with the particle surface. Black solid lines represent temperature contours ranging from –0.1 to –0.9.

Table 5.

Time-averaged particle Reynolds number ( $\widetilde {{\textit{Re}}_p}$ ) and Strouhal numbers ( ${\textit{St}}$ ) for particles at ${\textit{Ga}} = 40$ and $60$ . Here NA indicates that no dominant frequency could be identified.

Figure 18.

Spectrum of the fluctuation of horizontal velocity of particles for different ${\textit{Ri}}$ : (a) ${\textit{Ga}}=40$ ; (b) ${\textit{Ga}}=100$ .

Horizontal velocity oscillations exhibit monotonic amplification with increasing ${\textit{Ri}}$ as shown in figure 18, demonstrating direct correlation between the intensity of thermal buoyancy and path instability. Notably, ${\textit{Ga}}=100$ particles transition to chaotic motion occurring more immediately than ${\textit{Ga}}=40$ particles. The transitional stages (e.g. planar oscillations, spiralling trajectories) are negligible for particles with higher ${\textit{Ga}}$ . Furthermore, the spectra indicate that particles at lower ${\textit{Ga}}$ exhibit more distinct periodic or quasiperiodic behaviour. Accordingly, we also evaluated the Strouhal numbers for these cases. The time-averaged particle Reynolds number and Strouhal number for particles at ${\textit{Ga}} = 40$ and $60$ under various Richardson numbers are summarized in the table 5. As ${\textit{Ri}}$ increases, the Reynolds number increases while the Strouhal number exhibits a slight decrease. This trend is consistent with the observations reported by Kotouč et al. (Reference Kotouč, Bouchet and Dušek2009) in their §8.

As discussed in the above, particles with higher ${\textit{Ga}}$ exhibit reduced destabilization thresholds under identical ${\textit{Ri}}$ condition compared with the particles with lower ${\textit{Ga}}$ . To isolate gravitational effects, we perform comparative analysis of vertical ascent velocities and horizontal velocity spectrum under fixed ${\textit{Ri}}=\{0.25,0.5\}$ across varying Galileo numbers $(40 \leqslant {\textit{Ga}} \leqslant 120)$ , as shown in figure 19. We also examine several cases with the same classification but distinct dynamics, i.e. $({\textit{Ga}},{\textit{Ri}})=(40,0.5)$ , $(60,0.5)$ and $(100,0.25)$ . Although all fall within the spiralling regime in the regime map (figure 10), their vertical velocity (figure 19 a,b) reveal notable differences. The case of $({\textit{Ga}}, {\textit{Ri}}) = (40, 0.5)$ develops more slowly, with weak periodic oscillations in vertical velocity appearing only after $tU_C/L_C \approx 250$ , and with limited amplitude. By contrast, the $({\textit{Ga}}, {\textit{Ri}}) = (60, 0.5)$ case shows clearer periodicity and stronger deviations, while the $({\textit{Ga}}, {\textit{Ri}}) = (100, 0.25)$ case exhibits irregular, non-periodic oscillations, consistent with stronger vortex shedding and enhanced flow unsteadiness. The horizontal frequency spectra (figure 19 c,d) further support these observations, highlighting differences in dominant spectral peaks. With constant buoyancy forcing $({\textit{Ri}})$ , increasing ${\textit{Ga}}$ amplifies ascent velocities and modifying path instability progressions. Lower ${\textit{Ga}}$ particles exhibit periodic zigzagging or spiralling motions characteristic of buoyancy-modulated vortex shedding. While particles with higher ${\textit{Ga}}$ transition to chaotic ascent at reduced thermal buoyancy intensities. The stability of particle trajectories is governed by the synergistic interplay between gravitational forcing and thermal buoyancy, and the accelerated chaotic onset at elevated ${\textit{Ga}}$ suggests gravitational amplification of thermal convection effects through momentum advection.

Figure 19.

Vertical velocities (a,b) and horizontal velocity spectrum (c,d) under different Richardson number: (a,c) ${\textit{Ri}}=0.25$ ; (b,d) ${\textit{Ri}}=0.5$ .

We also examine the magnitude of the horizontal hydrodynamic force (not shown), which exhibits similar trends to that of the horizontal particle velocity. For a given Galileo number, the horizontal force increases with stronger thermal buoyancy (larger ${\textit{Ri}}$ ), and the force signal becomes more periodic at lower ${\textit{Ga}}$ . In contrast, for larger ${\textit{Ga}}$ , the particle experiences enhanced horizontal forcing accompanied by augmented chaotic motions.

3.3. Heat transfer coefficient

Complementary heat transfer analysis examines the time-averaged particle Reynolds numbers $(\widetilde {{\textit{Re}}_p})$ and corresponding Nusselt numbers $(\widetilde {\textit{Nu}})$ during the ascent of light particles $(\rho _r = 0.5)$ across varying ${\textit{Ri}}$ and ${\textit{Ga}}$ , as shown in figure 20. Nusselt number $\textit{Nu}$ , measuring the efficiency of heat transfer between fluid and particle, is defined as

(3.4) \begin{equation} \textit{Nu}(t) = \frac {\int _S \textit{Nu}_{loc} {\textrm d}S}{\pi D_p^2 } = \int _S(- \bar {\boldsymbol{\nabla }} \theta \boldsymbol{\cdot }\boldsymbol{n}) {\textrm d}\bar {S}. \end{equation}

Figure 20.

Time-averaged (a) particle Reynolds number $(\widetilde {{\textit{Re}}_p})$ and (b) corresponding Nusselt number $(\widetilde {\textit{Nu}})$ during the rising of light particles $(\rho _r=0.5)$ under varying thermal buoyancy $({\textit{Ri}})$ and Galileo number $({\textit{Ga}})$ . (c) Projection of the data onto the $(\widetilde {{\textit{Re}}_p}, \widetilde {\textit{Nu}})$ plane. (d) Three-dimensional representation of the coupled dependence of $\widetilde {\textit{Nu}}$ on $\widetilde {{\textit{Re}}_p}$ , ${\textit{Ga}}$ and ${\textit{Ri}}$ .

The dimensionless temperature gradient distribution $\bar {\boldsymbol{\nabla }} \theta$ on the particle surface is mainly influenced by convective transport and thermal conduction. For the isothermal case $({\textit{Ri}}=0)$ , the Nusselt numbers $(\textit{Nu})$ align closely with classical empirical correlations, such as the Ranz–Marshall equation (Rans & Marshall Reference Rans and Marshall1952), $\textit{Nu} = 2.0 + 0.6\textit{Re}^{0.5}{\textit{Pr}}^{0.33}(10 \lt \textit{Re} \lt 10^4,{\textit{Pr}} \gt 0.7)$ , validating the baseline of heat transfer. Figure 20(d) presents a 3-D diagram illustrating the relationship among ${{\textit{Ga}}}$ , $\widetilde {{\textit{Re}}_p}$ and $\widetilde {\textit{Nu}}$ , with figure 20(a–c) corresponding to its projections onto the respective parameter planes. As shown in figure 20(a), positive thermal buoyancy ( ${\textit{Ri}} \gt 0$ ) increases $\widetilde {{\textit{Re}}_p}$ relative to the isothermal case ( ${\textit{Ri}} = 0$ ), consistent with the recirculation-induced particle velocity enhancement mechanism discussed previously. In contrast, for negative thermal buoyancy ( ${\textit{Ri}} = -0.5$ ), $\widetilde {{\textit{Re}}_p}$ is significantly reduced, aligning with previous observations that stabilizing thermal buoyancy suppresses the average particle vertical velocity. At a fixed ${\textit{Ga}}$ , the Nusselt number $\widetilde {\textit{Nu}}$ generally increases with ${\textit{Ri}}$ (figure 20 b), reflecting the strengthened convective transport induced by thermal-buoyancy-enhanced recirculation. However, a different trend emerges when viewed at equivalent Reynolds numbers $\widetilde {{\textit{Re}}_p}$ , as in figure 20(c). For ${\textit{Ri}} \gt 0$ , particles with low ${\textit{Ga}}$ exhibit slightly lower $\widetilde {\textit{Nu}}$ than the empirical predictions at the same $\widetilde {{\textit{Re}}_p}$ , highlighting the complex interplay between thermal buoyancy and heat transfer. Conversely, for ${\textit{Ri}}\lt 0$ , the $\widetilde {\textit{Nu}}$ consistently exceeds empirical values at equivalent $\widetilde {{\textit{Re}}_p}$ .

When considering thermal-buoyancy effects with ${\textit{Ri}} \gt 0$ on heat transfer, there are two competing mechanisms: thermal buoyancy enhanced wake recirculation amplifies near surface flow velocities, compressing thermal boundary layers and augmenting forced convection, while heat accumulation in vortex cores may establish reverse thermal gradients that suppress conductive heat transfer. At low ${\textit{Ga}}$ , the introduction of positive thermal buoyancy $({\textit{Ri}} \gt 0)$ increases $\widetilde {{\textit{Re}}_p}$ while paradoxically reducing $\widetilde {\textit{Nu}}$ below isothermal values at equivalent $\widetilde {{\textit{Re}}_p}$ , as recirculation redistributes temperature through vortex-mediated flows. In contrast, higher ${\textit{Ga}}$ regimes exhibit chaotic rising that enhances heat transfer efficiency through turbulent mixing. For negative thermal buoyancy $({\textit{Ri}} \lt 0)$ , convective heat transfer is substantially suppressed due to reduced particle velocity. However, interestingly, these cases exhibit higher $\widetilde {\textit{Nu}}$ at equivalent ${\textit{Re}}_p$ than predicted by empirical models. This enhancement is attributed to the formation of a downward cold wake, as opposed to upward heated plumes for light particles with ${\textit{Ri}} \gt 0$ , which maintains stronger local temperature gradients near the particle surface. As a result, thermal conduction near the particle surface is more pronounced, becoming the dominant heat transfer mechanism under stabilizing thermal buoyancy conditions.

These observations indicate that $\widetilde {\textit{Nu}}$ is no longer a sole function of $\widetilde {{\textit{Re}}_p}$ under thermal buoyancy. To further explore this point, we compare two cases with similar particle Reynolds numbers $\widetilde {{\textit{Re}}_p}$ but differing ${\textit{Ga}}$ and ${\textit{Ri}}$ : Case A with $({\textit{Ga}},\ {\textit{Ri}}) = (120, 0)$ , and Case B with $({\textit{Ga}},\ {\textit{Ri}}) = (100,\ 0.5)$ , both yielding $\widetilde {{\textit{Re}}_p} \approx 130$ . Despite the matched Reynolds number, the wake dynamics vary significantly. Case A exhibits a relatively stable wake characterized by an axisymmetric vortex ring and a recirculation length of approximately $1.62$ , with $\widetilde {\textit{Nu}} \approx 8.54$ . In contrast, Case B displays a more extended and chaotic recirculation region ( $\widetilde {L_v} \approx 2.87$ ) due to enhanced vortex stretching and unsteadiness, yet its Nusselt number is slightly lower ( $\widetilde {\textit{Nu}} \approx 8.48$ ). This comparison suggests that enhanced recirculation may lead to heat accumulation in the wake, which suppresses conductive transport and reduces heat transfer. Therefore, $\widetilde {{\textit{Re}}_p}$ alone cannot uniquely determine the Nusselt numbers for particles with thermal buoyancy. Accurate characterization of heat transfer instead requires accounting for the coupled effects of ${\textit{Ga}}$ and ${\textit{Ri}}$ , as they collectively modulate wake structure and near-surface temperature gradients.