Kinematics and dynamics of freely rising spheroids at high Reynolds numbers

Abstract

We experimentally investigate the effect of geometrical anisotropy for buoyant spheroidal particles rising in a still fluid. All other parameters, such as the Galileo number (the ratio of gravitational to viscous forces) $Ga \approx 6000$, the ratio of the particle to fluid density $\varGamma \approx 0.53$ and the dimensionless moment of inertia $\boldsymbol{\mathsf{I}}^*= \boldsymbol{\mathsf{I}}_p/\boldsymbol{\mathsf{I}}_f$ (with $\boldsymbol{\mathsf{I}}_p$ being the moment of inertia of the particle and $\boldsymbol{\mathsf{I}}_f$ that of the fluid in an equivalent volume), are kept constant. The geometrical aspect ratio of the spheroids, $\chi$ , is varied systematically from $\chi = 0.2$ (oblate) to 5 (prolate). Based on tracking all degrees of particle motion, we identify six regimes characterised by distinct rise dynamics. Firstly, for $0.83 \le \chi \le 1.20$, increased rotational dynamics are observed and the particle flips over semi-regularly in a ‘tumbling’-like motion. Secondly, for oblate particles with $0.29 \le \chi \le 0.75$, planar regular ‘zig–zag’ motion is observed, where the drag coefficient is independent of $\chi$. Thirdly, for the most extreme oblate geometries ($\chi \le 0.25$), a ‘flutter’-like behaviour is found, characterised by precession of the oscillation plane and an increase in the drag coefficient. For prolate geometries, we observed two coexisting oscillation modes that contribute to complex trajectories: the first is related to oscillations of the pointing vector and the second corresponds to a motion perpendicular to the particle's symmetry axis. We identify a ‘longitudinal’ regime ($1.33 \le \chi \le 2.5$), where both modes are active and a different one, the ‘broadside’-regime ($3 \le \chi \le 4$), where only the second mode is present. Remarkably, for the most prolate particles ($\chi = 5$), we observe an entirely different ‘helical’ rise with completely unique features.

1. Introduction

A single particle rising or falling in a still fluid is an experiment of striking simplicity. Yet, the complexity arising from the interaction between the particle and fluid is astounding and has captivated a number of prominent scientists and engineers over the last centuries such as Leonardo da Vinci (see the review by Marusic & Broomhall 2020) and Isaac Newton (Newton 1999). Besides the fundamental scientific appeal, the problem is also relevant in numerous practical applications. This holds both in natural and industrial settings, where understanding and modelling of particle dynamics and wake structures can be of critical importance. One example is the use of particles in the chemical industry to mix a flow more efficiently enhancing reaction rates or heat transfer (Risso 2018). The problem is of further relevance to the sedimentation of granular naturally shaped particles in rivers and oceans (Lowe 1982; Meiburg & Kneller 2010), as well as to the precipitation of snow, hail and rain in the atmosphere (Byron et al. 2015; Gustavsson & Mehlig 2016; Gustavsson et al. 2017).

Most practical applications involve some interaction with strong background turbulence or particle-particle interactions (either directly through collisions or indirectly via the particle wakes). Overviews on these aspects of dispersed flows are provided by Toschi & Bodenschatz (2009), Balachandar & Eaton (2010) and Mathai, Lohse & Sun (2020). Voth & Soldati (2017) reviewed results related to particle geometrical anisotropy for light, neutrally buoyant and heavy particles in turbulent flows, specifically for small particles in turbulence. However, besides all these more complicated scenarios, the study of isolated particles in quiescent fluid is still very relevant (especially for buoyant particles), since the motion is often unaffected by the presence of turbulence or adjacent bodies (Magnaudet & Eames 2000; Risso 2018).

Directly motivated by practical problems, much of the previous work on the topic has dealt with particles of a higher density than that of the carrier fluid. In contrast, in the current work we will focus on light, rising, particles and on the influence of shape anisotropy on their rise behaviour. This expands the explored parameter space significantly and provides valuable insight on the role of the most important parameters governing the particle behaviour. These follow from the Newton–Euler equations applied to a submerged body, where the dimensionless control parameters are the density ratio $\varGamma \equiv {\rho _p}/{\rho _f}$, the dimensionless moment of inertia tensor $\boldsymbol{\mathsf{I}}^* \equiv {\boldsymbol{\mathsf{I}}_p} /{\boldsymbol{\mathsf{I}}_f}$, the particle Galileo number $Ga\equiv {\sqrt {|1-\varGamma | g l^3}}/{\nu }$ and indirectly also the particle geometry. Here, $\rho _{f}$ and $\rho _{p}$ are the densities of the fluid and particle, respectively, $g$ is the acceleration due to gravity, $l$ is a characteristic length scale of the geometry, $\nu$ is the kinematic viscosity of the fluid and $\boldsymbol{\mathsf{I}}_f$ and $\boldsymbol{\mathsf{I}}_p$ are the rotational moment of inertia tensors of the displaced fluid and the particle. Note that the Galileo number is used rather than the particle Reynolds number, $Re \equiv Vl/\nu$, because the relevant velocity scale, $V$, is not known a priori but is in fact an output parameter. In the definition of $Ga$ the buoyancy velocity, $V_b = \sqrt {|1-\varGamma |gl}$, is used instead. The main focus of the current work will be on the dependence on particle geometry.

The simplest case of a freely rising body is that of a spherical (isotropic) particle, which has been studied widely (see, e.g. Preukschat 1962; Murrow & Henry 1965; Karamanev & Nikolov 1992; Veldhuis, Biesheuvel & Lohse 2009; Horowitz & Williamson 2010; Auguste & Magnaudet 2018). For freely rising and falling spheres, the onset of path instabilities was characterised numerically by Jenny, Dušek & Bouchet (2004) and was found to occur between $Ga = 150$ and $Ga = 225$. It should be noted that the onset of path instabilities is only weakly dependent on $\varGamma$ and $\boldsymbol{\mathsf{I}}^*$. However, the influence of these parameters is significant for the complex dynamics and rise patterns at even higher $Ga$. This complexity arises from the laminar separation of the boundary layers past a blunt body, which leads to significant unsteadiness in the wake due to vortex shedding (Achenbach 1974). As a consequence, the variety of particle behaviours observed when varying $Ga$ or $\varGamma$ is rich, even for the simplest (i.e. isotropic) geometry (Veldhuis & Biesheuvel 2007; Horowitz & Williamson 2010).

Path oscillations of isotropic bodies are induced by horizontal asymmetries in the pressure field and are unaffected by particle orientation. Furthermore, it should be kept in mind that the rotation of the particle is coupled to translation, which can be modelled using the Kelvin–Kirchhoff equations (Ern et al. 2012; Mathai et al. 2020). Particle rotation can induce a Magnus lift force on the body causing additional pressure forcing which can even affect regime transitions (Mathai et al. 2018). Thus, in this manner, skin friction affects the horizontal oscillations indirectly. Therefore, the rotational moment of inertia of the body needs to be controlled when performing these types of experiments. This effect, combined with a potential dependence on $\varGamma$, could be one of the factors responsible for the large spread in drag coefficient data observed for spherical geometries at high $Re$ (Horowitz & Williamson 2010, figure 2), as previously suggested by (Mathai et al. 2018).

In practical situations, the particle's shape is almost never perfectly spherical. Therefore, studying the effect of anisotropic geometries is of significant relevance (Corey 1949; Dupleich 1949; Alger 1964; Dietrich 1982). Fernandes et al. (2007) showed that the onset of path instabilities is very similar for anisotropic bodies compared with spheres, albeit the critical Reynolds numbers are slightly altered. Beyond the transition point, however, the dynamics are fundamentally different since for non-spherical particles the centre of pressure, in general, does not coincide with the geometric centre of the particle. This creates an additional coupling between the fluctuating pressure in the wake and the rotational dynamics of the particle. Furthermore, the particle alignment relative to the incoming flow induces additional circulation around the body which results in a lift force contribution that is non-existent for spheres. Recent work by Sheikh et al. (2020) underlines the importance of fluid inertia for the alignment of particles settling in turbulence. To date, related studies have focused on the transitional dynamics at low $Ga$. A classification into regimes is reported for spheroids in the numerical work by Zhou, Chrust & Dušek (2017). Similarly, experimental investigations have been performed for both cylinders (Toupoint, Ern & Roig 2019) and disks (Fernandes et al. 2005, 2007, 2008; Zhong, Chen & Lee 2011; Lee et al. 2013; Auguste, Magnaudet & Fabre 2013; Zhong et al. 2013). For disks, a classification of the behaviour into regimes was first presented by Willmarth, Norman & Harvey (1964) and later extended by others (e.g. Field et al. 1997; Zhong et al. 2011).

It is the goal of this paper to elucidate how varying degrees of anisotropy affect the rise behaviour of light particles. Such an undertaking requires a precise definition and control of the ‘anisotropy’ in the particle shape – and evidently also a restriction of the infinitely many possible particle shapes. To achieve both, we investigate spheroids, ellipsoids of revolution, ranging from oblate (disks) to prolate (needles). Naturally, these shapes are well defined mathematically. Compared with, for instance, disks or cylinders, an additional benefit is that spheroids allow for a gradual transition from the isotropic sphere towards both prolate and oblate particles all the way to their extremes.

Among the specific aspects we want to address is the dependence of the drag coefficient on particle geometry. Both for spheres and anisotropic particles, it is known that the drag is dependent on $Re$. For low $Re$, the drag force is dominated by skin friction, this is called ‘Stokes drag regime’. For $Re \gtrapprox 1000$, the dominant contribution to the particle drag is from vertical asymmetry in the pressure distribution on the body's surface resulting from flow separation and the formation of a wake. Here, the overall drag becomes independent of $Re$ (Newton's drag regime). This has extensively been documented in studies on natural, settling, particles (Alger 1964; Dietrich 1982; Haider & Levenspiel 1989; Ganser 1993; Hölzer & Sommerfield 2008). In these papers the drag coefficient in Newton's regime is also shown to depend on parameters that characterise geometrical anisotropy. In Newton's regime a potential dependence of the drag coefficient on $\varGamma$ and on $\boldsymbol{\mathsf{I}}^*$ has not been investigated yet. This is one of the main aspects to which the present study aims to contribute.

Another point to consider is under what circumstances a particle can tumble, i.e. ‘flip over’. Related studies have characterised this transition in terms of $Re$ and a non-dimensionalised moment of inertia. Thus far, these works have predominantly focused on the dynamics of flat plates and strips (Smith 1971; Lugt 1983; Tanabe & Kaneko 1994; Belmonte, Eisenberg & Moses 1998; Mahadevan, Ryu & Samuel 1999; Andersen, Pesavento & Wang 2005a,b) or disks (Willmarth et al. 1964; Field et al. 1997; Zhong et al. 2011, 2013; Lee et al. 2013). Belmonte et al. (1998) identified a critical Froude number (defined as the ratio of characteristic times scale for downward motion and pendular oscillations) of $Fr \approx 0.67$, which governs the transition from flutter-to-tumble for quasi two-dimensional flat plates. Here, we will explore how the transition from flutter-to-tumble changes for non-slender geometries. Recently, Essmann et al. (2020) investigated the importance of added mass on the dynamics of ellipsoids with a set initial rotational and translational kinetic energy in both inviscid and viscous fluids.

Finally, we would like to point out that most of the work on anisotropic bodies has thus far focused on either heavy or close to neutrally buoyant particles i.e. $\varGamma \geq 1$. The value of $\varGamma$ in the present work is significantly lower. Our results can therefore shed light on the importance of the effects of density ratio and moment of inertia for anisotropic bodies. A lower particle mass will result in larger amplitude translations and rotations for the same pressure distribution around the geometry, which could result in different regimes for the same geometry and $Re$. Thus, a stronger coupling between fluid forcing and particle motion is expected. No experimental data (or otherwise) is available, especially for Newton's drag regime. Therefore, the study of light particles is of fundamental and general interest in understanding the mechanisms resulting in specific regimes at any density ratio.

2. Experimental method

2.1. Particle characteristics

All particles used in this study are spheroids. Their shape is defined by $4(X^2/h^2 + Y^2/d^2 + Z^2/d^2) = 1$, where $h$ and $d$ specify the full lengths of the major and minor axes and capital letters indicate the particle coordinate system, which is aligned with the three perpendicular planes of symmetry (see figure 1a,b). Thus, the geometry can be defined by the aspect ratio $\chi \equiv h/d$. For $\chi < 1$, the geometry is called oblate (disk shaped), $\chi = 1$ corresponds to a sphere and, for $\chi > 1$, the spheroids are prolate (needle shaped). The particle orientation can be expressed in terms of a pointing vector $\boldsymbol {\hat {p}}$, which by definition aligns with the axis of rotational symmetry (throughout this work,  $\hat {\ }$ will be used to denote unit vectors). Frequently, it is useful to consider the alignment of the particle with respect to the direction of gravity and with respect to the direction of motion of the body. These angles, which are schematically depicted in figure 1(c), will be denoted by $\theta _{\hat {g}}$ and $\theta _{\hat {v}}$, respectively. The rotation around $\boldsymbol {\hat {p}}$ is denoted by the rotation angle $\psi$.

Figure 1. Schematic of (a) prolate ($\chi >1$) and (b) oblate ($\chi < 1$) particles along with the relevant length scales and the pointing vectors. In both cases the grey shaded areas indicate the maximum cross-sectional area of the geometry. (c) The lab reference frame in which an arbitrary particle pointing vector and unit velocity vector are shown. The angles $\theta _{\hat {g}}$ and $\theta _{\hat {v}}$ represent the angle between the pointing vector and, respectively, the vertical $-\boldsymbol {\hat {g}}$ and the unit velocity vector $\boldsymbol {\hat {v}}$. (d) A picture of a selection of particles with varying aspect ratio. The particles have a pattern painted on their surface to aid with the orientation tracking.

Throughout all experiments presented here, only the aspect ratio ($\chi$) was varied in 23 steps from 0.2 to 5. All other particle characteristics, namely the volume, $V$, and the particle mass, $m_p$, were kept constant. Introducing the volume equivalent sphere diameter, $D \equiv \sqrt [3]{6V/{\rm \pi} }$, as the relevant length scale we obtain a buoyancy velocity $V_b = \sqrt {|1-\varGamma | gD}$ and the Galileo number

(2.1)\begin{equation} Ga \equiv \dfrac{\sqrt{\left|1-\varGamma\right|gD^3}}{\nu}. \end{equation}

The average value of $Ga$ is nominally constant across all particles and aspect ratios at approximately 6000 with a standard deviation of 122. The mean density ratio of the particles was kept as close to constant as possible with a mean $\langle \varGamma \rangle _{all} = 0.53$ and a standard deviation of 0.015. Here, $\langle \cdot \rangle _{all}$ is the average over all particles used in the experiments. The average volume equivalent diameter was $D = 19.90 \pm 0.3\ \textrm {mm}$. Finally, the internal structure of the particles is designed to mimic a particle produced from a uniform, homogeneous material with density $\rho _p$. This implies that the dimensionless moment of inertia tensor in the particle coordinate system in all cases is equal to $I^*_{ij} = \varGamma \delta _{ij}$, where $\delta _{ij}$ is the Kronecker delta. A full list of particle properties is provided in appendix A.

The rigid spheroidal particles were created using three-dimensional (3-D) printing on a RapidShape S30 SLA printer with a vertical resolution of $25\ \mathrm {\mu }\textrm {m}$ and a horizontal resolution of $21\ \mathrm {\mu }\textrm {m}$. A non-porous resin was used with a density of approximately $1130\ \textrm {kg}\,\textrm {m}^{-3}$. The particles shells were printed in two parts and glued together forming a watertight seal. The particles were smoothened by hand to remove any edges at the adhesive joint and layers due to the printing process. A base coat of white paint was applied. Next, the particle was masked and black paint was applied producing the patterns (see figure 1d) which facilitates the orientation tracking. Two types of patterns were used: a complex one identical to that employed by Mathai et al. (2016) for $0.83 \le \chi \le 1.20$ and a simpler one (as shown in figure 1d) for the others. The mass of the paint was found to be negligible. The primary dimensions of each particle and its mass were measured to obtain accurate values for the control parameters. No effect of surface roughness was found when testing different painting methods or unpainted particles. Surface roughness is often an important parameter when dealing with transitioning boundary layers, however, in the present work we suspect that the effect of roughness is limited as the Reynolds numbers encountered are too low to trigger a transition to turbulent boundary layers, which, for a fixed sphere, occur at substantially higher values of $Re \approx 3\times 10^5$ (Achenbach 1972).

2.2. Set-up and measurement procedure

The experiments were performed in the approximately 3 m high test section of the Twente water tunnel (TWT) depicted in figure 2. The temperature in the laboratory was kept constant at $20\,^{\circ }\textrm {C}$ and the water had ample time to equilibrate, therefore, the density and kinematic viscosity of the water are assumed to be constant at $998\ \textrm {kg}\,\textrm {m}^{-3}$ and $1.00\times 10^{-6}\ \textrm {m}^2\ \textrm {s}^{-1}$, respectively. The particles were released using a release mechanism 1.8 m ($90D$ for the particles used here) below the measurement domain in order to obtain a statistically steady state, as verified by means of the vertical acceleration statistics. The release mechanism consists of a water lock that allows the particle to be inserted without draining the tank. The particle is pushed to the release position inside a ‘basket’. The end of the pipe in which the basket rests is open at the top, such that the particle can be released by rotating the basket slowly. Between subsequent experiments the fluid was given 10 min to settle. No notable changes in particle behaviour were observed when extending the waiting period for up to an hour. Furthermore, we excluded any runs that were contaminated by bubbles on the surface of or around the particle. The weight of the particles was checked after every experiment in order to make sure no water had leaked into the particle shell.

Figure 2. (a) Schematic representation of the experimental set-up showing side and top views of the measurement section as well as the camera positions. (b) Three-dimensional graphic showing the full set-up including the release mechanism at the bottom of the TWT.

During the experiments, we tracked the position and orientation of the particles as they rose through the quiescent fluid. The lab coordinate frame is defined as shown in figure 2. The origin is located at the base of the measurement volume, 1.8 m above the release. The $z$-direction is defined upwards in the vertical direction. The translational and rotational tracking was performed by image analysis of recordings made using two stationary Photron AX200 cameras with $1024\times 1024$ resolution at 256 grey levels and a recording rate of 250 fps. This recording rate was found to be sufficient to resolve the particle dynamics. The two cameras were positioned perpendicular to one another, as shown in the top view in figure 2, to track the 3-D motion of the particles. The cameras were aligned with the axes of the global coordinate system along the $x$- and $y$-axes, perpendicular to the glass walls of the set-up to minimize optical distortion. In order to obtain longer particle tracks a second set of cameras with the same specifications was positioned above the first. The two sets were arranged such that their fields of view slightly overlapped in the vertical direction. The ‘overlap region’, where the particle could be seen by all four cameras simultaneously, was approximately 60 mm in height. All four cameras were positioned almost three metres away from the centre of the tunnel and outfitted with objectives with a 100 mm focal length. Each camera had approximately a $580\ \textrm {mm}\times 580\ \textrm {mm}$ field of view in the central plane of the tunnel resulting in a spatial resolution of around 0.560 mm per pixel. The difference in magnification between the individual cameras was less than 0.0025 mm per pixel. The particle tracks obtained using this set-up extend approximately 1 m in the vertical direction.

The tank was illuminated by eight continuous LED light sources. The lights were positioned such that the illumination is as homogeneous as possible while avoiding shadows in the particle images which would render the tracking inaccurate. As a backdrop, two grey PVC plates were used. The camera apertures are adjusted such that the background corresponds to a grey value of approximately 128 out of the 256 to maximise the contrast with the black and white patterns on the particles. We further ensured that the depth of focus was sufficiently large to cover the full test section.

The tracking algorithm is based on the work of Mathai et al. (2016) with only moderate modifications. Therefore, a complete description of the data processing method used to obtain the particle position and orientation is relegated to appendix B.

2.3. Data processing

The position of the centre of mass was extracted from the recorded frames. Consequently, the position data was smoothed in time using a convolution with a Gaussian kernel. Similarly, derivatives of the Gaussian kernel were used to obtain the first and second derivatives with respect to time (Mordant, Crawford & Bodenschatz 2004). For the position, a kernel with a standard deviation of 4 and window size of 13 frames was used. These values were increased to $(6,21)$ for the velocity and $(12,33)$ for the acceleration. These values were determined following the approach proposed by Mathai et al. (2016) to filter out the high frequency noise but not the much lower frequency particle position, velocity and accelerations. This was possible thanks to temporal oversampling such that the filtering leaves the particle dynamics and its statistics unaffected. Additionally, the orientation was extracted in terms of three Euler angles $\phi$, $\theta$ and $\psi$. For the remainder of this work, we will express the orientation in terms of the pointing vector $\boldsymbol {\hat {p}}$, i.e. the direction of the particles axis of rotational symmetry (see figure 1), and a rotation around $\boldsymbol {\hat {p}}$, which we call $\psi$. The orientation data was smoothed using the same parameters as the position data.

Further processing of the trajectories is required in order to extract properties such as the frequency and amplitude of oscillation. The simplest methods, such as fitting periodic functions to the data as employed in Fernandes et al. (2007), are not suited for the more chaotic trajectories encountered at higher $Ga$ (see figure 8 and supplementary movies available at https://doi.org/10.1017/jfm.2020.1104). It was therefore necessary to develop an alternate approach in order to determine properties consistently across the full range of aspect ratios considered. This will be described based on a sample trajectory as shown in figure 3(a–c) in the following.

Figure 3. Panels (a–c) show the same particle trajectory, for $\chi = 0.25$ ($Ga \approx 6263$ and $\varGamma \approx 0.53$), in various states of processing. Panel (a) shows the raw data in blue, and the period-averaged trajectory in red. Points of maximum and minimum amplitude are also shown (small black symbols). In (b,c) we show the drift corrected trajectory (dashed line), and the precession corrected path (solid line) from, respectively, the side and top. The schematic in (d) illustrates the precession correction based on the points of maximum amplitude, the symbols match those in (a). The points $\tilde {x}_n$ for the precession corrected trajectory are identical to the maxima in amplitude of the position vector in the horizontal plane.

As a starting point, the smoothed trajectories (blue line in figure 3a) and its derivatives were used. First, the autocorrelation functions of the horizontal particle velocity components in the lab coordinate frame were computed. By determining the interval between subsequent peaks of this function, a first estimate for the particle oscillation period was obtained. This procedure, however, does not yield the correct frequency in cases where the trajectory is also precessing, as illustrated in figure 3(d).

To overcome this issue, we defined a period-averaged trajectory (shown in red in figure 3a), where the window size of the moving average was based on the frequency obtained previously from the velocity autocorrelation. Subtracting the phase-averaged trajectory from the original one yields a drift corrected trajectory (dashed line in figure 3b,c). We defined new coordinates $\tilde {x},\tilde {y}$, for which the origin coincides with the position of the phase-averaged trajectory. Next, we searched for local maxima in the distance between the phase-averaged and the original trajectory (i.e. the norm of the position vector in the $\tilde {x},\tilde {y}$-frame). The magnitude of the distance at these points corresponds to the maximum amplitudes $a_n$, where $n$ is an index (see black lines in figure 3a). Based on the locations of the maxima, the trajectory was corrected for precession as shown by the solid lines in figure 3(b,c). This was done by point-wise rotating segments of the trajectory around the origin (corresponding to the period-averaged trajectory), such that the points of maximum amplitude lie on the $\tilde {x}$-axis, i.e. $\tilde {x}_n = a_n$, as schematically shown in figure 3(d). In this subfigure the designated half-period is rotated by $\zeta$/2 to end up across from the previous point of maximum amplitude, here $\zeta$ is the precession angle over a complete particle oscillation cycle.

Finally, the precession corrected velocity data in the $\tilde {x}$-direction was used to obtain a new estimate of the oscillation frequency and the mean was taken over all runs. This whole process is then repeated using this new frequency as the smoothing time scale for the phase averaging until the results converge (typically within 3–4 iterations). Statistics of the frequencies determined in this way form the basis for the discussion in the subsequent § 4.1 and appendix E.

2.4. Data set

In total the data set consists of 269 runs. For each run, the particle rises through the approximately 1 m high measurement domain. A minimum of nine runs have been obtained per aspect ratio. For each aspect ratio, a minimum of two particles were produced in order to rule out potential flaws in the production process. No aberrant behaviour was found that correlated to a single particle. Generally, the experiments show repeatable behaviour, however, in some cases there do appear to be multiple states that the particle motion can be in; transitioning from one to the other and then back akin to nonlinear dynamical systems. This appears to be a property of the particle dynamics and not a flaw of the experiments.

The temporal coupling between particle motion and orientation is more clearly evident in the 10 supplementary movies provided. These movies also provide a more intuitive understanding of the dynamic behaviour of the spheroids and we will refer back to them throughout this work.

3. Vertical motion: rise velocity, Reynolds number and drag coefficient

In this section we will focus on the particle's vertical motion, most importantly, their rise velocities. We will show that the emergent behaviour of the particles is strongly dependent on $\chi$. We classify this dependence into six distinct regimes of the particle motion. Their definitions, characteristics and crossovers will be further elaborated and supported as we study various properties of the rise patterns in more detail throughout this work.

3.1. Reynolds number

Following the work by Fernandes et al. (2007) on rising disks, we define the particle Reynolds number

(3.1)\begin{equation} Re = \dfrac{v_z d_A}{\nu}. \end{equation}

Here, $v_z$ is the vertical velocity and the length scale $d_A$ denotes the area equivalent disk diameter, which is based on the maximum cross-sectional area of the geometry ($A^+ (\chi )$) such that $d_A = \sqrt {4A^+ /{\rm \pi} }$ and, consequently,

(3.2)\begin{equation} d_A = \begin{cases} d, & \text{for} \ \chi < 1,\\ D, & \text{for} \ \chi = 1,\\ \sqrt{dh}, \quad & \text{for} \ \chi > 1.\end{cases} \end{equation}

Assuming a constant $\varGamma$ and a balance of buoyancy ($\sim d^2h$) and pressure drag ($\sim v_zd_A^2$) forces, $Re(\chi ) = \text {const.}$ implies the scalings $v_z\sim \sqrt {h}$ for oblate and $v_z\sim \sqrt {d}$ for prolate geometries, respectively. Note that $A^+$ corresponds to the grey shaded areas shown in figure 1(a,b). The choice of the length scale $d_A$ is motivated by the facts that, for blunt bodies, pressure drag dominates and that, according to potential flow theory, spheroids orient with the largest area perpendicular to the incoming flow (see Lamb 1932, article 124).

In figure 4 we show $\langle Re \rangle _n$ (black symbols), where the ensemble average $\langle \cdot \rangle _n$ denotes averaging over all data points obtained for a specific aspect ratio, i.e. an average over time and across runs. To complement the results for the average, we also consider two types of fluctuations in $Re$. The first kind relates to using the instantaneous value of $v_z$, as presented in (3.1). The corresponding results are indicative of the overall fluctuations in $Re$ and they are represented by the grey-shaded regions in figure 4. The different shadings relate to the quantiles of the distribution of data that lie within a certain region, e.g. for the 90 % region, 5 % of the instantaneous data had a higher and 5 % a lower $Re$ value than the bounds of this region. This method of visualizing the result will be used throughout this paper. Note that fluctuations in $Re$ are largely due to velocity fluctuations during an oscillation cycle. Additionally, we quantify how much the mean $Re$ varies over individual oscillations cycles of the particles. To this end, the error bars in the figure show the standard deviation of $\langle Re \rangle _p$, where $\langle \cdot \rangle _p$ indicates a moving averaged $Re$ over one oscillation period.

Figure 4. Reynolds number as a function of aspect ratio for $Ga \approx 6000$. The grey symbols show $\langle Re_D \rangle _n$ and black symbols indicate the alternative definition $\langle Re \rangle _n$, using the length scale based on the maximum cross-flow area. Note that $Re_D$ is a scalar multiple of the mean rise velocity, which is indicated on the right-hand side axis. The error bars represent the standard deviation of the phase-averaged fluctuations in $Re$. The grey shaded areas show the fraction of instantaneous data points that fall within these respective regions. For $\chi = 4$ and $5$, two values are shown; these correspond to the two modes occurring simultaneously as described in § 6.2, the data points marked with crosses mark the helical regime. The coloured regions indicate the different regimes that are defined based on the analysis of the particle kinematics.

The variation of $Re$ across different aspect ratios evident from figure 4 unveils the signature of distinct regimes of particle motion, which will be introduced in the following. These regimes are indicated by colours in the background of the figure.

Particles close to isotropic (0.83 $\leq \chi \leq$ 1.2) are observed to be in the ‘tumbling’ regime (green). Here, $\langle Re \rangle _n$ as a function of $\chi$ attains a global maximum for $\chi = 1$. On both, the oblate and the prolate sides, $\langle Re \rangle _n$ drops significantly once the particle becomes anisotropic. The characteristic feature of the tumbling regime is that the pointing vector can have any orientation and the particles flip over as will be shown in §§ 5.1 and 6.1. For even more oblate particles (lower $\chi$), we encounter the ‘zigzag’ regime (blue) for 0.29 $\leq \chi \leq$ 0.75. In this regime the motion of the particles is spiraling with varying degrees of eccentricity, ranging from nearly circular to almost planar orbits (see figure 20b and appendix E.2). The motion is very regular, as evidenced by the small error bars denoting the variation in the period-averaged Reynolds number. The most important property of the zigzag regime is the near constant value of $\langle Re \rangle _n$, indicating that $v_z \sim \sqrt {h}$. The most extreme oblate particles investigated here ($\chi \leq 0.25$) again display a different behaviour, which we refer to as the ‘flutter’ regime (purple). Particle motion in this regime resembles inverted falling leaves (Pesavento & Wang 2004). Horizontal motions become dominant and the vertical velocity varies more strongly throughout an oscillation cycle. The latter is evident from the wide spread in the instantaneous values of $Re$ while the period-to-period variation, $\langle Re \rangle _p$, remains small. As a consequence, $\langle Re \rangle _n$ decreases strongly with decreasing $\chi$ within the flutter regime.

The rise patterns of prolate particles with aspect ratios between approximately $1.25 \leq \chi \leq 2.5$ are characterised by strong oscillations of the pointing vector (the ‘long’ direction of these particles). However, unlike the tumbling prolate particles, spheroids in this regime do not ‘flip over’. Therefore, we call this the ‘longitudinal’ regime (yellow). Trajectories in this regime are more chaotic when compared with their oblate counterparts in the zigzag regime. This is related to the fact that with the symmetry axis oriented perpendicular to the velocity vector on average, there are two competing cross-flow length scales in the longitudinal regime as will be discussed in detail in § 4.2. In this regime $\langle Re \rangle _n$ also appears constant, thus, $v_z \sim \sqrt {d}$. At more extreme prolate aspect ratios ranging from $2.5 < \chi \leq 4.5$, we encounter the ‘broadside’ regime (orange). Here, the oscillations of the pointing vector almost completely vanish. The dynamics appear consistent with the forcing by vortex shedding on an infinite fixed cylinder. The resulting horizontal translation of the particle is almost exclusively in the broadside direction of the particle, giving the regime its name. Remarkably, the value of $\langle Re \rangle _n$ in this regime is consistently higher than in the longitudinal regime, even though $v_z \sim \sqrt {d}$ continues to hold. Finally, we observed another transition to what we call the ‘helical’ regime (red) for the most extreme prolate aspect ratio. Here, we find two rise patterns that coexist, the first being similar to the longitudinal mode indicated in the figure by the circular markers and the second, a near perfect helix represented by the crosses. The helical rise pattern was also found to exist for the particle of $\chi = 4$ and is in this case dependent on the release conditions. Not surprisingly, the values of $\langle Re \rangle$ differ significantly between the two modes. We will revisit this regime in § 6.2. It should also be noted that the onset of the helical pattern appears to be discontinuous, while all other regime transitions (as a function of $\chi$) are gradual.

It is instructive to also consider an alternate definition of the particle Reynolds number based on $D$,

(3.3)\begin{equation} \langle Re_D \rangle_n = \dfrac{\langle v_z \rangle_n D}{\nu}. \end{equation}

The resulting Reynolds numbers are also included in figure 4. Note that $\langle Re_D \rangle _n$ only varies due to changes in $\langle v_z \rangle _n$, the mean rise velocity (also provided for these data points in the same figure). Two interesting observations can be made by comparing results for $Re$ and $Re_D$. First, the rise velocity continuously decreases with increasing anisotropy everywhere but in the broadside regime where a secondary local maximum, i.e. a peak in the mean rise velocity, is observed. Secondly, it becomes obvious, that the approximately constant values of $Re$ in the zigzag, longitudinal and broadside regimes critically rely on rescaling with the appropriate cross-flow length scale.

3.2. Drag coefficient

A key finding from figure 4 is that neither $Re_D$ nor $Re$ are uniquely determined by $Ga$ but vary significantly with $\chi$. It is interesting to compare this result to the work by Fernandes et al. (2007), as is done in figure 5. This figure clearly shows that there is no discernible $\chi$ dependence for flat disks at lower $Re$ ($80 < Re < 330$) for $0.1 \leq \chi \leq 0.5$, while the spread with $\chi$ is significant in the data at large $Re$. This indicates that the effect of aspect ratio variations is dependent on $Ga$ or $Re$. Nevertheless, also the difference in geometry (disks versus spheroids) should be kept in mind when interpreting these results.

Figure 5. Particle Reynolds number versus Galileo number. The black markers show the data by Fernandes et al. (2007). The current results are indicated by the star shaped markers where the colour indicates the regimes as defined in figure 4. The inset shows a more detailed overview of the current data set along with lines of constant drag coefficient.

The ratio of $Ga$ to $Re$, i.e. the slope in figure 5 is related to the commonly used drag coefficient $C_d$. To demonstrate this, we assume that the particles have reached their terminal velocity (as is the case for all results considered in this study). This implies that the sum of time-averaged forces applied to the body in the vertical direction is equal to zero: $\langle \sum F_z \rangle _t = 0$. Here, the net driving force $\boldsymbol {F}_B = (\rho _f-\rho _p) V g \boldsymbol {\hat {z}}$ is balanced by the sum of all fluid forces. The fluid forces can be decomposed into two perpendicular components; the drag ($\boldsymbol {F}_d$) acting parallel to the direction of motion, and lift perpendicular to it. A convenient definition for rising or settling particles defines drag as the force balancing the net buoyancy force, such that the force balance reads as

(3.4)\begin{equation} \left\langle \sum F_z \right\rangle_t = \underbrace{(1-\varGamma) \rho_f \tfrac{1}{6}{\rm \pi} D^3 g}_{\boldsymbol{F}_B} - \underbrace{\tfrac{1}{8}{\rm \pi}\rho_f \langle v_z \rangle^2_t C_d d_A^2}_{\boldsymbol{F}_d} = 0. \end{equation}

Note that similar to (3.1) we use $d_A^2 \propto A^+(\chi )$ as the cross-flow area in the definition of $C_d$ in (3.4). Additionally, it is worth mentioning that adhering to common practice we use $\langle v_z \rangle ^2_n$ in the definition of $C_d$, instead of $\langle v_z ^2 \rangle _t$. While the latter appears more appropriate conceptually, it cannot simply be deduced from the terminal velocity and is therefore not readily available in practice. From (3.4), it follows that $C_d$ relates $Re$ and $Ga$ according to

(3.5)\begin{equation} C_d = \dfrac{4(1-\varGamma)D^3g}{3 \langle v_z \rangle^2_t d_A^2} = \dfrac{4}{3} \dfrac{Ga^2}{\langle Re \rangle_n^2}. \end{equation}

We include lines of constant $C_d$ in figure 5. Fernandes et al. (2007) noted that for their range of $Ga$, $C_d$ was approximately constant at 1.2 for $0.1 \le \chi \le 0.5$. In contrast, we find that the particle drag coefficient varies from close to 0.59 for a sphere to around 2.08 for the most oblate disk. In the following sections we will investigate the behaviour of the drag coefficient in further detail, looking for physical mechanisms explaining these observations.

The inset of figure 5 shows that $C_d \approx 1.2$ is also encountered for the present data in the zigzag regime (blue). And indeed also other features of the data in Fernandes et al. (2005, 2007) appear consistent with this regime. This suggests that, for higher $Ga$, regime transitions occur for lower levels of anisotropy. This hypothesis will be further explored in § 5.2.

Specifically for the sphere, for which there is plenty of reference data in the literature, we observe $C_d = 0.59$. This result lies in between the values reported in Horowitz & Williamson (2010) for what they call the ‘vertical’ ($C_d \approx 0.45$) and ‘zigzag’ ($C_d \approx 0.75$) regimes. However, our sphere exhibits oscillations similar to the ‘zigzag’ (appendix E.1) and the drag coefficient is very similar to the results by Preukschat (1962) ($C_d \approx 0.54$) for a very similar density ratio ($\varGamma \approx 0.6$). Both these experimental results are contrary to the idea by Horowitz & Williamson (2010) that low drag is associated with the ‘vertical’ and high drag for the ‘zigzag’ regime. Similar observations regarding the more subtle effects of $\varGamma$ and the importance of path oscillations on the drag of a sphere were made in the numerical work by Auguste & Magnaudet (2018).

As a side note, we mention that the inset in figure 5 shows the spread in $Ga$ for the particles used in the present experiments. Additionally, note that the variations in $Ga$ for different particles at the same $\chi$ are small (typically less than 2 %) (appendix A, table 1), therefore, the results are consolidated to one nominal $Ga$.

Table 1. Average particle properties per aspect ratio group $\langle \cdot \rangle _\chi$, we provide the measured particle aspect ratio $\chi '$, density ratio and the Galileo number along with the standard deviation per aspect ratio grouping $\sigma _\chi (Ga)$. The bottom two rows show the average and the standard deviation over all particles, i.e. not grouped by aspect ratio.

In figure 6 we plot $C_d$, as defined in (3.5), explicitly as a function of $\chi$. Note that according to (3.5) and with $Ga$ being constant, $C_d$ is directly related to $Re$. Hence, the regions of approximately constant $Re$ in figure 4 show up as near constant $C_d$ in figure 6 and also other aspects of the discussion around figure 4 apply. It is nevertheless useful to consider the drag coefficient separately. Firstly, because it is an important and widely used quantity in applications. Moreover, adopting the definition of $C_d$ has the added benefit of accounting for the slight variation in $Ga$ between different particles. But most importantly, studying the drag coefficient allows for direct comparison to related results in the literature and provides the opportunity to elucidate the cause for variations in the drag as a function of $\chi$.

Figure 6. Particle drag coefficient as a function of $\chi$ using two different definitions; $C_d$ (black circles) and $\langle C_d^* \rangle _n$ (blue circles) with the grey shaded area indicating the range of fluctuations in this quantity. The inset (blue circles) shows the mean projected area perpendicular to the direction of the instantaneous particle velocity over $A_{\chi = 1} = 1/4{\rm \pi} D^2$. The black lines indicate the minimum and maximum cross-flow area that is possible for a given aspect ratio.

3.3. Instantaneous drag coefficient

The first factor we consider that will affect the behaviour of $C_d$ with $\chi$ is the relevant cross-sectional area, which we expect to be close to $A^+(\chi )$ as was already argued in the context of $Re$ in § 3.1. We check this assumption in the inset of figure 6 by showing the mean of $A_{\perp \hat {v}}$, which denotes the area that is perpendicular to $\boldsymbol {\hat {v}}$, the instantaneous direction of motion and its probability distribution. In this figure we also show the maximum ($A^+$) and minimum ($A^-$) cross-sectional area as a function of $\chi$. All results are normalized by the cross-section of the sphere $A_{\chi = 1}$. It becomes apparent that, for almost all $\chi$, the $A_{\perp \hat {v}}$ follows $A^+$ very closely, with the notable exceptions of the flutter (purple), the helical (red) and the tumbling (green) regimes. For the latter however, the effect of this change is very small since the difference between $A^+$ and $A^-$ is small at $\chi \sim 1$. In general, these results confirm that even for very complex trajectories, $A^+$ is an appropriate choice for the cross-flow area in almost all cases.

One may further suspect that changes in the drag coefficient ($C_d$), which is defined using the vertical velocity only, are related to variations in the rise patterns. The drag might then be fairly uniform across $\chi$ when the velocity is considered along the path instead. To test this hypothesis, we define an instantaneous drag coefficient $C_d^*$, which is based on the instantaneous force balance of the particle along its trajectory

(3.6)\begin{equation} m_p\dfrac{\text{d}\boldsymbol{v}}{\text{d}t}\boldsymbol{\cdot}\boldsymbol{\hat{v}} = \boldsymbol{F}_B \boldsymbol{\cdot}\boldsymbol{\hat{v}} - \dfrac{1}{2}\rho_f\boldsymbol{v}^2A_{\perp \hat{v}}C_d^*. \end{equation}

Note that we retain the acceleration term, which only vanishes in the mean (when there is no horizontal drift) and we use the instantaneous cross-sectional area $A_{\perp \hat {v}}$ for consistency. Rewriting (3.6) to obtain an explicit expression for $C_d ^*$ yields

(3.7)\begin{equation} C_d^* = \dfrac{{\rm \pi} D^3\left((1-\varGamma)\boldsymbol{g} - \varGamma \dfrac{\text{d}\boldsymbol{v}}{\text{d}t}\right)\boldsymbol{\cdot}\boldsymbol{\hat{v}}} {3\boldsymbol{v}^2A_{\perp\hat{v}}}. \end{equation}

In this definition, $C_d^*$ accounts for all fluid forces including added mass in the direction of instantaneous motion. The ensemble-average $\langle C_d ^* \rangle _n$ is calculated for all aspect ratios and included in figure 6 (blue circles). We observe that in the ‘flutter’ regime using $C_d ^*$ indeed eliminates the strong increase in the mean drag coefficient that was observed for $C_d$ (black symbols) and we find that $\langle C_d ^* \rangle _n \approx 1.2$, which agrees with the neighbouring zigzag regime. This suggests that, when correcting for the reduced cross-flow area and increased path oscillations, the behaviour is similar between these two regimes. Note, however, that instantaneous values of $C_d^*$ fluctuate significantly in the ‘flutter’ regime as indicated by the grey shaded regions. Furthermore, we find that the instantaneous drag coefficient, $C_d^*$, is approximately single valued in the alignment $\boldsymbol {\hat {p}}(t)\boldsymbol {\cdot } \boldsymbol {\hat {v}}(t)$ (not shown), suggesting that the balance in (3.6) captures the essential dynamics. Besides in the flutter (purple) regime, we do not observe a significant difference between $C_d^*$ and $C_d$. This implies that differences between the path velocity and $\langle v_z \rangle$ are effectively compensated by changes in the alignment and the changes in the relevant cross-section associated with the latter. This observation indicates that using the maximum cross-flow area and the vertical velocity in fact are good approximations for spheroidal geometries, which is convenient as both quantities are easily obtained. We also conclude that, for all but the most oblate particles, variations in drag are not simply caused by the geometry of the trajectories but must have different origins, e.g. the structure of the wake.

On the basis of yet another definition of the drag coefficient (using $D$ as the length scale) a comparison of our results can be made to empirical drag models developed for heavy particles. This is presented in appendix D, where it can be seen that significant differences arise between the two, highlighting the importance of the particle density ratio $\varGamma$ as a control parameter.

3.4. Contribution of the added mass force on particle dynamics

Using our data, we can also single out the effect the added mass force (Maxey & Riley 1983; Magnaudet 1997; Limacher, Morton & Wood 2018) has on the particle motion. To this end, we consider the classical Kelvin–Kirchhoff equation for the translation of a rigid body in a fluid, i.e.

(3.8)\begin{equation} (m_p \boldsymbol{\mathsf{I}} + \boldsymbol{\mathsf{A}}) \boldsymbol{\cdot} \dfrac{\text{d} \boldsymbol{v}'}{\text{d} t} + \boldsymbol{\varOmega}' \times (m_p \boldsymbol{\mathsf{I}} + \boldsymbol{\mathsf{A}}) \boldsymbol{\cdot} \boldsymbol{v}' = \boldsymbol{F}'_\omega + \boldsymbol{F}'_B. \end{equation}

Here, $\boldsymbol{\mathsf{I}}$ is the identity tensor and $\boldsymbol{\mathsf{A}}$ is the added mass tensor. Furthermore, (3.8) is defined in a coordinate frame (indicated with superscript $'$) with the origin fixed in space but whose axes are parallel to the particle coordinate frame and rotate with the particle rotation rate $\boldsymbol {\varOmega }$. This choice ensures that the components of the added mass tensor $\boldsymbol{\mathsf{A}}$, which can be obtained analytically from potential flow solutions as derived in the work by Lamb (1932) (see appendix C for details), remain constant in time. The transformation from the lab coordinate frame to the instantaneous particle orientation is given by the rotation matrix $\boldsymbol{\mathsf{R}}(t)$. On the right-hand side of (3.8), the term $\boldsymbol {F}'_\omega$ represents all fluid forces on the particle (besides the added mass force) and can be decomposed into drag (acting along $\boldsymbol {\hat {v}}$) and lift (perpendicular to it): $\boldsymbol {F}'_\omega = \boldsymbol {F}'_d + \boldsymbol {F}'_l$. For completeness, results for the lift component are included in terms of a lift coefficient $\langle C_l^* \rangle _n = \langle 2\|\boldsymbol {F}_l\|/\rho _f \boldsymbol {v}^2 A_{\perp \hat {v}} \rangle _n$ in table 2 in appendix A.

Table 2. A tabulated version of the results presented throughout this work for different properties as a function of the aspect ratio. The two red highlighted lines at the bottom show the results for the helical rise rise pattern for the particles of aspect ratios $\chi =$ 4 and 5.

Typically, added mass contributions are lumped in with lift and drag, but here we aim to specifically consider their contributions to the particle dynamics. Remarkably, we find that the net contribution of the added mass terms cancels when taken along the path for all $\chi$, i.e.

(3.9)\begin{equation} \biggl\langle \Big(\underbrace{\boldsymbol{\mathsf{A}} \boldsymbol{\cdot} \dfrac{\text{d}\boldsymbol{v}'}{\text{d}t}}_{\boldsymbol{F}_{a1}'} + \underbrace{\boldsymbol{\varOmega}' \times [ \boldsymbol{\mathsf{A}} \boldsymbol{\cdot v}']}_{\boldsymbol{F}_{a2}'} \Big) \boldsymbol{\cdot} \boldsymbol{\hat{v}}' \biggl\rangle_n \approx 0. \end{equation}

This implies that, on average, the added mass does not contribute to the motion along the path under steady state conditions. Noting that also $( \boldsymbol {\varOmega }' \times [ m_p \boldsymbol {v}'] )\boldsymbol {\cdot } \boldsymbol {\hat {v}}' \equiv 0$, we thus find that the drag coefficient along the path, $\langle C_d^* \rangle _n$ as defined in (3.7), is in fact largely unaffected by added mass and inertial effects. Thus, $\langle C_d^* \rangle _n \approx \langle 2\|\boldsymbol {F}_d\|/\rho _f \boldsymbol {v}^2 A_{\perp \hat {v}} \rangle _n$.

Finally, we evaluate if there is a net contribution of the added mass terms to the vertical motion of the particle. In doing so, we consider a decomposition of the added mass force $\boldsymbol {F}_a = \boldsymbol {F}_{a1}+ \boldsymbol {F}_{a2}$, with the components defined in (3.9). In figure 7(a) we show the mean contribution of these terms in the vertical direction ($\hat {\boldsymbol {z}}$) normalized by the magnitude of the net driving force $\|\boldsymbol {F}_B\|$ as a function of $\chi$. For $\boldsymbol {F}_a$, we observe that for almost all spheroids, besides the most oblate ones, the net contribution remains small, accounting for less than $10\,\%$ of the net gravitational force. For the flutter regime, however, the contribution is much greater, reaching a value as high as 0.65 for $\chi = 0.2$. This effect is almost entirely due to $\boldsymbol {F}_{a1}$ as can be seen from the blue line in figure 7(a). Furthermore, there is a small but noticeable net vertical contribution of $\boldsymbol {F}_{a2}$ in the tumbling (green) regime, while this term remains insignificant for more anisotropic particles. For completeness, we also show $\boldsymbol {F}_{m2} \boldsymbol {\cdot }\boldsymbol {\hat {z}}/ \|\boldsymbol {F}_B\|$, where $\boldsymbol {F}_{m2} = \boldsymbol{\mathsf{R}}^T \boldsymbol {\cdot } ( \boldsymbol {\varOmega }' \times m_p \boldsymbol {v}' )$ in figure 7(a). This quantity is seen to be almost identical to $\boldsymbol {F}_{a2}$ for $\chi \approx 1$. This result is consistent with the facts that the added mass tensor for these particles is almost a scalar multiple of the identity tensor and that the translational added mass is approximately $0.5\rho _f$, which is close to $\rho _p$ used in our experiments. However, towards the fluttering regime we notice that the contribution of $\boldsymbol {F}_{m2}$ outgrows the added mass effect. This implies that the particle rotation plays a noticeable role in the vertical force balance of tumbling and fluttering particles.

Figure 7. (a) Ensemble-averaged added mass contributions as a function of $\chi$. (b,c) Scatter plots showing the instantaneous added mass force ($\boldsymbol {F}_{a1} + \boldsymbol {F}_{a2}$) along the $z$-direction and alignment between the direction of particle acceleration and the $z$-direction. In (b) four oblate and in (c) four prolate aspect ratios are shown.

In figure 7(b,c) we examine the combined vertical added mass contribution more closely to explain the asymmetry that causes a net vertical force. To this end, we present scatter plots of the normalized vertical component of the added mass force, $\boldsymbol {F}_a \boldsymbol {\cdot }\boldsymbol {\hat {z}}/ \|\boldsymbol {F}_B\|$, versus $\boldsymbol {\hat {a}} \boldsymbol {\cdot } \boldsymbol {\hat {z}}$, i.e. the alignment between the particle acceleration vector and the vertical.

Initially focusing on the results for $\chi = 0.75$ in figure 7(b), we observe that negative and positive excursions of the added mass force are more or less symmetrical, i.e. when the particle accelerates downwards $\boldsymbol {\hat {a}} \boldsymbol {\cdot } \boldsymbol {\hat {z}} < 0$, the added mass force is positive and when $\boldsymbol {\hat {a}} \boldsymbol {\cdot } \boldsymbol {\hat {z}} > 0$, it is negative and of comparable magnitude. As expected, this symmetry results in a zero net contribution of the added mass to the vertical added mass force as is shown in figure 7(a). The same applies to the results for prolate particles in figure 7(c).

Now we consider the remaining oblate particles, with more extreme anisotropy, shown in figure 7(b). For increasing anisotropy, the symmetry in $\boldsymbol {F}_a$ completely vanishes. Remarkably, even when $\boldsymbol {\hat {a}} \boldsymbol {\cdot } \boldsymbol {\hat {z}} < 0$, the direction of the added mass force is in the positive $z$-direction for a large portion of the data points. There are two effects at play that cause this asymmetry. First, an asymmetric distribution is produced if the particle orientation (relative to $\boldsymbol {\hat {a}} \boldsymbol {\cdot } \boldsymbol {\hat {z}}$) is different when the particle is accelerating versus decelerating. Second, it is important to note that $(\boldsymbol{\mathsf{R}}^{T}\boldsymbol {\cdot } (\boldsymbol{\mathsf{A}} \boldsymbol {\cdot } \boldsymbol{\mathsf{R}})) \boldsymbol {\cdot } \text {d}\boldsymbol {v}/\text {d}t$ is not parallel to $\text {d}\boldsymbol {v}/\text {d}t$. Combined with the fact that the particle acceleration in the horizontal direction becomes more and more dominant for increasing anisotropy, this can result in an upward added mass force even when the particle decelerates in the vertical direction. Both these effects can be observed in the supplementary movies showing the particle instantaneous accelerations and alignment along the path, for instance, when comparing the video for $\chi = 0.57$ to that for $\chi =0.25$.

The correlation between $\boldsymbol {\hat {a}} \boldsymbol {\cdot } \boldsymbol {\hat {z}}$ and $\boldsymbol {F}_{a} \boldsymbol {\cdot } \boldsymbol {\hat {z}}/\|\boldsymbol {F}_{B} \|$ translates into the net force in the vertical direction shown in figure 7(a). For the most oblate particle, $\langle \boldsymbol {F}_{a} \boldsymbol {\cdot }\boldsymbol {\hat {z}} \rangle _n$ amounts to more than half of the buoyancy force, but the effect quickly decreases with increasing $\chi$. Entirely consistent with the observations made in figure 7(b); $\langle \boldsymbol {F}_{a} \boldsymbol {\cdot }\boldsymbol {\hat {z}} \rangle _n$ is essentially zero for $\chi \gtrapprox 0.5$ or all non-tumbling cases. We observe that when added mass contributes significantly $\langle \boldsymbol {F}_{a} \boldsymbol {\cdot }\boldsymbol {\hat {z}} \rangle _n >0$, which implies that the added mass acts along the rise direction, thus leading to an increased $\langle v_z \rangle _n$, compared with a hypothetical case without added mass.

4. Horizontal periodic motion of the particles

In this section we will focus on the motion in the horizontal plane, which is induced by unsteady motions in the wake of the rising particles. We plot horizontal projections of representative particle trajectories for 18 different aspect ratios in figure 8 in order to illustrate the richness of the patterns of motion. The data is smoothed and corrected for drift.

Figure 8. Each subfigure shows a representative trajectory for a specific $\chi$ as seen from the top with the particle moving from the black circle to the cross. Particle position is indicated by the dots, the colour indicates the degree of alignment between the pointing vector and the velocity vector. The solid grey lines are shown every 0.02 s (5 frames) and indicate the horizontal projection of the pointing vector, with the dot at the end indicating the upward direction. The coloured dot in the bottom right of each subfigure indicates the assigned regime.

In general, the patterns for the oblate particles (figure 8a–i) are regular and periodic, but even there the differences with varying $\chi$, e.g. in amplitude or eccentricity, are significant. For prolate particles (figure 8k–t), the trajectories are a lot more irregular compared with their oblate counterparts. The shape of the trajectories appears similar between the longitudinal (yellow) and the broadside (orange) regimes, but the amplitude of the excursions varies significantly between these two. What clearly stands out against the seemingly chaotic tracks observed otherwise in this regime is how stable and close to circular the helical modes for $\chi = 4$ and $\chi =5$ (figure 8r,t) are, this in addition to the near constant particle alignment. It should be noted that the existence of two rather distinct rise patterns is not limited to the most prolate particles. Albeit less pronounced, another example for which two different rise patterns are encountered can be seen in figure 8(d) for $\chi = 0.33$. Here, an almost planar oscillation switches to a close to circular one as the particle rises. However, in this case the transition is reversible. The distinct patterns therefore appear to be transient states of the particle dynamics in this instance. Such behaviour is similar to that observed for rising bubbles that transition from planar to helical motion (Mougin & Magnaudet 2006; Shew & Pinton 2006). We will revisit the interesting case of the most prolate spheroids in § 6.2. First, we proceed to describe and classify the horizontal motions in terms of frequency ($f$) in the next section. Results for other basic properties (the amplitude, $a$, and eccentricity, $\eta$, of the oscillations) are presented in appendix E to keep the main text succinct.

4.1. Frequency analysis

We present results for the oscillation frequency of the particles (see figure 9) in terms of the dimensionless Strouhal number given by

(4.1)\begin{equation} St \equiv \dfrac{f\kern0.7pt D}{V_b}. \end{equation}

This definition makes use of the buoyancy velocity scale $V_b$ and the volume equivalent sphere diameter $D$, which are both independent of $\chi$. The frequency $f$ is the one obtained using the approach detailed in § 2.3.

Figure 9. Oscillation frequency as a function of particle aspect ratio presented in terms of the Strouhal number using the volume equivalent sphere diameter and the buoyancy velocity scale. The black filled circles show the Strouhal number using the frequency of the precession corrected trajectory. For prolate particles, $St$ is also shown using the frequencies $f_{\perp }$ (blue) and $f_{\parallel }$ (red) of the two different oscillation modes identified in § 4.2. The grey lines at the top show the current results using the $\chi$ dependence as proposed in the literature. The inset of the figure depicts (for prolate particles) a measure of the kinetic energy in the two respective oscillation modes normalized by the vertical kinetic energy. Showing that, except for $\chi = 5$, the oscillation associated with $v_{\perp }$ is dominant.

Figure 9 shows that $St \approx 0.12$ and this value is almost constant for the full range of aspect ratios investigated here. Data from Horowitz & Williamson (2010) for rising spheres show $St \approx 0.1$, which is within the scatter of our experiments, albeit slightly lower than the mean observed here. The mean over all $\chi$ is $St = 0.117$, corresponding to an oscillation frequency of 1.77 Hz with a standard deviation of 0.1 Hz. This result is very surprising and indicates that the frequency is neither dependent on $\chi$ nor on $Re$ in the present case. The constant frequency for varying $\chi$ is different from related findings in literature. Fernandes et al. (2007) found that for oblate disks, the relevant length scale at low $Ga$ is $(dh)^{1/2}$ with a weak dependence on $Re$. This scaling is derived from work on tumbling cards by Mahadevan et al. (1999) and later by Jones & Shelley (2005), giving an empirical dependence of the Strouhal number on $Re$ (or equivalently on $C_d$) and $\chi$ according to $St^* = St(2/C_d)^{1/2}\chi ^{1/2} = fd/\sqrt {(\varGamma -1)gd}$. The latter is found to collapse all their data in the range $80 < Re < 330$ and $0.1<\chi <0.5$. In the work by Toupoint et al. (2019) similar results are found for $St$ for heavy prolate cylinders (which behave similar to our longitudinal regime): $St^{**} =f\sqrt {hd}/\sqrt {(\varGamma -1)gd}$. Both definitions, $St^*$ (for oblate) and $St^{**}$ (for prolate particles), clearly indicate a dependence on $\chi$. For comparison, we show both definitions in figure 9 as grey lines with triangular markers. For both oblate and prolate particles, the use of $St^*$ and $St^{**}$ does not collapse the results obtained here. It should be mentioned however that, for prolate particles, both expressions lie within the spread in the observed frequencies (as evidenced by the error bars). Therefore, no definitive conclusions can be drawn regarding the scaling even though $D$ appears to be the more appropriate choice based on the present results. However, the fact that the scaling observed in Fernandes et al. (2007) does not apply here is not entirely surprising, since the results by Mahadevan et al. (1999) and Jones & Shelley (2005) are strictly valid for heavy particles $\varGamma = {O}(10^3)$ only. While they do agree with findings at $\varGamma \approx 1$ (Fernandes et al. 2007), they fail for the present case where $\varGamma <1$.

In interpreting the apparently different frequency scalings, it is important to note that for both Fernandes et al. (2007) and Toupoint et al. (2019), the density ratios ($\varGamma \approx 1$ and $\varGamma \approx 1.16$, respectively) are significantly higher than that considered here ($\varGamma \approx 0.53$). It therefore appears that the stronger coupling to the fluid motions are a likely explanation for the different behaviour observed in the present case. This is in line with conclusions drawn for spherical bodies by Karamanev & Nikolov (1992), Horowitz & Williamson (2010), Mathai et al. (2018) and Auguste & Magnaudet (2018) who found regime changes related to changes in the density ratio, and the rotational moment of inertia. Details, especially on how the apparent insensitivity to geometry variations at low density ratios emerges, remain nevertheless unclear at this point.

4.2. Prolate spheroids; two interacting modes

Next, we disentangle the complexity of the motion of the prolate particles, which leads to the rather irregular patterns in the trajectories displayed in figure 8(k–q,s). The general tendency of particles to rise with their largest cross-section perpendicular to gravity exposes an ellipsoidal cross-section to the flow in the case of prolate spheroids. This is in contrast to oblate particles, where this area is circular. Hence, instead of only one length scale for oblate particles (the diameter), two unique length scales, i.e. the long ($h$) and the short axis ($d$) of the ellipsoidal cross-section, are relevant for the periodic motions of prolate spheroids. In this section we will show that this gives rise to two distinct oscillation modes that coexist and interact with each other producing the observed dynamics. Each of these modes will be shown to have a distinct frequency.

For the analysis, we consider projections of the velocity $\boldsymbol {v}$ and of the pointing vector $\boldsymbol {\hat {p}}$ onto the horizontal plane, as indicated in figure 10(a). We decompose the horizontal velocity into a component along the pointing vector ($v_{\parallel }$) and a component perpendicular to it ($v_{\perp }$). Without loss of generality, we define the perpendicular direction as $\boldsymbol {\hat {n}}= \boldsymbol {\hat {z}}\times \boldsymbol {\hat {p}} / \| \boldsymbol {\hat {z}} \times \boldsymbol {\hat {p}} \|$, so that the velocity $v_{\perp }$ attains both positive and negative values. This decomposition results in two periodic velocity signals that individually describe the motion of the particle in a plane spanned by $\boldsymbol {\hat {p}}$ and $\boldsymbol {\hat {z}}$ ($v_{\parallel }$) and perpendicular to it ($v_{\perp }$), as illustrated in figure 10(a,b). On this basis, it is now possible to discern the frequencies $f_\perp$ and $f_\parallel$ of these different motions by considering the primary peak of the autocorrelation of $v_{\perp }$ and $v_{\parallel }$, respectively. Note that no evidence for significant interactions of the two modes at frequencies corresponding to $f_{\perp } + f_{\parallel }$ was observed. Results at $\chi = 2$ are displayed in figure 10(c) in terms of a probability density function (PDF) of frequencies corresponding to individual periods.

Figure 10. (a) Decomposition of the velocity vector in the horizontal plane into a component along ($v_\parallel$) and perpendicular ($v_\perp$) to the pointing vector. (b) Graphic showing the oscillation modes associated with the decomposition of the velocity; the top sketch shows the longitudinal and the bottom sketch the broadside oscillations for the same particle. (c) Probability density function of the frequencies obtained for aspect ratio $\chi = 2$. Two distinct frequency peaks are visible which are associated with the two modes of the geometry.

It is clearly evident that the periodic motion for $f_\perp$ and $f_\parallel$ are distinctly different for this particle. Furthermore, we included the PDF of frequencies corresponding to oscillations of $\theta _{\hat {g}}$, i.e. the angle between the pointing vector and the $z$-axis (see figure 10b) as yellow symbols in figure 10(c). We see that the frequency of particle alignment is identical to that of $v_{\parallel }$. This is similar to oblate bodies where these frequencies also match and are related to the same flow mechanism of vortex shedding inducing both horizontal translation and a periodic oscillation in the particle orientation. The oscillation in the $v_{\perp }$ direction, on the other hand, more resembles the pattern induced by vortex shedding behind an infinitely long cylinder. The latter might also induce a rotation around the $\boldsymbol {\hat {p}}$ ($\psi$), which remains inconsequential to the particle alignment however.

The mean of $f_\perp$ and $f_\parallel$ is non-dimensionalised using (4.1) and plotted in figure 9 for all prolate non-tumbling particles. This figure confirms that the frequencies corresponding to motions in different planes are indeed different at all $\chi$ considered. When rescaled according to (4.1), both $f_\perp$ and $f_\parallel$ individually also exhibit a significant dependence on $\chi$, making it even more surprising that this normalization works well without the decomposition (the black data points).

For $v_{\perp }$, we see the Strouhal number increasing with $\chi$. This is consistent with the notion that $d$ is the relevant length scale in this case, since $d$ decreases with $\chi$, while $D$, as employed in (4.1), is constant. In the case of $v_{\parallel }$ we observe a general decrease of the Strouhal number for increasing anisotropy in the longitudinal regime, which is consistent with the increase of the length scale $h$. This suggests that, for the oscillation perpendicular to the pointing vector, the length scale $d$ is appropriate and, similarly, $h$ is characteristic for the oscillation parallel to $\hat {p}$. This yields a more constant value of $St_{\perp }$ and $St_{\parallel }$ in the longitudinal regime; however, the behaviour in the broadside regime (orange) is not consistent for $f_\parallel$, this could be explained by the fact that the oscillations in this direction are almost non-existent in this regime.

Apart form their frequencies, it is also relevant to consider how the two modes differ in energy. To this end, we show a measure for the kinetic energy in the respective oscillation modes in the inset of figure 9. Here, the values are normalized by the mean kinetic energy in the vertical according to $\langle \tilde {v}^2 \rangle _t / \langle v_z^2 \rangle _t$, where $\tilde {v}$ stands for either $v_{\perp }$ or $v_{\parallel }$. In general, the energy in the oscillation along $\boldsymbol {\hat {n}}$ is larger than the one in the plane containing $\boldsymbol {\hat {p}}$, such that the vortex-shedding mechanism that induces oscillations normal to $\boldsymbol {\hat {p}}$ is the dominant one for most particles. It is also important to note that for both, $v_{\perp }$ or $v_{\parallel }$, the energy of the horizontal motions decreases to very low values towards the broadside regime. Additionally, the pointing vector oscillations are minute in this regime, yet measurable. Finally, the situation changes for $\chi = 5$, where the energy in the horizontal motion picks up significantly and $v_\parallel$ becomes the stronger mode, in § 6.2 this will be associated to the transition to the helical rising pattern.

A final point concerns the question whether the oscillations associated with $v_\perp$ are connected to rotations around the pointing vector, $\psi$. In the numerical work by Namkoong, Yoo & Choi (2008) rising and falling cylinders were studied at particle Reynolds numbers ranging from 66 to 185. They found that translations of the particle in the $\boldsymbol {\hat {n}}$-direction are indeed related to body rotations around its symmetry axis inducing a lifting force on the body, i.e. a Magnus lift force. In our experiments we do not observe this direct connection as the rotation $\textrm {d}\psi /\textrm {d}t$ remains weak and is non-periodic as well as only weakly correlated to $v_\perp$. A similar observation was made by Toupoint et al. (2019), but was hampered by their limited accuracy in resolving rotations around $\boldsymbol {\hat {p}}$. Here, this accuracy is $\pm 0.5^\circ$, which puts this finding on stronger footing. We therefore conclude that strong path oscillations in the perpendicular direction are not always induced by a rotation induced Magnus lift force and, thus, are not necessarily coupled to the particle's rotational dynamics.

5. Statistics of the particle orientation

In this section we will consider particle orientations and alignment. For anisotropic bodies, the interaction between the particle's orientation and its direction of translation is a crucial aspect that determines the particle kinematics and dynamics. We will see that this effect is strongly related to the onset of regime changes.

5.1. Preferential orientation and alignment

To capture the mean orientation of the particles around which the pointing vector oscillates, we consider the angle with respect to the vertical direction as defined in figure 1(c) using the definitions

(5.1)\begin{equation} \bar{\theta}_{\hat{g}} =\left\{ \begin{matrix} \bar{\theta}_{{hor}} = \sin^{-1}\left(\left\|\left\{ \langle\hat{p}_x\rangle_n, \langle \hat{p}_y \rangle_n\right\}\right\|\right), & (\chi < 1),\\ \bar{\theta}_{\hat{z}} = \cos^{-1}\left(\left| \left\langle \hat{p}_z\right\rangle_n\right|\right), & (\chi > 1), \end{matrix} \right. \end{equation}

for oblate and prolate particles, respectively (with overlap in the tumbling regime). Our results in figure 11 confirm that the expected mean orientation of $\boldsymbol {\hat {p}}$ is pointing along $\hat {z}$, $\bar {\theta }_{{hor}} = 0$ for oblate, and perpendicular to it ($\bar {\theta }_{\hat {z}} = 90^\circ$) for prolate spheroids.

Figure 11. Statistics of particle pointing vector alignment, on the left $y$-axis in degrees and on the right in terms of the dot products. The black and grey markers show the ensemble-averaged absolute alignments of the pointing vectors with, respectively, the direction of gravity and that of the instantaneous particle motion. The shaded regions show the distribution of $\theta _{\hat {g}}$. Additionally, the statistics of the alignment of the reference orientation are shown indicated by the triangular mark. To observe the convergence of these quantities towards $0^\circ$ and $90^\circ$, respectively, would require more statistics for the tumbling regime.

It is important to note, however, that the mean orientation might never actually be encountered during the rise of the particles. This applies in particular for oblate objects, where the pointing vector tends to orbit around the vertical direction without ever perfectly aligning with it. It is therefore insightful to consider additional orientation statistics including the mean alignment of the pointing vector with the particle motion ($\boldsymbol {\hat {v}}$), but also the mean offset from the reference angle $\bar { \theta }_{\hat {g}}$. To this end, we introduce the angles

(5.2)\begin{gather} \theta_{\hat{g}} = \cos^{-1}{(| \boldsymbol{\hat{p}}(t) \boldsymbol{\cdot} \boldsymbol{\hat{g}} |)}, \end{gather}

(5.3)\begin{gather}\theta_{\hat{v}} = \cos^{-1}{(| \boldsymbol{\hat{p}}(t) \boldsymbol{\cdot} \boldsymbol{\hat{v}}(t) |)}, \end{gather}

as a measure of the instantaneous alignments. For convenience, we also use $| \boldsymbol {\hat {p}}(t) \boldsymbol {\cdot } \boldsymbol {\hat {g}} |$ and $| \boldsymbol {\hat {p}}(t) \boldsymbol {\cdot } \boldsymbol {\hat {v}}(t) |$ without conversion to an angle and, therefore, these quantities are also included in figure 11. The absolute values are taken because of symmetry in the particle geometry, which renders the positive direction of $\boldsymbol {\hat {p}}$ arbitrary. The angles $\theta _{\hat {g}}$ and $\theta _{\hat {v}}$ are schematically shown in figure 1(c). Slight differences arise in the interpretation of the statistics of $\theta _{\hat {g}}$ for $\chi <1$ and $\chi >1$ since $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}$ always changes sign along the trajectory for prolate particles. In essence, however, $\theta _{\hat {g}}$ remains a measure of the actual mean particle orientation encountered in both cases.

In figure 11 the ensemble-averaged alignments $\langle | \boldsymbol {\hat {p}}(t) \boldsymbol {\cdot } \boldsymbol {\hat {g}} |\rangle _n$ and $\langle | \boldsymbol {\hat {p}}(t) \boldsymbol {\cdot } \boldsymbol {\hat {v}}(t) |\rangle _n$ are shown as black circles and grey crosses, respectively. For $|\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}|$, we also show the probability density distribution contours indicating the range of angles covered during an oscillation cycle. The error bars for $|\boldsymbol {\hat {p} }\boldsymbol {\cdot } \boldsymbol {\hat {v}}|$ provide similar information.

The synopsis of the results in figure 11 substantiates the previously defined regimes. Firstly, the tumbling (green) regime is evident from the fact that $|\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}|$ spans the full range from 0 to 1, which is the case for $0.83 \le \chi \le 1.20$. It is surprising to observe that the statistical distribution of $|\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}|$ is almost the same throughout this regime even though the dynamics are very different for $\chi \neq 1$. In particular, the non-spherical tumbling particles show a much more dynamic rotational behaviour. A stronger dependence on $\chi$ is seen in the statistics of $\langle | \boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}} |\rangle _n$ in the tumbling regime, where a more gradual transition between oblate and prolate is evident.

For the special case of the sphere one expects random orientation for the pointing vector direction. This corresponds to $\langle |\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}|\rangle _n = 0.5$ or $\theta _{\hat {g}} = 60^\circ$. Our result slightly deviates from this value, namely we observe $\langle |\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}|\rangle _n = 0.55$. The difference is however within one standard deviation of the statistical convergence of the data, as determined from considering random pointing vector orientations.

The mean alignment of $\boldsymbol {\hat {p}}$ with gravity is very strong and consistent throughout both non-tumbling oblate regimes ($\chi \le 0.75$). In the zigzag (blue) regime we observe a gradual decrease in $\langle \theta _{\hat {g}}\rangle _n$ with decreasing $\chi$. At the same time, also the distribution of $\theta _{\hat {g}}$ becomes narrower, indicating lower fluctuation levels, which is most likely associated with an increase in (added) rotational inertia. For the same range of aspect ratios, $\langle \theta _{\hat {v}} \rangle _n$ is seen to steadily increase with decreasing $\chi$. The main reason for this trend is the increased phase delay between the oscillations of $\boldsymbol {\hat {p}}$ and $\boldsymbol {\hat {v}}$ as we will discuss in more detail in § 5.2.

In the flutter (purple) regime, $\langle \theta _{\hat {g}}\rangle _n$ rises only slightly compared with the zigzag regime. However, $\langle \theta _{\hat {v}}\rangle _n$ increases considerably below $\chi =0.29$ which clearly distinguishes the two regimes. The angle $\langle \theta _{\hat {v}}\rangle _n$ reaches a mean value as high as $60^\circ$ at $\chi = 0.20$. This implies that the particle is predominantly moving ‘sideways’ in the direction of its smallest cross-section (consistent with the statistics of $A_{\perp v}$ reported in the inset of figure 6). Note also that the particles never get close to $\theta _{\hat {v}} \approx 0$ in the flutter regime.

On the prolate side, the different regimes are distinctly noticeable in both the behaviour of $\langle \theta _{\hat {g}}\rangle _n$ and of $\langle \theta _{\hat {v}}\rangle _n$. In the longitudinal (yellow) regime we observe large oscillations in the orientation of $\boldsymbol {\hat {p}}$ relative to $\boldsymbol {\hat {g}}$ as well as to $\boldsymbol {\hat {v}}$. The amplitude of this pointing vector oscillation decreases with increasing $\chi$, as can be seen from the angle $\langle \theta _{\hat {g}}\rangle _n$ which increases and the spread in the data which decreases, drawing closer to the reference orientation of $90^\circ$. When transitioning to the broadside (orange) regime, the averages $\langle \theta _{\hat {g}}\rangle _n$ and $\langle \theta _{\hat {v}}\rangle _n$ increase, while the fluctuations decrease considerably, such that the maximum deviation of the pointing vector from the horizontal becomes less than $10^\circ$.

At $\chi = 5$, data is shown for the rise patterns corresponding to the longitudinal (yellow) state. In this state the inclination $\langle \boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}} \rangle _n$ is identical to that in the yellow regime, however, the value of $\langle \boldsymbol { \hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}} \rangle _n$ is significantly higher (i.e. $\langle \theta _{\hat {v}} \rangle _n$ is lower). This implies that the particle is more aligned with the $\boldsymbol {\hat {v}}$-direction which will be shown to play an important role in triggering the transition to the helical mode in § 6.2. The particle alignment is also visualized in the supplementary movies, showing the coupling between orientation and the particle motion.

5.2. Particle reorientation

Within the Newton drag regime, where pressure induced drag and lift are due to a deflection of the flow, it is immediately obvious that an inclination of the particle's symmetry axis will lead to a lift force. The statistics of particle alignment provide some insight into this facet but lack details on the temporal sequence. The latter is important since it is the temporal dependence that leads, for example, to the net contribution of the added mass force observed in § 3.4. To complete the analysis, we will therefore examine the phase difference $\Delta \phi$ between the periodic motion of the $\boldsymbol {\hat {p}}$- and $\boldsymbol {\hat {v}}$-directions in this section, which is defined such that a positive value of $\Delta \phi$ corresponds to a lagging pointing vector.

To illustrate this, we show examples of trajectories and alignments for five different oblate particles in figure 12. Let us first consider $\chi = 0.75$. In this case, the rotation of the particle and the direction of motion are in phase. This means that when the particle's horizontal velocity is maximal, the pointing vector is also at a maximum deflection from its reference orientation. And when the velocity is minimal the pointing vector direction is also close to the vertical. Note that this does not necessarily mean that $|\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}}|$ remain close to 1 at all times. This can be seen close to the zero crossings, where the horizontal particle position $\tilde {x} = 0$, here the particle orientation ‘overshoots’ resulting in a slight dip in $|\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}}|$. This overshooting is related to the tumbling motion, which is described in § 6.1. With increasing anisotropy, the orientation starts lagging behind the direction of motion more and more (see figure 12: $\chi \le 0.29$ ), which corresponds to a positive phase lag $\Delta \phi$. This causes stronger fluctuations and a general decrease in the alignment between particle orientation and path velocity as $\chi$ decreases in figure 12.

Figure 12. Representative precession corrected trajectories. Trajectories are colour-coded with the instantaneous alignment $|\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}}|$. The grey lines are shown every 5 frames (0.02 s) and indicate the major axis of the projected geometry, which is always perpendicular $\boldsymbol { \hat {p} }$.

To extract quantitative information about the phase delay, we have to employ slightly different approaches for oblate and prolate particles to account for the different pointing vector alignments. In the oblate case we consider the cross-correlation of the periodic signals $\boldsymbol {\hat {p} }_{\tilde {x}}$ and $\boldsymbol {\hat {v}}_{\tilde {x}}$, i.e. the horizontal components of the pointing and of the velocity vectors in the precession corrected frame of reference. For prolate particles, where the pointing vector is oscillating around the horizontal plane, we use a cross-correlation between $\hat {v}_\parallel$ (as defined in figure 10b) and $-\boldsymbol {\hat {p}}_z$. Note that here we focus exclusively on the mode associated with $\hat {v}_\parallel$, since only this one occurs at the same frequency as the pointing vector oscillations (see § 4.2 and, in particular, figure 10c).

Based on the peak of the cross-correlations defined in this way, we determine the phase lag. The corresponding results are shown in figure 13 in terms of the phase angle $\Delta \phi$. In this figure we observe an interesting similarity between oblate and prolate geometries. Close to spherical, the phase lag is essentially zero for both cases. Towards more extreme anisotropy however, $\boldsymbol {\hat {p}}$ increasingly lags behind $\boldsymbol {\hat {v}}$. For oblate particles, we observe that the regime change from zigzag (blue) to flutter (purple) coincides with a levelling-off of the phase lag at around $\Delta \phi \approx 100^\circ$. The flattening out of the $\Delta \phi$-curve also coincides with the onset of a strong precession at $\chi = 0.2$. On the prolate side, the regime transition from longitudinal (yellow) to broadside (orange) oscillations is marked by the onset of a phase lag. The transition to the helical regime occurs at a value of $\Delta \phi$ similar to the one at which oblate particles level off. Yet, in the prolate case, $\Delta \phi$ appears to keep increasing towards more extreme anisotropy, which might be one of the factors that promote the transition to the helical pattern as described in § 6.2.

Figure 13. Phase angle between velocity and pointing vector oscillations at different aspect ratios. The results for the current data set ($Ga \approx 6000$) are represented by the black symbols, while the grey dashed lines are fits from the work by Fernandes et al. (2005) for disks at lower $Ga$.

It is further noteworthy that, for both prolate and oblate geometries, the intermediate phase lags, where the behaviour transitions from $0^\circ$ lag to $100^\circ$, are associated with minima in the amplitude of the particle oscillations (see appendix E). For prolate particles, oscillations of the pointing vector also almost vanish in this regime. It remains unclear, however, to what extent these trends are connected.

The observed trend in the phase lag with $\chi$ for oblate particles is in line with results by Fernandes et al. (2005) for oblate disks at low Galileo number. Their fitted results are also shown in figure 13. We observe that, for increasing $Ga$, the curve for $\Delta \phi$ shifts towards smaller particle anisotropy. Intuitively this makes sense: for larger $Ga$, the pressure forcing on the particle surface becomes larger in magnitude. Neglecting changes in the rotational inertia, thus, the time scale of the particle response becomes shorter with increasing $Ga$, causing the observed $Ga$ trend in figure 13.

In § 3.1 it was noted that the drag coefficient in our ‘zigzag’ (blue) regime collapsed to the data from Fernandes et al. (2007). It appears from figure 13 that the experiments performed in the cited paper all would fall in the blue regime, thus, potentially explaining the collapse of $C_d$ for their experiments, even though their particle anisotropy is much larger, going as low as $\chi = 0.1$.

6. Tumbling and helical rise – a closer look

In this last section we will have a closer look at the two most outstanding regimes: tumbling particles with shapes close to spherical (§ 6.1) and the helical pattern observed for the most prolate spheroids (§ 6.2).

6.1. Tumbling regime: Froude number dependence

In the current data there is a gradual transition from a fluttering to a tumbling state. This regime change has been extensively studied mostly for thin bodies, where the particle length scale in the flow direction is much smaller than the one perpendicular to it. It is common (see, e.g. Dupleich 1949; Willmarth et al. 1964; Smith 1971; Lugt 1983) to characterise the phase space and transition between regimes in terms of a dimensionless moment of inertia $\tilde {I} = I/I_{ref}$, where $I_{ref}$ is the moment of inertia of the cross-section of a reference geometry, often chosen to be a cylinder. This approach works for particles that are quasi-two-dimensional, i.e. whose cross-section does not vary along the axis around which they tumble. The flutter-to-tumble transition received additional attention more recently in experiments by Belmonte et al. (1998), Zhong et al. (2011) and Wang et al. (2013), in numerical work (Andersen et al. 2005b; Lau, Wei-Xi & Chun-Xiao 2018) or through a combination of both in Andersen et al. (2005a). Nevertheless, in all these works, the influence of the geometry, especially of a finite and varying thickness in the flow direction, has not been considered systematically yet. Therefore, the current experiments can provide additional insight by applying this analysis to particles with a significant thickness in the direction of the flow.

The approach using a dimensionless moment of inertia is not well suited for our particles, since there is no universal logical choice for $I_{ref}$ due to the variations of the cross-section along the axis of rotation. Therefore, we follow the approach by Belmonte et al. (1998), and define a dimensionless Froude number

(6.1)\begin{equation} Fr \equiv \dfrac{\tau_p}{\tau_r} = \sqrt{\dfrac{\varGamma d^*}{L}}, \end{equation}

as the ratio of two intrinsic particle time scales $\tau _p$, which is the time scale for a buoyant pendulum, and $\tau _r$, the settling time scale. These are defined by

(6.2a,b)\begin{equation} \tau_p = \sqrt{\dfrac{ L \varGamma}{(1-\varGamma)g}} \quad\text{and} \quad\tau_r = \dfrac{L}{v_0}. \end{equation}

Here, $v_0$ is a velocity scale based on the balance between particle drag force and buoyancy given by $v_0 \sim \sqrt {(1-\varGamma ) g d^*}$ as derived from (3.5) and $d^*$ is the particle length scale in the direction parallel to the flow, i.e. $d^* = h$ for $\chi < 1$ and $d^* = d$ for $\chi \geq$ 1. In ((6.1) and (6.2a,b)) we additionally introduce the length scale $L$, which is defined as $L = d$ for $\chi \leq$ 1 and $L = h$ for $\chi >$ 1 and characterises the largest length scale perpendicular to the mean flow direction. With these definitions, (6.1) simplifies to $Fr = \sqrt {\varGamma \chi }$ for oblate particles and $Fr = \sqrt {\varGamma \chi ^{-1}}$ for prolate ones.

The Froude number as function of aspect ratio is shown in figure 14(a) for the tumbling regime as well as for the two neighbouring ones. In this figure we mark particles that tumble with open and those that do not with filled symbols. We find that at least sporadic tumbling is observed for all particles with $Fr \gtrapprox 0.64$, whereas particles with lower values of $Fr$ never tumble. Hence, the critical $Fr$ appears to be the same regardless of whether the shape is oblate or prolate. Moreover, the critical value of $Fr$ is also similar to that obtained by Belmonte et al. (1998), where the flutter-to-tumble transition for quasi-two-dimensional strips was observed to occur at $Fr = 0.67 \pm 0.05$. This is somewhat surprising and may be coincidental since the mechanism of flipping is actually different for strips and spheroids, as we will discuss later in this section.

Figure 14. Particle Froude number as calculated using the scaling approach by Belmonte et al. (1998) for particles close to the transition to tumbling. The symbols indicate the observed behaviour and particle type and the dashed horizontal line shows the estimated critical $Fr$ value above which tumbling occurs. The $\chi$ dependencies of $Fr$ indicated by the black lines follow from Belmonte et al. (1998) (see text for details). Note also the log-scaling on the $\chi$-axis.

To analyse the transition to tumbling in more detail, we consider the temporal evolution of the particle inclination in terms of $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}$ for which an example is shown in figure 15(a). From this it becomes clear that the tumbling of the particle, indicated by the change in sign of $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}$, does not occur for every single oscillation, which is unlike the behaviour of flat quasi-two-dimensional strips observed in Belmonte et al. (1998) at similar $Re$ ($3\times 10^3 < Re < 4 \cdot 10^4$). A similar intermittent fluttering and tumbling behaviour was encountered in the ‘apparently chaotic’ regime reported by Andersen et al. (2005a,b). For the present data, a gradual increase in amplitude in the particle inclination, $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}$, leading up to flipping events can be observed in figure 15(a). It appears that the particle has to build up enough rotational momentum before it can eventually flip over. Only once the amplitude is large enough, after a couple of oscillations, a flip happens and $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}$ changes sign. The direction of the flip (clockwise or counterclockwise) appeared to be random and both appeared equally probable.

Figure 15. (a) Particle alignment $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {g}}$ over time for a single experiment for $\chi = 0.83$. (b) Schematic representation of the tumbling mechanisms for spheroids and flat plates. On the left a section of one of the trajectories for $\chi = 0.83$ is shown and on the right the rise pattern of a strip (taken from Belmonte et al. 1998) is displayed very close to the tumbling transition (but not actually tumbling). The orientation of the particles largest cross-section is indicated by the black lines. The points of the trajectory at which the particle flips (or would flip) are indicated by the large arrows and the label I marks the same instances in (a,b).

However, differences between the tumbling motion of two-dimensional (2-D) strips and 3-D spheroids are not restricted to the frequency of the tumble. Also the way in which the spheroids tumble differs significantly from the mechanism that is observed for thin bodies, which also renders the tumbling trajectories fundamentally different. For thin plates, the tumbling motion is associated with oblique trajectories that exhibit a very strong mean drift (Belmonte et al. 1998; Pesavento & Wang 2004; Andersen et al. 2005b), whereas we find that for spheroids the trajectory remains vertical. The two distinct mechanisms of tumbling are exemplified by the two particles trajectories in figure 15(b); the left one corresponds to a spheroid in the present work and the one on the right shows the trajectory of a thin plate close to tumbling as reported in Belmonte et al. (1998). The key difference is the point of the trajectory at which the flip happens. For spheroids, the flip occurs when the path oscillation is at a minimum, whereas for flat plates, the tumbling occurs at the maximum path amplitude. The difference is related to the phase difference $\Delta \phi$. For a flat plate, $\Delta \phi$ is close to $90^\circ$, thus, the pointing vector is horizontal at points of maximum amplitude, whereas $\Delta \phi = 0$ in our case results in a flip at the centre. The difference is most likely related to the mechanism that causes the particle to rotate. For flat plates, the particle rotates primarily due to a torque induced by pressure forces. However, such torque is minimal for spheroids close to spherical ($\chi \sim 1$). In this case, the main mechanism is therefore the asymmetry in the skin friction caused by vortex shedding. Additionally, rotational added mass plays an important role for flat plates while its effect is negligible for particles close to spherical.

As mentioned, this results in markedly distinct trajectories from those of flat plates. The oscillation frequencies and amplitudes observed for tumbling spheroids are very close to those of the perfect sphere. However, the trajectories differ between tumbling spheroids and the sphere, as can be seen in figure 8(h–l). The former show more sudden changes in direction which can be related to the strong rotations and changing alignment for these particles. The rise velocity and drag are also clearly affected by this.

6.2. The odd one out: The helical regime

All cases discussed so far have in common that the pointing vector performs an oscillatory motion around the reference orientation for which the largest area of the geometry is perpendicular to gravity. Remarkably, however, we found that, for $\chi \geq 4$, also a distinctly different rise pattern is possible. In what we refer to as the ‘helical’ rise pattern, it was found that the particles enter an orbit in which both angles $\theta _{\hat {g}}$ and $\theta _{\hat {v}}$ remain approximately constant in time. In the following we provide details on the helical orbit and on how the particles at $\chi = 4$ and $\chi = 5$ transition to this pattern. All properties of this regime are tabulated along with the other results from this paper in table 2 of appendix A.

The difference between the non-helical (figure 8q,s) and helical patterns (figure 8r,t) is stunningly clear from the particle trajectories alone. With $\langle \eta \rangle _n = 0.99$ for both $\chi = 4$ and $\chi = 5$, the eccentricity (as defined in § E.2) of the helical trajectories is very close to that of a perfect circle for which $\eta = 1$. Furthermore, the helical path is almost perfectly periodic, which is in complete contrast to all other observed prolate modes that generally display very chaotic and complex motion.

The 3-D representation in figure 16(a–c) shows all the different rise modes possible for the particles with $\chi \geq 4$ in the form of a render. For $\chi = 4$ (figure 16a), the broadside (orange) and helical (red) regimes coexist. In this case, the regime appears to be determined exclusively by the particles initial release since in all experiments performed we have never observed a regime transition during the rise. In case the particle is released with minimal disturbance and $\boldsymbol {\hat {p} } \boldsymbol {\cdot } \boldsymbol {\hat {g}}$ is close to 0, the particle tends to go into the broadside regime. However, if the release is less controlled or the particle is inclined at release, it was found to gain significant velocity along the direction of the pointing vector and enter the helical pattern.

Figure 16. Particle trajectories obtained from the experimental measurements with superimposed renders of the particles at their instantaneous orientation plotted every 0.1 s. (a) The two coexisting rise patterns for $\chi = 4$ are shown. These correspond to the broadside (orange) regime on the left where the particle rises almost vertically with minimal oscillations of the pointing vector and the helical (red) regime on the right. (b) The two regimes for $\chi = 5$: on the left is the longitudinal (yellow) regime showing large amplitude oscillations in the pointing vector oscillations, on the right its helical (red) regime state. (c) Two cases for $\chi = 5$ where a transition from the longitudinal regime to the helical one is observed (marked by arrows). (d) For the two transition cases in (c) we show the alignment $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}}$ as a function of time, clearly illustrating the transition as well as the distinct behaviour in the two regimes. The coloured arrows correspond to the marked points in (c).

The dependence on the release conditions is similar at $\chi = 5$. However, when not in the helical mode the particle at this aspect ratio displays a pattern that is reminiscent of the longitudinal regime rather than the broadside one with large oscillations of the pointing vector inclination (see figure 16b). Contrary to the case at $\chi = 4$, we identified transitions from the longitudinal to the helical regime in two experiments at $\chi = 5$ (but never the reverse transition). The transitional trajectories are shown in figure 16(c) with the approximate point of transition indicated by the arrows.

The gradual transition is also visible from the corresponding plot of $\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}}$ in figure 16(d). Here, we initially see the large changes in particle inclination characterising the longitudinal regime. As time advances these fluctuations slowly decrease and the particle inclination subsequently remains largely constant at a value close to $|\boldsymbol {\hat {p}} \boldsymbol {\cdot } \boldsymbol {\hat {v}}| = 1$. This reduction in angular fluctuations is accompanied by a gradual growth of the particle horizontal velocity $v_{\parallel }$ (not shown), which may help to stabilize the particle at a non-zero inclination.

It should also be remembered that the fact that $\boldsymbol {\hat {p}}$ comes close to aligning with $\boldsymbol {\hat {v}}$ is related to the increase of the phase lag $\Delta \phi$ with increasing $\chi$. Unlike for oblate particles, $\Delta \phi$ was seen to increase beyond $\Delta \phi \approx 100^\circ$ for prolate spheroids (see § 5.2). Based on the above observations, we suspect that the longitudinal regime is not stable at $\chi = 5$ and that all particles may eventually transition to the helical regime at this aspect ratio. Due to the small oscillations of the pointing vector alignment in the helical regime it is deemed unlikely that the particle will transition back to the longitudinal one without outside forcing.

To characterise the helical regime in more detail, we plot the amplitude ($\langle a/D \rangle _r$), the oscillation frequency ($\langle\, f\rangle _r$) and the rise velocity (($\langle v_z \rangle _r$) as a function of $\langle \theta _{\hat {v}}\rangle _r$ in figure 17 for both $\chi = 4$ and $\chi = 5$. Note that the averaging here, $\langle \cdot \rangle _r$, is performed over individual runs to highlight run-to-run variations. Also included in the plot are the standard deviations around this mean for all quantities in the form of error bars. Furthermore, for some runs a rotation around $\boldsymbol {\hat {p}}$ was observed, while in other cases it was absent. Somewhat surprisingly, such rotations appeared to have no significant effect on the particle motion.

Figure 17. Properties of the helical trajectories for $\chi = 4$ and $\chi = 5$ as a function of the angle $\langle \theta _{\hat {v}}\rangle _r$: (a) amplitude $\langle a/D \rangle _r$, (b) oscillation frequency $\langle\, f\rangle _r$ and (c) rise velocity $\langle v_z \rangle _r$. All averages are per run and error bars correspond to $\pm$ one standard deviation. Three different particles were used at $\chi =5$ and data corresponding to the same particle are marked by different shades of red.

The first thing to note about $\langle a/D \rangle _t$ in figure 17(a) is that the amplitude of all helical trajectories is much larger than for any other regime where $\langle a/D \rangle _n < 1$ consistently (see figure 20a). Moreover, the standard deviations for $\langle a/D \rangle _r$ but also for $\langle \theta _{\hat {v}}\rangle _r$ are minute (average of the per run standard deviation: 0.0210 $a/D$ for $\chi = 4$ and 0.0195 $a/D$ for $\chi = 5$), demonstrating how consistent the trajectories are. Another interesting observation is that there is quite some spread in the data between runs, and that the amplitude appears to be inversely correlated to the angle $\langle \theta _{\hat {v}}\rangle _r$. This trend is consistent for both aspect ratios considered.

Also the frequency (see figure 17b) varies very little over the course of a single run. However, we do note a spread in the data between subsequent runs. The frequencies of the helical patterns are substantially lower than for all other regimes for which the average frequency is 1.77 Hz. Also $\langle\, f\rangle _r$ is seen to decrease with increasing $\langle \theta _{\hat {v}}\rangle _r$. However, there is an offset between the data at $\chi = 4$ and $\chi = 5$.

Finally, we show the mean rise velocity in figure 17(c). Similar to the amplitude and frequency, we observe nearly constant velocities within individual runs with small standard deviations. Although, again there is substantial variation between the runs that is inversely correlated with $\langle \theta _{\hat {v}}\rangle _r$. The mean rise velocity of the particles for $\chi = 4$ is lower in the helical regime ($0.265\ \textrm {m}\,\textrm {s}^{-1}$) compared with the broadside velocity ($0.334\ \textrm {m}\,\textrm {s}^{-1}$). For $\chi = 5$, on the other hand, the mean vertical velocity is higher in the helical regime ($0.297\ \textrm {m}\,\textrm {s}^{-1}$) than the longitudinal rise velocity ($0.255\ \textrm {m}\,\textrm {s}^{-1}$).

The observation that the particle inclination is so constant suggests that for this regime there is no periodic vortex shedding which would cause variations in rise velocity or particle orientation. We would therefore expect a continuous wake behind the particle with minimal discrete vortex-shedding events similar to that for slender ‘hydrodynamic’ objects and it will be interesting to study this in the future.

To rule out that the helical pattern was caused by some imbalance of the particle caused by inaccuracies in the manufacturing, we tested three different particles at $\chi = 5$. All three of them entered similar helical patterns upon release. Results corresponding to each individual particle are marked in figure 17, as indicated by the different shades of red, and reveal that the parameters of the helical orbit do vary somewhat per particle. However, also for $\chi = 4$, where only a single particle was used in this regime, there are still significant differences between the runs. It appears therefore that the inclination angle $\theta _{\hat {v}}$ is very sensitive to the particle geometry and mass distribution, but potentially also to the release conditions. The sensitivity to particle properties is interesting since it opens up possibilities to design for a specific desired behaviour. It should be noted that in the non-helical rise pattern the dynamics and kinematics of the different particles at $\chi = 5$ were indistinguishable from one another.

7. Conclusion and discussion

The main outcome of this study is a classification of the rich behaviour of freely rising spheroidal particles across a wide range of particle aspect ratios, at high Galileo number $Ga \approx 6000$.

The spheroidal geometry allows for a gradual transition from spherical to (slightly) anisotropic particles. Our data reveals that even small levels of anisotropy in the shape have a profound impact on the particle dynamics. This is reflected in significantly reduced particle rise velocities and substantially increased rotational dynamics once $\chi$ deviates from 1. The latter can cause the particles to flip over, resulting in a tumbling motion for $0.83 \leq \chi \leq 1.20$. Interestingly, the bounds of the tumbling regime here are governed by the same time scales ratio, captured by the Froude number $Fr$, as the ‘flutter-to-tumble’ transition for flat plates. However, the precise nature of the tumbling motion differs in several aspects between spheroids and flat plates. For spheroids, (i) the tumbling does not result in oblique trajectories, (ii) the torque inducing rotation is primarily due to skin friction and not pressure induced, and (iii) the tumble occurs at a different time during the oscillation cycle.

For non-tumbling oblate spheroids ($\chi \le 0.75$), period-to-period variations in the oscillation cycle are lower. In particular, the erratic tumbling motion gives way to regular oscillations of the pointing and the velocity vectors with a phase shift $\Delta \phi (\chi )$. This more regular motion is due to the fact that there is only a single transverse length scale, namely the diameter $d$ of the circular cross-section exposed to the flow. We found that for the ‘zigzag’ regime (blue), ranging from $0.29 \lessapprox \chi \lessapprox 0.75$, the average cross-flow area is almost exactly the maximum cross-sectional area of the particle. This results in a constant drag coefficient within this regime of $C_d \approx 1.2$ ($v_z \sim \sqrt {h}$), identical to that found for disks at $Ga = 100-300$ by Fernandes et al. (2007), but unlike these authors we do not observe a constant velocity along the path.

At the most extreme oblate aspect ratios considered here ($0.2 \leq \chi \lessapprox 0.29$), the direction of particle motion is no longer primarily perpendicular to the largest cross-section of the geometry. Instead, the particles perform periodic leave-like oscillations that lead to significant horizontal movement but also significant fluctuations in the vertical velocity. As a consequence, $C_d$ is seen to increase in this ‘flutter’ regime. Remarkably, the added mass force is found to significantly contribute to the vertical particle motion, even when terminal velocity has been reached. This effect amounts to more than $0.6\|\boldsymbol {F}_B\|$ for $\chi = 0.2$. In contrast, we found that along the path the added mass force does not result in a net drag contribution for any $\chi$ here.

We further observe a trend from planar (at $\chi = 0.25$) to strongly precessing oscillations ($\chi = 0.2$) within the flutter regime, which is very similar to findings by Zhong et al. (2011) for thin flat disks (from ‘planar-zigzag’ to ‘transitional’ in their wording). Our transition occurs around the same $Re$ and for a similar dimensionless moment of inertia as in their data when plotted in the phase space according to Field et al. (1997). This analogy invites some speculation regarding a potential ‘spiral’ regime for even more oblate particles (see Zhong et al. 2011, figure 5). The existence of such a regime would give rise to a rather pleasing symmetry with helical trajectories for both, extremely prolate and oblate spheroids.

Generally, the trajectories for prolate particles appear more chaotic compared with their oblate counterparts. We were able to relate this to the fact that the mean cross-flow area of the geometry is no longer circular (as is the case for oblate particles) but ellipsoidal. This results in two transverse length scales that grow increasingly more distinct with increasing anisotropy. Our analysis revealed that these two length scales give rise to two different modes of oscillation, one in the plane containing the pointing vector and one perpendicular to it. Those oscillations occur at distinctively different frequencies and their interaction gives rise to the very complex and seemingly chaotic trajectories.

For $1.2 \lessapprox \chi \lessapprox 2.5$, we define the ‘longitudinal’ regime which is characterised by large amplitude pointing vector oscillations. These oscillations disappear almost completely for $2.5 \lessapprox \chi \lessapprox 4.5$ , we refer to this range as the ‘broadside’ regime. Both of these regimes have a near constant drag coefficient of $C_d = 0.86$ and $C_d =0.67$, respectively, i.e. $v_z \sim \sqrt {d}$. The fact that the drag actually decreases for increasing anisotropy is a remarkable result that is unique for the ‘broadside’ regime. This gives rise to a secondary local maximum in the particle rise velocity at $\chi \neq 1$.

Arguably, the most surprising observation was made for the most extreme prolate aspect ratio, $\chi = 5$. For this case, we found that two completely different rise patterns coexist: one is an oscillation similar to the ‘longitudinal’ regime but with a larger phase difference $\Delta \phi \approx 110$, the other an almost perfectly helical trajectory. The helical rise pattern is triggered by an alignment between the velocity direction and the pointing vector and this transition occurs naturally at $\chi =5$. The same helical pattern can also be triggered for $\chi = 4$, but it has to be forced via a large inclination at the particle release in this case. The helical pattern truly stands out among all other trajectories observed since its the only case where the particle inclination is fixed and not oscillating around the orientation of maximum cross-flow area.

It is tempting to draw a parallel between the present observations and the helical trajectories reported for bubbles (Lindt 1972; Clift, Grace & Weber 1978). Mougin & Magnaudet (2006) showed numerically that deformed bubbles transition from planar type oscillations to a helical path and this was confirmed by Shew & Pinton (2006) in a simplified model based on the Kelvin–Kirchhoff equations. Obviously there are differences compared with the present work. This concerns not only the boundary conditions at the particle/bubble surface, but also the alignment. For uncontaminated rising bubbles, the major axis is aligned with the radial direction, i.e. normal to the path (Veldhuis 2007, chap. 7), whereas in our case the pointing vector is directed tangentially to the spiraling motion.

On the oblate side, it is possible to get insight into the $Ga$-variation of the regime boundaries by comparing the present results to those of Fernandes et al. (2005, 2007). Across a wide range of different properties (Reynolds number, amplitude, eccentricity, phase difference), it appears consistently that transitions appear at lower levels of anisotropy (i.e. larger $\chi$ for oblate particles), the higher $Ga$. This is also in agreement with the trend over the (small) range of $Ga$ considered in Fernandes et al. (2005) for the phase lag.

We find that results for $C_d$ can be divided into three approximately constant values: $C_d\approx 1.22$ for the ‘zigzag’ regime, $C_d \approx 0.83$ for the ‘longitudinal’ and ‘tumbling’ regimes and $C_d \approx 0.67$ for the ‘broadside’ regime. These values seem related to the relevant particle dimension setting the size of the vortex-shedding system. For the ‘zigzag’ regime, this is the major axis $d$ (oblate), which results in a large vortex system and, thus, the highest drag coefficient. For the ‘broadside’ motions, the vortex shedding is associated with only the minor axis $d$ (prolate) resulting in the lowest drag. Finally, for the intermediate ‘tumbling’ and ‘longitudinal’ regimes, we observe oscillations associated with both particle dimensions $d$ and $h$ (major and minor axes). Correspondingly, the resulting drag is intermediate between the values obtained for the other cases.

Another important aspect of the current work is that it sheds light on the effect of particle (rotational) inertia. In the present case, $\varGamma \approx 0.53$ and, thus, also the relative (rotational) inertia is smaller. This implies that the particle reacts more swiftly and stronger to the fluid forcing. We observe these effects indirectly by comparing our results to those of settling particles. From this we find that the drag coefficient of particles with similar geometry is much lower for heavy particles compared with the current data set (with the notable exception of the ‘broadside’ regime). Keep in mind here, however, that the settling data is compiled across a wide range of geometries with no distinction made between oblate and prolate cases. Therefore, this does not necessarily imply that prolate particles rise faster than they settle.

Furthermore, we observe that the particle frequency of oscillation is largely independent of $\chi$. This is contrary to findings by Fernandes et al. (2007) who identified a dependence on the cross-flow length scale $d$ at lower $Ga$. This difference is most likely related to the lower density ratio ($\varGamma$) and rotational moment of inertia ($I^*$) in our case. Further research will be required to shed light on how the oscillation frequency depends on these parameters.

As a final remark, we would like to point out that no significant rotations around the particle symmetry axis were observed for any aspect ratio in this work.