2.1 Problem definition

We use direct numerical simulation to explore the dynamics of light spheres throughout the range of solid-to-fluid density ratios smaller than unity in regimes where viscous effects, although small, are still significant. The sphere is considered to have a homogeneous density $\unicode[STIX]{x1D70C}_{s}$ , diameter $d$ , volume ${\mathcal{V}}=(1/6)\unicode[STIX]{x03C0}d^{3}$ , moment of inertia $I_{s}=(1/10)\unicode[STIX]{x1D70C}_{s}{\mathcal{V}}d^{2}$ and rises in an infinite body of Newtonian fluid with uniform density $\unicode[STIX]{x1D70C}$ and kinematic viscosity $\unicode[STIX]{x1D708}$ . The problem depends on two control parameters, namely the body-to-fluid density ratio $\overline{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}_{s}/\unicode[STIX]{x1D70C}$ and the Archimedes number $Ar=V_{g}d/\unicode[STIX]{x1D708}$ which is a Reynolds number based on the gravitational velocity $V_{g}=\{|\overline{\unicode[STIX]{x1D70C}}-1|gd\}^{1/2}$ ( $g$ denotes gravity). It is also frequently useful to refer to a Reynolds number $Re_{m}=V_{m}d/\unicode[STIX]{x1D708}$ based on the average vertical velocity $V_{m}$ measured throughout the sphere ascent (or during some specific stages), since this velocity is usually the one determined in experiments. Periodic paths and sphere oscillations with frequency $f$ are then characterized by a reduced frequency or Strouhal number $St=fd/V_{m}$ , and this is the definition used throughout the paper.

In what follows we explore the dynamics of spheres with $\overline{\unicode[STIX]{x1D70C}}\leqslant 0.99$ and $125\leqslant Ar\leqslant 700$ , which, as will be shown later, approximately corresponds to $150\leqslant Re_{m}\leqslant 1.2\,10^{3}$ for $O(1)$ sphere-to-fluid density ratios. Available experiments have of course covered a much wider range of Reynolds number ( $Re_{m}\leqslant 1.5\times 10^{4}$ approximately), and most of them focused on values of $Re_{m}$ larger than those achieved here. Nevertheless extensive results are available at $Re_{m}=450$ , where a variety of paths was encountered (Horowitz & Williamson Reference Horowitz and Williamson2010); detailed experiments with light ( $\overline{\unicode[STIX]{x1D70C}}\approx 0.55$ ) and very light ( $\overline{\unicode[STIX]{x1D70C}}=0.02$ ) spheres were also carried out for $Re_{m}\lesssim 400$ (Veldhuis & Biesheuvel Reference Veldhuis and Biesheuvel2007) and $Re_{m}=O(10^{3})$ (Veldhuis et al. Reference Veldhuis, Biesheuvel and Lohse2009), respectively. The Reynolds number range spanned by present computations makes comparison with these data sets possible.

2.2 Computational strategy

The starting point of our computational strategy follows the general approach described in Mougin & Magnaudet (Reference Mougin and Magnaudet2002a ) and used in Mougin & Magnaudet (Reference Mougin and Magnaudet2002b ), Auguste (Reference Auguste2010) and Auguste et al. (Reference Auguste, Magnaudet and Fabre2013) to study path instability of oblate spheroidal bubbles and falling disks, respectively. In brief, the Navier–Stokes equations governing the absolute fluid velocity, $\boldsymbol{U}$ , are written in a system of axes translating and rotating with the sphere. As shown in Mougin & Magnaudet (Reference Mougin and Magnaudet2002a ), the corresponding form of these equations (assuming $\unicode[STIX]{x1D70C}=1$ ) is

with boundary conditions $\boldsymbol{U}=\boldsymbol{V}+\unicode[STIX]{x1D734}\times \boldsymbol{r}$ at the sphere surface ( $\Vert \boldsymbol{r}\Vert =d/2$ ) and $\boldsymbol{U}\rightarrow \mathbf{0}$ for $\Vert \boldsymbol{r}\Vert \rightarrow \infty$ , $\boldsymbol{r}=\mathbf{0}$ corresponding to the sphere centre. The reader is referred to the above reference for a discussion on the computational advantages of this formulation compared to others. Equations (2.1)–(2.2) are integrated in time by keeping the translational sphere velocity, $\boldsymbol{V}$ , and rotation rate, $\unicode[STIX]{x1D734}$ , constant during each time step. This integration is achieved using the finite volume in-house JADIM code described in various papers, e.g. Calmet & Magnaudet (Reference Calmet and Magnaudet1997). This code makes use of a third-order Runge–Kutta/Crank–Nicolson algorithm to advance the velocity field. A Poisson equation is solved at the end of each time step $[t,t+\unicode[STIX]{x0394}t]$ to enforce the incompressibility condition, yielding the absolute fluid velocity field $\boldsymbol{U}_{QS}(t+\unicode[STIX]{x0394}t)$ and the pressure field $P_{QS}(t+\unicode[STIX]{x0394}t/2)$ , both of which are second-order accurate in time. The ‘quasi-static’ hydrodynamic force, $\boldsymbol{F}_{QS}(t+\unicode[STIX]{x0394}t/2)$ , and torque, $\unicode[STIX]{x1D71E}_{QS}(t+\unicode[STIX]{x0394}t/2)$ , acting on the sphere at time $t+\unicode[STIX]{x0394}t/2$ are straightforwardly obtained by integrating the corresponding stresses over the sphere surface. Then $\boldsymbol{V}$ and $\unicode[STIX]{x1D734}$ are updated as $\boldsymbol{V}(t+\unicode[STIX]{x0394}t)=\boldsymbol{V}(t)+\unicode[STIX]{x0394}\boldsymbol{V}(t+\unicode[STIX]{x0394}t)$ and $\unicode[STIX]{x1D734}(t+\unicode[STIX]{x0394}t)=\unicode[STIX]{x1D734}(t)+\unicode[STIX]{x0394}\unicode[STIX]{x1D734}(t+\unicode[STIX]{x0394}t)$ , respectively, by solving the appropriate form of the generalized Kirchhoff–Kelvin equations expressing Newton’s second law for a rigid body moving in an incompressible viscous fluid (Mougin & Magnaudet Reference Mougin and Magnaudet2002a ). These equations are integrated in time with a third-order Runge–Kutta algorithm that guarantees the stability of the ( $\boldsymbol{V},\unicode[STIX]{x1D734}$ ) solution. After this update, the strategy employed by Mougin & Magnaudet (Reference Mougin and Magnaudet2002b ) and Auguste et al. (Reference Auguste, Magnaudet and Fabre2013) to obtain the final fluid velocity field, $\boldsymbol{U}$ , and pressure field, $P$ , at time $t+\unicode[STIX]{x0394}t$ simply consisted in changing ( $\boldsymbol{U}_{QS},P_{QS}$ ) into ( $\boldsymbol{U}=\boldsymbol{U}_{QS}+\boldsymbol{u}^{\prime },P=P_{QS}+p^{\prime }$ ) by computing the irrotational velocity correction $\boldsymbol{u}^{\prime }$ satisfying the updated no-penetration condition at the body surface, namely

where $\boldsymbol{n}$ denotes the unit normal to the surface (the last term in the right-hand side of (2.3) vanishes in the case of a sphere since $\boldsymbol{r}=(d/2)\boldsymbol{n}$ ). However an important improvement had to be introduced to deal with light spheres for two reasons. First, the sphere moment of inertia vanishes in the limit $\overline{\unicode[STIX]{x1D70C}}\rightarrow 0$ . Second, in the inviscid limit $Ar\rightarrow \infty$ , a sphere does not displace any fluid when rotating, owing to its point-symmetric nature. Hence, an infinitely light sphere rotating in a fluid has zero overall rotational inertia and the variations of its rotation rate are governed by viscous effects. The leading contribution to these effects is due to the unsteady diffusion of vorticity within the Stokes layer created by the update of the tangential fluid velocity at the sphere surface, i.e.

As changes $\unicode[STIX]{x0394}\boldsymbol{V}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D734}$ go to zero with the time step $\unicode[STIX]{x0394}t$ , it may be shown (Mougin & Magnaudet Reference Mougin and Magnaudet2002a ) that the correction $\boldsymbol{u}^{\prime }$ obeys the unsteady Stokes equation

subject to (2.3) and (2.4) and to the initial condition $\boldsymbol{u}^{\prime }(t)=\mathbf{0}$ at the beginning of the time step $[t,t+\unicode[STIX]{x0394}t]$ . The viscous term in (2.5a ) is zero outside the Stokes layer, making the corresponding solution reduce in that outer region to the irrotational flow field computed in the original approach of Mougin & Magnaudet (Reference Mougin and Magnaudet2002a ). The pressure field associated with this irrotational contribution is responsible for the added-mass force resulting from the instantaneous acceleration of the sphere. In contrast, the solution to (2.5) carries a non-zero vorticity within the Stokes layer, the consequences of which were neglected in the above approach. The viscous stresses resulting from $(\boldsymbol{u}^{\prime },p^{\prime })$ yield leading contributions to the force and torque analogous to the so-called Basset–Boussinesq ‘history’ force. Since the sphere translational and rotational accelerations, $\unicode[STIX]{x0394}\boldsymbol{V}/\unicode[STIX]{x0394}t$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D734}/\unicode[STIX]{x0394}t$ , may be considered constant during the time interval $[t,t+\unicode[STIX]{x0394}t]$ , these contributions are of $O(\unicode[STIX]{x0394}t^{1/2})\unicode[STIX]{x0394}\boldsymbol{V}/\unicode[STIX]{x0394}t$ and $O(\unicode[STIX]{x0394}t^{1/2})\unicode[STIX]{x0394}\unicode[STIX]{x1D734}/\unicode[STIX]{x0394}t$ , respectively. Actually the ‘history’ torque also comprises a contribution of $O(\unicode[STIX]{x0394}t)\unicode[STIX]{x0394}\unicode[STIX]{x1D734}/\unicode[STIX]{x0394}t$ (Feuillebois & Lasek Reference Feuillebois and Lasek1978; Gatignol Reference Gatignol1983; Zhang & Stone Reference Zhang and Stone1998) but the latter is negligibly small compared to the former when $\unicode[STIX]{x0394}t\rightarrow 0$ . The solution to (2.5)–(2.4) also yields Stokes-like additional contributions to the force and torque, namely $-3\unicode[STIX]{x03C0}\unicode[STIX]{x1D708}d\unicode[STIX]{x0394}\boldsymbol{V}$ and $-\unicode[STIX]{x03C0}\unicode[STIX]{x1D708}d^{3}\unicode[STIX]{x0394}\unicode[STIX]{x1D734}$ , respectively. Neglecting the latter three terms (this simplification is discussed in appendix A), the Kirchhoff–Kelvin equations at time $t+\unicode[STIX]{x0394}t/2$ reduce to

where $m$ is the mass of fluid that could be enclosed in the sphere volume, and $\unicode[STIX]{x1D6FF}$ is defined as $\unicode[STIX]{x1D6FF}=\{\unicode[STIX]{x1D708}\unicode[STIX]{x0394}t/(\unicode[STIX]{x03C0}d^{2})\}^{1/2}$ . This quantity may be thought of as the dimensionless thickness of the Stokes layer resulting from the no-slip condition (2.4) after a period of time $\unicode[STIX]{x0394}t$ . The pre-factors $B_{V}$ and $B_{\unicode[STIX]{x1D6FA}}$ involved in the unsteady viscous contributions are known analytically: according to Feuillebois & Lasek (Reference Feuillebois and Lasek1978), Gatignol (Reference Gatignol1983) and Zhang & Stone (Reference Zhang and Stone1998), one has $B_{V}=18,B_{\unicode[STIX]{x1D6FA}}=2$ . Thanks to the non-zero value of $B_{\unicode[STIX]{x1D6FA}}$ , there is no longer a singularity in (2.7) when $\overline{\unicode[STIX]{x1D70C}}\rightarrow 0$ , so that the sphere dynamics can be properly computed whatever the body-to-fluid density ratio.

Once $\unicode[STIX]{x0394}\boldsymbol{V}(t+\unicode[STIX]{x0394}t)$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D734}(t+\unicode[STIX]{x0394}t)$ have been determined thanks to (2.6)–(2.7), (2.5) is solved with a two-step method. Using a Crank–Nicolson algorithm, we first compute a predictor velocity correction, $\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}(t+\unicode[STIX]{x0394}t)$ , satisfying $[\unicode[STIX]{x2202}_{t}-\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}]\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}=\mathbf{0}$ together with boundary conditions $\boldsymbol{n}\times \boldsymbol{u}^{\unicode[STIX]{x1D6FF}}=\boldsymbol{n}\times \{\unicode[STIX]{x0394}\boldsymbol{V}+\unicode[STIX]{x0394}\unicode[STIX]{x1D734}\times \boldsymbol{r}\}$ and $\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}\boldsymbol{\cdot }\boldsymbol{n}=0$ at the sphere surface and to the initial condition $\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}(t)=\mathbf{0}$ . Then we obtain the divergence-free velocity field $\boldsymbol{u}^{\prime }(t+\unicode[STIX]{x0394}t)$ through a projection step: introducing an auxiliary potential $\unicode[STIX]{x1D719}^{\prime }$ and setting $\boldsymbol{u}^{\prime }=\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}-\unicode[STIX]{x0394}t\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}^{\prime }$ , we solve the Poisson equation $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}^{\prime }=\unicode[STIX]{x0394}t^{-1}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}$ subject to the Neumann boundary condition $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}^{\prime }=-\unicode[STIX]{x0394}t^{-1}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x0394}\boldsymbol{V}$ at the sphere surface and obtain the pressure $p^{\prime }(t+\unicode[STIX]{x0394}t/2)$ as $p^{\prime }=\unicode[STIX]{x1D719}^{\prime }-\unicode[STIX]{x1D708}\unicode[STIX]{x0394}t\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}^{\prime }$ . Various tests carried out with prescribed values of $\unicode[STIX]{x0394}\boldsymbol{V}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D734}$ proved that the surface stresses resulting from the correction $(\boldsymbol{u}^{\prime },p^{\prime })$ yield the expected added-mass force $(m/2)\unicode[STIX]{x0394}\boldsymbol{V}/\unicode[STIX]{x0394}t$ in (2.6) and the ‘history’ contributions to the force and torque, $B_{V}m\unicode[STIX]{x1D6FF}\unicode[STIX]{x0394}\boldsymbol{V}/\unicode[STIX]{x0394}t$ and $B_{\unicode[STIX]{x1D6FA}}md^{2}\unicode[STIX]{x1D6FF}\unicode[STIX]{x0394}\unicode[STIX]{x1D734}/\unicode[STIX]{x0394}t$ in (2.6) and (2.7), respectively, with an excellent accuracy. In appendix A we prove that the above fluid/body coupling algorithm is second-order time accurate whatever $\overline{\unicode[STIX]{x1D70C}}$ , i.e. the truncation errors in the body translational and rotational velocities $\boldsymbol{V}$ and $\unicode[STIX]{x1D734}$ and in the velocity and pressure corrections $\boldsymbol{u}$ and $p^{\prime }$ are all of $O(\unicode[STIX]{x0394}t^{2})$ . Since the truncation error resulting from the third-order Runge–Kutta/Crank–Nicolson algorithm employed to integrate (2.1), (2.2) is also of $O(\unicode[STIX]{x0394}t^{2})$ (Calmet & Magnaudet Reference Calmet and Magnaudet1997), the fluid velocity $\boldsymbol{U}=\boldsymbol{U}_{QS}+\boldsymbol{u}^{\prime }$ and the pressure $P=P_{QS}+p^{\prime }$ are also second-order time accurate. Hence all steps involved in the present computational approach are consistent with second-order time accuracy.

Figure 1.

Paths of two spheres with the same $Ar$ but different $\overline{\unicode[STIX]{x1D70C}}$ computed with (green curve) or without (red curve) the $\unicode[STIX]{x1D6FF}$ -terms in (2.6)–(2.7) and the contribution $\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}$ in $\boldsymbol{u}^{\prime }$ . (a)  $\overline{\unicode[STIX]{x1D70C}}=0.33$ , (b) $\overline{\unicode[STIX]{x1D70C}}=0.01$ ; the horizontal ( $x$ ) and vertical ( $z$ ) axes are normalized with the sphere diameter, $d$ .

A typical example of the influence of viscous effects induced by the Stokes layer on the path evolution is shown in figure 1. The standard fluid/body coupling algorithm (Mougin & Magnaudet Reference Mougin and Magnaudet2002a ) and the modified version including the $\unicode[STIX]{x1D6FF}$ -terms in (2.6), (2.7) and the viscous contribution $\boldsymbol{u}^{\unicode[STIX]{x1D6FF}}$ in $\boldsymbol{u}^{\prime }$ produce virtually indiscernible results for the sphere with $\overline{\unicode[STIX]{x1D70C}}=0.33$ . In contrast, the two algorithms predict dramatically different evolutions in the case of the very light sphere with $\overline{\unicode[STIX]{x1D70C}}=0.01$ . This is fully in line with the analysis reported above: variations of the rotation rate in the former case are essentially governed by the sphere moment of inertia, while those in the latter case are controlled by the viscous torque associated with the unsteady diffusion of vorticity across the Stokes layer. This example clearly establishes that the path of very light spheres cannot be predicted without taking into account these subtle ‘history’ effects.

It is probably worth pointing out that the numerical strategy described in this section may easily be extended to bodies with more complex geometries, even though no closed-form expression is available for viscous (or even inviscid) effects due to the body acceleration. For this purpose, before running the complete problem, one simply needs to impose successively prescribed uniform translational or rotational accelerations to the body in each direction of space, record the evolutions of all force and torque components during a short time period and analyse them by using the methodology described in Rivero, Magnaudet & Fabre (Reference Rivero, Magnaudet and Fabre1991). This provides once and for all the whole set of coefficients involved in the counterpart of (2.6), (2.7).

2.3 Computational protocol

The fluid problem is spatially discretized on a spherical grid. Extensive tests were carried out to determine the grid characteristics, the goal being to guarantee grid-independent paths (which requires grid-independent loads on the sphere) while keeping the computational cost acceptable. The loads on the sphere are especially sensitive to the resolution of both the boundary layer and the near wake, and to the disturbances that may be generated by vortices leaving the computational domain through its outer boundary. Moreover, as the sphere freely moves, the wake may take any orientation with respect to the grid axis, which makes effects of the anisotropy inherent to the use of an axisymmetric grid of particular importance in the present problem. Last, since instabilities may act on short time scales, controlling the time step also requires specific attention. For these reasons, we very carefully considered the above five aspects. In appendix B, we detail the determination of the various grid parameters, discuss the condition that controls the time step, and summarize the results of several test cases.

Simulations to be described below were generally started from rest; nevertheless, to examine the influence of initial conditions, we also frequently ran computations for a given set of parameters by using fields $(\boldsymbol{U},P,\boldsymbol{V},\unicode[STIX]{x1D734})$ from a previous run corresponding to neighbouring values of $Ar$ and $\overline{\unicode[STIX]{x1D70C}}$ . When starting from rest, a small disturbance was added during several time steps in the right-hand side of (2.6) to trigger the instability. Similarly, the stability of planar non-vertical paths was examined by introducing momentarily a small disturbance in the right-hand side of (2.6) and/or (2.7). In each run, the sphere was followed over a vertical distance ranging typically from $250d$ to $700d$ . To decide when a run was stopped, we combined different criteria, depending on the nature of the regime under consideration. Close to one of the first bifurcations of the system, we monitored the growth rate of the path deviations and stopped the computation once saturation was reached. In most cases we rather considered the various components of the force and torque acting on the sphere, and examined their variations. In particular, in periodic or quasi-periodic regimes, we stopped the computation when these changes were negligibly small (i.e. of the order of the computational accuracy expected on these dynamical quantities) over typically 10 periods. With the grid defined above, the CPU time of a typical run ranged from $500$ to $1000$ hours in sequential mode on a single-core PC.

3.1 Preliminary comment

Before we discuss computational results, it is of interest to remark that some physical insight can be obtained directly from (2.6), (2.7) regarding the influence of the solid-to-fluid density ratio on the rise regimes. Indeed, the left-hand side of these equations involves three different types of contributions: those related to the body inertia, weighted by $\overline{\unicode[STIX]{x1D70C}}$ , those due to the fluid inertia, which in the present case reduce to the added-mass coefficient $1/2$ in (2.6), and those resulting from the permanent existence of a Stokes layer around the sphere, due to the continual variations of its translational and rotational velocities. In the framework of the algorithm described above, the magnitude of the latter scales with the parameter $\unicode[STIX]{x1D6FF}$ because the characteristic time involved in the time-advancement procedure is $\unicode[STIX]{x0394}t$ . However, from a physical point of view, the relevant time scale is the gravitational time $t_{g}=d/V_{g}=\sqrt{d/((1-\overline{\unicode[STIX]{x1D70C}})g)}$ . Hence the influence of the Stokes layer on the body dynamics is measured by the dimensionless parameter $\unicode[STIX]{x0394}=(\unicode[STIX]{x1D708}t_{g}/d^{2})^{1/2}=Ar^{-1/2}$ . Based on the comparison of the three characteristic orders of magnitude $\overline{\unicode[STIX]{x1D70C}},1/2$ and $Ar^{-1/2}$ , one may anticipate three different classes of dynamical regimes, keeping in mind that $Ar^{-1/2}\ll 1$ : those of ‘heavy’ spheres corresponding to $\overline{\unicode[STIX]{x1D70C}}\gg 1/2$ (here presumably restricted to $\overline{\unicode[STIX]{x1D70C}}\lesssim 1$ ), those of spheres with intermediate density ratios such that $(1/2)\gg \overline{\unicode[STIX]{x1D70C}}\gg Ar^{-1/2}$ , and those of very light spheres with $\overline{\unicode[STIX]{x1D70C}}\ll Ar^{-1/2}$ (here $0.04\leqslant Ar^{-1/2}\leqslant 0.08$ ). In the first regime, possible changes in the translational and rotational sphere velocities are entirely driven by quasi-steady hydrodynamic stresses and body inertia. In contrast, in the second case, fluid acceleration plays a dominant role in the translational dynamics, while changes in the rotation rate are still controlled by the sphere moment of inertia. Last, in the regime $\overline{\unicode[STIX]{x1D70C}}\ll Ar^{-1/2}$ , the fluid also controls the rotational dynamics through the viscous processes involved in the Stokes layer. Of course intermediate regimes are expected to take place for $\overline{\unicode[STIX]{x1D70C}}\approx 1/2$ and $\overline{\unicode[STIX]{x1D70C}}\approx Ar^{-1/2}$ , where the body and fluid both drive the translational and rotational dynamics, respectively. The findings described in the rest of this paper will confirm the above analysis, especially the strikingly distinct behaviour of ‘heavy’ and ‘very light’ spheres in several ranges of $Ar$ . Note that the dividing density ratio $\overline{\unicode[STIX]{x1D70C}}=1$ separating falling and rising spheres does not play any role in the above classification. Hence spheres with $\overline{\unicode[STIX]{x1D70C}}\gtrsim 1$ (which are not considered here) and with $\overline{\unicode[STIX]{x1D70C}}\lesssim 1$ are expected to behave similarly. This is indeed what was observed by ZD. Obviously the above classification only makes sense in unsteady regimes, since the sphere inertia does not play any role as far as the translational and rotational velocities do not change over time.

Figure 2.

State diagram of the first rise regimes in the range $140\leqslant Ar\leqslant 180$ . Grey: steady vertical (SV); pale orange: steady oblique (SO); dark orange: periodic oblique (PO); purple: intermittent (IT); green: large-amplitude zigzagging (ZZ).

3.2 First non-vertical regimes

Figure 2 gathers the various regimes of rise observed up to $Ar=180$ . Starting from left, i.e. from low $Ar$ , the first encountered regime corresponds of course to a steady vertical (SV) rise associated with an axisymmetric flow field. Increasing $Ar$ , this base state is succeeded by a steady oblique (SO) regime which was first identified numerically by JDB before it was unambiguously observed experimentally by Veldhuis & Biesheuvel (Reference Veldhuis and Biesheuvel2007) and Horowitz & Williamson (Reference Horowitz and Williamson2010). Somewhat later, the SV–SO transition was shown to correspond to a stationary (pitchfork) bifurcation by Fabre, Tchoufag & Magnaudet (Reference Fabre, Tchoufag and Magnaudet2012). The same study determined that the corresponding threshold takes place at $Ar=155.6$ , irrespective of the sphere-to-fluid density ratio. The recent computations reported in ZD confirmed this threshold and its remarkable independence with respect to $\overline{\unicode[STIX]{x1D70C}}$ which is also in line with the qualitative analysis developed above. Here, since we are essentially interested in time-dependent regimes, we did not attempt to recover precisely this threshold throughout the whole range of $\overline{\unicode[STIX]{x1D70C}}$ ; hence in figure 2, the vertical line separating the SV and SO regimes is merely based on the prediction of Fabre et al. (Reference Fabre, Tchoufag and Magnaudet2012). Nevertheless some specific runs allowed us to observe that the SV–SO threshold stands in the range $155<Ar<156$ for small $\overline{\unicode[STIX]{x1D70C}}$ as well as for $\overline{\unicode[STIX]{x1D70C}}=0.99$ , in line with this prediction; we also checked that the system is still in the SV regime for $Ar=150$ whatever $\overline{\unicode[STIX]{x1D70C}}$ , whereas it is always in the SO regime when $Ar=162.5$ . Moreover, as figure 3 reveals, the sphere rotation rate and drift angle are still independent of $\overline{\unicode[STIX]{x1D70C}}$ for $Ar=162.5$ (i.e. they only depend on the distance to the threshold), the drift angle with respect to the $z$ -axis being approximately $4.3^{\circ }$ (throughout the paper, $z$ denotes the vertical ascending coordinate and $x$ and $y$ lie in the horizontal plane; all three coordinates are normalized by $d$ ). On the same figure, one may also notice that the vertical position at which the SV–SO transition takes place slightly increases as $\overline{\unicode[STIX]{x1D70C}}$ decreases (compare the green and orange curves). However, one has to keep in mind that the terminal rise velocity resulting from the balance between drag and buoyancy forces roughly varies as $(1-\overline{\unicode[STIX]{x1D70C}})^{1/2}$ since the flow is dominated by inertia effects ( $Ar\gg 1$ ). Hence, the time at which the transition takes place is actually approximately one order of magnitude smaller for the lightest spheres than for the heaviest ones. This is consistent with the fact that the smaller the sphere inertia, the larger the rate of change of the sphere velocity and rotation rate in response to a given force or torque disturbance, hence the faster the transition.

Figure 3.

Paths for various solid-to-fluid density ratios at $Ar=162.5$ (SO regime). From top to bottom: $\overline{\unicode[STIX]{x1D70C}}=0.001$ (red, almost hidden under the green curve), $0.01$ (green), $0.05$ (blue), $0.1$ (purple), $0.25$ (cyan), $0.5$ (yellow), $0.75$ (black), $0.99$ (orange).

Increasing $Ar$ , the SO regime is quickly replaced by a periodic oscillating (PO) regime. The threshold of the corresponding transition increases with $\overline{\unicode[STIX]{x1D70C}}$ , from $164.4$ for $\overline{\unicode[STIX]{x1D70C}}\approx 0$ to $173.3$ for $\overline{\unicode[STIX]{x1D70C}}\approx 1$ . These values are slightly smaller than those reported by ZD (167.2 and 178.5, respectively). The PO regime is only observed up to $Ar=167.3$ when $\overline{\unicode[STIX]{x1D70C}}\approx 0$ ; this upper limit increases with $\overline{\unicode[STIX]{x1D70C}}$ (mostly for $\overline{\unicode[STIX]{x1D70C}}>0.5$ ), up to $Ar=187.7$ when $\overline{\unicode[STIX]{x1D70C}}\approx 1$ . Whatever $\overline{\unicode[STIX]{x1D70C}}$ , the dimensionless frequency (or Strouhal number) is approximately 0.045 at the SO–PO threshold and decreases with the distance to the threshold. This finding is in good agreement with the results of ZD that indicate values of $St$ close to 0.05 at the SO–PO threshold throughout the density ratio range $0\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 2$ . Compared with usual Strouhal numbers associated with vortex shedding in bluff-body wakes, the values involved here are typically 3–4 times smaller, revealing a low-frequency oscillation. To better understand this feature, it is important to realize that the sphere average Reynolds number, $Re_{m}$ , is in the range [220, 255] throughout the range of $Ar$ within which the PO regime is observed. Hence, no vortex shedding would take place in the wake, were the sphere forced to rise in straight line, since the onset of vortex shedding past a sphere held fixed in a uniform stream corresponds to a critical Reynolds number of approximately 275 (Ghidersa & Dušek Reference Ghidersa and Dušek2000). Therefore the observed low-frequency path oscillation unambiguously results from the couplings of the sphere dynamics and the surrounding flow. This is similar to what is for instance observed with disks which, in some ranges of $\overline{\unicode[STIX]{x1D70C}}$ , start fluttering at a critical Reynolds number much lower than the threshold corresponding to the loss of axisymmetry in the wake of a disk held fixed in a uniform flow (Auguste et al. Reference Auguste, Magnaudet and Fabre2013). Figure 4 shows a typical PO path, together with the associated wake evolution observed through the $\unicode[STIX]{x1D706}_{2}$ criterion (Jeong & Hussain Reference Jeong and Hussain1995). The path inclination oscillates from $2.6^{\circ }$ to $7.3^{\circ }$ , with an average value of approximately $4.5^{\circ }$ . As the correspondence between the path and wake reveals, the streamwise vortex pair is virtually absent when the path is almost vertical (stage 1), and the larger the inclination the more intense this vortex pair, with a maximum during stages 2 and 4. This is no surprise since streamwise vortices are responsible for the lift force acting on the sphere, which is the cause of the bending of its path. Similar to the SO regime, the sign of the streamwise vorticity remains constant within each thread, which yields the non-zero average horizontal drift of the trajectory.

Figure 4.

Characteristics of the PO regime for $\overline{\unicode[STIX]{x1D70C}}=0.5$ and $Ar=172$ . (a) Sphere trajectory; (b) wake evolution during about $3/4$ of a period; the wake is visualized using the $\unicode[STIX]{x1D706}_{2}$ criterion (Jeong & Hussain Reference Jeong and Hussain1995), with iso- $\unicode[STIX]{x1D706}_{2}$ surfaces coloured by the local value of the vertical component, $\unicode[STIX]{x1D714}_{z}$ , of the vorticity normalized by $V_{g}/d$ ( $-1.94\leqslant \unicode[STIX]{x1D714}_{z}\leqslant 1.94$ , $\unicode[STIX]{x1D714}_{z}\approx 0$ in green regions).

Figure 5.

Paths for various solid-to-fluid density ratios at $Ar=175$ . Same colour code as in figure 3.

Throughout the whole range of $\overline{\unicode[STIX]{x1D70C}}$ , the PO regime is succeeded by an intermittent regime within which the sphere ‘hesitates’ more or less often to perform periodic planar zigzags. That is, it drifts laterally for a while until its path becomes nearly vertical. At this point, it either starts to drift in the opposite direction, then performing a full zigzag, or does not succeed to do so, then keeping on drifting in the same direction, possibly in a different vertical plane. Indeed, since ‘hesitations’ take place during the time lapse when the wake is almost axisymmetric (as in stage 1 of figure 4), the system is momentarily free to select a new symmetry plane. This succession of ‘successes’ and ‘failures’ in performing complete zigzags occurs randomly. Such rare events are known to be characteristic of intermittency in dynamical systems and transition to turbulence (Manneville & Pomeau Reference Manneville and Pomeau1979; Pomeau & Manneville Reference Pomeau and Manneville1980; Bergé, Pomeau & Vidal Reference Bergé, Pomeau and Vidal1984). Hence this new regime, which we refer to as IT, marks the onset of chaos in the present dynamical system whatever $\overline{\unicode[STIX]{x1D70C}}$ (its chaotic nature will be established later). The IT regime covers only a narrow parameter range, from 2 to 10 $Ar$ -units depending on $\overline{\unicode[STIX]{x1D70C}}$ , the narrowest part being found in the range $0.3\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 0.75$ . The styles of path observed in this intermittent regime and the influence of the sphere-to-fluid density ratio on the PO–IT transition may be appreciated in figure 5 which corresponds to $Ar=175$ . Starting from $\overline{\unicode[STIX]{x1D70C}}\approx 1$ , the two heaviest spheres are still in the PO regime, the amplitude of oscillations increasing when $\overline{\unicode[STIX]{x1D70C}}$ decreases. For $\overline{\unicode[STIX]{x1D70C}}=0.5$ , the sphere first describes several complete zigzags until its horizontal velocity fails to reverse. Then it starts to drift along an inclined path reminiscent of those observed in the PO regime. In the range $0.25\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 0.05$ , the paths consist essentially of periodic zigzags. However a ‘failure’ still occurs after the sphere has travelled a vertical distance of $O(350d)$ for $\overline{\unicode[STIX]{x1D70C}}=0.1$ ; no such feature is observed for $\overline{\unicode[STIX]{x1D70C}}=0.25$ and $0.05$ throughout the entire ascent, although there is no guarantee that it would not occur at a higher vertical position; one may also notice that a significant average drift exists in the latter case. The lightest two spheres exhibit a somewhat different behaviour: throughout their ascent, their paths follow complete zigzags, but in the late stage one of these zigzags is asymmetric, i.e. the sphere spends a longer time drifting in one direction than in the other, making the average path suddenly exhibit a slight lateral shift. This may be interpreted as a weak form of the above defect: the horizontal velocity does not change sign in due time to maintain the strict periodicity of the zigzag; however, it succeeds to do so with some limited delay, keeping the appearance of the path close to that of a pure zigzag (as will be shown in figure 10, slightly increasing the Archimedes number up to $Ar=178$ yields a strictly periodic zigzagging path for $\overline{\unicode[STIX]{x1D70C}}=0.01$ ). Overall, what can be concluded from this sequence is that the smaller $\overline{\unicode[STIX]{x1D70C}}$ , the less intermittent the sphere dynamics at this particular value of $Ar$ . In contrast, as soon as the sphere trajectory is essentially made of zigzags, the amplitude of the latter is almost independent of $\overline{\unicode[STIX]{x1D70C}}$ and is about one sphere diameter.

3.3 Periodic and intermittent regimes up to $Ar=300$

As figure 6(a) shows, the situation changes significantly when $Ar$ is increased up to $187.5$ . The heaviest sphere is still in the PO regime with an average drift angle similar to that previously found ( $4.2^{\circ }$ ); however path oscillations are now highly asymmetric. All spheres corresponding to $0.1\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 0.75$ exhibit regular planar zigzagging paths, demonstrating the existence of a well-defined periodic ZZ (zigzagging) regime. The amplitude of the corresponding lateral excursions is still of $O(d)$ and slightly increases when $\overline{\unicode[STIX]{x1D70C}}$ decreases; one may also notice a weak average drift showing no clear trend with $\overline{\unicode[STIX]{x1D70C}}$ (this drift is almost zero for $\overline{\unicode[STIX]{x1D70C}}=0.5$ ). As $\overline{\unicode[STIX]{x1D70C}}$ is decreased, the zigzags become more and more nonlinear, higher harmonics contributing to make their extremities sharper. For $\overline{\unicode[STIX]{x1D70C}}=0.1$ , a slight lateral shift of the uppermost zigzag may be discerned, suggesting that the corresponding conditions are close to the transition between a strictly periodic regime and an intermittent one. This is confirmed by the evolution of the paths observed with the lightest three spheres which exhibit erratic, albeit still planar, paths (the ratio of the out-of-plane to in-plane lateral displacements is ${\leqslant}5\times 10^{-5}$ ). The various styles of path encountered in figure 6(a) provide a clear illustration of the three distinct behaviours predicted on the basis of the qualitative analysis presented in § 3.1, with one PO path for $\overline{\unicode[STIX]{x1D70C}}\lesssim 1$ , several intermittent paths for $\overline{\unicode[STIX]{x1D70C}}\lesssim Ar^{-1/2}$ (here $Ar^{-1/2}=0.073$ ), and periodic ZZ paths in between.

Figure 6.

Paths for various solid-to-fluid density ratios at (a) $Ar=187.5$ and (b) $Ar=200$ : $\overline{\unicode[STIX]{x1D70C}}=0.001$ (red), $0.01$ (green), $0.05$ (blue), $0.1$ (purple), $0.175$ (grey), $0.25$ (cyan), $0.5$ (yellow), $0.75$ (black), $0.99$ (orange).

Increasing further $Ar$ up to $200$ leads to a simpler situation, with essentially two styles of periodic paths (figure 6 b). Spheres with $\overline{\unicode[STIX]{x1D70C}}\geqslant 0.175$ are found to follow well-defined zigzagging trajectories; their characteristics and variations with $\overline{\unicode[STIX]{x1D70C}}$ are similar to those described above. In contrast, very light spheres with $\overline{\unicode[STIX]{x1D70C}}\leqslant 0.05$ follow a new type of periodic path combining characteristic features of PO and ZZ paths. More precisely, the corresponding trajectories exhibit a clear lateral drift (about $3^{\circ }$ for the lightest sphere, then decreasing with increasing $\overline{\unicode[STIX]{x1D70C}}$ ) but the oscillations are much more pronounced than in the PO regime, the horizontal velocity being able to change sign twice during each period, as in the ZZ regime. These combined characteristics make this regime specific and we refer to it as the ZO regime. The corresponding trajectories are planar, with maximum local deviation angles within horizontal planes of approximately $0.05^{\circ }$ throughout the sphere ascent. The Strouhal number stands in the range 0.04–0.045, slightly increasing with both $\overline{\unicode[STIX]{x1D70C}}$ and $Ar$ . A periodic zigzagging path with a non-zero mean drift was also observed in ZD in the same range of parameters, i.e. small $\overline{\unicode[STIX]{x1D70C}}$ and $Ar=O(200)$ (see their figure 13). However the corresponding path exhibits two crucial differences with those characterizing the present ZO regime: the corresponding Strouhal number is approximately $0.1$ and the drift is perpendicular to the plane of the zigzagging motion, making the sphere rise in a non-vertical plane. No such feature was detected here in this region of the $(Ar,\overline{\unicode[STIX]{x1D70C}})$ plane, implying that the two codes provide qualitatively different predictions in that range of parameters.

Figure 7.

Characteristics of the ZZ regime for $\overline{\unicode[STIX]{x1D70C}}=0.5$ and $Ar=212.5$ . (a) Sphere trajectory; (b) wake evolution during $3/4$ of a period visualized with the $\unicode[STIX]{x1D706}_{2}$ criterion ( $-5.0\leqslant \unicode[STIX]{x1D714}_{z}\leqslant 5.0$ , see figure 4 for caption).

Figures 7 and 8, taken at the slightly larger Archimedes number $Ar=212.5$ , display a zoom of the typical sphere trajectory and the correspondence between path oscillations and wake dynamics in the ZZ and ZO regimes, respectively. As could be anticipated from figure 6, the sphere dynamics in the ZZ regime generally involves several harmonics of the fundamental frequency, so that the path is not a pure sinusoid; as shown by JDB, symmetry conditions imply that only odd harmonics are present. The simulations indicate that the dominant frequency, which governs the reversal of the horizontal velocity, is an increasing function of both $Ar$ and $\overline{\unicode[STIX]{x1D70C}}$ , with $St$ ranging from $0.016$ to $0.036$ ( $St\approx 0.02$ in the case displayed in figure 7). As shown in figure 7(b), two distinct shedding sequences are involved. That is, two pairs of counter-rotating vortices with opposite signs of the streamwise vorticity are successively shed in the wake during a zigzag half-period but one pair is much stronger than the other, keeping the path close to a sinusoid. The generation of the dominant pair starts just after the sphere has left one extremity of the zigzag (stages 1 and 3 in figure 7). This pair grows gradually until it forms a horseshoe vortex that detaches from the sphere approximately when it passes through the inflection point, i.e. midway between stages 1–2 or 3–4 (the remnants of this primary horseshoe vortex are still visible in the bottom part of stages 2 and 4 in figure 7(b)). The secondary pair, in which the streamwise vorticity within a given thread has changed sign compared to the primary pair, starts growing right after this shedding. It reaches its maximum magnitude approximately between the inflection point and the next extremity of the zigzag, i.e. in stages 2 and 4, and detaches from the sphere short before the latter reaches that extremity. A similar sequence was first observed experimentally by Horowitz & Williamson (Reference Horowitz and Williamson2010) with freely rising spheres and was referred to as the ‘4R’ shedding mode. It was also detected with nearly spheroidal rising bubbles in the recent computational study of Cano-Lozano et al. (Reference Cano-Lozano, Martínez-Bazán, Magnaudet and Tchoufag2016). In both cases, the body was found to follow a zigzagging path, so that this shedding mode appears to be intimately associated with this specific path geometry. However, for other sets of control parameters but still in the zigzagging regime, a ‘2R’ shedding mode with only one horseshoe vortex shed during each half-period of the path was also found in both studies. The mechanism that selects the ‘4R’ mode appears to be connected to the existence of oscillations in the rise velocity: when these oscillations (with frequency $2St$ , owing to obvious symmetry constraints) have a significant magnitude (typically 5 %), the ‘4R’ mode emerges, whereas the usual ‘2R’ mode is observed when they remain small (typically 1 % or less).

Figure 8.

Characteristics of the ZO regime for $\overline{\unicode[STIX]{x1D70C}}=0.01$ and $Ar=212.5$ . (a) Sphere trajectory; (b) wake evolution during $3/4$ of a period visualized with the $\unicode[STIX]{x1D706}_{2}$ criterion ( $-4.34\leqslant \unicode[STIX]{x1D714}_{z}\leqslant 4.34$ , see figure 4 for caption).

Figure 8 displays the same information for the ZO path observed, still for $Ar=212.5$ , when the sphere-to-fluid density ratio is decreased to $\overline{\unicode[STIX]{x1D70C}}=0.01$ . The zigzags are now strongly asymmetric and the average path is inclined by approximately $3.6^{\circ }$ with respect to the vertical. The Strouhal number is approximately $0.046$ , which is significantly larger than in the ZZ regime. Again, the evolution of the corresponding wake structure reveals that four pairs of counter-rotating vortices are shed downstream during a period. However, in contrast to ZZ paths, the sign of the streamwise vorticity changes every time a new pair appears. Hence, considering the thread located on a given half-space with respect to the plane of rise, the successive stages 1–4 in figure 7 correspond to, say, a $+--+$ sequence, whereas those in figure 8 correspond to a $+-+-$ one. Moreover, one pair (instead of two in the ZZ regime) is much stronger than the other three (respectively two) and is responsible for the lateral drift (this pair is the one that starts growing when the horizontal velocity reverses almost symmetrically, i.e. in stage 1 of figure 8).

Figure 9.

Paths for various solid-to-fluid density ratios at $Ar=237.5$ ; (a) projection onto the vertical $(x,z)$ plane, (b) projection onto the horizontal plane. Same colour code as in figure 6.

Figure 10.

Velocity diagrams built on the vertical $(V_{v})$ and horizontal $(V_{h})$ sphere velocities (both normalized with $V_{g}$ ), obtained for spheres with $\overline{\unicode[STIX]{x1D70C}}=0.01$ in the various regimes encountered in the range $150\leqslant Ar\leqslant 237.5$ . Diagrams corresponding to the three successive intermittent (IT) regimes observed for $Ar=172$ (a), $187.5$ (b) and $237.5$ (c) are compared with those found in the successive stationary and periodic regimes, namely SV ( $Ar=150$ ), SO ( $Ar=162.5$ ), PO ( $Ar=166$ ), ZZ ( $Ar=178$ ) and ZO ( $Ar=200$ ). The colour scale refers to that used in the state diagrams of figures 2 and 13; a horizontal flip has been applied to the ZO curve to make the figure clearer.

We observe the ZO regime typically in the range $200\leqslant Ar\leqslant 250$ for density ratios that may be up to $0.175$ . The corresponding ‘island’ in the $(Ar,\overline{\unicode[STIX]{x1D70C}})$ diagram (see figure 13 below) has a limited area and looks like a thin stripe, the right edge of which makes a positive angle with respect to the $Ar$ -axis. Indeed, increasing $Ar$ within the ZO regime, the path of the lightest sphere ( $\overline{\unicode[STIX]{x1D70C}}=0.001$ ) starts to exhibit some intermittency for $Ar=212.5$ after the sphere has risen a vertical distance of $O(150d)$ . Other very light spheres switch similarly from a ZO path to an intermittent path at a critical $Ar$ increasing with $\overline{\unicode[STIX]{x1D70C}}$ . Thus the lightest sphere still exhibiting a periodic ZO path when $Ar=237.5$ corresponds to $\overline{\unicode[STIX]{x1D70C}}=0.05$ . At the same time, the maximum $\overline{\unicode[STIX]{x1D70C}}$ at which a ZO path is observed also increases with $Ar$ , from $0.05$ for $Ar=200$ to $0.175$ for $Ar=250$ . Consequently, as figure 9(a) shows, only spheres with $\overline{\unicode[STIX]{x1D70C}}=0.05$ and $\overline{\unicode[STIX]{x1D70C}}=0.1$ still exhibit a periodic ZO path at $Ar=237.5$ . Figure 9(b) indicates that these paths are still planar. The same conclusion holds for the intermittent paths observed with the lightest two spheres, as well as for the periodic ZZ paths of the two heaviest ones. In contrast, the two intermediate spheres corresponding to $\overline{\unicode[STIX]{x1D70C}}=0.25$ and $\overline{\unicode[STIX]{x1D70C}}=0.5$ develop a strikingly different behaviour. They both follow an intermittent large-amplitude zigzagging path. However, every time a ‘hesitation’ of the zigzag happens, a new symmetry plane is selected for the next, almost periodic, zigzagging sequence. As a result, these paths lose their planarity after the first defect, i.e. for $z\gtrsim 100d$ and $150d$ , respectively. At lower values of $Ar$ , the maximum out-of-plane path deviations we measured were always less than $10^{-2}d$ , and most frequently less than $10^{-5}d$ . Here in contrast the ratio of the maximum horizontal excursions in two perpendicular directions is of $O(1)$ . Hence $Ar=237.5$ marks the threshold beyond which intermittent paths start to exhibit a clear three-dimensionality. Interestingly, since $V_{m}\approx 1.45V_{g}$ at this value of $Ar$ , the average Reynolds number $Re_{m}$ is close to $350$ , which is approximately the critical value beyond which three-dimensionality is known to exist in the wake of a sphere held fixed in a uniform stream (Mittal Reference Mittal1999).

Figure 11.

Paths for various solid-to-fluid density ratios at $Ar=275$ ; (a) projection onto the vertical $(x,z)$ plane, (b) projection onto the horizontal plane (with a zoom for $\overline{\unicode[STIX]{x1D70C}}=0.5,0.75$ and $0.99$ ). Same colour code as in figure 6.

The chaotic or periodic nature of the various regimes encountered with very light spheres ( $\overline{\unicode[STIX]{x1D70C}}=0.01$ ) up to $Ar=237.5$ is made apparent in figure 10, using phase portraits built on the two components of the body velocity (the transients corresponding to the initial vertical rise were removed so as to focus on the ulterior evolution). Increasing the Archimedes number from $Ar=150$ , the system successively stands (see figures 2 and 13) in the SV, SO, PO, IT, ZZ, IT, ZO and IT regimes (in our nomenclature, the IT label gathers all intermittent regimes, even though they sometimes occupy separate areas in the $(Ar,\,\overline{\unicode[STIX]{x1D70C}})$ plane). Curves associated with the three successive intermittent regimes tend to fill densely a region of the $(V_{v},V_{h})$ plane, demonstrating the chaotic nature of the corresponding paths. The rate at which this filling operates clearly increases with $Ar$ (the three graphs were obtained by considering the velocity evolution over a similar height), highlighting the fact that the system becomes more and more chaotic as $Ar$ increases. The phase portrait corresponding to the ZZ regime displays a slight asymmetry in the region corresponding to the maxima of the lateral sphere excursion ( $V_{h}\approx 0$ ), underlining the influence of higher harmonics on the sphere kinematics. In the ZO regime, a long transient precedes the fully periodic evolution. Then, unlike the Lissajous-type curve characterizing the ZZ regime, the phase portrait consists of a single, highly asymmetric, loop. This is because the vertical velocity takes two different values at the two locations where the horizontal velocity changes sign, which is not unlikely given the strong asymmetry of the path around points 1 and 3 in figure 8. Note also the sharp bump associated with the maximum of $V_{v}$ , which corresponds to the region in between points 3 and 4 in figure 8 where the path remains almost vertical during a significant time lapse.

Figure 12.

Characteristics of the ZZ $_{2}$ regime for $\overline{\unicode[STIX]{x1D70C}}=0.5$ and $Ar=300$ . (a) Sphere trajectory; (b) wake evolution during $3/4$ of a period visualized with the $\unicode[STIX]{x1D706}_{2}$ criterion ( $-4.73\leqslant \unicode[STIX]{x1D714}_{z}\leqslant 4.73$ , see figure 4 for caption).

In the range $250\lesssim Ar\lesssim 300$ , a new regime is observed for large enough solid-to-fluid density ratios, mostly $\overline{\unicode[STIX]{x1D70C}}\gtrsim 0.5$ . Figure 11 provides typical examples of that regime by displaying the various paths for $Ar=275$ . All spheres with $\overline{\unicode[STIX]{x1D70C}}<0.5$ follow intermittent paths, some of them being fully three-dimensional (see below for the definition of the corresponding IT regime). Nevertheless, in the range $0.05\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 0.25$ , long sequences reminiscent of the ZO regime are still noticed. As figure 11(b) shows, several of the corresponding trajectories are strongly three-dimensional but some others remain planar. The three heaviest spheres, corresponding to $\overline{\unicode[STIX]{x1D70C}}\geqslant 0.5$ , start by following a well-defined large-amplitude planar zigzagging trajectory. However, this path becomes unstable at some point and the zigzags disappear. After a transitional stage, a new style of path emerges in the form of another zigzagging motion with characteristics distinctly different from those of the ZZ regime. This new regime, which we refer to as ZZ $_{2}$ , is characterized by small-amplitude high-frequency lateral sphere displacements: in the example of figure 11, the crest-to-crest amplitude is typically approximately 20 % of the sphere diameter and the Strouhal number is close to $0.1$ . For other $(Ar,\,\overline{\unicode[STIX]{x1D70C}})$ sets, the Strouhal number stays in the range 0.10–0.11 and the crest-to-crest amplitude is always less than $0.3d$ . Hence, compared to the ZZ regime, the amplitude of the zigzags is typically one order of magnitude smaller, whereas their frequency is three to five times larger. Several additional features of the ZZ $_{2}$ regime are made apparent in figure 11. First, unlike the ZZ and ZO paths, the ZZ $_{2}$ path is found to emerge after a (very) long transient (up to ${\approx}330d$ for $\overline{\unicode[STIX]{x1D70C}}=0.75$ ). Next, it may be noticed that in this regime the sphere drifts slowly but consistently throughout its ascent, which implies that its wake delivers a non-zero average lift force. Last, as the inset in figure 11(b) reveals, ZZ $_{2}$ paths are generally non-planar because the plane within which the zigzags take place slowly rotates.

Figure 13.

State diagram of rise regimes in the intermediate range $180\leqslant Ar\leqslant 320$ . Dark orange: periodic oblique (PO); purple: intermittent (IT); green: large-amplitude zigzagging (ZZ); red: zigzagging oblique (ZO); dark green: small-amplitude zigzagging (ZZ $_{2}$ ); blue: three-dimensional chaotic (3DC).

Figure 12 shows the correspondence between path oscillations and wake dynamics in the case $\overline{\unicode[STIX]{x1D70C}}=0.5,Ar=300$ ; for this set of parameters, the amplitude of the out-of-plane horizontal motions is less than 4 % of that of in-plane excursions so that the path is nearly two-dimensional. The streamwise vorticity in a given trailing vortex is seen to keep the same sign in stages 1–4 and 2–3, respectively. Hence this quantity only changes sign every time the sphere crosses the inflection point of the zigzag, so that only two vortex pairs are released downstream during a ZZ $_{2}$ period, instead of four in the ZZ regime. It must be noticed that ZD also observed a rise/fall regime characterized by small-amplitude high-frequency oscillations with $St\approx 0.1$ for $250\lesssim Ar\lesssim 300$ and $0.5\lesssim \overline{\unicode[STIX]{x1D70C}}\lesssim 2.0$ . Hence, for rising spheres, there is a perfect match between the domains where regimes involving such fast oscillations are observed in both studies. However, ZD found the corresponding paths to be strictly planar and vertical on average, which led them to term this regime ‘vertical oscillating’. Hence, here again, some important structural characteristics of the zigzagging regime observed in the present study differ from those reported by ZD.

Figure 13 summarizes the succession of regimes encountered in the range $180\leqslant Ar\leqslant 300$ . It underlines the particularly strong influence of the sphere-to-fluid density ratio on the style of path in that range. The extent of the planar large-amplitude ZZ regime is quite remarkable. In particular, it is found to dominate the scene for spheres with $\overline{\unicode[STIX]{x1D70C}}\gtrsim 0.1$ throughout the range $195\lesssim Ar\lesssim 220$ . This clearly differs from the findings of ZD who did not detect any large-amplitude zigzag for sphere-to-fluid density ratios larger than $0.5$ , nor for $Ar>200$ . Up to $Ar\approx 260$ , most situations for very light spheres, say $\overline{\unicode[STIX]{x1D70C}}\lesssim 0.1$ , correspond to intermittent regimes; note also the significant tongue of IT regimes in the range $220\lesssim Ar\lesssim 250$ which separates the ZZ and ZZ $_{2}$ domains up to $\overline{\unicode[STIX]{x1D70C}}\approx 0.6$ . In all cases, these intermittent regimes are associated with planar paths below $Ar\approx 235$ and with three-dimensional paths beyond that limit. Nevertheless, we do not classify the latter as fully three-dimensional chaotic (3DC), although both are of the same nature from a dynamical system viewpoint, because the observed three-dimensional intermittent paths still essentially exhibit the ZZ style, with only one or two defects throughout the entire path (see figure 9). This trend stops beyond a certain $Ar$ , and this is where we located the qualitative IT/3DC border in figure 13. We did the same with the ‘almost ZO’ intermittent paths found for $\overline{\unicode[STIX]{x1D70C}}\lesssim 0.175$ beyond the ZO regime.

Figure 14.

Two three-dimensional paths in the 3DC regime at $Ar=350$ . (a) $\overline{\unicode[STIX]{x1D70C}}=0.99$ , (b) $\overline{\unicode[STIX]{x1D70C}}=0.05$ . The horizontal projection of the path is seen in the plane $z=0$ .

3.4 Beyond $Ar=300$

The ZZ $_{2}$ regime ends in the range $260\lesssim Ar\lesssim 300$ , depending on $\overline{\unicode[STIX]{x1D70C}}$ . It is succeeded by the regime we refer to as 3DC, characterized by chaotic, fully three-dimensional paths in which only few traces of the previous periodic regime subsist. For $\overline{\unicode[STIX]{x1D70C}}\lesssim 0.3$ , the same happens beyond the three-dimensional intermittent regime defined above.

Figure 14 displays two trajectories characteristic of the 3DC regime at $Ar=350$ , for $\overline{\unicode[STIX]{x1D70C}}=0.99$ and $\overline{\unicode[STIX]{x1D70C}}=0.05$ , respectively. The total lateral excursion of the sphere throughout its ascent is approximately $8d$ and $5d$ , respectively. This is typical of what we observed in that regime, although significantly smaller values ( ${\approx}2d$ ) were noticed in several cases. No clear trend regarding variations of this amplitude as a function of $\overline{\unicode[STIX]{x1D70C}}$ and $Ar$ emerged from the simulations throughout the explored range corresponding to the 3DC regime. In contrast, figure 14 reveals a striking difference between the dynamics of ‘light’ and ‘heavy’ spheres, the former repeatedly experiencing large lateral deviations over short periods of time. In other words, the horizontal fluctuations of the instantaneous sphere velocity are found to be much larger for $\overline{\unicode[STIX]{x1D70C}}=0.05$ than for $\overline{\unicode[STIX]{x1D70C}}=0.99$ . This trend subsists until the upper limit of our simulations, i.e. $Ar=700$ . It is in line with the findings reported by ZD who, for $Ar=350$ , noticed root-mean-square horizontal velocity fluctuations approximately twice as large for $\overline{\unicode[STIX]{x1D70C}}\approx 0$ than for $\overline{\unicode[STIX]{x1D70C}}\lesssim 1$ . This increase of horizontal velocity fluctuations for light spheres has a direct consequence on the rise velocity: since a larger part of the potential energy available in the system is used to ‘feed’ horizontal motions, less energy is injected in the vertical direction, resulting in a reduction of the average rise velocity. This is why, in the example of figure 14, the rise velocity of the lighter sphere is smaller by approximately 6 % than that of the ‘heavy’ one.

Figure 15.

Paths for various solid-to-fluid density ratios at $Ar=400$ ; (a) projection onto the vertical $(x,z)$ plane, (b) projection onto the horizontal plane. Same colour code as in figure 6.

The specificity of very light spheres takes a spectacular turn for $Ar\geqslant 400$ , as figure 15 shows for $Ar=400$ . Spheres corresponding to $\overline{\unicode[STIX]{x1D70C}}\geqslant 0.25$ follow fully chaotic trajectories. Decreasing $\overline{\unicode[STIX]{x1D70C}}$ , the motion of the sphere with $\overline{\unicode[STIX]{x1D70C}}=0.1$ is also found to be chaotic, but its projection onto the vertical $(x,z)$ plane reveals a more ordered structure, with in particular the emergence of a well-defined dominant frequency. The lightest two spheres confirm this evolution toward an ordered state as $\overline{\unicode[STIX]{x1D70C}}\rightarrow 0$ . Indeed, after some transient, they are seen to follow perfectly periodic paths corresponding to circular helices. Thus, a new periodic regime, hereinafter referred to as SP (spiralling), emerges within the ‘ocean’ of the 3DC regime. Figure 16 displays the characteristics of the corresponding spiralling path in the case $\overline{\unicode[STIX]{x1D70C}}=0.01$ . The diameter of the helix is about $1.3d$ and its pitch is $13.5d$ (i.e. $St\approx 0.074$ ), so that the path is inclined by about $17^{\circ }$ with respect to the vertical. Not surprisingly, the wake structure is similar to that encountered with spiralling bubbles (Mougin & Magnaudet Reference Mougin and Magnaudet2002b ; Cano-Lozano et al. Reference Cano-Lozano, Martínez-Bazán, Magnaudet and Tchoufag2016). That is, the streamwise vorticity is concentrated within two threads which wrap up around one another, the vorticity contained in each of them keeping a constant sign all along the path. This wake structure is stationary in a system of axes rotating with the sphere, thus delivering constant force and torque. Note the maximum magnitude of the vertical vorticity in the wake, which is typically three times larger than in the various zigzagging regimes encountered in § 3.3, suggesting a strong influence of the wake on the overall sphere dynamics. Note also the slight difference between the positive and negative extrema, which determines the clockwise rotation.

Figure 16.

The spiralling (SP) regime for $\overline{\unicode[STIX]{x1D70C}}=0.01$ and $Ar=400$ . (a) Three-dimensional path (the horizontal projection of the path is displayed in the plane $z=0$ ); (b) wake structure visualized with the $\unicode[STIX]{x1D706}_{2}$ criterion ( $-16.0\leqslant \unicode[STIX]{x1D714}_{z}\leqslant 14.5$ , see figure 4 for caption).

Using the refined grid described in § 2.3, we carried out a coarser exploration of the $(Ar,\overline{\unicode[STIX]{x1D70C}})$ plane from $Ar=400$ to $Ar=700$ , especially to determine whether the SP regime identified at $Ar=400$ exists over a significant parameter range or is rather an exception. To save CPU resources, only cases corresponding to $Ar=500,600$ and $700$ were simulated and the solid-to-fluid density ratio $\overline{\unicode[STIX]{x1D70C}}=0.001$ was omitted. Figure 17 shows how the path of a sphere with $\overline{\unicode[STIX]{x1D70C}}=0.05$ varies from $Ar=375$ until $Ar=700$ . Premises of the SP regime are already present at $Ar=375$ but no well-defined behaviour could be reached at that value of $Ar$ (the same conclusion holds for both smaller and larger $\overline{\unicode[STIX]{x1D70C}}$ ). In contrast, all paths obtained for $Ar\geqslant 400$ are seen to exhibit a clear spiralling structure. Most of the corresponding helices have a significant ellipticity, typically $1.5\pm 0.1:1$ , depending on $Ar$ . Some of them keep a vertical axis (i.e. the sphere rises vertically on average) but others do not. In particular, paths corresponding to $Ar=500$ and $Ar=600$ contain a significant horizontal drift which may either keep a constant direction ( $Ar=600$ ) or slowly rotate; in the latter case, the time-averaged path eventually describes an almost closed loop after a large number of periods of the primary helix (more than $25$ for $Ar=500$ ). Such paths, where a ‘slow’ periodic drift superimposes onto a ‘fast’ helical motion, are reminiscent of three-dimensional paths observed under certain conditions with falling disks in the fluttering or tumbling regimes (Auguste et al. Reference Auguste, Magnaudet and Fabre2013). To appreciate how the domain of existence of the SP regime varies with $Ar$ , figure 18 shows the trajectories obtained for the various solid-to-fluid density ratios at $Ar=700$ . All paths are fully chaotic and three-dimensional for spheres with $\overline{\unicode[STIX]{x1D70C}}>0.25$ . In contrast, after a long chaotic transient extending up to $z\approx 250$ , the sphere with $\overline{\unicode[STIX]{x1D70C}}=0.25$ starts to follow a clear, although not perfectly periodic, helical path. The remaining fluctuations discernible in the horizontal motion result in some meandering of the helix axis, and in some time-dependent variations of its ellipticity. As $\overline{\unicode[STIX]{x1D70C}}$ decreases, so does the relative magnitude of these fluctuations. The helical motion develops much more rapidly for $0.05\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 0.175$ and the corresponding paths are close to perfectly circular helices. Surprisingly, as figure 18(b) reveals, the sphere with $\overline{\unicode[STIX]{x1D70C}}=0.01$ follows an almost planar zigzagging path with only a weak lateral drift over a long time, rather than a spiralling trajectory. Not unlikely, during that time, the wake structure is similar to that of figure 7, i.e. the shedding mode is of the ‘4R’ type, instead of the twisted vortex geometry displayed in figure 16. However the path geometry slowly evolves toward an elliptical helix with, in the last stages of the simulation, an aspect ratio somewhat less than 3 : 1. All dynamical characteristics recorded in that peculiar case are similar to those found throughout the SP regime: the Strouhal number is $0.083$ (which corresponds to a maximum inclination of $24^{\circ }$ with respect to the vertical), and the crest-to-crest amplitude of lateral displacements is ${\approx}1.9d$ , similar to the diameter of the helix observed for $\overline{\unicode[STIX]{x1D70C}}=0.05$ ; in the next section it will be shown that the drag coefficient also behaves similarly. Hence there is no doubt that this case belongs to the SP regime, although it exhibits a long transient with a specific path geometry (interestingly, $Ar^{-1/2}\approx 0.038$ here, which according to § 3.1, makes this case with $\overline{\unicode[STIX]{x1D70C}}=0.01$ the only one in which the sphere’s angular dynamics is controlled by viscous effects at $Ar=700$ ).

Figure 17.

Path of a sphere with $\overline{\unicode[STIX]{x1D70C}}=0.05$ from $Ar=375$ to $Ar=700$ ; (a) projection onto the vertical $(x,z)$ plane, (b) projection onto the horizontal plane. Red: $Ar=375$ , green: $Ar=400$ , blue: $Ar=500$ , purple: $Ar=600$ , cyan: $Ar=700$ .

Figure 18.

Paths for various solid-to-fluid density ratios at $Ar=700$ ; (a) projection onto the vertical $(x,z)$ plane, (b) projection onto the horizontal plane. Same colour code as in figure 6.

Figure 19.

State diagram of rise regimes in the range $320\leqslant Ar\leqslant 700$ . Blue: three-dimensional chaotic (3DC); cyan: spiralling (SP).

The extent of the SP regime observed in present simulations is specified in figure 19 (the boundary between the 3DC and SP regimes was drawn without distinguishing between strictly periodic and ‘almost’ periodic helical paths). As this figure reveals, the maximum solid-to-fluid density ratio at which helical paths are observed sharply increases with $Ar$ , from $0.05$ for $Ar=400$ to nearly $0.25$ for $Ar=700$ . Throughout the SP regime, the diameter of the helices does not vary much and is always of the order of the sphere diameter. In contrast, for a given $\overline{\unicode[STIX]{x1D70C}}$ , their pitch decreases as $Ar$ increases, yielding increasing values of $St$ (e.g. for $\overline{\unicode[STIX]{x1D70C}}=0.1$ , from $St=0.058$ at $Ar=500$ to $St=0.085$ at $Ar=700$ ). Similarly, for a given $Ar$ , the smaller $\overline{\unicode[STIX]{x1D70C}}$ , the shorter the pitch (i.e. the larger the inclination of the path with respect to the vertical), with for instance at $Ar=700$ , $St=0.075$ for $\overline{\unicode[STIX]{x1D70C}}=0.25$ and $St=0.085$ for $\overline{\unicode[STIX]{x1D70C}}=0.1$ and $\overline{\unicode[STIX]{x1D70C}}=0.05$ . Circular helices with characteristics close to those discussed above were also identified by ZD in the range $375\leqslant Ar\leqslant 500$ for massless spheres ( $\overline{\unicode[STIX]{x1D70C}}=0$ ), as well as in the case $\overline{\unicode[STIX]{x1D70C}}=0.1,Ar=500$ ( $Ar=500$ was the upper limit considered in ZD’s simulations). Experiments reported in Veldhuis et al. (Reference Veldhuis, Biesheuvel and Lohse2009) also evidenced this type of trajectories for $\overline{\unicode[STIX]{x1D70C}}=0.02$ in the range $615\leqslant Ar\leqslant 875$ , the corresponding paths exhibiting a significant ellipticity (typically from 1.5 : 1 to 2 : 1). With $Ar=615$ (their figure 2(a)), the Strouhal number was found to be $0.078$ , which compares very well with present findings which, at $Ar=600$ , indicate $St=0.080$ for both $\overline{\unicode[STIX]{x1D70C}}=0.05$ and $\overline{\unicode[STIX]{x1D70C}}=0.01$ . Interestingly, nearly perfect zigzags with Strouhal numbers about $0.085$ and crest-to-crest amplitudes about $2d$ were observed in the same study, when $Ar$ was increased up to $O(10^{3})$ values (see cases $d$ , $e$ in their figures 2 and 3). These characteristics are essentially those found in the specific case $Ar=700,\overline{\unicode[STIX]{x1D70C}}=0.01$ described above. Hence the whole set of observations presented above closely agree with the phenomenology reported in that experimental study.

5.1 Overview of rise regimes

Figure 22 obtained by merging figures 2, 13 and 19, displays all rise regimes observed throughout this study. It helps make a global assessment of the influence of the sphere-to-fluid density ratio on the flow and path transitions up to $Ar=700$ . Note that in some cases, paths on which this regime map is based might not reflect the ultimate sphere trajectory (e.g. the late intermittency observed for $\overline{\unicode[STIX]{x1D70C}}=0.1$ in figure 5 or the late switching to the ZZ $_{2}$ regime for $\overline{\unicode[STIX]{x1D70C}}=0.99$ in figure 11). Although they only represent very long transients in such cases, they are relevant for comparison with available laboratory observations which were all carried out in tanks of limited height, with height-to-sphere diameter ratios of some hundreds (400 and 500 for the smallest spheres in Veldhuis & Biesheuvel (Reference Veldhuis and Biesheuvel2007) and Horowitz & Williamson (Reference Horowitz and Williamson2010), respectively).

Increasing $Ar$ from an $O(10^{2})$ value, all spheres successively pass through the SV, SO and PO regimes, and the latter is always succeeded by an intermittent regime spanning a narrow $Ar$ range. Then the ‘life’ of spheres becomes dramatically dependent on $\overline{\unicode[STIX]{x1D70C}}$ . Basically, those with $\overline{\unicode[STIX]{x1D70C}}\gtrsim 0.6$ reach the 3DC state through a succession of two periodic regimes corresponding to zigzagging paths with markedly distinct characteristics, planar large-amplitude low-frequency ZZ paths being succeeded by three-dimensional small-amplitude high-frequency ZZ $_{2}$ paths. Spheres with intermediate density ratios, say $0.2\lesssim \overline{\unicode[STIX]{x1D70C}}\lesssim 0.6$ , follow somewhat more complicated routes with an intermediate intermittent regime (with planar paths up to $Ar\approx 235$ and three-dimensional paths beyond that threshold) separating the ZZ and ZZ $_{2}$ regimes if $\overline{\unicode[STIX]{x1D70C}}\gtrsim 0.3$ , and a direct transition from the three-dimensional region of this IT regime to the 3DC state without any ZZ $_{2}$ range if $\overline{\unicode[STIX]{x1D70C}}\lesssim 0.3$ . Obviously the behaviour of very light spheres (say with $\overline{\unicode[STIX]{x1D70C}}\lesssim 0.2$ ) is by far the richest. Increasing $Ar$ from the first IT stripe ending at $Ar\approx 175$ , one first encounters a ZZ range which becomes very narrow as $\overline{\unicode[STIX]{x1D70C}}\rightarrow 0$ . It is followed by a transition region, the detailed shape of which strongly depends on $\overline{\unicode[STIX]{x1D70C}}$ , dominated by planar oblique zigzagging paths. These paths develop, first intermittently, then in a uniform manner for slightly larger $Ar$ , yielding the ZO periodic regime. Intermittent two- or three-dimensional oblique zigzags succeed this specific periodic regime (depending on the highest $Ar$ at which it stops for the considered $\overline{\unicode[STIX]{x1D70C}}$ ) before the system reaches the 3DC state. Last, a new ordered state emerges in the form of the SP regime characterized by paths with large-amplitude horizontal displacements associated with a Strouhal number of approximately $0.08$ . Paths in that regime are spiralling, with circular or elliptical horizontal projections (with ellipticity ratios up to 3 : 1), sometimes reached after a long nearly planar zigzagging transient.

5.2 Comparison with available results

Some partial comparison between the above findings and available computational and experimental results was performed in the course of §§ 3 and 4. However a more global discussion is in order to identify the points of consensus and disagreement that have been reached.

Starting with computations, to the best of our knowledge, only JDB and ZD dealt in detail with this problem. Present findings are in good agreement with those reported in these two studies up to $Ar\approx 175$ for $\overline{\unicode[STIX]{x1D70C}}\approx 0$ and $Ar\approx 190$ for $\overline{\unicode[STIX]{x1D70C}}\approx 1$ , which corresponds to the upper limit of the periodic oscillating (PO) regime. Only some minor differences were noticed regarding the $Ar$ threshold at the SO–PO transition. All three studies also conclude that the onset of chaos in the present system takes place at the upper bound of the PO regime. In the opposite limit, say $Ar>300$ , our results and those of ZD also yield converging conclusions (JDB’s investigation was limited to $Ar=350$ ). That is, both studies predict that all spheres follow chaotic three-dimensional paths, except those corresponding to very small $\overline{\unicode[STIX]{x1D70C}}$ which exhibit spiralling paths for large enough $Ar$ . Here also the thresholds are only slightly different: ZD detected the onset of the SP regime at $Ar=375$ for $\overline{\unicode[STIX]{x1D70C}}=0.0$ while we only observe it for $Ar\geqslant 400$ . No comparison can be performed regarding the extent of the SP subdomain, since ZD did not consider solid-to-fluid density ratios intermediate between $\overline{\unicode[STIX]{x1D70C}}=0.0$ and $\overline{\unicode[STIX]{x1D70C}}=0.2$ , except for $Ar=500$ where they also observed a helical path for $\overline{\unicode[STIX]{x1D70C}}=0.1$ . At $Ar=500$ , they reported drag coefficients of $0.67$ and $0.59$ for $\overline{\unicode[STIX]{x1D70C}}=0.0$ and $\overline{\unicode[STIX]{x1D70C}}=0.1$ , respectively. These values compare well with present findings, namely $\langle C_{D}\rangle =0.654$ and $0.598$ for $\overline{\unicode[STIX]{x1D70C}}=0.01$ and $\overline{\unicode[STIX]{x1D70C}}=0.1$ , respectively.

The situation is more complex in the intermediate $Ar$ range, i.e. in between the PO and 3DC regimes. Roughly speaking, ZD found that a vast majority of the $(Ar,\,\overline{\unicode[STIX]{x1D70C}})$ combinations in that range corresponds to chaotic regimes, with a few exceptions associated with periodic paths. Most of the latter belong to two disconnected subdomains, namely that of large-amplitude planar zigzags extending up to $\overline{\unicode[STIX]{x1D70C}}\approx 0.5$ for $180\lesssim Ar\lesssim 200$ , and that of ‘vertical oscillating’ paths extending beyond $\overline{\unicode[STIX]{x1D70C}}\approx 0.5$ in the range $250\lesssim Ar\lesssim 300$ . This contrasts with the large subdomains corresponding to periodic regimes in figure 22, especially for $\overline{\unicode[STIX]{x1D70C}}\geqslant 0.35$ where the planar ZZ regime and the three-dimensional ZZ $_{2}$ regime dominate the scene. Even for lower $\overline{\unicode[STIX]{x1D70C}}$ , the ZZ regime still spans a substantial range of $Ar$ down to $\overline{\unicode[STIX]{x1D70C}}\approx 0.15$ and the ZO regime represents a well-identified, albeit small, island of periodicity for low $\overline{\unicode[STIX]{x1D70C}}$ , which contrasts with ZD’s results in which only two isolated points were found to correspond to a ‘fast oblique zigzagging’ regime for $\overline{\unicode[STIX]{x1D70C}}\leqslant 0.2$ . At this stage, it must be mentioned that in this intermediate $Ar$ range, the regime map provided by ZD for $\overline{\unicode[STIX]{x1D70C}}<1$ is quite different from JDB’s original map, although both studies were carried out with the same code. Some of these differences may be attributed to the longer physical times reached in ZD’s simulations but the origin of others is less clear and presumably results from small changes between the two versions of the code. The most remarkable of these differences is the extent of the domain corresponding to the ZZ regime, which turns out to be much larger in JDB, where this type of path was observed up to $\overline{\unicode[STIX]{x1D70C}}\approx 1$ within the range $200\lesssim Ar\lesssim 215$ . Even for lighter spheres, the $Ar$ stripe corresponding to the ZZ regime is much broader in JDB: for $\overline{\unicode[STIX]{x1D70C}}=0.2$ it spans the range $175\leqslant Ar\lesssim 200$ , while it is restricted to the narrow interval $175\lesssim Ar\lesssim 182$ in ZD. Considering figure 22, it is clear that our results are in much better agreement with those of JDB in that respect; nevertheless, the ZZ regime subsists up to $Ar\approx 250$ when $\overline{\unicode[STIX]{x1D70C}}\gtrsim 0.6$ in our computations, while JDB reported three-dimensional chaotic paths beyond $Ar\approx 215$ in that range of $\overline{\unicode[STIX]{x1D70C}}$ . Another difference between JDB’s and ZD’s conclusions stands in the area of coexistence of the 3DC regime with the ‘vertical oscillating’ regime, which JDB found to take place within a wide region covering approximately $60$ $Ar$ -units for $\overline{\unicode[STIX]{x1D70C}}=0.0$ and broadening up to $100$ $Ar$ -units for $\overline{\unicode[STIX]{x1D70C}}=1.0$ , while it was only observed at $Ar=250$ and $260$ for $\overline{\unicode[STIX]{x1D70C}}=1.0$ by ZD: by running simulations over much longer times, the latter study could exclude most cases of bi-stability identified or suspected by JDB. Despite several attempts to disturb ZZ $_{2}$ paths or to modify initial conditions in the range where they are seen to exist in figure 22, we never observed any indication of bi-stability close to the ZZ $_{2}$ –3DC boundary; in contrast we noticed some sporadic switches to the ZZ regime close to the ZZ–ZZ $_{2}$ boundary, which might correspond to the small area of bi-stability detected by ZD. Last but not least, we stress again that the periodic ZZ $_{2}$ and ZO regimes observed here exhibit structural differences with, respectively, the ‘vertical oscillating’ and ‘fast oblique zigzagging’ regimes described in ZD, although their domains of existence coincide (for ZZ $_{2}$ ) or suggest a close connection (for ZO). Indeed we found the former to exhibit significant three-dimensional features (the only exception is the path corresponding to $\overline{\unicode[STIX]{x1D70C}}=0.99$ in figure 11) and a non-zero in-plane lateral drift, while ZD insisted on the ‘remarkably accurate planarity and verticality’ of the ‘vertical oscillating’ regime. Conversely, ZO paths stand in strictly vertical planes and correspond to low-frequency oscillations (with $St\approx 0.04-0.045$ ), while the ‘fast oblique zigzagging’ paths observed by ZD develop in inclined planes and are associated with high-frequency oscillations with $St\approx 0.1$ .

Elucidating the precise origin of the above differences is out of the scope of the present paper. However some qualitative clues can easily be obtained. First, some evolutions were probably not computed over the same amount of physical time, which could yield a different classification in cases where the ‘final’ regime occurs very late. In some other cases, a path may have been considered periodic in one study while it is actually only ‘almost perfectly’ periodic and was classified differently in the other. However most of the observed discrepancies cannot be resolved by these arguments: e.g. the planar large-amplitude zigzagging path observed in figure 9 over nearly $500d$ for $\overline{\unicode[STIX]{x1D70C}}=0.99$ and $Ar=237.5$ which clearly falls into the ZZ regime and was found to be chaotic, both by JDB and ZD. Consequently, such unsolved discrepancies question the sensitivity of the physical problem to numerical details. To fully appreciate this aspect, it is worth reiterating that the same two codes were used to explore the first non-vertical regimes of freely rising or falling disks (Auguste et al. Reference Auguste, Magnaudet and Fabre2013; Chrust, Bouchet & Dušek Reference Chrust, Bouchet and Dušek2013) and cylinders with various aspect ratios (Auguste Reference Auguste2010; Chrust, Bouchet & Dušek Reference Chrust, Bouchet and Dušek2014). In those situations, they were found to provide very close state diagrams, thresholds and characteristics of non-vertical paths, despite their many differences. Indeed, instead of the technical options summarized in § 2.2, the code used by JDB and ZD combines spectral elements and Fourier decomposition to discretize the Navier–Stokes equation, imposes zero-traction and zero-pressure boundary conditions on the domain outlet, and solves Newton’s equations using an implicit scheme (Jenny & Dušek Reference Jenny and Dušek2004). That two codes employing so different methods provide very similar predictions in the aforementioned problem represents a remarkable cross-validation. Hence one has to conclude that something specific makes computational predictions more delicate in the sphere case. There is little doubt that the crux of the problem is related to the point-symmetric nature of the sphere, which gives more freedom to the development of translational and rotational velocity disturbances in all directions of space than the axial symmetry of disks and cylinders. Any spatial discretization breaks this point symmetry, as discussed in appendix B for the spherical grid employed here. This is why we devoted a great deal of effort to minimize spurious grid-induced anisotropy effects. ZD and JDB employed a cylindrical grid that translates (but does not rotate) with the sphere. Influence of this grid geometry on the details of the evolution of light free spheres is difficult to appreciate, although convincing comparisons with a fixed-grid approach based on an immersed boundary technique were reported up to $Ar=250$ for falling spheres with $\overline{\unicode[STIX]{x1D70C}}=1.5$ (Uhlmann & Dušek Reference Uhlmann and Dušek2014). We believe that much of the differences observed above find their root in the distinct tiny anisotropy effects induced by these two different grid geometries.

Comparison with experiments was already examined in some detail in § 4 regarding the drag coefficient. Apart from that global parameter, few studies considered the path of freely rising spheres in detail. In Veldhuis & Biesheuvel (Reference Veldhuis and Biesheuvel2007), experiments were carried out with hollow spheres corresponding to $\overline{\unicode[STIX]{x1D70C}}\approx 0.56$ and $Ar=166,\,198$ and $218$ respectively (with uncertainties between $15$ and $20$ $Ar$ -units). According to figure 2, the first of them lies in the SO regime, close to the SO/PO border. Figure 13 predicts that the second lies in ZZ regime and the last one stands at the border of the ZZ and IT regimes. All paths recorded in that study exhibit non-negligible random fluctuations, including those of spheres with $Ar=166$ which are indeed inclined on average as expected from the SO regime, but display a significant curvature (their figure 5(d)). In the other two cases, the path was found to be only ‘occasionally zigzagging’, which clearly corresponds to an intermittent regime. Nevertheless, comparison of their figure 9(c,d) indicates that random fluctuations were much less intense with $Ar=198$ ; this was confirmed by ZD who computed specifically that case and found the corresponding path to be ‘almost exactly periodic’. As the authors themselves pointed out, the density distribution of the hollow spheres they used was certainly not perfectly spherically symmetric, and they attributed the randomness found in most of the recorded paths to this imperfect mass distribution. This is indeed highly plausible as JDB demonstrated by performing some computations with the sphere centre of mass shifted by $0.025d$ from its geometric centre. In particular, with $\overline{\unicode[STIX]{x1D70C}}=0.5$ and $Ar=200$ , they observed that the path quickly became chaotic instead of exhibiting periodic zigzags. To conclude with results from Veldhuis & Biesheuvel (Reference Veldhuis and Biesheuvel2007), one may only say that the issue of imperfect mass distribution prevents a one-to-one comparison with present predictions but experimental findings are not incompatible with them, since inclined and zigzagging paths clearly emerge as the dominant features in the first two cases and we found the third one to correspond to transitional conditions. Paths reported by Veldhuis et al. (Reference Veldhuis, Biesheuvel and Lohse2009) all concern expanded polystyrene spheres ( $\overline{\unicode[STIX]{x1D70C}}\approx 0.02$ ) with Archimedes numbers ranging roughly from $600$ to $1300$ . As pointed out in § 3.4, these paths look like elliptical spirals up to $Ar\approx 900$ , and are close to planar zigzags at larger $Ar$ . In all cases the Strouhal number is in the range $0.08\pm 0.005$ and the diameter or crest-to-crest amplitude of horizontal motions ranges from $1.5d$ to $2d$ . Moreover, as discussed in § 4, the drag coefficients were found to be nearly constant, irrespective of the precise path geometry (notwithstanding the 20 % uncertainty due to the background noise, which the authors attributed to the inhomogeneity of the spheres inherent to the process (air-assisted expansion) used to produce them). The corresponding average drag coefficient, $\langle C_{D}\rangle \approx 0.85$ , is 15 % higher than the largest values determined here but remains compatible with them, given this significant uncertainty. Hence, although limited to $\overline{\unicode[STIX]{x1D70C}}\approx 0.02$ , these experimental observations clearly support the existence of a SP regime with characteristics essentially similar to those predicted by present computations.

The last and most complete set of experimental observations available to date is provided by the extensive study of Horowitz & Williamson (Reference Horowitz and Williamson2010). As mentioned in § 1, they observed a dramatic split, according to the sphere-to-fluid density ratio, beyond $Ar\approx 195$ ( $Re_{m}\approx 260$ ). At larger $Ar$ , all spheres corresponding to $\overline{\unicode[STIX]{x1D70C}}\leqslant 0.36$ described almost planar large-amplitude periodic zigzags with horizontal displacements of the order of the sphere diameter and a Strouhal number about $0.06$ (see their figure 9 corresponding to $Re_{m}=450$ ). In contrast, spheres with $\overline{\unicode[STIX]{x1D70C}}>0.36$ went on rising along an oblique path up to $Ar\approx 370$ ( $Re_{m}\approx 600$ ), and followed an ‘intermittent oblique’ path beyond that point, up to $Ar\approx 970$ ( $Re_{m}\approx 1.5\times 10^{3}$ ). These findings are at deep variance with the regime map resulting from our computations, as well as with those of JDB and ZD. According to figure 13, no sphere with $\overline{\unicode[STIX]{x1D70C}}<1$ follows a steady or oscillating oblique path beyond $Ar\approx 190$ . In figure 22, nothing specific is seen to happen for sphere-to-fluid density ratios in the range 0.3–0.4 whatever $Ar$ . The large-amplitude zigzagging ZZ regime is observed both above and below that range of $\overline{\unicode[STIX]{x1D70C}}$ , and the corresponding Strouhal number is typically twice as small as that reported for ‘vibrating’ spheres in these experiments; hence, although both of them involve planar large-amplitude zigzags, the present ZZ regime and the ‘vibrating’ regime are entirely distinct. Here, spheres with $0.3\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 0.4$ stand in the 3DC regime for $Ar\gtrsim 250$ to $280$ , depending on $\overline{\unicode[STIX]{x1D70C}}$ . They never recover a periodic behaviour at higher $Ar$ , at least up to the upper limit reached in our computations. As is obvious from figure 21, their drag coefficient follows the approximate $Re_{m}^{-1/3}$ trend typical of fixed spheres, instead of exhibiting the $Re_{m}$ -independent behaviour observed in the above study. Values of $\langle C_{D}\rangle$ close to $0.75$ are only reached in the SP regime the bounds of which correspond to much higher $Ar$ and much smaller $\overline{\unicode[STIX]{x1D70C}}$ than those of the domain of ‘vibrating’ spheres determined by Horowitz & Williamson (Reference Horowitz and Williamson2010). The characteristics of the peculiar spiralling path observed within that regime for $\overline{\unicode[STIX]{x1D70C}}=0.01$ and $Ar=700$ , i.e. $Re_{m}\approx 950$ , may be considered compatible with those of these ‘vibrating’ spheres. However this does not resolve the above major discrepancies between the two state diagrams, nor does it give a clue to understand their origin.

5.3 Summary

To conclude this work, it sounds useful to summarize its main findings, together with the essential points of consensus now reached on the basis of all three available computational studies and experimental data, mostly those provided by Veldhuis & Biesheuvel (Reference Veldhuis and Biesheuvel2007) and Veldhuis et al. (Reference Veldhuis, Biesheuvel and Lohse2009):

1. (i) Starting from the steady vertical regime and increasing $Ar$ , all spheres encounter two successive bifurcations that yield steady oblique (SO) and periodic oblique (PO) paths, respectively.

2. (ii) The onset of chaos takes place at the third bifurcation which corresponds to the upper bound of the periodic oscillating regime, in line with the conclusions of JDB and ZD.

3. (iii) Beyond an intermittent regime that only spans some $Ar$ -units, all rising spheres enter a large-amplitude planar zigzagging (ZZ) regime. This regime subsists up to $Ar\approx 250$ for $\overline{\unicode[STIX]{x1D70C}}\gtrsim 0.6$ but only up to $Ar\approx 220$ for $0.15\lesssim \overline{\unicode[STIX]{x1D70C}}\lesssim 0.6$ . It only spans a few $Ar$ -units for lower $\overline{\unicode[STIX]{x1D70C}}$ . The zigzagging regime exhibits $O(d)$ lateral displacements and Strouhal numbers in the range 0.016–0.036, increasing with both $Ar$ and $\overline{\unicode[STIX]{x1D70C}}$ . Vortex shedding is characterized by the ‘4R’ structure evidenced by Horowitz & Williamson (Reference Horowitz and Williamson2010), with four pairs of vortices shed over one complete zigzag, two of which are much stronger that the other two.

4. (iv) The next steps on the route to the three-dimensional chaotic regime depend on $\overline{\unicode[STIX]{x1D70C}}$ . Spheres with $\overline{\unicode[STIX]{x1D70C}}\gtrsim 0.6$ directly switch to a three-dimensional periodic regime characterized by fast ( $St\approx 0.1$ ) small-amplitude ( $O(0.1d)$ ) zigzags with two vortex pairs of equal strength shed in the wake during one period. Spheres with $0.3\lesssim \overline{\unicode[STIX]{x1D70C}}\lesssim 0.6$ first exhibit an intermittent behaviour until $Ar\approx 250$ before they switch to the ZZ $_{2}$ regime. The latter disappears at $\overline{\unicode[STIX]{x1D70C}}\approx 0.3$ , so that, after the intermittent range $220\lesssim Ar\lesssim 250$ , spheres with $0.2\lesssim \overline{\unicode[STIX]{x1D70C}}\lesssim 0.3$ directly switch to the three-dimensional chaotic regime. Very light spheres with $\overline{\unicode[STIX]{x1D70C}}\lesssim 0.2$ exhibit a specific behaviour. From the end of the ZZ regime at $Ar\approx 180$ , until the onset of the three-dimensional chaotic regime in the range $250\lesssim Ar\lesssim 270$ , they develop large-amplitude oblique zigzags. This regime is mostly intermittent, except within a small region (PO regime) where paths are fully periodic with a significant ( ${\approx}4^{\circ }$ ) non-zero inclination and a Strouhal number of approximately $0.045$ ; these periodic paths are planar and develop within a vertical plane. The corresponding shedding mode involves four vortex pairs per period, but one is much stronger than the other three and is responsible for the average drift of the path.

5. (v) Whatever $\overline{\unicode[STIX]{x1D70C}}$ , the onset of three-dimensionality is found to correspond to $Ar\approx 235$ , i.e. $Re_{m}\approx 350$ , in all intermittent regimes (except within the distinct stripe separating the PO and ZZ regimes); only the ZZ and ZO periodic regimes remain planar somewhat beyond that threshold. The three-dimensional nature of ZZ $_{2}$ paths, which are only observed beyond $Ar\gtrsim 260$ , is consistent with this threshold.

6. (vi) Spheres follow three-dimensional fully chaotic (3DC) paths when $Ar$ exceeds a threshold value in the range 250–300, depending on $\overline{\unicode[STIX]{x1D70C}}$ . In that regime, the path of ‘light’ spheres undergoes much faster fluctuations than that of ‘heavy’ ones, resulting in a somewhat larger drag coefficient of the former.

7. (vii) In line with ZD’s observations, a specific spiralling (SP) periodic regime emerges from the 3DC regime beyond $Ar\approx 400$ for very light spheres. The maximum relative sphere density at which this regime is found is an increasing function of $Ar$ and is about $0.25$ at $Ar=700$ . Paths are characterized by circular or elliptical horizontal projections with an $O(d)$ average diameter and Strouhal numbers ranging approximately from $0.06$ to $0.085$ , increasing with $Ar$ and decreasing with $\overline{\unicode[STIX]{x1D70C}}$ . Similar to what is observed past spiralling bubbles, the corresponding wake structure is made of two disconnected intertwined vorticity filaments when the path is close to a circular helix, but is of the ‘4R’ type when the path ellipticity is large.

8. (viii) The drag of light enough spheres starts to exceed that of fixed spheres when the path becomes fully chaotic and three-dimensional. However one has to wait until the SP regime for the drag to exhibit clear non-standard variations with the rise Reynolds number. Then, these variations are characterized by an increase of $\langle C_{D}\rangle$ with the Reynolds number beyond a certain $Re_{m}$ ; the smaller $\overline{\unicode[STIX]{x1D70C}}$ the larger this increase. The maximum drag coefficients reached in the present study are about $0.75$ . This is consistent with (although somewhat smaller than) observations reported by Veldhuis et al. (Reference Veldhuis, Biesheuvel and Lohse2009) for spheres corresponding to $\overline{\unicode[STIX]{x1D70C}}=0.02$ and Reynolds numbers of $O(10^{3})$ , which follow spiralling or zigzagging paths with characteristics close to those determined here. This is also consistent with the drag of ‘vibrating’ spheres reported by Horowitz & Williamson (Reference Horowitz and Williamson2010), but the extent of the $(Ar,\overline{\unicode[STIX]{x1D70C}})$ domain where this phenomenology was observed in these experiments deeply differs from that of the SP regime revealed by present computations.

9. (ix) Comparison of present results with the available literature evidences the fact that the point-symmetric nature of the sphere makes the first bifurcations and transitional stages of the body+fluid system extremely sensitive to small spurious anisotropy effects, in both experiments and computations. Experiments have to face sphere inhomogeneities (which tend to make the body centre of mass distinct from its geometric centre), but computational grids also induce subtle anisotropy effects which may influence the predictions.

Finally, from a technical point of view, this study provided an opportunity to develop an original numerical algorithm summarized in § 2.2 to account for the instantaneous viscous force and torque resulting from the body acceleration. Following the steps provided in that section, we believe that this strategy which results in a second-order time accurate scheme, may be useful to other groups interested in an accurate tracking of the free motion of bodies with an arbitrary shape.