4.1. Steady vertical motion at Ga = 150 (Re_T ≈ 186)

Figure 5 shows time histories of the velocities for spheres with different densities at Ga = 150, indicating that all cases eventually attained a steady state. The vertical velocities of the spheres (figure 5a) reached a similar value ($|U_z^\ast |\, \approx \,1.24$) for all the sphere-to-fluid density ratios considered, resulting in the same terminal Reynolds number of Re_T ≈ 186. This is reasonable because our system completely maps its reverse case (i.e. uniform flow past a fixed sphere) at steady state irrespective of the density of the spheres.

Figure 5. Time history of the velocities for spheres with different densities at Ga = 150: (a) vertical component ($|{U_z^\ast } |$), (b) horizontal component ($U_x^\ast $). The terminal Reynolds number is Re_T ≈ 186 for all particle densities considered.

We also include the result of Deloze et al. (Reference Deloze, Hoarau and Dušek2012) in figure 5(a), that is, the terminal velocity of a heavy sphere ($|U_z^\ast |\, \approx \,1.26$) in a 5d circular tube. These two values were very close to each other. In addition, the terminal velocity of an unconfined sphere at Ga = 150, obtained using (1.4), is approximately 1.3, slightly larger than that of a confined sphere. Thus, it can be inferred that at low Ga the influence of the tube walls is limited when the tube sizes are larger than 5d. However, the particle density affects the motion of the spheres in terms of the transient nature of both $U_z^\ast $ (figure 5a) and $U_x^\ast $ (figure 5b).

In figure 6, we show iso-surfaces of the vorticity (figure 6a,b), where λ ₂ = −0.01, and vertical fluid velocity (figure 6c,d), where $U_z^\ast=- 0.1$ for ρ_s/ρ = 2 and $U_z^\ast= 0.1$ for ρ_s/ρ = 0.1. Both results, which were obtained when the motion of the spheres reached a steady state, showed axial symmetrical wakes that were downwards for the heavy sphere (figure 6a,c) and upwards for the light sphere (figure 6b,d). Our study indicates that the spheres in the square tube will lose their steadiness beyond a critical Galileo number (Ga_crit ≈ 157), along with loss of the axial symmetry of the wake seen at Ga < Ga_crit.

Figure 6. Simulation results for a heavy sphere (a,c) with ρ_s/ρ = 2 and a light one (b,d) with ρ_s/ρ = 0.1 at steady state. Note that (a,b) show iso-surfaces of the vorticity (λ ₂ = −0.01), and (c,d) show iso-surfaces of the vertical fluid velocity where $u_z^\ast = {-}0.1$ for the heavy sphere and $u_z^\ast= 0.1$ for the light sphere, respectively. The light sphere (d) seems to be orange in colour because of the transparent iso-surface of the fluid velocity (the same below).

4.2. Onset of unsteadiness at Ga_crit ≈ 157 (Re_T ≈ 196)

Our simulations indicate that the steadiness of spheres confined to the square tube is lost at Ga ≈ 157 for all sphere-to-fluid density ratios considered (0.1 ≤ ρ_s/ρ ≤ 2.3). To demonstrate this, we provide trajectories of both a heavy sphere with ρ_s/ρ = 2.3 (figure 7a,b) and a light sphere with ρ_s/ρ = 0.1 (figure 7c,d) at Ga = 156 (figure 7a,c) and Ga = 157 (figure 7b,d). The blue circles in each panel represent the initial positions of the spheres (the same below). Small but finite oscillations are seen for both spheres at in figure 7(b,d), indicating a transition to unsteadiness occurs at Ga_crit ≈ 157, corresponding to a terminal Reynolds number of Re_T ≈ 196. This critical value is consistent with those of both an unconfined sphere (Ga_crit = 155 for ρ_s/ρ = 0 and Ga_crit = 160 for ρ_s/ρ → ∞, Jenny et al. Reference Jenny, Dušek and Bouchet2004) and a confined one in a 5d circular tube (155 < Ga_crit < 160, Deloze et al. Reference Deloze, Hoarau and Dušek2012). Later studies (Fabre, Tchoufag & Magnaudet Reference Fabre, Tchoufag and Magnaudet2012; Zhou & Dušek Reference Zhou and Dušek2015) reported a critical Ga falling between 155 and 156. This further confirmed our results when we considered the effects of the tube confinement. More importantly, the present bifurcation is Hopf type owing to wall repulsion (Deloze et al. Reference Deloze, Hoarau and Dušek2012), which is different from the case of an unconfined sphere. To further illustrate the role of the wall repulsion in the transition, we show contours of the 2-D vorticity in an (x, y) plane for the light sphere (ρ_s/ρ = 0.1) at Ga = 150–160 in figure 8. All the results were obtained when the sphere reached a statistically steady state. Note that the 2-D vorticity is given by $\omega _z^\ast= {(\boldsymbol{\nabla } \times \boldsymbol{u})_z}\ast d/{U_g}$ .

Figure 7. Trajectories of a heavy sphere with ρ_s/ρ = 2.3 (a,b) and a light one ρ_s/ρ = 0.1 (c, d) at Ga = 156 (a,c) and Ga = 157 (b,d), respectively. Note that both spheres are seen to lose steadiness at Ga = 157. The corresponding Reynolds number is Re_T ≈ 196. The blue circle in each panel denotes the initial position of the spheres (the same below). In this study we normalized the particle trajectory through x* = x/d, y* = y/d, and z* = (z − z ₀)/d. Note that $(x_0^\ast ,y_0^\ast ,z_0^\ast ) = ( - 1,0,0)$ for all spheres.

Figure 8. Instantaneous contours of the 2-D vorticity ($\omega _z^\ast $) in the (x, y) plane for the light sphere (ρ_s/ρ = 0.1) at different Galileo numbers: (a) Ga = 150, (b) Ga = 156, (c) Ga = 157, (d) Ga = 158 and (e) Ga = 160. Note that the 2-D vorticity is given by $\omega _z^\ast= {(\boldsymbol{\nabla } \times \boldsymbol{u})_z}\ast d/{U_g}$.

Similar vorticity distribution is seen near each tube wall at Ga < Ga_crit (figure 8a,b). Therefore, the sphere which is falling along the tube axis does not experience any lateral force from the walls. As Ga increases (Ga ≥ Ga_crit) the magnitude of the vorticity grows and the vortex becomes unstable (see the upper and bottom walls in figure 8c−e), reflecting an unsteady motion of the light sphere in the (x, y) plane.

Figure 7 also shows that the sphere-to-fluid density ratio has a significant influence on the oscillatory pattern of the spheres at near-critical Galileo numbers. It appears that the heavy sphere oscillated in the central plane (figure 7b). However, the light sphere eventually rose along a helical path (figure 7d). To address the effects of ρ_s/ρ we projected the particle trajectories onto the (x, y) plane for different spheres at Ga = 158 (figure 9). We chose Ga = 158, a Galileo number slightly higher than Ga_crit, to provide a better view of the particle motion. Interestingly, light spheres tend to undergo helical motion, with their path projections being circular at small ρ_s/ρ (figure 9a), which transform into an ellipse at large ρ_s/ρ (figure 9b). In contrast, the heavy spheres tended to oscillate near the central (y* = 0) plane, where they were released (figure 9c,d). This tendency was stronger for heavier spheres. The behaviours shown in figure 9 are reasonable because lighter spheres with smaller inertia are more sensitive to solid walls, and accordingly, are more likely to alter their directions of motion. As Ga increases, heavy spheres may exhibit helical trajectories, as discussed in the following section (§ 4.3).

Figure 9. Top views (x–y projections) of particle trajectories for spheres with different densities at Ga = 158: (a) ρ_s/ρ = 0.1, (b) ρ_s/ρ = 0.5, (c) ρ_s/ρ = 1.5 and (d) ρ_s/ρ = 2.3.

4.3. Helical motion of heavy spheres (ρ_s/ρ > 1)

Similar to a circular tube, heavy spheres are widely observed undergoing helical motion in a square tube at Ga > Ga_crit. Figure 10 shows falling trajectories of a sphere (ρ_s/ρ = 2) at different Galileo numbers (Ga = 165, 170, 185 and 195), providing a clear view of the helical motion.

Figure 10. Falling trajectories of a heavy sphere (ρ_s/ρ = 2) showing RHM at (a) Ga = 165, (b) Ga = 170 and (c) Ga = 185, and AHM at (d) Ga = 195.

In all cases, the spheres exhibited a zigzagging-like movement once they were released from rest before giving way to the helical motion. However, a significant difference is clearly observed in the result of Ga = 195 (figure 10d) as compared with the others (figure 10a−c): the sphere changes its helical direction in an alternate manner while falling in the tube. This is unique because heavy spheres generally circle in the same direction when confined to a circular tube (Deloze et al. Reference Deloze, Hoarau and Dušek2012). Heavy spheres confined to a square tube also underwent RHM at low Ga. The RHM is here referred to as a type of helical motion with unchanged helical direction and a constant helical size (figure 10a−c). Figure 11 provides a clear view of this, which shows the respective top views of figure 10(a−d). Accordingly, we use the term ‘alternate helical motion’ to represent another type of helical motion (figure 10d) which is extensively seen for heavy spheres in the square tube as well (figure 4). For instance, our study indicates that the sphere (ρ_s/ρ = 2) undergoes AHM for the Galileo number ranging between 190 and 205, with its helical direction changing more frequently at a higher Ga (figure 4). Our regime map indicates the heavy spheres exhibit RHM for Ga approximately ranging from 160 to190. This range is surprisingly consistent with that of the steady oblique regime for an unconfined heavy sphere (Jenny et al. Reference Jenny, Dušek and Bouchet2004; Zhou & Dušek Reference Zhou and Dušek2015; Raaghav et al. Reference Raaghav, Poelma and Breugem2022).

Figure 11. Top views of figure 10(a−d). The arrows denote the circling directions of the sphere.

Figure 11 shows that the 2-D projections of the particle trajectories are nearly circular at low Ga (figure 11a,b). The limit cycle expands outwards as Ga increases, suggesting a stronger helical motion of the sphere at a higher Ga. Further increasing Ga, however, results in the tendency of the sphere to move towards the tube corners because of the confinement. This resulted in an elliptical limit cycle (figure 11c). The situation is quite different for the AHM (figure 11d): no limit cycle exists. Owing to the change in the helical direction, the sphere alternately moves towards the corners located in either of the diagonal planes, leading to a random and homogeneous distribution of the particle trajectory in its 2-D projection (figure 11d).

On the other hand, the rotational angles of the sphere (figure 12) provide a better view on the difference between RHM and AHM. We here define θ_x, θ_y and θ_z as the angles of the sphere rotating around the x-, y- and z-axes, as shown in figure 12(a). It is seen that the rotational angle θ_z increases monotonically after the initial transients die out for Ga = 185 (RHM), representing a positive rotation of the sphere around the z-axis (figure 12a). For Ga = 195 (AHM), however, the change in θ_z indicates the sphere was rotating around the z-axis in a positive or negative direction alternately (figure 12a). Figure 12(b) shows other rotational angles (θ_x and θ_y) of the sphere at Ga = 185. Note that t ₁ ~ t ₃ correspond to those in figure 11(c). Both components suggest that the sphere was ‘rolling’ forward as it circled in the counterclockwise direction.

Figure 12. Time history of the rotational angles of the same sphere as in figure 10: (a) rotational angles around the z-axis ($\theta _z^\ast $) at Ga = 185 and 195, (b) rotational angles around the x-($\theta _x^\ast $) and y-axes ($\theta _y^\ast $) at Ga = 185. The angles were normalized through $\theta _x^\ast= {\theta _x}/{\rm \pi}$, $\theta _y^\ast= {\theta _y}/{\rm \pi}$ and $\theta _z^\ast= {\theta _z}/{\rm \pi}$. The initial value of each of θ_x, θ_y and θ_z was zero, i.e. $\theta _x^\ast (0) = 0$, $\theta _y^\ast (0) = 0$ and $\theta _z^\ast (0) = 0$. Note that t ₁ ~ t ₃ in (b) correspond to those in figure 11(c).

Figure 13 shows time history of the velocities ($U_z^\ast $ and $U_h^\ast $) of the same sphere as in figure 10 and instantaneous vortical structures of the wake for three Galileo numbers, that is, Ga = 165, 185 and 195. It is clear that the oscillations in both components of the particle velocity are significantly suppressed after the transition from ZM to RHM (figure 13a,b). In contrast, the AHM of the sphere exhibited the opposite trend (figure 13c). The double-threaded vortical structures in the wake, which are obtained when the sphere is closest to the right wall, are clearly seen for all three Ga values in figure 13(d). The sphere helical motion is a direct result of the interaction between the double-threaded vortex and the walls, as Deloze et al. (Reference Deloze, Hoarau and Dušek2012) stated. In addition, the principal threads became longer and more bent as Ga increased, suggesting stronger oscillatory motion of the sphere at higher Ga (figure 13d). However, no visible vortex shedding was observed even at Ga = 195.

Figure 13. Time history of the velocities (top: vertical, bottom: horizontal) for the same sphere as in figure 10 at (a) Ga = 165 (RHM), (b) Ga = 185 (RHM) and (c) Ga = 195 (AHM), and (d) instantaneous vortical structures (λ ₂ = −0.005) in the wake of the sphere at the three Galileo numbers.

To further elucidate the behaviour of AHM, we provide a close-up view of the change in the circling (helical) direction of a heavy sphere (ρ_s/ρ = 2.3) falling at Ga = 205 (figure 14). The solid and dashed lines represent a specific time period during which the sphere circles in the counterclockwise and clockwise directions (figure 14a), respectively. The opposite was true for the sphere during another time period (figure 14b). When this sphere undergoes AHM, it comes closer to the tube walls every time it circles in the tube (note the solid lines in figure 14). As a result, the wall repulsion becomes stronger at smaller sphere–wall separations resulting from the increased pressure and flow shear rate between the sphere and tube walls. At a time of, say t* = t ₂ (as in figure 14a), the sphere could be ‘bounced’ back by a tube wall once the repulsion is strong enough. Then, the sphere runs towards the opposite wall in a nearly straight line following the ‘bounce back’ (say t* = t ₂ ~ t ₃ as in figure 14a). A second ‘bounce back’ occurs when the sphere approaches the opposite wall where a mirror-reflection-like path is formed, leading to a change in the direction of circling (say t* = t ₃ ~ t ₅ as in figure 14a). A similar pattern is shown in figure 14(b). It should be remarked that here we use ‘bounce back’ or ‘bounced’ just to emphasize a very strong sphere–wall interaction. The sphere did not collide with the tube walls, as indicated by the simulations.

Figure 14. Top views of the trajectories for a heavy sphere (ρ_s/ρ = 2.3) falling at Ga = 205 during two time periods: (a) t* ≈ 1220–1435, and (b) t* ≈ 1614–1714, respectively. The arrows denote the circling direction of the sphere, indicating the AHM.

The strong interactions between a heavy sphere and the tube walls may also be reflected by the wake features generated by the falling sphere. Figure 15 shows instantaneous vortical structures (principally displaying the double-threaded character) at t ₁ ~ t ₅ (see figure 14a), during which the sphere changes its helical direction. There is no vortex shedding at t* = t ₁ before the first ‘bounce back’ takes place (figure 15a). By contrast, an intense vortex shedding is seen at both t* = t ₂ and t* = t ₃ when the sphere is successively ‘bounced’ back by the left and right walls (figure 15b,c). The shape of the wake is reminiscent of the oscillating oblique motion of an unconfined sphere at similar Galileo numbers (Jenny et al. Reference Jenny, Dušek and Bouchet2004; Uhlmann & Dušek Reference Uhlmann and Dušek2014; Zhou & Dušek Reference Zhou and Dušek2015; Raaghav et al. Reference Raaghav, Poelma and Breugem2022). The major difference is that the double-threaded vortex detaches from the unconfined sphere regularly. After the helical direction of the sphere changed, the double-threaded vortex was quickly restored until the next cycle began (figure 15d,e).

Figure 15. Instantaneous iso-surfaces of the vorticity (λ ₂ = −0.005) for the same sphere (AHM) as in figure 14 at different times: (a) t* = t ₁, (b) t* = t ₂, (c) t* = t ₃, (d) t* = t ₄ and (e) t* = t ₅. Note that t ₁ ~ t ₅ are indicated in figure 14(a).

Further increasing Galileo number (Ga > 205) may result in another transition (from AHM) to the ZM of heavy spheres (ρ_s/ρ ≥ 1.8). This is discussed in § 4.5. After the zigzagging regime, it is interesting to reveal that these confined spheres will fall in a helical path (RHM) again (figure 4), due to a significant change in the wake of the spheres. To illustrate this we show instantaneous iso-surfaces of the vorticity in the wake of a heavy sphere (ρ_s/ρ = 2.2) at Ga = 230 (the highest Ga considered) in figure 16. The vortical structures exhibiting a chain of vortex loops (so-called hairpin vortices) in the wake of the sphere were significantly different from those at low Ga (figure 13d or figure 15) in terms of wake curvature and distortion. We believe these structures are relevant to the helical/rotating regime reported for heavy spheres (ρ_s/ρ = 3.19 and 3.9) at Re_T > 270 (Raaghav et al. Reference Raaghav, Poelma and Breugem2022). Figure 16 also illustrates how the vortex loops rolling down from the sphere connect and evolve as the vorticity convects and develops downstream. The general feature is very similar to those revealed for flow over a fixed sphere (Re = 300, Johnson & Patel Reference Johnson and Patel1999; Re = 350, Mittal Reference Mittal1999). This is reasonable because the terminal Reynolds number of our case (figure 16) is Re_T ≈ 300. In addition, the sphere oscillated much faster at Ga = 230 compared with the RHM of the same sphere at low Ga. Furthermore, the normalized period (T*) was 80 at Ga = 190, which declined to 69 at Ga = 230.

Figure 16. Instantaneous iso-surfaces of the vorticity (λ ₂ = −0.01) for a heavy sphere (ρ_s/ρ = 2.2) at Ga = 230 during one period (T* ≈ 69): (a) t* = 0, (b) t* ≈ 0.25T*, (c) t* ≈ 0.45T* and (d) t* ≈ 0.85T*. Note that the sphere undergoes the RHM.

4.4. Helical motion of light spheres (ρ_s/ρ < 1)

Light spheres exhibited helical motion at a lower Galileo number (near Ga_crit) than their counterparts (figure 4). Figure 17 shows trajectories of a light sphere (ρ_s/ρ = 0.2) for the Galileo number ranging from 160 to 195. Similar to the case of a heavy sphere (figure 10), a zigzagging-like ascending path quickly gives way to a helical path at each Ga for a light sphere (figure 17). Top views provide a better understanding of this phenomenon (figure 18). A circular-shaped limit cycle is eventually formed for each case except Ga = 190 and 195 (figure 18d,e), suggesting the sphere undergoes the RHM at Ga ≤ 180. The limit cycle expands rapidly outwards as Ga increases from 160 to 180. However, it appears that the helical path of the sphere becomes less smooth at Ga = 190 (figure 17d) compared with those at Ga ≤ 180. This tendency is more significant at Ga = 195 (figure 17e). Here, we refer to this pattern as JHM because its 2-D projections show clear jaggedness (figure 18). Our study indicates that the JHM occurs only for very light spheres (ρ_s/ρ ≤ 0.3), as shown in figure 4.

Figure 17. Rising trajectories of a light sphere (ρ_s/ρ = 0.2) showing the RHM at (a) Ga = 160, (b) Ga = 170, (c) Ga = 180, and showing the JHM at (d) Ga = 190 and (e) Ga = 195.

Figure 18. Top views of figure 17(a−e).

As stated above, the AHM and JHM exist for heavy and light spheres, respectively, at high Ga (Ga > 185). Figures 10(d) and 17(d,e) show that they differ from each other in an essential manner because of the significant difference in inertia between the heavy and light spheres. For a heavy sphere exhibiting the AHM, wall repulsion has a slow effect on the change in the direction of its translational movement until the sphere approaches the wall (figure 14). By contrast, light spheres are more sensitive to wall repulsion. In the case of intermediate Galileo numbers, the wall repulsive effect is more significant for the motion of a lighter sphere. Therefore, the limit cycles of a light sphere were generally more circular than those of a heavy sphere (figures 11 and 18). In particular, for a light sphere exhibiting JHM at a high Ga, the 2-D projection of its trajectory has a circular-like shape, suggesting that it is rising helically with a similar helical size and shape (figure 18d,e). However, the heavy spheres exhibiting AHM showed significant variations in both helical sizes and shapes (figure 11d). More importantly, we reveal that vortex shedding occurs when a light sphere undergoes JHM, which is not a necessary phenomenon for AHM, as seen in figure 19, which shows instantaneous vortical structures corresponding to figure 17. A double-threaded vortex was clearly observed at Ga = 170 and 180 (figure 19b,c), whereas it was invisible at Ga = 160 (figure 19a). At a higher Ga, when the light sphere rises in a helical path with jaggedness, shedding of the double-threaded vortex is demonstrated (figure 19d,e).

Figure 19. Instantaneous iso-surfaces of the vorticity (λ ₂ = −0.01) for the same sphere (ρ_s/ρ = 0.2) as in figure 17. Each result is chosen at the time when the sphere moves to the rightmost position.

To cast more light on this issue, we show iso-surfaces of the vertical vorticity ($\omega _y^\ast= \pm 0.5$) at different Ga in figure 20(a−e), which corresponds to figure 19(a−e). Note that $\omega _y^\ast= {(\boldsymbol{\nabla } \times \boldsymbol{u})_y}\ast d/{U_g}$. From $\omega _y^\ast $ one could observe the vortex features of the wake from a different aspect. For instance, the length of the sphere wake is clearly increasing with Ga, and the wake becomes wavy for JHM, suggesting stronger unsteadiness in the motion of the sphere (figure 20d,e). More importantly, it seems that the tube walls have a negligible effect on the sphere at Ga = 160 (figure 20a). At Ga ≥ 170, the wall effect becomes increasingly important as Ga increases, as indicated by the vorticity on the right wall shown in figure 20(b−e). This could not be revealed from λ ₂ (figure 19).

Figure 20. Iso-surfaces of the vertical vorticity ($\omega _y^\ast= \pm 0.5$) corresponding to figure 19(a−e). Note that $\omega _y^\ast= {(\boldsymbol{\nabla } \times \boldsymbol{u})_y}\ast d/{U_g}$.

The jaggedness of the trajectories, which we believe is attributed to the added mass effect, became more noticeable for the spheres with lower densities. To illustrate this, we provide simulation results for light spheres with different densities (ρ_s/ρ = 0.1, 0.3 and 0.5) at Ga = 195 in figure 21. We obtain a clear view of JHM in figure 21(a,b) while considering RHM in figure 21(c). For spheres with very small ρ_s/ρ (say ρ_s/ρ ≤ 0.3), the added mass forces become dominant over the wall repulsion as the spheres are significantly accelerating or decelerating, which is the case of figure 21(a,b). Under the action of the added mass forces, the spheres may not run along a smooth path until the wall repulsion forces dominate and cause an abrupt change in the direction of movement. To clarify this, we present top views of the trajectory of the lightest sphere studied (ρ_s/ρ = 0.1) at two specific time periods in figure 21(d). The 2-D projections are shaped like a regular pentagon (figure 21d). Between the two neighbouring corners of the pentagon, the sphere runs along a nearly straight path, resulting from a strong added mass force. At each corner where the sphere approaches a tube wall, it responds promptly to wall repulsion by abruptly changing its direction of movement. In addition, the vortical structures (figure 22) provide evidence of an intense vortex shedding as this sphere rises helically. Compared with figure 19(e), the vortices in the wake of the lighter sphere are more twisted and bent.

Figure 21. Rising trajectories of a light sphere with different densities at Ga = 195: (a) ρ_s/ρ = 0.1 (JHM), (b) ρ_s/ρ = 0.3 (JHM), (c) ρ_s/ρ = 0.5 (RHM). (d) Top views of (a) at two specific time periods.

Figure 22. Instantaneous iso-surfaces of the vorticity (λ ₂ = −0.01) for a light sphere (ρ_s/ρ = 0.1) at Ga = 195 (JHM): (a) t* = t ₁, (b) t* = t ₂, (c) t* = t ₃ and (d) t* = t ₄. Note that t ₁ ~ t ₄ are indicated in figure 21(d).

Note that the present simulations were performed under conditions that the initial particle positions were fixed at $(x_0^\ast ,y_0^\ast ,z_0^\ast ) = ( - 1,0,12)$ and $(x_0^\ast ,y_0^\ast ,z_0^\ast ) = ( - 1,0,25)$, for heavy and light spheres, respectively (figure 2). To check the effects of the initial position we varied both $x_0^\ast $ and $y_0^\ast $ while keeping z ₀* unchanged and show the results in figure 23 for ρ_s/ρ = 0.2 and Ga = 180 (RHM). The trajectory in figure 23(a) is equivalent to that in figure 18(c). It seems that the calculations based on different initial positions will converge to a similar solution (i.e. a helical path with the same size) after the initial transients die out (figure 23).

Figure 23. Top views of the trajectories for a rising sphere (ρ_s/ρ = 0.2) at Ga = 180 (RHM): effects of the initial position, i.e. ($x_0^\ast $, $y_0^\ast $). Note that the trajectory resulting from (−1, 0) is equivalent to that in figure 18(c).

4.5. Zigzagging motion of spheres

Here, we discuss the ZM of heavy and light spheres confined to a square tube. Previous studies (Jenny et al. Reference Jenny, Dušek and Bouchet2004; Horowitz & Williamson Reference Horowitz and Williamson2010; Zhou & Dušek Reference Zhou and Dušek2015; Auguste & Magnaudet Reference Auguste and Magnaudet2018) have reported that only light spheres may rise in a zigzagging path in the absence of confinement, although a significant disagreement on the critical value of ρ_s/ρ exists among these studies. However, ZM has not been reported for heavy spheres confined to a circular tube (Yu et al. Reference Yu, Phan-Thien and Tanner2004; Deloze et al. Reference Deloze, Hoarau and Dušek2012). Surprisingly, here, we report that heavy spheres (ρ_s/ρ ≥ 1.8) confined to a square tube may fall in a zigzagging trajectory for Ga ranging from 210 to 225 (210 ≤ Ga ≤ 225), as shown in figure 4.

Figure 24 shows the trajectories of a heavy sphere (ρ_s/ρ = 1.9, 2, 2.2 and 2.3) freely falling at Ga = 220. It is seen that the sphere vibrates in a perfect zigzagging path in the square tube for ρ_s/ρ ≥ 2 (figure 24b−d). However, the lighter sphere (ρ_s/ρ = 1.9) exhibited a zigzagging-like path (figure 24a). We intend to include this motion because it displays a helical path with oscillations along the x direction that are much stronger than those along the y direction (figure 24a), which is different from RHM. Our study indicates that such zigzagging-like motion exists narrowly for 1.8 ≤ ρ_s/ρ ≤ 2 at similar Galileo numbers. Similar to figure 23, we show effects of the initial particle position on the ZM of a heavy sphere (ρ_s/ρ = 2.3 and Ga = 220) in figure 25. It seems that the sphere underwent ZM only if it was released from a central plane (figure 25a,b). Otherwise (figure 25c−e), we observed helical paths of alternate circling direction (AHM).

Figure 24. Falling trajectories of a heavy sphere with different densities at Ga = 220 showing the ZM for (a) ρ_s/ρ = 1.9, (b) ρ_s/ρ = 2, (c) ρ_s/ρ = 2.2 and (d) ρ_s/ρ = 2.3.

Figure 25. Top views of the trajectories a falling sphere (ρ_s/ρ = 2.3) at Ga = 220 (ZM): effects of the initial position. Note that (a) is the top view of figure 24(d).

We show instantaneous vortical structures of the wake of a heavy sphere (ρ_s/ρ = 2.3) during one oscillation period at Ga = 220 in figure 26. Compared with those at Ga = 205, where the same sphere undergoes AHM (figure 15), vortex shedding is much more noticeable (figure 26). The double-threaded vortices detach from the sphere in succession, form an oblique vortex street downstream, and show no twist because of the sphere's oscillating motion in a fixed plane (the central plane). Another significant difference is that once a double-threaded vortex was generated behind the wavering sphere, it quickly evolved into a vortex loop which helps to prevent the sphere from leaving the central plane and maintain a planar symmetry of the wake. This tendency is more noticeable for light spheres (see also figure 29). However, this was not the case for helically falling spheres (figure 15). In particular, a vortex ring eventually developed downstream as the sphere approaches the wall (figure 26a,c).

Figure 26. Instantaneous iso-surfaces of the vorticity (λ ₂ = −0.01) for a heavy sphere (ρ_s/ρ = 2.3) at Ga = 220 (ZM) during one period (T* ≈ 69): (a) t* = 0, (b) t* ≈ 0.25T*, (c) t* = 0.5T* and (d) t* = 0.75T*.

Similar to the unconfined case, light spheres may exhibit a zigzagging trajectory when confined to a square tube. Our simulations indicate that the ZM is seen for light spheres with ρ_s/ρ ≤ 0.7 when the Galileo number falls within [200, 220], as shown in figure 4. This suggests that the transition from helical to ZM occurs earlier for light spheres than for their counterparts. Note that our critical sphere-to-fluid density ratio (≈0.7) is larger than the values obtained by Horowitz & Williamson (Reference Horowitz and Williamson2010) and Zhou & Dušek (Reference Zhou and Dušek2015), who reported 0.36 (at low Re_T) or 0.61(at high Re_T) and approximately 0.5, respectively. In contrast, some studies (Jenny et al. Reference Jenny, Dušek and Bouchet2004; Auguste & Magnaudet Reference Auguste and Magnaudet2018) revealed that all light spheres (ρ_s/ρ < 1) may oscillate in a zigzagging manner at certain ranges of Ga. This discrepancy warrants further investigation. However, we did not intend to focus on this issue. Our simulations are consistent with those of one of these studies in terms of both the critical ρ_s/ρ and the relevant range of Ga that is 175−215, as revealed by Jenny et al. (Reference Jenny, Dušek and Bouchet2004).

Figure 27 shows oscillating trajectories experienced by the light spheres with different densities (ρ_s/ρ = 0.3–0.7) at Ga = 210. All the results reflect a zigzagging pattern. Moreover, it appears that the ZM is more regular at a larger ρ_s/ρ (figure 27c–e). In other words, the sphere may occasionally deviate from the central plane where it is released from rest if it is much lighter than the fluid (figure 27a). Similar to figure 12, we show time history of the rotational angles of a zigzagging sphere (ρ_s/ρ = 0.7) in figure 28. The trajectory of the sphere is given in figure 27(e) along with notations ‘t ₁ ~ t ₃’ corresponding to those in figure 28. It is seen that the sphere only rotated around the y-axis (θ_y) as it was rising in a wavy manner. From t = t ₁ to t = t ₂, the sphere has a positive rotation direction (figures 28e) while moving from the leftmost position to the rightmost one (figure 27e). We see the opposite from t = t ₂ to t = t ₃. Similar to figure 12(b), the sphere appeared to ‘roll’ forward as it vibrated in the tube.

Figure 27. Rising trajectories of a light sphere with different densities at Ga = 210 showing the ZM: (a) ρ_s/ρ = 0.3, (b) ρ_s/ρ = 0.4, (c) ρ_s/ρ = 0.5, (d) ρ_s/ρ = 0.6 and (e) ρ_s/ρ = 0.7.

Figure 28. Time history of the angles of a light sphere (ρ_s/ρ = 0.7) exhibiting the ZM at Ga = 210. Note that t ₁ ~ t ₃ correspond to those in figure 27(e).

Figure 29. Instantaneous iso-surfaces of the vorticity (λ = −0.01) in the wake of the same light sphere as in figure 27 (Ga = 210).

Figure 29(a‒e) shows instantaneous vortical structures in the wake of the same sphere, as shown in figure 27(a‒e), indicating a more intense motion for the lighter sphere. A vortex loop (hairpin vortex) quickly formed as it extended from each sphere, which is different from the case of the helical motion at low Ga (figure 19) where a pair of threads was generally observed. The vortex shedding became more noticeable as ρ_s/ρ decreases. At least four pairs of the double-threaded vortices were detached from the sphere of ρ_s/ρ = 0.3, along with the appearance of a vortex ring downstream (figure 29a). By contrast, one could hardly see the occurrence of the vortex shedding at ρ_s/ρ = 0.7 (figure 29e). It is also clearly shown that the two sides of the vortex rings did not reconnect fully (figure 27a–c), as revealed by Horowitz & Williamson (Reference Horowitz and Williamson2010). It is worth stating that the vortical structures revealed here (e.g. figures 15, 16, 19, 26 and 29), stemming from the motion of a single sphere, might be of use for understanding the clustering behaviour of settling particles at high Reynolds numbers (Kajishima Reference Kajishima2004; Daniel et al. Reference Daniel, Ecke, Subramanian and Koch2009; Uhlmann & Doychev Reference Uhlmann and Doychev2014) because the wake pattern plays a key role in how particles interact with each other. For instance, sedimenting spheres at Ga = 178 may form a ‘tube-like’ cluster where the particle wakes exhibit complex tangles of vortices extending from particle to particle, as Uhlmann & Doychev (Reference Uhlmann and Doychev2014) noticed. By contrast, for those spheres outside the cluster region the double-threaded vortices were seen in their wakes (Uhlmann & Doychev Reference Uhlmann and Doychev2014).

The ZM may become irregular for a very light sphere owing to the small inertia of the sphere. To illustrate this, we show trajectories of light spheres with ρ_s/ρ ≤ 0.3 in figure 30(a−c), along with their top views in figure 30(d). Compared with figure 27, these results reflect at least three types of irregular ZM in terms of the oscillating orientations of the spheres. For example, the lightest sphere studied (ρ_s/ρ = 0.1) is alternating between the diagonal planes as it rises in the tube (figure 30a). This is rational for a sphere driven by a large buoyant force if one recognizes that the diagonal planes allow a larger space than the central plane. The top panel of figure 30(d) provides a clearer view. For ρ_s/ρ = 0.2 (figure 30b), the sphere exhibits a zigzagging path in one of the diagonal planes. This time the sphere no longer shifts to the other diagonal plane. A further increase ρ_s/ρ leads to the fact that spheres are more likely to stay in the central plane, where they are initially released. Figure 30(c) confirms this tendency. Note that these irregular patterns result directly from the confined condition of our square tube and, therefore, should not appear for an unconfined sphere.

Figure 30. Trajectories of a light sphere showing irregular ZM: (a) ρ_s/ρ = 0.1 (Ga = 200), (b) ρ_s/ρ = 0.2 (Ga = 205), (c) ρ_s/ρ = 0.3 (Ga = 205) and (d) top views of (a‒c).

4.6. Period of oscillation and drag coefficient

In this section, we focus on the oscillation periods and drag coefficients the spheres. A comparison was also performed to illustrate the differences in the dynamical features exhibited by the heavy and light spheres.

Figure 31(a,b) summarizes the normalized oscillation periods (i.e. the inverses of the Strouhal number, T* = 1/St) for the spheres undergoing RHM and ZM, respectively. It can be seen that both heavy and light spheres oscillate much faster when they display a zigzag path. This is reasonable, because when a sphere moves helically in a tube, it is more likely to experience a wall effect. More significantly, it appears that the light spheres oscillate at a similar speed at a fixed Ga, both helically and waveringly. For their counterparts, however, the period increases rapidly as ρ_s/ρ increases (figure 31). Note that we normalized the periods through the time scale T_g, that is, T* = T/T_g. However, a smaller time scale (T_g) was obtained for a heavier sphere when ρ_s/ρ > 1 and for a lighter sphere when ρ_s/ρ < 1 [see (2.7)]. Consequently, T_g varies according to ρ_s/ρ. To illustrate how fast a sphere circles around the tube axis, we normalized the periods at Ga = 190 using the same time scale (i.e. T_g at ρ_s/ρ = 0.8 or 1.2) and show them ($T_r^\ast $) in figure 32. It is interesting to find that the value of $T_r^\ast $ is ‘rising’ with ρ_s/ρ for the rising spheres while the opposite is seen for the falling spheres, suggesting an increase in the particle inertia results in a different dynamical response from a heavy sphere and a light one (figure 32a). In addition, the period ($T_r^\ast $) is rapidly increasing as ρ_s/ρ approaches one (the neutrally buoyant case) from both sides.

Figure 31. Oscillation periods for both heavy (ρ_s/ρ > 1) and light (ρ_s/ρ < 1) spheres undergoing: (a) the RHM (Ga = 170 and 190), and (b) ZM (Ga = 210 and 220), respectively.

Figure 32. Oscillation periods ($T_r^\ast $) normalized by the time scale at ρ_s/ρ = 0.8 (or 1.2) for Ga = 190 depending on (a) ρ_s/ρ and (b) |ρ_s/ρ − 1|, respectively.

However, we show the relationship between the period and |ρ_s/ρ − 1| at Ga = 190 in figure 32(b), which may provide a better understanding of the effects of the sphere-to-fluid density ratio. A similar dependence of $T_r^\ast $ on the accelerating factor (|ρ_s/ρ − 1|) is seen for both types of the spheres. A sphere falling or rising in the tube circles the tube axis faster as |ρ_s/ρ − 1| increases (figure 32b). Surprisingly, a light sphere always achieves a higher circling speed than a heavy sphere when it has the same driving force (i.e. |ρ_s/ρ − 1|). This may be primarily owing to the confined conditions of the square tube. To address this, we show limit cycles of the helical trajectories projected onto the (x, y) plane for heavy spheres (figure 33) and light spheres (figure 34) at Ga = 190. It can be clearly observed that each heavy sphere exhibits an ellipse-shaped limit cycle that is more elliptical at larger ρ_s/ρ (figure 33). However, the limit cycles remained nearly circular for all the light spheres (figure 34). This is due to the difference in particle inertia between the heavy and light spheres, which respond differently to wall repulsion, as mentioned above. We can expect a stronger wall effect for a heavy sphere because it is more likely to approach the tube walls as it runs along an ellipse-shaped path. Consequently, the movement of the heavy spheres was hindered by the tube walls more noticeably than their counterparts.

Figure 33. Top views of helical trajectory for heavy spheres with different densities at Ga = 190: (a) ρ_s/ρ = 1.2, (b) ρ_s/ρ = 1.5, (c) ρ_s/ρ = 1.7 and (d) ρ_s/ρ = 1.8. The limit cycle is indicated by a thick line for each case (the same below).

Figure 34. Top views of particle trajectory for light spheres with different densities at Ga = 190: (a) ρ_s/ρ = 0.25, (b) ρ_s/ρ = 0.5, (c) ρ_s/ρ = 0.6 and (d) ρ_s/ρ = 0.8.

Figure 35 provides a direct view of the wall ‘hindering’ effect by displaying the 2-D vorticity ($\omega _z^\ast $) for a light (ρ_s/ρ = 0.8) and heavy (ρ_s/ρ = 0.8) spheres when they move close to the right and upper walls, respectively. Note that both spheres have the same driving force (|ρ_s/ρ − 1|). The sphere–wall interaction, characterized by the magnitude of $\omega _z^\ast $, is more significant for the heavy sphere (figure 35b,d) than for the light one (figure 35a,c) whenever they approach a tube wall, suggesting a stronger ‘hindering’ effect experienced by the heavy sphere. For instance, the heavy sphere attained a lower speed (≈0.13U_g) as compared with the light one (≈0.14U_g) when both spheres moved close to the right wall (figure 35a,b). In particular, both spheres showed a similar speed (≈0.14U_g) despite a larger distance from the upper wall seen for the heavy sphere (figure 35c,d). This leads to a smaller period (or a faster circling motion) for light spheres (figure 34b).

Figure 35. Instantaneous contours of the 2-D vorticity ($\omega _z^\ast $) for a light (a,c: ρ_s/ρ = 0.8) and heavy (b,d: ρ_s/ρ = 1.2) spheres, respectively. Note that (a,b) represent results when the spheres are closest to the right wall while (c,d) represent those when the spheres are closest to the upper wall.

It is known that both the terminal velocity and drag coefficient of an unconfined sphere are considerably dependent on the Galileo number. Similarly, we verified the dependence of Re_T and C_D of our confined sphere on Ga in figure 36(a,b). The drag coefficient was computed through the following expression:

(4.1)\begin{equation}{C_D} = \left|{\frac{{{\rho_s}}}{\rho } - 1} \right|\frac{{2gd}}{{3U_T^2}}.\end{equation}

The solutions of the empirical formulations, that is (1.4) and (1.5) from Cabrera-Booman et al. (Reference Cabrera-Booman, Plihon and Bourgoin2024) are also provided for reference. The general features of our results (both Re_T and C_D) agree with those of previous studies (Cabrera-Booman et al. Reference Cabrera-Booman, Plihon and Bourgoin2024). A lower Re_T and larger C_D were observed for our spheres because of the wall effects, which is rational. In addition, at steady state, similar terminal velocities were obtained for both the confined and unconfined spheres (figure 36a), irrespective of ρ_s/ρ. The onset of flow unsteadiness led to an increasing discrepancy between our results and the solutions in the literature as Ga increased (figure 36).

Figure 36. Dependence of (a) the terminal Reynolds number and (b) drag coefficients on the Galileo number. The relationship between Re_T and Ga and that between C_D and Ga from Cabrera-Booman et al. (Reference Cabrera-Booman, Plihon and Bourgoin2024) are indicated for reference.

As far as we know, the influence of the sphere-to-fluid density ratio on the drag coefficient is an arguable issue. It can be inferred from (1.5), C_D is independent of ρ_s/ρ. This is reasonable because (1.5) was developed based on the C_D−Re_T relationship for a fixed sphere (Brown & Lawler Reference Brown and Lawler2003). Interestingly, an experimental study on the free motion of spheres (ρ_s/ρ = 0.87, 1.12, 3.19 and 3.9) at 100 < Ga < 700 reflected a similar behaviour (Raaghav et al. Reference Raaghav, Poelma and Breugem2022): their drag coefficients nearly coincided with those of a fixed sphere. However, previous studies (Karamanev & Nikolov Reference Karamanev and Nikolov1992; Horowitz & Williamson Reference Horowitz and Williamson2010; Auguste & Magnaudet Reference Auguste and Magnaudet2018) have revealed that a very light sphere (${\rho _s}/\rho \ll 1$) with vibrations experiences a significantly larger drag coefficient than a fixed sphere. In contrast, Cabrera-Booman et al. (Reference Cabrera-Booman, Plihon and Bourgoin2024) reported a systematic effect of the sphere-to-fluid density ratio on C_D of heavy spheres. In other words, the drag coefficients of very heavy spheres (large ρ_s/ρ) were considerably smaller than those of fixed spheres. The opposite is true for nearly neutrally buoyant spheres (Cabrera-Booman et al. Reference Cabrera-Booman, Plihon and Bourgoin2024). Things were different for our confined spheres. From figure 36(b) it seems that the effect of ρ_s/ρ on C_D cannot be neglected. Small but visible deviations in the drag coefficients are seen between the result of ρ_s/ρ = 2.3 and the other ones (0.1 ≤ ρ_s/ρ ≤ 2). The C_D of our heaviest sphere studied (ρ_s/ρ = 2.3) shows an opposite trend as compared with the observations of Cabrera-Booman et al. (Reference Cabrera-Booman, Plihon and Bourgoin2024). The primary reason for this is the same as mentioned above. The tube walls had a stronger ‘hindering’ effect on the heaviest sphere than on the other spheres.

Finally, it should be stated that the transition scenario illustrated in figure 4 only applies to the confined conditions of a 5d square tube. The results essentially differ from those seen for an unconfined sphere. So, the tube width-to-sphere diameter ratio (L/d) is certainly a key factor in determining the motion patterns of the spheres as well as the boundaries between different regimes. Due to high costs of the computations, we do not provide a quantitative analysis on the effect of L/d. For a larger tube, the irregular helical motions could only occur for lighter spheres (AHM) and heavier spheres (JHM), respectively, due to requirement of a stronger sphere–wall interaction. Both motions will switch to the regular one (RHM) if the tube width is larger than a threshold. The ZM of heavy spheres is sensitive to the confined conditions (see figure 25). As a result, the zigzagging path may disappear quickly for heavy spheres when L/d increases.