Freely rising or falling of a sphere in a square tube at intermediate Reynolds numbers

Abstract

The motion of a sphere freely rising or falling in a 5d (d is the diameter of the sphere) square tube was numerically studied for the sphere-to-fluid density ratio ranging from 0.1 to 2.3 (0.1 ≤ ρ_s/ρ ≤ 2.3, ρ_s is the density of spheres and ρ the fluid density) and Galileo number from 140 to 230 (140 ≤ Ga ≤ 230). We report that Hopf bifurcation occurs at Ga_crit ≈ 157, where both the heavy and light spheres lose stability. The helical motion is widely seen for all spheres at Ga > 160 resulting from a double-threaded vortex interacting with the tube walls, which becomes irregular at Ga ≥ 190 where heavy spheres act differently from their counterparts; that is, heavy spheres change their helical directions alternately while light spheres exhibit helical trajectories with jaggedness in connection with the shedding of the double-threaded vortices. This is because of the difference in inertia between the heavy and light spheres. We also checked the oscillation periods for the helical motion of the spheres. They show opposite variations with ρ_s/ρ for the two types of spheres. Light spheres (ρ_s/ρ ≤ 0.7) reach a zigzagging regime at Ga ≥ 200 where a vortex loop (hairpin-like vortical structure) is formed which may develop into a vortex ring downstream at small ρ_s/ρ. This might be the first time a transition from the helical motion to the zigzagging motion for heavy spheres (ρ_s/ρ ≥ 1.8) has been reported. Finally, we examined the dependence of both the terminal Reynolds number and the drag coefficient of the spheres on the Galileo number.

1. Introduction

The free motion of solid particles in air or water is a common phenomenon in nature and industry. Heavy particles generally fall, while light particles rise, owing to competition between the effects of gravity and buoyancy. One may think that this is a simple problem to understand and can be mathematically solved. However, existing studies indicate that this common phenomenon reflects a complex and rich dynamics when the inertia of both particles and fluids becomes significant (Ern et al. 2012). Owing to its importance from both theoretical and practical perspectives, most relevant studies have been performed numerically or experimentally at moderate Reynolds numbers. For the free motion of a single sphere in an unbounded domain, much effort (Natarajan & Acrivos 1993; Mittal 1999; Ormières & Provansal 1999; Ghidersa & Dušek 2000; Tomboulides & Orszag 2000; Fabre, Auguste & Magnaudet 2008; Chrust, Goujon-Durand & Wesfreid 2013) was devoted to the reverse case in the early years, that is, the case of uniform flow past a fixed sphere. Accordingly, the transition scenario for a fixed sphere is well understood in terms of the onset of instability, loss of flow planarity and vortex dynamics. For instance, it is widely believed that the flow loses its axial symmetry at a Reynolds number (Re) of approximately 212 (Natarajan & Acrivos 1993; Johnson & Patel 1999; Tomboulides & Orszag 2000; Chrust et al. 2013), in association with the formation of a double-threaded vortex in the wake that exhibits planar symmetry. The second transition occurred at Re = 270–280, representing the Hopf bifurcation (Natarajan & Acrivos 1993; Johnson & Patel 1999; Ghidersa & Dušek 2000). As Re increased, the shedding of the double-threaded vortex occurred and a hairpin-like vortex formed downstream of the sphere. At a high Reynolds number (Re ≈ 420, Sakamoto & Haniu 1995; Re = 350–375, Mittal 1999; Re ≈ 500, Tomboulides & Orszag 2000; Re ≈ 375, Chrust et al. 2013) the flow falls into chaos, leading to disappearance of the planarity.

Things become more complicated when a sphere can freely fall or rise in a fluid. As previously reported (Karamanev & Nikolov 1992; Karamanev, Chavarie & Mayer 1996), light spheres may experience a significantly larger drag than the standard drag curve (Turton & Levenspiel 1986). Horowitz & Williamson (2010) indicated that the increasing drag coefficient of a light sphere (ρ_s/ρ < 0.36, ρ_s is the density of spheres and ρ the fluid density) could be nearly twice that of a fixed sphere because of vibration at Re > 260. However, Auguste & Magnaudet (2018) argued that a deviation from the standard drag could only be observed for light spheres (ρ_s/ρ ≤ 0.1) undergoing a helical motion at high Reynolds numbers (Re > 350). This contradicts the findings of an experimental study by Veldhuis, Biesheuvel & Lohse (2009), who claimed that a very light sphere (ρ_s/ρ = 0.02) experiences enhanced drag as it rises helically or along a zigzag path. A recent experiment work (Will & Krug 2021a) revealed the existence of a helical regime for light spheres with ρ_s/ρ ≤ 0.42 where the drag is significantly increased. They believed that a change in the wake structure leads to the higher drag. Further attention should be paid to addressing the disagreement in the ranges of both ρ_s/ρ and Re where the deviation from the standard drag curve occurs and how large the deviation is.

Another concern is the wake-induced motion patterns of the spheres. Jenny, Dušek & Bouchet (2004) presented a pioneering study on the motion and dynamics of spheres over a broad region of sphere-to-fluid density ratio and Galileo number, which is defined as follows:

(1.1)\begin{equation}Ga = \sqrt {\left|{\frac{{{\rho_s}}}{\rho } - 1} \right|\frac{{g{d^3}}}{{{\nu ^2}}}} ,\end{equation}

where d is the diameter of the spheres, g is the gravitational acceleration and ν is the kinematic viscosity. A variety of motions were revealed for a sphere under the action of gravity: vertical steady, oblique steady, oblique oscillating, zigzagging and chaotic states. More significantly, Jenny et al. (2004) established a regime diagram in the parametric space (Ga, ρ_s/ρ), which provided a clear view of the transition scenario and bifurcations of the motion of falling or rising spheres. They reported that a zigzagging motion existed for all light spheres (ρ_s/ρ < 1) for Ga values ranging between 175 and 215. A subsequent experimental study (Veldhuis & Biesheuvel 2007) confirmed ‘the main features’ of the motion patterns revealed by Jenny et al. (2004), despite the small discrepancy observed in the frequency of the oscillating regime. Using the same numerical scheme as Jenny et al. (2004), Zhou & Dušek (2015) performed a systematic study on the same problem and provided an updated regime map that, in principle, resembled the previous one. In contrast, several new patterns (helical/rotating and vertical periodic patterns) were reported in the latter. The regimes for the falling of heavy spheres (Zhou & Dušek 2015) agree well with those of a recent experimental study (Cabrera-Booman, Plihon & Bourgoin 2024). A new (vertical periodic) pattern was confirmed by Cabrera-Booman et al. (2024), who performed experiments on settling spheres with density ratio close to unity (ρ_s/ρ = 1.1). However, Zhou & Dušek (2015) claimed light spheres with ρ_s/ρ ≤ 0.5 may undergo the zigzagging motion. This density ratio range considerably differs from that of Jenny et al. (2004) as well as those reported by several other studies (ρ_s/ρ ≤ 0.36 at Re = 260–1550 and ρ_s/ρ ≤ 0.61 at Re = 1550–15 000, Horowitz & Williamson 2010; ρ_s/ρ < 1, Auguste & Magnaudet 2018; Raaghav, Poelma & Breugem 2022). Raaghav et al. (2022) reported that a moderately light sphere (ρ_s/ρ ≈ 0.87) exhibits a zigzagging motion at Ga values ranging between 190 and 250. The reason for the diversity in the critical ρ_s/ρ remains unclear. It is worth mentioning here that Cabrera-Booman et al. (2024) derived formulations of the drag coefficient (C_D) and terminal Reynolds number (Re_T = U_Td/ν, U_T is the terminal velocity of freely falling or rising spheres) as a function of the Galileo number, based on the relationship between C_D and Ga

(1.2)\begin{equation}{C_D} = \frac{4}{3}{\left( {\frac{{Ga}}{{R{e_T}}}} \right)^2}.\end{equation}

To solve (1.2), Cabrera-Booman et al. (2024) adopted the drag coefficient of a fixed sphere proposed by Brown & Lawler (2003)

(1.3)\begin{equation}{C_D} = \frac{{24}}{{R{e_T}}}(1 + 0.15\,Re_T^{0.681}) + \frac{{0.407}}{{1 + \dfrac{{8710}}{{R{e_T}}}}}.\end{equation}

The following formulations were obtained (Cabrera-Booman et al. 2024):

(1.4)\begin{gather}R{e_T}(Ga) = \frac{{G{a^2}(22.5 + G{a^{1.364}})}}{{0.0258G{a^{2.6973}} + 2.81G{a^{2.0306}} + 18G{a^{1.364}} + 405}},\end{gather}

(1.5)\begin{gather}{C_D}({Ga} )= \frac{4}{3}{\left( {\frac{{Ga(22.5 + G{a^{1.364}})}}{{0.0258G{a^{2.6973}} + 2.81G{a^{2.0306}} + 18G{a^{1.364}} + 405}}} \right)^{ - 2}}.\end{gather}

Equations (1.4) and (1.5) allow one to directly assess the terminal velocity (or terminal Reynolds number) and drag coefficient of falling or rising spheres, respectively, because Ga is an input parameter instead of an output parameter, such as Re_T. Broadly speaking, existing studies show reasonable agreement on the motion regimes of a falling or rising sphere at low Ga and moderate ρ_s/ρ. As Ga increases, discrepancies emerge among these studies, especially for very heavy and light spheres. For instance, both studies (Jenny et al. 2004; Zhou & Dušek 2015) revealed an oblique oscillating motion with high frequency for heavy spheres (ρ_s/ρ > 2.5), but Raaghav et al. (2022) did not observe such a regime for spheres with ρ_s/ρ = 3.19 and 3.9 at relevant Ga. By contrast, the recent study (Cabrera-Booman et al. 2024) confirmed this regime through experiments on spheres with ρ_s/ρ = 7.9. Furthermore, according to Zhou & Dušek (2015), a vertical periodic pattern occurred at Ga > 250, which was not reported in the experimental work of Veldhuis & Biesheuvel (2007). The most noticeable disagreement may exist between the study by Horowitz & Williamson (2010) and other studies in terms of the range of Ga or Re for the zigzagging motion of light spheres and the motion patterns of heavy spheres. Numerical (Zhou & Dušek 2015) and experimental (Will & Krug 2021a; Raaghav et al. 2022) work revealed that a light sphere will rise in a helical path at high Ga. However, Horowitz & Williamson (2010) did not report this motion. Substantial efforts are required to shed light on the gravity/buoyancy-driven motion of spheres in the broader parametric space of (Ga, ρ_s/ρ) to provide a comprehensive and indisputable regime map.

In contrast, little is known about confined spheres under similar flow conditions. The presence of solid walls raises several issues which are absent for an unconfined sphere, such as the inertial migration (Segrè & Silberberg 1961; Matas, Morris & Guazzelli 2004; Yu, Phan-Thien & Tanner 2004; Yang et al. 2005), wake-induced lift (Bagchi & Balachandar 2002; Zeng, Balachandar & Fischer 2005; Feng, Gatewood & Michaelides 2021; Shi et al. 2021) and the suppression of flow instability (Zeng et al. 2005; Zhao et al. 2016). Yu et al. (2004) used a distributed Lagrange multiplier method to simulate the falling of a heavy sphere (ρ_s/ρ = 1.5) in a 5d circular tube at several Reynolds numbers, i.e. Re = 20, 100, 200, 300 and 400 (corresponding to Re_T ≈ 8, 84, 197, 306 and 424). These Reynolds numbers represent four different regimes of flow past a fixed sphere (Natarajan & Acrivos 1993; Johnson & Patel 1999; Tomboulides & Orszag 2000): unseparated flow (Re = 20), axisymmetric flow (Re = 100 and 200), periodic oscillatory flow (Re = 300) and chaotic flow (Re = 400). The confined sphere was seen to fall in a helical path at Re ≥ 200 (Yu et al. 2004), which essentially differs from the unconfined case for which one could only observe the helical motion of heavy spheres (ρ_s/ρ > 1.7) at a bi-stable region (212 < Ga < 240) where the chaotic motion and a high-frequency oscillating oblique motion coexist. As indicated by Deloze, Hoarau & Dušek (2012), spheres confined to a tube are quickly subjected to a wall repulsive interaction and tend to fall helically. Furthermore, Yu et al. (2004) checked the role of rotation on the lateral motion of the sphere (i.e. the Magnus effect). Note that Will & Krug (2021b) also addressed the significance of rotational dynamics in the motion of falling and rising spheres. Yu et al. (2004) provided a careful analysis on the vortex dynamics in the wake of the sphere and characterized the formation and evolution of a hairpin-like vortex at Re = 400. Similarly, Deloze et al. (2012) focused on the motion transitions of a heavy sphere (mostly at ρ_s/ρ = 2) within the same tube. The critical Galileo number, at which the flow loses its axial symmetry and instead attains a planer symmetry, was found to be lie somewhere between 155 and 160 (Deloze et al. 2012), slightly higher than that of an unconfined sphere (Ga_crit ≈ 155, Fabre et al. 2008). This bifurcation is of the Hopf type because the sphere begins to oscillate in the tube and quickly falls into the helical regime as Ga increases (Deloze et al. 2012). They predicted a possible jump in the vortex shedding frequency at Ga = 250 as ρ_s/ρ increases from 3 to 5. In the case of an unconfined sphere, this occurs at ρ_s/ρ values between 2 and 3, as Jenny et al. (2004) indicated. Despite considerable progress made in the past, the motion regimes and transition scenarios of falling or rising spheres confined to tubes are not well understood. Because of the limited number of computations, neither study (Yu et al. 2004; Deloze et al. 2012) established a regime map for the motion of a confined sphere in a relatively broad parametric space (Ga, ρ_s/ρ), unlike in the unconfined case. Furthermore, light spheres have not been considered in the literature. Considering what we know about an unconfined sphere, one may ask whether a confined sphere can run in a zigzag path or how heavy spheres behave differently from their counterparts in a tube. The wall effect on the drag coefficient is also a concern. In view of this, we performed a considerable number of direct numerical simulations of a sphere freely falling or rising in a 5d square tube within the ranges of 140 ≤ Ga ≤ 230 and 0.1 ≤ ρ_s/ρ ≤ 2.3, hoping to present a comprehensive understanding. To the best of our knowledge, there have been no similar report on a square tube. Furthermore, we can expect more complex dynamic features from a sphere released in a square tube as compared with a circular tube.

The remainder of this paper is organized as follows. In § 2, we briefly introduce the three-dimensional (3-D) lattice Boltzmann method (LBM) along with the boundary conditions. The problem description and key parameters are also included. Code validation is presented in § 3. In § 4, we provide the simulation results and a discussion on the critical Galileo number, helical motion and zigzagging motion of both heavy and light spheres. In § 4 we show a regime map of sphere motion in the Ga − ρ_s/ρ space. The oscillation periods and drag coefficients of both the heavy and light spheres are summarized. In § 5, the concluding remarks are provided.

2. Numerical method and problem description

2.1. Lattice Boltzmann method

The fluid motion was solved using a 3-D lattice Boltzmann method. The discrete lattice Boltzmann equations of a single-relaxation-time model are expressed as (Qian 1992)

(2.1)\begin{equation}{f_i}(\boldsymbol{x} + {\boldsymbol{e}_i}\Delta t,t + \Delta t) - {f_i}(\boldsymbol{x},t) =- \frac{1}{\tau }(\,{f_i}(\boldsymbol{x},t) - {f_i}^{(eq)}(\boldsymbol{x},t)),\end{equation}

where f_i(x, t) is the distribution function for the microscopic velocity e_i in the ith direction, $f_i^{(eq)}(\boldsymbol{x},t)$ is the equilibrium distribution function, Δt is the time step of the simulation, τ is the relaxation time and w_i are weights related to the lattice model. The fluid density ρ and velocity u are determined by the distribution function

(2.2a,b)\begin{equation}\rho = \sum\limits_i {{f_i}} ,\quad \rho \boldsymbol{u} = \sum\limits_i {{f_i}{\boldsymbol{e}_i}} .\end{equation}

For the D3Q19 (19 speeds in three dimensions) lattice model used in this study, the discrete velocity vectors (shown in figure 1) are given by

(2.3) \begin{equation}{\boldsymbol{e}_i} = \left\{ {\begin{array}{*{20}{@{}ll}} {(0,0,0)c}&{i = 0}\\ {({\pm} 1,0,0)c,(0, \pm 1,0)c,(0,0, \pm 1)c}&{i = 1\textrm{--6}}\\ {({\pm} 1, \pm 1,0)c,({\pm} 1,0, \pm 1)c,(0, \pm 1, \pm 1)c}&{i = 7\unicode{x2013} 18} \end{array}} \right.,\end{equation}

where c = Δx/Δt, and Δx is the lattice spacing. The speed of sound (c_s) was determined as ${c_s} = c/\sqrt 3 $.

Figure 1. Discrete velocity vectors of D3Q19 lattice model (e₀ ~ e₁₈).

Following Qian (1992), the equilibrium distribution function is chosen as follows:

(2.4)\begin{equation}{f^{(eq)}}(\boldsymbol{x},t) = \rho {w_i}\left[ {1 + \frac{{3({\boldsymbol{e}_i}\boldsymbol{\cdot }\boldsymbol{u})}}{{{c^2}}} + \frac{{9{{({\boldsymbol{e}_i}\boldsymbol{\cdot }\boldsymbol{u})}^2}}}{{2{c^4}}} - \frac{{3(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u})}}{{{c^2}}}} \right],\end{equation}

where the weights are set to be w ₀ = 1/3, w _1–6 = 1/18 and w _7–18 = 1/36.

By performing Taylor and Chapman–Enskog expansions, the macroscopic mass and momentum equations in the low-Mach-number limit can be recovered from (2.1)

(2.5)\begin{gather}\frac{{\partial \rho }}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}) = 0,\end{gather}

(2.6)\begin{gather}\frac{{\partial (\rho \boldsymbol{u})}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{uu}) =- \boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }[\rho \nu (\boldsymbol{\nabla }\boldsymbol{u} + {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}})].\end{gather}

A linear relation between the kinematic viscosity and the relaxation time is achieved, i.e. $\nu = c_s^2\mathrm{\Delta }t$ (τ − 0.5). Lattice units are commonly used instead of physical units in lattice Boltzmann simulations. For example, Δx and Δt represent 1 lu (lattice unit) and 1 ts (time step), respectively. This resulted in c = 1 and ${c_s} = 1/\sqrt 3 $, both in a lu/ts unit. Using these lattice units one could develop LBM codes in a simple manner.

2.2. Problem

In this work, we aim to provide a comprehensive study of the motion and wake patterns of a sphere falling or rising in an infinite square tube. The spheres are assumed to be uniform in density. As shown in figure 2, the sphere diameter and density are denoted as d and ρ_s, respectively, leading to the mass and moment of inertia of the sphere being M = ${\rm \pi} $ρ_sd ³/6 and I = Md ²/10, respectively. The square tube is filled with a fluid of density ρ and kinematic viscosity ν. The width of the tube was fixed at L = 5d. There were three reasons for choosing 5d: first, the effect of the tube walls is expected to be considerably large at this size. We may be able to reveal primary regimes and transitions of moving spheres confined to a square tube at intermediate Reynolds numbers. Second, several studies (Yu et al. 2004; Deloze et al. 2012) focused on the falling of a heavy sphere in a 5d circular tube. Therefore, we could possibly make a reasonable comparison between our study and previous analyses. Finally, the not-too-large tube size allowed us to perform 3-D simulations at an affordable computational cost. Similar to most of the relevant studies, both the sphere-to-fluid density ratio (ρ_s/ρ) and the Galileo number (Ga) are adopted as the control parameters for our system (figure 2), with ρ_s/ρ ranging from 0.1 to 2.3, and Ga from 140 to 230. These ranges allowed us to provide a fundamental understanding of the transition scenario of the gravity/buoyancy-driven motion of both light and heavy spheres.

Figure 2. Schematic diagram of the present problem. A heavy sphere (a: ρ_s/ρ > 1) or a light one (b: ρ_s/ρ < 1) is freely falling or rising in an infinite square tube which is filled with a fluid with density ρ and kinematic viscosity ν. Note that d and ρ_s denote the diameter and density of the spheres, respectively. The width of the tube was fixed at 5d, i.e. L = 5d. We used a moving grid scheme to make sure the spheres stay at a vertical position having distances of 12d and 25d from the upstream and the downstream, i.e. H_u = 12d and H_d = 25d, respectively. This results in a total height of the domain of 37d, i.e. H = H_u + H_d = 37d.

The sphere was initially placed in the central plane (y = 0 plane) at a distance d from the tube axis. The vertical position of the sphere was 12d apart from its upstream (H_u = 12d) and 25d from its downstream (H_d = 25d). Consequently, the height of the computational domain was H = 37d. In the simulations, a moving grid scheme was applied to ensure that the sphere maintained a distance of 25d from the downstream. In doing so, the fluid field and the position of the sphere are shifted upward/downward by one lattice spacing once the sphere falls/rises below a height of 12d from upstream. No-slip boundary conditions were imposed on all the tube walls. The normal derivative of the velocity is zero downstream, and the velocity upstream is zero. In accordance with most of previous studies, we chose U_g as the velocity scale, which is given by

(2.7)\begin{equation}{U_g} = \sqrt {\left|{\frac{{{\rho_s}}}{\rho } - 1} \right|gd} .\end{equation}

The time scale was T_g = d/U_g. Note that this velocity scale could be obtained from a force balance of gravity, drag and buoyancy for a given drag coefficient. For convenience, we use U_x, U_y and U_z to denote the three components of the velocity of the spheres, which are consistent with the coordinate system shown in figure 2. These components were normalized through U_g, that is, $U_x^\ast= {U_x}/{U_g}$, $U_y^\ast= {U_y}/{U_g}$ and $U_z^\ast= {U_z}/{U_g}$. Some parameters were fixed as follows: d = 32, ρ = 1 and g = 9.8 × 10⁻⁴. For all cases considered the initial positions (x ₀, y ₀, z ₀) of the heavy and light spheres are (−d, 0, 12d) and (−d, 0, 25d), respectively, unless otherwise stated.

2.3. Moving boundary conditions

To ensure a no-slip boundary condition on a moving surface, special treatments are usually required to redistribute the distribution functions (f_i) adjacent to the surface. In this study, we adopted an improved bounce-back scheme proposed by Lallemand & Luo (2003) to address the boundary conditions of moving spheres. This scheme, which originates from the bounce-back rule and is based on a second-order interpolation method, is easily implemented for 3-D problems and was successfully applied in our previous study (Nie & Lin 2020). Readers are referred to Lallemand & Luo (2003) for more details.

Once all distribution functions at the fluid nodes surrounding a sphere are known, the hydrodynamic force and torque experienced by the sphere can be obtained through a momentum exchange scheme (Mei et al. 2002; Lallemand & Luo 2003). In addition, to account for the effects of a sphere moving in a fluid, the method proposed by Aidun, Lu & Ding (1998) was used to calculate the added force and torque due to the covered and uncovered fluid nodes. Using the net force and torque, the motion of the spheres was determined by solving Newton's equations. Note that Peng et al. (2016) investigated the numerical performance in LBM simulations of particulate flows and carefully checked the impact of both the momentum exchange and interpolated bounce-back schemes used here. These implementations may lead to fluctuations in calculating the hydrodynamic forces of a moving particle, as Peng et al. (2016) revealed. However, one could still obtain generally accurate results by introducing Aidun's correction process (Peng et al. 2016), as demonstrated by our previous studies (Nie & Lin 2019; Nie & Lin 2020; Nie et al. 2023).

3. Validation

To validate our code and check the grid resolution, we performed simulations of a heavy sphere (ρ_s/ρ = 1.5) settling in an unbounded domain (7d × 7d × 25d) at two Galileo numbers, Ga = 144 and Ga = 178, representing patterns of steady vertical motion and steady oblique motion (Jenny et al. 2004), respectively. Periodic boundary conditions were implemented in both lateral directions.

In table 1, we show the terminal velocities of the sphere (horizontal: $U_h^\ast $ and vertical: $U_z^\ast $) obtained from different grid resolutions, i.e. d = 16, 24, 32 and 40, and compare our results with previous values from Uhlmann & Dušek (2014). Note that the horizontal velocity of the sphere is determined through $U_h^\ast= {(U_x^{{\ast} 2} + U_y^{{\ast} 2})^{0.5}}$, which is zero at Ga = 144 (steady vertical motion). Good agreement was observed in the terminal velocities ($U_h^\ast $ and $U_z^\ast $) for both Galileo numbers. Note that the literature (Uhlmann & Dušek 2014) provided two values for the terminal velocity at each Ga, which were obtained from a small lateral domain and a large one. We chose the latter for comparison and revealed that the relative errors were below 2% for both Ga. The results given by (1.4) are also included in table 1, showing no apparent discrepancy with our calculations ($U_z^\ast $). In addition, the values obtained from d = 32 are almost the same as those from a finer grid (d = 40), which makes d = 32 a good choice for the study.

Table 1. Terminal velocities of a heavy sphere (ρ_s/ρ = 1.5) at steady state for two Ga numbers.

^a Note that the horizontal velocity of the sphere is determined through $U_h^\ast= {(U_x^{{\ast} 2} + U_y^{{\ast} 2})^{0.5}}$.

Figure 3 further shows iso-surfaces of the vertical fluid velocity (figure 3a,c), where $U_z^\ast=- 0.1$, as well as those of the vorticity (figure 3b,d), where λ ₂ = −0.01, at Ga = 144 (figure 3a,b) and Ga = 178 (figure 3c,d). Here, λ ₂ is the second largest eigenvalue of the tensor S² + $\boldsymbol{\varOmega}$² (where S and $\boldsymbol{\varOmega}$ are the symmetrical and anti-symmetrical components of the velocity gradient tensor, i.e. $\boldsymbol{\nabla }\boldsymbol{u}$, respectively), represented by λ ₂ (Jeong & Hussain 1995), to identify the vortical structures in the wake of falling or rising spheres. The flow was perfectly axisymmetric at Ga = 144 (figure 3a,b). In contrast, an oblique wake was observed at Ga = 178 (figure 3c,d). The axial symmetry is replaced by a planar symmetry (Johnson & Patel 1999; Mittal 1999; Ghidersa & Dušek 2000; Chrust et al. 2013). Figure 3(d) provides a clear view of the double-threaded vortical structure generated by the falling sphere in accordance with previous studies (Uhlmann & Dušek 2014; Zhou & Dušek 2015). Both threads were steady and off centred, producing a horizontal force on the sphere. This is how the sphere underwent an oblique falling motion at Ga = 178. Note that the double-threaded wake was first reported by Magarvey & Bishop (1961), who used dye visualization to display the wake generated by liquid drops falling in water. In the range of 210 < Re < 270, they observed that the flow became non-axisymmetric and the dye was released from the wake in two parallel threads (Magarvey & Bishop 1961). The double-threaded vortical structure was also observed in the wake of flow past a fixed sphere (Natarajan & Acrivos 1993; Johnson & Patel 1999; Tomboulides & Orszag 2000) at Re ≈ 212 where a transition from the axisymmetric wake to non-axisymmetric wake occurs.

Figure 3. Simulation results at (a,b) Ga = 144 and (c,d) Ga = 178, respectively. Note that (a,c) show the iso-surfaces of the vertical fluid velocity where $u_z^\ast=- 0.1$ and (b,c) show the iso-surfaces of the vorticity where λ ₂ = −0.01. We use the second largest eigenvalue of the tensor S² + $\boldsymbol{\varOmega}$² (where S and $\boldsymbol{\varOmega}$ are the symmetrical and anti-symmetrical components of the velocity gradient tensor, i.e. $\boldsymbol{\nabla }\boldsymbol{u}$, respectively), represented by λ ₂ (Jeong & Hussain 1995), to identify the vortical structures in the wake of a falling or rising sphere (the same below). Note that it seems the heavy sphere is purple in colour because the iso-surface results were made 60% transparent. The sphere is actually red in colour as indicated in figure 2 (the same below).

4. Results

To provide an overview of the study, we show the motion patterns of the spheres confined to a 5d square tube in figure 4, which also illustrates the transition scenario between different regimes. Roughly speaking, the spheres (0.1 ≤ ρ_s/ρ ≤ 2.3) exhibit four kinds of motion when Ga ranges from 140 to 230 (140 ≤ Ga ≤ 230), i.e. the vertical motion, helical motion (regular helical motion, alternate helical motion and jagged helical motion), zigzagging motion (ZM) and chaotic motion.

Figure 4. Motion patterns of the spheres rising or falling in a 5d square tube in the parametric space of (Ga, ρ_s/ρ) for Ga ranging between 140 and 240, and ρ_s/ρ from 0.1 to 2.3.

Both heavy and light spheres will lose steadiness at a similar Galileo number (Ga_crit ≈ 157), which is discussed in § 4.2. Before that (Ga < Ga_crit) the spheres were seen to fall or rise vertically in the tube, exhibiting a steady vertical motion for all ρ_s/ρ studied (§ 4.1). A transition regime (transitional oscillatory motion, TOM) is seen for the heavy spheres after the onset of unsteadiness, which is not the case for the light ones. The TOM features small oscillations in an arbitrary plane (Deloze et al. 2012).

Helical motion was widely observed in both the heavy and light spheres at Ga > 160. The regular helical motion (RHM) is characterized by its unchanged helical direction and a constant helical size. A type of irregular helical motion (alternate helical motion, AHM) occurs at Ga ≥ 190 for all heavy spheres studied, which will be discussed in § 4.3. For spheres undergoing AHM, they change the helical direction in an alternate way as sedimenting in the tube. By contrast, very light spheres (ρ_s/ρ ≤ 0.3) undergo another type of irregular helical motion (jagged helical motion, JHM) at Ga ≥ 190 where the vortex shedding occurs (§ 4.4). Jagged helical motion here refers to the motion of the spheres exhibiting trajectories with visible jaggedness (or corners). As Ga increases the helical motion may give way to the ZM for both heavy (ρ_s/ρ ≥ 1.8 and Ga ≥ 210) and light (ρ_s/ρ ≤ 0.7 and Ga ≥ 200) spheres (figure 4). The relevant discussion is given in § 4.5. This is because of the confined conditions of the square tube. A transition from zigzagging to chaotic motion has been extensively observed for light spheres (Ga > 220). However, heavy spheres (ρ_s/ρ ≥ 1.8) are falling helically again after going through the zigzagging regime (Ga > 225), as seen in figure 4.

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. (2012) 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. 2004) and a confined one in a 5d circular tube (155 < Ga_crit < 160, Deloze et al. 2012). Later studies (Fabre, Tchoufag & Magnaudet 2012; Zhou & Dušek 2015) 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. 2012), 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. 2012). 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. 2004; Zhou & Dušek 2015; Raaghav et al. 2022).

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. (2012) 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. 2004; Uhlmann & Dušek 2014; Zhou & Dušek 2015; Raaghav et al. 2022). 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. 2022). 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 1999; Re = 350, Mittal 1999). 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. 2004; Horowitz & Williamson 2010; Zhou & Dušek 2015; Auguste & Magnaudet 2018) 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. 2004; Deloze et al. 2012). 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 (2010) and Zhou & Dušek (2015), 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. 2004; Auguste & Magnaudet 2018) 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. (2004).

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 (2010). 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 2004; Daniel et al. 2009; Uhlmann & Doychev 2014) 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 (2014) noticed. By contrast, for those spheres outside the cluster region the double-threaded vortices were seen in their wakes (Uhlmann & Doychev 2014).

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. (2024) 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. 2024). 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. (2024) 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 2003). 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. 2022): their drag coefficients nearly coincided with those of a fixed sphere. However, previous studies (Karamanev & Nikolov 1992; Horowitz & Williamson 2010; Auguste & Magnaudet 2018) 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. (2024) 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. 2024). 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. (2024). 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.

5. Conclusions

We simulated the freely falling or rising of a sphere confined to a 5d square tube using a 3-D lattice Boltzmann method combined with an improved bounce-back scheme. The sphere-to-fluid density ratio (ρ_s/ρ) and the Galileo number (Ga) were chosen as the control parameters with ρ_s/ρ ranging from 0.1 to 2.3 and Ga from 140 to 230. For spheres confined to a square tube, our simulations showed different behaviours compared with the cases of unconfined spheres or spheres confined to a circular tube. The conclusions are summarized as follows:

1. (i) The onset of unsteadiness occurred at Ga_crit ≈ 157 for both the heavy (ρ_s/ρ > 1) and light spheres (ρ_s/ρ < 1), representing a Hopf bifurcation. This G_crit is consistent with those of both an unconfined sphere (Jenny et al. 2004; Zhou & Dušek 2015) and a sphere confined to a 5d circular tube (Deloze et al. 2012). Near G_crit light spheres exhibited a helical path, whereas the heavy spheres oscillated near the central plane of the tube, where they were initially released from rest. Increasing Ga slightly (Ga > 160) also led to helical motion of the heavy spheres.

2. (ii) Both heavy and light spheres were widely observed to move helically within the ranges of ρ_s/ρ and Ga studied because of the double-threaded vortex in the wake of the spheres interacting with the tube walls. As Ga increases the helical motion of spheres will become irregular (Ga ≥ 190), with heavy spheres acting differently from their counterparts. The trajectories of light spheres (ρ_s/ρ ≤ 0.3) show jaggedness which is associated with the added mass effect. The double-threaded vortices begin to shed when the light spheres exhibit a JHM. In contrast, the heavy spheres alternately change their helical directions as they fall. These different behaviours were caused by the difference in inertia between the heavy and light spheres. More significantly, at a fixed Ga the oscillation period of the light spheres increases rapidly as ρ_s/ρ increases. However, the opposite was observed for the heavy spheres. This suggests a very slow circling motion of the nearly neutrally buoyant spheres. In addition, a light sphere always circles faster than a heavy one even though they have the same driving force (|ρ_s/ρ − 1|). This may be owing to the fact that a heavy sphere generally falls in an ellipse-shaped path (from top view) while a light one rises in a circle-shaped path (from top view), leading to a stronger wall effect on the hindering of the heavy sphere.

3. (iii) Interestingly, a transition from the helical motion to the ZM is seen for both heavy (ρ_s/ρ ≥ 1.8 and Ga ≥ 210) and light (ρ_s/ρ ≤ 0.7 and Ga ≥ 200) spheres. To the best of our knowledge, there are no previous reports on the ZM of a falling sphere. This certainly results from the confined situations of a square tube. A hairpin-like vortical structure (vortex loop) is formed downstream in the wake of a sphere that waveringly falls or rises at high Ga. In particular, for the ZM of light spheres, our study shows reasonable agreement with a pioneering study (Jenny et al. 2004) in terms of both the ranges of Ga and ρ_s/ρ. After passing through the zigzagging regime, the light spheres underwent chaotic motion at Ga > 220. However, the heavy spheres return to their helical regime (Ga > 225) due to the emergence of a chain of vortex loops.

4. (iv) The dependences of both the terminal Reynolds number (Re_T) and drag coefficient (C_D) on Ga show trends similar to those in the case of an unconfined sphere. It appears that the sphere-to-fluid density ratio had a small effect on both Re_T and C_D except for ρ_s/ρ = 2.3. Because of the stronger wall effect, the heaviest sphere considered in this study experienced a considerably larger drag.

The findings shed light on the near-wall motion and wake patterns of both heavy and light spheres at intermediate Reynolds numbers, and may contribute to a fundamental understanding of the transport mechanism of particulate flows confined to a square tube.