2.1. Governing equations and basic assumptions

We assume that all flow characteristics are independent of the azimuthal position around the rotation axis, but the local velocity has a non-zero azimuthal component. We further assume that the sphere rotates with the prescribed angular velocity of the container, which lies along the $z$-axis. This assumption is satisfied rigorously when the flow exhibits a perfect fore–aft symmetry with respect to the sphere equator, which is achieved in the limit ${R}o=0$. Nevertheless, we also carried out additional computations covering the whole range of flow conditions of interest here with the torque-free condition. In § 4.2, it will be shown that switching from one condition to the other has a negligible influence on the drag as long as ${R}o\lesssim 1$, and only a modest influence at higher ${R}o$, yielding relative drag differences of less than $10\,\%$. Assuming the flow to be incompressible and the fluid to be Newtonian, with density $\rho$ and kinematic viscosity $\nu$, the continuity and Navier–Stokes equations expressed in the reference frame rotating and translating with the sphere read

(2.1a,b)\begin{align} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{U}=0,\quad\partial_t \boldsymbol{U}+\boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{U} =-\rho^{-1}\, \boldsymbol{\nabla} P+ \nu\,\boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{U} + \boldsymbol{\nabla}\boldsymbol{U}^\text{T}) -2 {\varOmega}\boldsymbol{e}_z\times \boldsymbol{U}, \end{align}

with $\boldsymbol {U}$ the velocity field, $P$ the modified pressure including the centrifugal contribution, $\varOmega$ the imposed rotation rate, and $\boldsymbol {e}_z$ the unit vector in the $z$-direction.

2.2. Computational aspects

The computations are carried out with the second-order in-house finite volume code JADIM developed at IMFT. The spatial discretization of the velocity and pressure fields is performed on a staggered grid. Time integration of (2.1a,b) is achieved by combining a third-order Runge–Kutta scheme for advective and Coriolis terms with a semi-implicit Crank–Nicolson scheme for viscous terms. Incompressibility is satisfied at the end of each time step through a projection method. The accuracy of the complete time integration scheme is second order (Calmet & Magnaudet Reference Calmet and Magnaudet1997).

The boundary conditions are summarized in figure 2. A uniform velocity $U_\infty \boldsymbol {e}_z$ is imposed on the upstream and lateral boundaries. Since the reference frame translates and rotates with the sphere, the no-slip condition $\boldsymbol {U} = \boldsymbol {0}$ is enforced at the sphere surface, while on the flow axis the velocity components obey

(2.2a,b)\begin{equation} {{{{\boldsymbol{e}}_z\times\boldsymbol{U}}={\boldsymbol{0}}\quad \text{and} \quad {\boldsymbol{e}}_z\times\boldsymbol{\nabla}({\boldsymbol{U}}\boldsymbol{\cdot}{\boldsymbol{e}}_z)={\boldsymbol{0}}.}} \end{equation}

Hence only the axial velocity is non-zero on the axis, and its normal derivative vanishes there. Finally, the non-reflecting condition described by Magnaudet, Rivero & Fabre (Reference Magnaudet, Rivero and Fabre1995) is used on the downstream boundary. In short, the first (second) normal derivative of the tangential (normal) velocity component is set to zero on this boundary, together with the second-order cross-derivative of the pressure.

Figure 2. Sketch of the computational domain and boundary conditions (not to scale).

Since the axial sphere motion generates non-zero components of the Coriolis force, inertial waves take place when the Rossby number is small enough. These waves, whose wavelength is proportional to $U_\infty /\varOmega$, are emitted by the sphere and propagate downstream and outwards. Therefore, they are not ‘naturally’ evacuated from the computational domain. To prevent their energy from accumulating near the outer boundary, we add a sponge layer that damps them progressively without creating any reflection within the domain. For this purpose, we use the Rayleigh damping technique (Slinn & Riley Reference Slinn and Riley1998) which consists in replacing the exact velocity field $\boldsymbol {U}$ in this layer with the damped surrogate $\boldsymbol {U}^*$ defined as

(2.3)\begin{equation} \boldsymbol{U}^*= \boldsymbol{U} - \alpha \left( \boldsymbol{U} - \boldsymbol{U}_0 \right), \end{equation}

where $\boldsymbol {U}_0$ is some reference velocity reached by the flow close to the boundary, and $\alpha \in [0, 1]$ is a damping parameter. We select $\boldsymbol {U}_0 = U_{\infty } \boldsymbol {e}_z$ and $\alpha (\zeta ) = \exp [ - 6.125 ( \zeta / L_{sl} )^2 ]$, with $L_{sl}$ the thickness of the sponge layer, and $\zeta$ the local distance from the relevant outer boundary. We choose $L_{sl}$ such that at least five cells stand in the sponge layer, which was found sufficient to damp efficiently the inertial waves while limiting the thickness of this specific region within which the numerical solution is unphysical. The quality of the solutions provided by the present code in association with the above sponge layer technique may be appreciated in the work of Zhang, Mercier & Magnaudet (Reference Zhang, Mercier and Magnaudet2019) in the context of internal waves radiated by a sphere settling in a stratified fluid. A sketch of the computational domain specifying the treatment applied to each boundary is shown in figure 2.

In JADIM, the Navier–Stokes equations (2.1a,b) are expressed in a system of generalized orthogonal curvilinear coordinates. This makes it possible to use a variety of orthogonal boundary-fitted grids, as discussed by Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995). The detailed form of the governing equations expressed in this general coordinate system is also provided in Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995). Examples of solutions produced by this code associated with boundary-fitted grids for flows past spherical or spheroidal bodies, including in transitional or unstable regimes, may be found for instance in the works of Magnaudet & Mougin (Reference Magnaudet and Mougin2007) and Auguste & Magnaudet (Reference Auguste and Magnaudet2018). Here, following Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995), we employ an orthogonal grid built on the streamlines and iso-potential lines of the potential flow past a circular cylinder (figure 3). Accuracy of the solutions returned by the code on this type of grid may be appreciated in works such as Magnaudet et al. (Reference Magnaudet, Rivero and Fabre1995), Legendre & Magnaudet (Reference Legendre and Magnaudet1998) and Legendre, Magnaudet & Mougin (Reference Legendre, Magnaudet and Mougin2003), in which predictions for the forces acting on a spherical bubble in various two- and three-dimensional flow configurations are shown to compare very well with theoretical predictions in the limits of both low and high Reynolds number.

Figure 3. Computational grid (pressure nodes stand at the vertices). $(a)$ Close-up view in the sphere vicinity. $(b)$ Upper semi-domain $z\leq 0$ with $\mathcal {L}=180$ and $\mathcal {L}_\sigma =60$. (The sphere stands at the bottom right corner; 2 out of 3 cells have been removed in each direction for better visibility.)

With the above choice, the grid is nearly spherical in the sphere's vicinity (except close to the poles) and becomes gradually cylindrical as the distance to the sphere centre increases. Such a grid is particularly suitable for capturing efficiently not only the boundary layer surrounding the sphere, but also the wake and the near-axis upstream region even at very large distances from the body. As will become apparent later, such far-field regions are of particular importance in the present problem and could hardly be captured with a spherical grid. The grid is non-uniform close to the sphere and becomes uniform far from it. Uniformity in the far field allows the thickness of the sponge layer to be properly controlled. The use of very thin cells along the sphere surface, and the slow geometrical increase of the cell thickness as the distance to the sphere increases, allow the ‘inertial’ boundary layer (whose dimensionless thickness scales as ${R}e^{-1/2}$) and/or the Ekman boundary layer (scaling as ${T}a^{-1/2}$) to be captured properly throughout the considered range of parameters. Details on the grid design and a sensitivity study to some of the grid parameters are reported in Appendix A. When not stated otherwise, the half-length and radius of the computational domain (measured from the sphere centre and normalized by the sphere radius $a$) are given by $\mathcal {L} \times \mathcal {L}_\sigma =180 \times 60$, and the spatial discretization makes use of $314\times 96$ cells. Nevertheless, following the discussion of § 1, a detailed assessment of the influence of the axial confinement on the flow characteristics and the drag force is carried out in Appendix B, with $\mathcal {L}$ varied from $40$ to ${\approx }10^3$. Results of this sensitivity study are used in §§ 3 and 4 in the low-${R}o$ regime. Since the flow is expected to be invariant along the azimuthal direction over most of the conditions considered in Maxworthy's experiments, we opted for an axisymmetric resolution. Obviously, this simplification makes a parametric study much less expensive than a fully three-dimensional resolution. Nevertheless, it calls for some caution when the Reynolds number is large (typically ${R}e\gtrsim 100$) since the flow is known to be three-dimensional in that range in the absence of rotation. We will come back to this issue at the beginning of § 3. Starting from the uniform initial condition $\boldsymbol {U} = U_{\infty } \boldsymbol {e}_z$ throughout the flow domain, the computational time required to reach a converged stationary axisymmetric solution is approximately $2$ h on a standard single-processor workstation. The solution is considered converged when the relative time variation of the drag becomes less than $0.1\,\%$ over $5\times 10^4$ time steps.

2.3. Preliminary test

We first compare the local stress distribution at the sphere surface predicted with the above numerical set-up at small but non-zero Reynolds number with those obtained by Tanzosh & Stone (Reference Tanzosh and Stone1994), who made use of a boundary integral method in the creeping-flow limit. For this purpose we define the stress tensor $\boldsymbol {T}=-P\boldsymbol {I}+ \rho \nu (\boldsymbol {\nabla } \boldsymbol {U} + \boldsymbol {\nabla } \boldsymbol {U}^{\text {T}})$ (with $\boldsymbol {I}$ denoting the Kronecker delta), and the surface traction $\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {T}|_{r=a} =F_r{\boldsymbol {e}}_r+F_\theta {\boldsymbol {e}}_\theta +F_\phi {\boldsymbol {e}}_\phi$, with $\boldsymbol {n}\equiv {\boldsymbol {e}}_r$ the unit normal to the sphere pointing into the fluid, and $({\boldsymbol {e}}_r,{\boldsymbol {e}}_\theta, {\boldsymbol {e}}_\phi )$ the radial, polar and azimuthal unit vectors corresponding to the $(r,\theta, \phi )$ spherical coordinate system, with $r=0$ at the sphere centre, and $\theta =0$ (resp. ${\rm \pi}$) at the upstream (resp. downstream) pole. The $\theta$ variations of the three components of the surface traction are displayed in figure 4. The agreement is very good for the two tangential components, $F_\theta$ and $F_\phi$, although the values of the Taylor number in present simulations differ slightly from those of Tanzosh & Stone (Reference Tanzosh and Stone1994). The $F_r$ distributions also look similar, but differences growing from the equator to the poles and reaching approximately $10\,\%$ close to the latter may be noticed. In particular, while the $F_r$ distributions reported by Tanzosh & Stone (Reference Tanzosh and Stone1994) display an exact fore–aft antisymmetry (imposed by the ${R}o=0$ assumption), those provided by present results do not. This is obviously due to finite Reynolds number effects. These effects manifest themselves essentially on $F_r$ because this component of the traction reduces to the surface pressure, since the normal viscous stress vanishes on the sphere surface, owing to the combination of continuity and no-slip conditions. In contrast, only viscous stresses are involved in $F_\theta$ and $F_\phi$. Therefore, these traction components are less directly influenced by finite-${R}e$ effects, although a slight fore–aft asymmetry may be noticed in the central part of the distributions corresponding to ${R}e=5$, most notably on $F_\phi$.

Figure 4. Variations of the three components of the surface traction, normalized by $\rho \nu U_{\infty } /a$, versus the polar angle $\theta$. Blue solid and dashed lines display present results for $({R}e = 5,\ {T}a = 55.5)$ and $({R}e = 2,\ {T}a = 445)$, respectively; black solid and dashed lines display the zero-${R}o$ results of Tanzosh & Stone (Reference Tanzosh and Stone1994) for ${T}a=50$ and ${T}a=500$, respectively.

3.1. Preliminary comments

We now examine the salient features of the flow fields provided by the simulations in the parameter range ${T}a \in [20, 450]$, ${R}o \in [10^{-2}, 10]$ (which makes the Reynolds number vary in the range ${R}e \in [5, 300]$). By covering this range, we are in a position to compare numerical predictions with the full set of experimental data reported by Maxworthy (Reference Maxworthy1970). However, it must be stressed again that present results were all obtained in axisymmetric simulations, although it is known that for large enough Rossby numbers the flow is already three-dimensional at Reynolds numbers less than the upper limit considered here. Indeed, in the absence of rotation (${R}o\rightarrow \infty$), it is established that axial symmetry in the wake of a translating sphere breaks down at ${R}e\approx 105$ (Natarajan & Acrivos Reference Natarajan and Acrivos1993; Johnson & Patel Reference Johnson and Patel1999). In the presence of moderate rotation effects (${R}o=2$), Wang et al. (Reference Wang, Lu and Zhuang2004) showed that the flow past the sphere is still axisymmetric at ${R}e=100$, but is three-dimensional and unsteady at ${R}e=250$. However, since the governing equation for the vorticity becomes linear in the limit ${R}o\rightarrow 0$ (see the explicit form of the $z$-component of this equation below), no wake instability, hence no vortex shedding, can take place in this limit no matter how large the Reynolds number is. Consequently, it is expected that the lower ${R}o$ is, the higher the critical Reynolds number for the onset of three-dimensional effects becomes. A closely related increase of the critical ${R}e$ below which the wake remains stable was reported by Machicoane et al. (Reference Machicoane, Labarre, Voisin, Moisy and Cortet2018) with a circular cylinder towed perpendicularly to the axis of a rapidly rotating container under conditions ${R}o\lesssim 10$. More precisely, it was found that the cylinder's wake remains steady provided that ${R}e\lesssim 550/{R}o$, to be compared with ${R}e\leq 23.5$ in the limit ${R}o\rightarrow \infty$. Hence, considering that the constraints imposed on the flow in the low-${R}o$ limit delay drastically the transition to three-dimensionality in the sphere's wake, we expect present axisymmetric results to remain valid up to the maximum considered Reynolds number (${R}e=300$) at low enough Rossby number, typically ${R}o\lesssim 1$. (Unfortunately, how the critical ${R}e$ varies precisely with ${R}o$ is currently unknown.) Results corresponding to Reynolds numbers larger than $105$ and ${R}o\gtrsim 1$ require some more caution. However, even in that range, the influence of three-dimensional, possibly unsteady, effects on the drag is still limited up to ${R}e\approx 200$. For instance, in a non-rotating flow, the time-averaged drag at ${R}e=150$ is only $4\,\%$ larger than that predicted by constraining the flow to remain axisymmetric, and the relative amplitude of the drag oscillations is less than $1\,\%$ (Tomboulides & Orszag Reference Tomboulides and Orszag2000). Consequently, the comparison of present predictions for the drag with experimental data in the same range (we carried out a series of runs at ${R}e=167$) remains relevant. Only the few predictions corresponding to ${R}e=300$ and ${R}o > 1$ may really suffer from the fact that three-dimensional effects – which yield a chaotic but not yet turbulent regime in the wake at this Reynolds number in a non-rotating flow (Tomboulides & Orszag Reference Tomboulides and Orszag2000; Poon et al. Reference Poon, Ooi, Giacobello, Iaccarino and Chung2014) – are ignored in the present investigation.

3.2. General features

From now on, we analyse the flow field using the cylindrical coordinates $\sigma, \phi, z$, with $\sigma$ the cylindrical radius ($\sigma =0$ on the rotation axis), $\phi$ the azimuthal angle, and $z$ the axial distance from the sphere centre ($z<0$ upstream of the sphere, $z>0$ downstream of it). The flow past the sphere is presented in figures 5 and 6 for various values of ${T}a$ and ${R}e$. Figure 5 allows us to appreciate the details of the flow structure close to the sphere, while figure 6 makes use of a compression along the vertical axis at the higher two ${T}a$ values to display the entire recirculation regions. Figure 5 evidences the vertical and radial growth of the Taylor columns as the Taylor number is increased, and for a given ${T}a$, as ${R}o$ is decreased by decreasing ${R}e$. In line with previous observations, the fluid is seen to rotate more slowly (resp. faster) than the container in the upstream (resp. downstream) column. This feature may be rationalized by considering the governing equation for the axial vorticity $\omega _z=\partial _\sigma U_\phi$, namely

(3.1)\begin{equation} \partial_t\omega_z+U_\sigma\,\partial_\sigma\omega_z+U_z\,\partial_z\omega_z-\omega_z\, \partial_zU_z-\omega_\sigma\,\partial_\sigma U_z=2\varOmega\,\partial_zU_z+\nu\,\nabla^2_{\sigma,z}\omega_z, \end{equation}

where $\nabla ^2_{\sigma,z}$ stands for the two-dimensional Laplacian operator, and $\partial _\sigma$ and $\partial _z$ denote the partial derivatives with respect to the cylindrical coordinates $\sigma$ and $z$, respectively. Noting that $\omega _\sigma =-\partial _zU_\phi$ and $\partial _zU_\phi |_{\sigma \ll a}\approx \sigma (\partial _z\omega _z)|_{\sigma =0}$, the vortex tilting term $-\omega _\sigma \,\partial _\sigma U_z$ may be approximated as $\sigma \,\partial _z\omega _z|_{\sigma =0}\,\partial _\sigma U_z$ near the rotation axis. Therefore, all but one terms in (3.1) involve $\omega _z$ or its derivatives, which allows us to conclude that non-zero values of $\omega _z$ may be created only by the source term $2\varOmega \,\partial _zU_z$. Moving towards positive $z$ along the generatrix $\sigma =a$, the no-slip condition at the sphere surface forces the flow to decelerate ahead of the equatorial plane, implying $\partial _zU_z< 0$ for $z<0$. Conversely, the flow must accelerate downstream of the equatorial plane, yielding $\partial _zU_z>0$ for $z>0$. Therefore, starting from rest, negative (resp. positive) values of $\omega _z$ are generated in the upper (resp. lower) part of the cylindrical region $\sigma \leq a$. Normalizing velocities, distances and time by $U_\infty$, $a$ and $a/U_\infty$, respectively, and denoting provisionally normalized quantities with an overbar, the non-dimensional form of (3.1) reads ${R}o\,\overline {{lhs}}=2\,\partial _{\bar {z}}\bar {U}_z+{T}a^{-1}\,\nabla ^2_{\bar {\sigma }, \bar {z}}\bar {\omega }_z$, where $lhs$ stands for the left-hand side of (3.1). Now, the source term is $O(1)$, while the transport/stretching and diffusion terms are $O({R}o)$ and $O({T}a^{-1})$, respectively. The steady-state distribution of $\bar {\omega }_z$ depends on the relative intensity of advection/stretching and viscous diffusion at a given ${T}a$, hence on ${R}o$ (or equivalently ${R}e$). Considering panels (a–c) in both figures 5 and 6 for instance, the angular swirl $U_\phi /\sigma$ (which reduces to $\omega _z$ in the vicinity of the axis) is seen to approach a fore–aft symmetric distribution at ${R}e=8.9$ (figures 5a and 6$a$), and to become increasingly asymmetric as the Reynolds number increases, with $\omega _z\approx 0$ upstream of the sphere at ${R}e=167$ (figures 5c and 6$c$). In the latter case, advective effects are strong enough to reduce the flow region where $\omega _z$ exhibits significant values to a slender cylindrical zone in the wake.

Figure 5. Flow structure past the sphere in the parameter space $({R}e, {T}a)$. The flow is from top to bottom. The left half of each plot (red–blue scale) presents the angular velocity $U_{\phi } / \sigma$ (scaled by $U_\infty / a$); the right half (scale of greens) presents the velocity magnitude $\|\boldsymbol {U}\|$ (scaled by $U_{\infty }$) and streamlines in the sphere reference frame.

Figure 6. Same as figure 5 with the vertical axis compressed by a factor of $2$ (resp. $5$) for ${T}a = 117$ (resp. ${T}a = 445$) to capture the variations of the vertical extent of the recirculation regions.

Tanzosh & Stone (Reference Tanzosh and Stone1994) established that the recirculation regions appear at ${T}a \approx 50$ in the zero-${R}o$ limit. Figures 6(a,d), which correspond to a fairly low Reynolds number, support this prediction qualitatively, as the former (${T}a =23.2$) reveals no recirculation, while the latter (${T}a =117$) does. No upstream recirculation is found for ${T}a = 117$ and ${R}e=167$, i.e. ${R}o=1.43$ (figure 6$f$), which suggests that the condition required for an upstream recirculation region to be present is actually ${T}a \gtrsim 50$ and ${R}o\lesssim 1$. Note that in the three panels corresponding to Rossby numbers larger than unity (figure 6b,c,f), the spatial distribution of the angular velocity upstream of the sphere differs deeply from the columnar structure observed in all other cases. In figure 6$(\,f)$, the distribution downstream of the sphere looks also specific, with two well-separated maxima located on both sides of a tiny standing eddy detached from the body. The flow structure in figure 6$(c)$ (${R}o=7.2$) is similar to that observed in a non-rotating case, with a large standing eddy extending downstream of the sphere. In contrast, no such structure is present in figure 6$(b)$ (${R}o=2.25$), indicating that rotation is now controlling the flow structure in the near wake. Therefore, it may be concluded that rotation effects start to manifest themselves when the Rossby number is below some units, typically ${R}o\lesssim 5$. A similar transition is observed with particles settling in a linearly stratified fluid, the Froude number based on the Brunt–Väisälä frequency then playing the role of the Rossby number (Magnaudet & Mercier Reference Magnaudet and Mercier2020; Torres et al. Reference Torres, Hanazaki, Ochoa, Castillo and van Woert2000). When ${R}o<1$, the flow becomes more and more one-dimensional as the Taylor number increases, with the Taylor column extending far upstream and downstream of the body (figure 6g–i). In the same plots, the radius of the upstream column is seen to decrease with the Rossby number, the Stewartson layer getting closer to the surface of the fluid cylinder circumscribing the sphere as ${R}o\rightarrow 0$. At the largest ${T}a$, the upstream recirculation bubble extends more than 15 radii upstream of the sphere (figure 6$g$) and is shifted ahead of it by 2 radii. At such large ${T}a$, the flow experiences strong variations along the sphere circumference, within the thin Ekman layer that surrounds it. For instance, the fluid velocity near the equator is approximately $2.6$ times larger than the settling/rise velocity.

3.3. Near-body velocity distributions

Figure 7 shows how the three velocity components vary under different flow conditions with the radial position in four successive planes perpendicular to the axis, from $z=0$ (equatorial plane) to $z=-10$, a plane standing within the recirculation region in the low-${R}o$ limit. (Lengths and velocities are considered dimensionless throughout this section, being normalized by $a$ and $U_\infty$, respectively.) Disregarding provisionally the equatorial slice, one of the most significant features common to the three components is their large radial variation across the Stewartson layer standing around the mean position $\sigma =1$ and bounding externally the Taylor column. The peak values $U_z\approx 1.4$, $U_\sigma \approx 0.02$, $U_\phi /\sigma \approx -1.1$ reached by the three components in the plane $z=-2$ within this layer in the case ${R}o=4.5\times 10^{-3}$ agree well with the predictions of Tanzosh & Stone (Reference Tanzosh and Stone1994) for ${R}o=0$. Still with ${R}o=4.5\times 10^{-3}$, the near-axis plug-like profile of the axial velocity at $z=-2$ (figure 7b), with near-zero values up to $\sigma \approx 0.6$, is typical of the nearly geostrophic region. Moving upstream, $U_z$ is seen to take small negative values from the axis to $\sigma \approx 0.4$ (figure 7c,d), which gives an estimate of the radius of the recirculation region. In contrast, $U_z$ keeps significant positive values at any $z$ down to $\sigma=0$ in the most inertial case (${R}o=1.43$), which confirms the intuition that no nearly geostrophic or recirculation region exists under such conditions. Intermediate cases with $0.043\leq {R}o\leq 0.44$ (all with ${T}a=117$) exhibit a nearly geostrophic behaviour up to $\sigma \approx 0.3$ in the plane $z=-2$ (figure 7b). In contrast, the axial velocity keeps significant positive values down to the axis at $z=-10$ in these cases, showing that this plan stands beyond the tip of the recirculation region whatever the Rossby number for ${T}a=O(10^2)$.

Figure 7. Radial slices of the velocity field in several $z=const.$ planes ahead of the sphere: (a–d) axial component $U_z$; (e–h) radial component $U_\sigma$; (i–l) angular swirl $U_\phi /\sigma$. Velocities and distances are normalized by $U_\infty$ and $a$, respectively. Lines become lighter as the Rossby number increases: black indicates ${R}o = 4.5\times 10^{-3}$, ${R}e=2$ (${T}a = 445$); dark purple indicates ${R}o = 4.3\times 10^{-2}$, ${R}e=5$; dark blue indicates ${R}o = 0.138$, ${R}e=16.1$; medium blue indicates ${R}o = 0.444$, ${R}e=51.9$; pale blue indicates ${R}o = 1.43$, ${R}e=167$ (${T}a = 117$ in the latter four cases).

Returning to the case ${R}o=4.5\times 10^{-3}$, the near-axis profile of the angular swirl is seen to flatten gradually as the distance to the sphere increases, with on-axis values of $|U_\phi |/\sigma$ increasing from $0.6$ at $z=-2$ to $1.1$ at $z=-10$, approximately (figures 7j–l). Again, these findings are consistent with those of Tanzosh & Stone (Reference Tanzosh and Stone1994). Since $U_\phi /\sigma \approx \omega _z$ near the axis, the reason for this gradual increase and final plug-like profile may be understood by using (3.1). When ${R}o\rightarrow 0$, axial variations of $\omega _z$ ahead of the sphere can arise only through the non-zero source term resulting from the weak axial variations of $U_z$. Radial variations of $\omega _z$ being negligible near the axis, one then has $2\,{T}a\,\partial _zU_z\approx -\partial _{zz}\omega _z$. Thus viscous diffusion is seen to induce a non-zero curvature in the axial profile of $\omega _z$. The axial velocity increasing from small negative values in the recirculation region to near-zero values in the nearly geostrophic region, the axial gradient $\partial _zU_z$ is positive, yielding $\partial _{zz}\omega _z<0$. Moreover, at a given radial location $\sigma \neq 0$, $U_\phi$ increases from negative values upstream of the sphere to zero at its surface, while it remains null along the axis. Therefore, $\partial _\sigma (\partial _zU_\phi )$ is positive, implying $\partial _z\omega _z>0$ at the sphere surface. Combining the above two inequalities leads to the conclusion that $\partial _{z}\omega _z$ is necessarily positive (and larger than its surface value) ahead of the sphere, which translates into an increase of the angular swirl (in absolute value) as $|z|$ increases, in line with the behaviour observed in figures 7(j–l). The argument still holds up to $z=-5$ for the two intermediate cases with ${T}a=117$ and ${R}o<0.2$. However, the plane $z=-10$ stands beyond the recirculation region in these cases, as the significant positive values of the axial velocity ($U_z\approx 0.1$) confirm. Hence $\partial _zU_z$ is negative and quite large beyond $z=-5$. This makes $\partial _{zz}\omega _z$ positive and significantly larger than in the zone closer to the sphere, leading to $\partial _{z}\omega _z<0$ beyond the recirculation region, and therefore to a reduction of the angular swirl as $|z|$ increases.

Symmetry arguments imply that the radial and azimuthal velocity components must both vanish on the equatorial plane at ${R}o=0$. Hence their non-zero values in that plane (figures 7e,i) give insight into the strength of advective effects. Since these effects tend to enhance the amount of fluid transported from the upstream Taylor column to the downstream column through the Ekman layer, the main features of the $U_\sigma$ and $U_\phi$ near-surface distributions at $z=0$ when ${R}o$ is non-zero are expected to resemble those found slightly above the equatorial plane in the zero-${R}o$ limit. The $O(1)$ values of the axial velocity in the median part of the Stewartson layer upstream of the sphere (figure 7b), combined with its large peak values in the Ekman layer at $z=0$ (figure 7a), result in a positive $\partial _zU_z$ upstream of the equatorial plane at radial positions $\sigma \approx 1$. Continuity combined with the no-slip condition at the sphere surface then implies $U_\sigma <0$ for $\sigma \gtrsim 1$ above the equatorial plane. This is why, in the presence of finite inertial effects, one expects $U_\sigma$ to be negative near the sphere surface in that plane. This is indeed the case as long as the near-surface peak of $U_z$ subsists (figure 7e). More specifically, the magnitude of the (negative) peak value of $U_\sigma$ and $U_\phi /\sigma$ within the Ekman layer is seen to increase strongly with the Rossby number as long as ${R}o$ is less than unity. The peak shifts away from the sphere surface as ${R}o$ increases, and its magnitude at ${R}o=0.44$ is close to $0.2$ and $0.6$ for the (inward) radial velocity and angular swirl, respectively (figures 7e,i). The strongly inertial case corresponding to ${R}o=1.43$, ${R}e=167$ behaves differently, with especially $U_\sigma$ first taking positive values within the part of the boundary layer closest to the sphere surface (figure 7e). Within the Ekman layer, the axial velocity reaches a maximum close to $2.4$ at ${R}o=4.5\times 10^{-3}$ (figure 7a). This value is in line with the findings of Tanzosh & Stone (Reference Tanzosh and Stone1994), who reported a maximum of $2.25$ at ${R}o=0$. The large positive values of $U_z$ within the Ekman layer play a pivotal role in the overall dynamics of the flow in the low-${R}o$ regime, as they control directly the amount of fluid transported from the upstream Taylor column to the downstream column. Inertial effects are found to change the ${R}o=0$ picture dramatically, lowering the $U_z$ maximum to $1.7$ at ${R}o=0.44$, which, taking the unit freestream velocity as reference, corresponds to a $50\,\%$ reduction of the peak. Putting the findings observed in the equatorial plane on the three velocity components together, it appears that inertial effects deeply modify the local flow structure within the Ekman layer, which may be expected to have direct consequences on the stress distribution at the sphere surface, hence on the drag.

3.4. Extent of the upstream recirculation region

Figure 8 shows how the extent of the upstream recirculation region, $\ell _{s}$, varies as a function of ${R}o$ for various ${T}a$. We define $\ell _{s}$ as the distance (normalized by the sphere radius $a$) from the sphere centre to the farthest upstream location where the axial velocity changes sign on the rotation axis. Strictly speaking, as figure 1 shows, the recirculation region stands in between the two locations where the axial velocity changes sign, the one closest to the sphere defining the tip of the nearly geostrophic region. However, we follow Maxworthy (Reference Maxworthy1970), who, using dye visualizations, focused on the location of the tip of the recirculation region. In line with his observations, $\ell _{s}$ reaches a plateau when ${T}a$ is kept fixed and ${R}o\rightarrow 0$. Conversely, $\ell _{s}$ vanishes when ${R}o \rightarrow 1$, implying that no recirculation region exists for ${R}o>1$ (see also figures 6 and 7). For a fixed ${R}o$, $\ell _{s}$ decreases as ${T}a$ is decreased, down to a critical Taylor number close to $50$, below which the recirculation region disappears. As figure 8(b) shows, the simulations recover the prediction $\ell _{s} \approx 0.052\, {T}a$ (Tanzosh & Stone Reference Tanzosh and Stone1994) in the range $50 \lesssim {T}a \lesssim 200$, ${R}o \lesssim 5\times 10^{-2}$. However, as the symbols in the upper left corner of figure 8(a) reveal, the numerical results deviate from this prediction as well as from Maxworthy's data for the largest value of ${T}a$ considered here, i.e. ${T}a=445$. For instance, we find $\ell _{s} \approx 21$ for ${R}o = 2\times 10^{-2}$, which is significantly less than the value $\ell _{s} =23$ reported by Tanzosh & Stone (Reference Tanzosh and Stone1994) in the zero-${R}o$ limit. We attributed this discrepancy to axial confinement effects, a track already suggested by Ungarish & Vedensky (Reference Ungarish and Vedensky1995). To check this hypothesis, we increased the length of the computational domain from $\mathcal {L}=180$ to $\mathcal {L}\approx 10^3$ along the lines discussed in Appendix B, where a detailed analysis of the sensitivity of the recirculation length and the drag force to these effects is presented. As the star symbols in figure 8(b) show, the recirculation length obtained with this much longer domain is in excellent agreement with the zero-${R}o$ prediction. This is a clear indication that the characteristics of the recirculation region, and more generally those of the Taylor column, are extremely sensitive to axial confinement effects, even in containers with $\mathcal {L}=O(10^2)$. Indeed, although the tip of the recirculation region stands far away from the top and bottom ends of the domain, the tip of the Taylor columns interacts directly with them when ${T}a$ is large and ${R}o\rightarrow 0$, and the corresponding blocking effect is sufficient to alter the characteristics of the various zones of the flow located much closer to the body.

Figure 8. Extent of the upstream recirculation region: $(a)$ $\ell _{s}$ versus ${R}o$ for various values of ${T}a$; $(b)$ same with $\ell _{s}$ normalized by ${T}a$. Symbols $\lozenge, \star$ correspond to simulations with $\mathcal {L}=180$ and ${\approx }10^3$, respectively; $\triangleright$ corresponds to Maxworthy's experiments (Reference Maxworthy1970); $\triangleleft$ corresponds to boundary-integral simulations at ${R}o=0$ (Tanzosh & Stone Reference Tanzosh and Stone1994). The solid black line indicates the zero-${R}o$ limit $\ell _{s}=0.052\, {T}a$ (Tanzosh & Stone Reference Tanzosh and Stone1994). Dark blue, red, green and yellow symbols refer to ${T}a=55,117,193,445$, respectively.

3.5. Inertial wave pattern

To finish with the characterization of the flow field, it is of interest to look at the dominant feature of the flow outside the Taylor column, namely the inertial wave field radiated by the sphere. The generation of such waves by bodies moving in a rotating fluid, or rotating topographies subject to a transverse flow, is well documented (Greenspan Reference Greenspan1968). Taylor (Reference Taylor1922) predicted the existence of these waves and discovered that they exhibit an anisotropic dispersion property, with a radian frequency $\omega _{I0}$ obeying the orientation-dependent dispersion relation $\omega _{I0} = \pm 2 \varOmega \cos \psi$, with $\psi$ the angle between the wavevector ${\boldsymbol {k}}$ and the rotation axis. He also pointed out that, remarkably, this dispersion relation holds irrespective of the wave amplitude. However, unlike internal waves in a stably stratified fluid, pure inertial waves take place in a homogeneous fluid, which makes their experimental observation more difficult (Pritchard Reference Pritchard1969). For this reason, Taylor could not observe the waves whose existence he had predicted. Nevertheless, by releasing a light sphere on the axis of a rotating cylinder, he could visualize the existence of the resting column of fluid that was later named after him. Much later, waves generated by a pulsating, oscillating or transversely moving circular cylinder in a rotating tank could be visualized by Machicoane et al. (Reference Machicoane, Cortet, Voisin and Moisy2015, Reference Machicoane, Labarre, Voisin, Moisy and Cortet2018) using particle image velocimetry. We are not aware of similar experimental observations in the configuration considered here. The numerical investigation of Wang et al. (Reference Wang, Lu and Zhuang2004) provides some streamline maps for ${R}o=2$ and ${T}a=50$ and $125$, from which the wave pattern in a moderately rotating flow may be inferred.

In the present configuration, the waves are radiated by the sphere moving relatively to the undisturbed fluid with velocity $-U_\infty \boldsymbol{e}_z$. Therefore, in the reference frame attached to the body, the radian frequency has to be corrected from the corresponding Doppler shift and becomes $\omega _I=\omega _{I0}+U_\infty {\boldsymbol {k}}\boldsymbol {\cdot }\boldsymbol{e}_z$, so that $\omega _I$ obeys (Lighthill Reference Lighthill1967; Whitham Reference Whitham1974)

(3.2)\begin{equation} \omega_I= 2({\rm \pi} U_\infty/\lambda\pm\varOmega) \cos\psi. \end{equation}

Equation (3.2) indicates that the relative displacement of the body along the rotation axis allows the existence of axisymmetric standing waves with wavelength $\lambda ={\rm \pi} U_\infty /\varOmega$, i.e. $\varLambda \equiv \lambda /a={\rm \pi} \,{R}o$, as predicted by Taylor (Reference Taylor1922). In figure 9, we use the pressure disturbance field to visualize the wave pattern at three different values of the Rossby number. Figure 10 shows that the wavelength determined by seeking the minimum distance separating two successive crests (yellow arrows in figure 9) agrees closely with Taylor's theoretical prediction.The details of the wave field, i.e. the spatial distribution of $\psi$, are dictated by the no-penetration condition at the body surface and are influenced by viscous effects, especially those controlling the Ekman layer (Cheng & Johnson Reference Cheng and Johnson1982; Johnson Reference Johnson1982).

Figure 9. Visualization of the inertial wave pattern for several ${R}o$: (a) $Re=51.9$, $Ro=2.24$; (b) $Re=167$, $Ro=1.43$; (c) $Re=167$, $Ro=0.38$. The relative flow with respect to the sphere is from left to right. Colours refer to the amplitude of pressure variations (scaled by $\frac {1}{2}\rho U_\infty ^2$).

Figure 10. Variations of the wavelength with respect to the Rossby number: dots indicate simulations, and the line indicates the theoretical prediction $\varLambda ={\rm \pi} \,{R}o$.

Energy is radiated by the waves with the Doppler-shifted group velocity ${\boldsymbol {c}}_g=\partial _{\boldsymbol {k}}\omega _{I}={\boldsymbol {c}}_{g0}+U_\infty {\boldsymbol {e}}_z$, with ${\boldsymbol {c}}_{g0}=\partial _{\boldsymbol {k}}\omega _{I0}$. The axial and radial components of ${\boldsymbol {c}}_g$ are $c_{gz}=U_\infty -{\rm \pi} ^{-1}\varOmega \lambda \sin ^2\psi$ and $c_{g\sigma }=(2{\rm \pi} )^{-1}\varOmega \lambda \sin 2\psi$, respectively. Consequently, standing waves with $\lambda ={\rm \pi} U_\infty /\varOmega$ have axial and radial group velocities $c_{gz}=U_\infty \cos ^2\psi$ and $c_{g\sigma }=\frac {1}{2}U_\infty \sin 2\psi$, respectively, and their energy propagates along straight rays $\sigma =(z-z_0)\tan \psi$, with $z_0$ the origin of the ray on the rotation axis. According to figure 9, the angle $\psi$ increases from approximately ${\rm \pi} /3$ downstream of the sphere to values close to ${\rm \pi} /2$ upstream. Therefore, the axial and radial components of the energy flux are positive everywhere, i.e. they are directed downstream and outwards, respectively. Moreover, they decrease continuously as one moves upstream, and eventually vanish when the wave crests become parallel to the axis ($\psi ={\rm \pi} /2$). Therefore, far upstream of the sphere, the wave energy does not propagate at all, i.e. it just travels with the sphere. Comparing figures 9(a–c) indicates that the lower ${R}o$, the more the wave fronts are parallel to the rotation axis at a given position upstream of the sphere. Therefore, the lower ${R}o$, the shorter the upstream position at which the wave energy stops propagating.

Examination of the whole set of computational results reveals the presence of inertial waves with characteristics similar to those discussed above for Rossby numbers in the range $0.2 \lesssim {R}o \lesssim 5$. These two limits result from totally distinct reasons. At ${R}o= 5$, the wavelength is approximately one-third of the radius of the computational domain (and even half that size if the sponge layer is not considered). Hence the outer cylindrical boundary affects the distribution of the disturbances radiated by the body at larger ${R}o$, preventing the formation of standing waves. Conversely, viscous effects are responsible for the disappearance of waves for ${R}o\lesssim 0.2$. Indeed, a disturbance with wavevector ${\boldsymbol {k}}$ is damped at a rate $-\nu \,\|{\boldsymbol {k}}\|^2$. Hence the ratio $r_v$ between the viscous force and the restoring Coriolis force acting on the disturbance is $2{\rm \pi} ^2\nu /(\varOmega \lambda ^2)=2{\rm \pi} ^2/(\varLambda ^2\,{T}a)$, which for $\varLambda ={\rm \pi} {R}o$ yields $r_v=2({T}a\,{R}o^2)^{-1}$. Consequently, the lower ${R}o$, the larger $r_v$ at a given ${T}a$, with for instance $r_v=1$ for ${R}o=0.1$ and ${T}a=200$.

4.1. Drag

Figure 11 presents the drag coefficient obtained through a direct integration of the surface traction defined in § 2.3 over the sphere. Figure 11$(a)$ shows the compensated drag coefficient $C_D\,{R}e/12$ as a function of the Reynolds number for various ${T}a$, while figure 11$(b)$ shows $C_D$ as a function of the Rossby number for various ${R}e$. The standard drag curve for a sphere translating in a quiescent fluid, based on the empirical correlation $C_D({R}e)=12\,{R}e^{-1}\,(1+0.241\,{R}e^{0.687})$ (Schiller & Naumann Reference Schiller and Naumann1933), and the inviscid prediction (1.1), are also shown as references. For reasons discussed in § 3.1, only the few numerical predictions corresponding to ${R}e=300$ and ${R}o>1$ (i.e. the three rightmost lozenges located below the dotted line in figure 11a) are expected to be altered significantly by the absence of three-dimensional effects in the computed solutions.

Figure 11. Drag coefficient versus the Reynolds and Rossby numbers: (a) compensated drag coefficient $C_D Re/12$ versus ${R}e$; (b) $C_D$ versus ${R}o$. Symbols $\blacklozenge, \bigstar$ correspond to simulations with $\mathcal {L}=180$ and $\approx 10^3$, respectively; $\circ$ corresponds to Maxworthy's experiments (Reference Maxworthy1970). In (a), the solid line is the standard drag curve for a sphere translating in a fluid at rest; the vertical dashed line corresponds to the threshold beyond which the wake is three-dimensional in the limit ${R}o\rightarrow \infty$; and the dotted line is a guide for the eye separating data belonging to the range ${R}o>1$ (below the line) from those for which ${R}o<1$ (above the line). The solid line in $(b)$ corresponds to the inviscid prediction (1.1). Symbols are coloured according to the value of ${T}a$ in $(a)$, and ${R}e$ in $(b)$.

At low to moderate Reynolds number, say ${R}e\lesssim 50$, the drag is significantly larger than predicted by the above correlation, highlighting the influence of the rigid-body rotation. Present results agree well with those of Maxworthy (Reference Maxworthy1970) up to ${T}a\approx 80$, i.e. in the range where rotation effects are moderate. In contrast, they clearly deviate from experimental data for larger ${T}a$, predicting a lower drag. The lower ${R}e$, the larger the deviation at a given ${T}a$, the relative difference between the two values exceeding $40\,\%$ at ${R}e=5$ for the highest ${T}a$. Conversely, the larger ${T}a$, the larger the Reynolds number at which the deviation starts. Thus numerical predictions and experimental data still agree for large enough ${R}e$ when ${T}a$ is large. In the $C_D$ versus ${R}o$ representation of figure 11(b), numerical predictions are seen to fall within the somewhat scattered interval of experimental values for ${R}o \gtrsim 2\times 10^{-1}$. In contrast, below this threshold, the numerical series departs from the experimental one, and the departure increases as ${R}o$ decreases. It may be noticed that all numerical data obtained for ${R}o\lesssim 0.3$ stand beyond the inviscid prediction (1.1). As these data correspond to Reynolds numbers less than $200$, viscous effects are likely to be responsible for the observed difference. This will be confirmed later.

In figure 11(a), the drag is observed to be lower than predicted by the standard law at large enough ${R}e$ and low enough ${T}a$, say ${R}e\gtrsim 70$ and ${T}a\lesssim 80$ (consider the last three purple lozenges and the very last dark green lozenge at the bottom right). This is in line with Maxworthy's experimental findings as the dots confirm, the associated Rossby number being such that ${R}o\gtrsim 1$ throughout this regime. The same behaviour was observed numerically by Rao & Sekhar (Reference Rao and Sekhar1995) and Sahoo et al. (Reference Sahoo, Sarkar, Sivakumar and Sekhar2021). As non-axisymmetric effects not accounted for in present simulations are known to increase the drag in a non-rotating flow, one might suspect their absence to be at the origin of the low numerical drag values found in the high-${R}e$ range. However axisymmetry in the sphere's wake breaks down only at ${R}e\approx 105$ when ${R}o\rightarrow \infty$ (dashed line in figure 11a), and the computed drag at ${R}e=93$ and ${T}a=23.2$ (most left lozenge below the solid line in the figure) is $15\,\%$ lower than predicted by the standard drag law. Similarly, for the same ${T}a$ but ${R}e=167$, the drag is $20\,\%$ smaller than expected on the basis of the standard drag law, whereas the axisymmetric prediction is known to underestimate the drag by only $4\,\%$ in the limit ${R}o\rightarrow \infty$ in that ${R}e$ range (see the discussion in § 3.1). Therefore, one can conclude that the difference observed in the figure is really the result of rotation effects and has the same physical origin as that revealed by Maxworthy's experimental data. In the corresponding ${R}e$ range, a large standing eddy is present behind the sphere. Figure 12 shows how the rigid-body rotation alters the size of this eddy, together with the pressure and axial velocity distributions. Even a modest level of rotation (the lower half of each plot corresponds to ${R}o=5.4$) is seen to reduce significantly the negative axial velocity within the eddy, increasing the pressure in its core. Therefore, compared with the case of a sphere translating in a fluid at rest, the overall pressure difference between the front and rear stagnation points is reduced, lowering the pressure drag. This effect is significant, as the drag may be reduced by $10$–$20\,\%$ with respect to the standard law in the range $100 \lesssim {R}e \lesssim 300$. Maxworthy (Reference Maxworthy1970) argued that the pressure increase in the core of the standing eddy is due to the fact that ‘the outward flow of rotating fluid over this [recirculation] bubble causes it to rotate at a rate less than the applied value’. However, present results contradict this explanation. For instance, figure 6(c) for ${R}e=167$, ${T}a=23.2$ shows that the angular velocity within the standing eddy is larger than the applied rotation rate, and the corresponding drag (penultimate purple lozenge in the bottom right corner of figure 11a) stands $15\,\%$ below the standard drag curve. Actually, the origin of the drag reduction may be understood by considering the governing equation for the azimuthal vorticity, $\omega _\phi =\partial _zU_\sigma -\partial _\sigma U_z$. Rotation enters the $\omega _\phi$ balance through the source term $2\varOmega \,\partial _zU_\phi$, similar to that involved in (3.1), but with $\partial _zU_z$ replaced with $\partial _zU_\phi$. In the vicinity of the axis, this source term virtually equals $2\varOmega \sigma \,\partial _z\omega _z$. As discussed in § 3.2, $\omega _z$ increases downstream of the sphere with the distance to the rear stagnation point (as figure 6c confirms). Therefore, this source term is positive within the standing eddy, bringing a positive variation in $\omega _\phi$ compared to the non-rotating configuration. Near the axis, $\omega _\phi \approx -\partial _\sigma U_z$, so that this change in $\omega _\phi$ translates into an increase in $U_z$ as $\sigma \rightarrow 0$, i.e. a positive variation of the axial velocity as the rotation axis is approached. Hence when $U_z$ is negative in the limit ${R}o\rightarrow \infty$, finite rotation effects decrease its magnitude, leading to an increase in the local pressure, from which the observed drag reduction ensues.

Figure 12. Influence of the rigid-body rotation on the velocity and pressure in the sphere vicinity at high Reynolds number (${R}e=300$): (a) axial velocity (scaled by $U_{\infty }$); (b) pressure (scaled by $\frac {1}{2}\rho U^2_{\infty }$). The upper half of each plot corresponds to a non-rotating flow (${R}o=\infty$); the lower half corresponds to ${R}o = 5.4$. The relative flow with respect to the sphere is from left to right.

In figure 13, the experimentally and numerically determined drag coefficients are plotted versus the semi-empirical prediction (1.3) suggested by Tanzosh & Stone (Reference Tanzosh and Stone1994). This prediction is expected to be valid at arbitrary Taylor number. In contrast, it is supposed to apply only as long as the Reynolds number is low, given the range of validity of (1.2) from which the first two terms of (1.3) are borrowed. Numerical results are seen to be in excellent agreement with (1.3) as long as $C_D\gtrsim 6$, provided that the computational domain is long enough. Indeed, following the conclusions of Appendix B, results corresponding to $C_D\geq 10^2$ were obtained using the extended domain with a half-length $\mathcal {L}\approx 10^3$, while those corresponding to lower $C_D$ were obtained on the standard domain with $\mathcal {L}=180$. Note that the predictions provided by the two domains match properly throughout the intermediate range $8\lesssim C_D\lesssim 10^2$. In stark contrast with present results, Maxworthy's 1970 experimental data stand beyond the prediction (1.3) as soon as $C_D>15$, the difference being up to $50\,\%$ for large $C_D$. The empirical extrapolation (1.5) supposed to account for axial confinement effects in his device (which had $\mathcal {L}\approx 80$ or $120$ depending on the sphere size) brings only a marginal improvement, leaving a $40\,\%$ over-prediction for $C_D=O(10^3)$. When $C_D$ is low enough, typically $C_D\lesssim 6$, (1.3) is found to overestimate the drag. This regime corresponds to large Reynolds numbers and moderate Rossby numbers. These are the conditions under which the above drag reduction mechanism related to the rotation-induced shortening of the standing eddy operates. The drag modification resulting from this mechanism is obviously not included in the low-${R}e$ asymptotic result (1.2), making the simple drag law (1.3) inaccurate in this regime. Conversely, (1.3) is found to hold even in the moderate-to-large Reynolds number regime provided that the Rossby number is somewhat lower than unity. For instance, setting ${R}o=0.5$ and ${R}e=100$ yields $C_D\approx 7.9$, which is close to the lower limit of validity of (1.3) according to figure 13. This leads to the conclusion that this simple semi-empirical prediction is actually valid well beyond the low-${R}e$ regime within which it is in principle supposed to hold.

Figure 13. Comparison of the measured drag coefficient with the semi-empirical prediction (1.3). The red solid line indicates the prediction (1.3). Symbols $\blacklozenge, \bigstar$ correspond to present simulations with $\mathcal {L}=180$ and $\approx 10^3$, respectively; $\bullet$ corresponds to Maxworthy's experimental data (Reference Maxworthy1970). The dashed line indicates the ‘corrected’ experimental law (1.5). The red and blue colour bars refer to the numerical and experimental data, respectively. In each series, symbols are coloured according to the value of ${T}a$, darkening as ${T}a$ increases.

That present numerical results agree closely with the semi-empirical prediction (1.3) over two decades of $C_D$ (hence ${R}o$) proves that axial confinement effects are responsible for the long-standing but previously unresolved disagreement between experimental results and theoretical models. The restored agreement obtained by considering stationary axisymmetric solutions of the Navier–Stokes equations also rules out the possibility that the problem could be due to non-axisymmetric or unsteady effects, as was suggested previously (Minkov et al. Reference Minkov, Ungarish and Israeli2002). Although Maxworthy rightly identified the origin of the problem, the correction (1.5) that he proposed was biased because it was in a good part based on an extrapolation of his previous data obtained in a much shorter device with $5\lesssim \mathcal {L} \lesssim 10.5$, depending on the sphere size (Maxworthy Reference Maxworthy1968). This extrapolation was not appropriate because end effects in short and long containers do not involve the same mechanisms at all, and therefore do not influence the drag in the same way. In short containers, the direct interaction of the Taylor column with the Ekman layers present along the end walls controls the drag to leading order, making $C_D$ depend on viscosity as (1.4) shows. This is not the case in long containers, in which the end walls produce only a (mostly inviscid) blocking effect that slightly compresses the Taylor column. This difference induces dramatic consequences on the flow structure in the vicinity of the body. For instance, Ungarish & Vedensky (Reference Ungarish and Vedensky1995) showed that, for a thin disc, the upstream recirculation exists only if the ratio $\delta =\mathcal {L}/ {T}a$ is larger than $0.08$. Maxworthy's 1968 experiments were carried out at very large Taylor numbers, ${T}a\geq 2.5\times 10^3$, so that the corresponding data all correspond to $\delta \leq 4\times 10^{-3}$, a regime in which the flow within the Taylor column has little to do with that sketched in figure 1 (for which $\delta =0.93$). Because of these structural differences, there was little chance that an extrapolation mixing two fundamentally different regimes could work.

It is also worth noting that axial confinement effects in Maxworthy's 1970 experiments were actually more severe than can be expected on the basis on the container-to-particle size ratios $\mathcal {L}\approx 80$ and $\mathcal {L}\approx 120$. Indeed, since the drag was obtained by determining the time of flight of rising particles between two sets of horizontal lines, these particles were closer to the bottom wall when the stopwatch was unlocked, and closer to the top wall when it was stopped. Some quantitative details are missing in Maxworthy's description of the experimental protocol. Nevertheless, it may be hypothesized reasonably that the two sets of lines were close to the bottom and upper ends of the ‘viewing box’ that surrounded the middle part of the cylindrical rotating container. With this, it may be estimated that the container length available downstream of the sphere varied over time in the range $48\leq \mathcal {L}\leq 112$ for the large particles with which the large-${T}a$, low-${R}o$ conditions were achieved. Obviously, the length available upstream of the particle followed opposite time variations. Therefore, the actual container-to-particle size ratio that determines the strength of confinement effects rather stood in the range $50\lesssim \mathcal {L}\lesssim 80$ (grey bars in figures 18 and 19). In contrast, the axial confinement does not vary over time in present computations, since the sphere stays midway between the two end ‘walls’ throughout a run. In Appendix B, we examine in two low-${R}o$ cases how the drag varies as the length of the computational domain is increased. Based on these variations, we determined the fit (B1) predicting the artificial drag increase induced by axial confinement effects. This fit may be useful to design or interpret future experiments, although some caution is required, given the differences between the experimental and numerical set-ups.

Figure 14 summarizes the various ‘regimes’ encountered in present simulations in the parameter space $({R}e, {R}o)$, the shaded area sketching the range covered by Maxworthy's 1970 experiments. Three main regions may be identified. Beyond the solid line, inertial effects dominate over those induced by the rigid-body rotation, making the drag depart from the semi-empirical prediction (1.3). Below this line, numerical predictions are in good agreement with (1.3), provided that the computational domain is long enough. This constraint is fulfilled with $\mathcal {L}=180$ in between the solid and dashed lines. Confinement effects become more severe below the latter, i.e. for ${R}o < 0.125$ when ${T}a > 150$, and we had to use the extended domain with $\mathcal {L}\approx 10^3$ to get rid of these effects in that range. Although figure 13 shows that the agreement with (1.3) extends up to the highest Reynolds number considered in the simulations (${R}e=300$) provided that ${R}o$ is low enough, it must again be stressed that the actual flow is no longer axisymmetric at such Reynolds numbers when rotation effects are moderate or low, as the vertical dotted line in figure 14, which corresponds to the transition to three-dimensionality in the limit ${R}o\rightarrow \infty$, reminds us. However, as discussed in § 3.1, three-dimensional effects affect the drag only marginally for ${R}e\approx 150$ in the non-rotating limit. This is why present results in that range (penultimate vertical series of lozenges in figure 14) are still relevant for a comparison with experimental data, and only results corresponding to the three lozenges with the green outline in the rightmost series (${R}e=300$) are expected to be significantly modified by non-axisymmetric effects.

Figure 14. Regimes covered by the present simulations in the parameter space $({R}e, {R}o)$: orange $\blacklozenge$ and black $\bigstar$ correspond to simulations performed on domains with $\mathcal {L} = 180$ and $\mathcal {L}\approx 10^3$, respectively. The shaded area indicates the parameter range covered by Maxworthy's 1970 experiments. The solid line corresponds to the prediction (1.3) for $C_D\approx 6$. The dashed line indicates the limit below which confinement effects are observed with $\mathcal {L}=180$. The dotted line indicates the threshold beyond which the wake is three-dimensional in the limit ${R}o\rightarrow \infty$. Lozenges with a green outline correspond to conditions (${R}e=300$, ${R}o>1$) under which the drag is suspected to be affected significantly by three-dimensional effects.

4.2. Torque

As stated in § 2.1, present computations were carried out by imposing that the sphere rotates at the same rate as the undisturbed flow. Therefore, it experiences a non-zero torque and it is of interest to examine how this torque varies with the flow parameters. Making use of the definitions introduced in § 2.3, especially the spherical coordinate system whose origin stands at the sphere centre, the $z$-component of the torque is

(4.1)\begin{equation} M_z={\boldsymbol{e}}_z\boldsymbol{\cdot}\int_\mathcal{S}a{\boldsymbol{e}}_r\times({\boldsymbol{e}}_r \boldsymbol{\cdot}{\boldsymbol{T}}|_{r=a})\,{\rm d}S, \end{equation}

where $\mathcal {S}$ denotes the sphere surface. Expanding the surface traction ${\boldsymbol {e}}_r\boldsymbol {\cdot }{\boldsymbol {T}}|_{r=a}$ component-wise, (4.1) is found to reduce to

(4.2) \begin{equation} M_z=-\rho\nu a\int_{\mathcal{S}}({\boldsymbol{e}}_z \boldsymbol{\cdot}{\boldsymbol{e}}_{\theta})\,\partial_r U_\phi \,{\rm d}S=2{\rm \pi}\rho\nu a^3 \int_0^{\rm \pi}\left.\sin^2\theta\,(\partial_rU_\phi)\right|_{r=a} {\rm d}\theta. \end{equation}

Noting that $\sin ^2\theta$ is symmetric with respect to the equatorial plane $\theta ={\rm \pi} /2$, it is relevant to expand $U_\phi$ into a component that shares this property, i.e. an even function of $z$, and a component that is antisymmetric with respect to the equatorial plane, i.e. an odd function of $z$ (formally, this could be achieved via Fourier transform, for instance). Only the even component of $U_\phi$ contributes to $M_z$. As such a component results from the downstream advection of the negative upstream axial vorticity (see figures 5a,b,e,f,h,i), $M_z$ is expected to be negative, which in the case of a torque-free sphere would make it rotate slower than the undisturbed fluid. Thus it is appropriate to introduce a torque coefficient $C_T$, related to the axial torque through $M_z=-({{\rm \pi} }/{2})C_Ta^3\rho U_\infty ^2$. Based on this definition and on the above remark, one has

(4.3) \begin{equation} C_T=-4\,{R}e^{-1}\int_0^{\rm \pi}\left.\sin^2\theta\,(\partial_{r/a}U_\phi^e)\right|_{r/a=1} {\rm d}\theta, \end{equation}

where $U_\phi ^e$ stands for the dimensionless even contribution to $U_\phi$.

Outside the boundary layer, the dimensionless thickness of which is denoted $\delta _B$, one has $U_\phi ^e\approx (\sigma /a)\omega _z^e$, where $\omega _z^e$ stands for the (dimensionless) even component of $\omega _z$. Since $U_\phi ^e=0$ at the sphere surface, and $\sigma /a\approx 1$ in the equatorial region, $(\partial _{r/a}U_\phi ^e)|_{r/a=1}\sim \omega _z^e/\delta _B$, so that $C_T\sim {R}e^{-1}\,\omega _z^e/\delta _B$. To determine the scaling laws obeyed by $C_T$, one must consider the governing equation (3.1) for $\omega _z$, keeping in mind that $\delta _B\sim {T}a^{-1/2}$ if ${R}o\ll 1$ and ${T}a\gg 1$, while $\delta _B\sim {R}e^{-1/2}$ in the opposite limit ${R}o\gg 1$, ${R}e\gg 1$. In the latter regime, making use of the near-axis approximation discussed in § 3.2 for the tilting term $-\omega _\sigma \,\partial _\sigma U_z$, the $\omega _z$ balance outside the boundary layer reduces at leading order to $U_\sigma \,\partial _\sigma \omega _z +(U_z+\sigma \,\partial _\sigma U_z)\partial _z\omega _z-\omega _z\partial _zU_z\approx 2\varOmega \,\partial _zU_z$. Since the axial velocity at radial positions $\sigma /a\approx 1$ decreases (resp. increases) as $z$ increases upstream (resp. downstream) of the sphere, the leading-order contribution to $\partial _zU_z$ is an odd function of $z$. Hence the flow past the sphere is dominated by an even component in $U_z$ and, owing to continuity, an odd component in $U_\sigma$. Then, according to the above form of the $\omega _z$ balance, it turns out that the leading contribution to $\omega _z$ is even with respect to $z$. Effects of the Coriolis force do not put a severe restriction on the variations of the flow field in the $z$-direction in that regime. Therefore, outside the boundary layer, $\partial _z \sim a^{-1}\sim \partial _\sigma$, $U_z\sim U_\infty$, $U_\sigma \sim \delta _BU_\infty$ and the $\omega _z$ balance implies $U_\infty \omega _z/a\sim \varOmega U_\infty /a$, i.e. $\omega _z\sim \varOmega$, so that $\omega _z^e\sim {R}o^{-1}$. Hence $\omega _z^e/\delta _B\sim {R}e^{1/2}\,{R}o^{-1}$ and

(4.4)\begin{equation} \left. C_T\right|_{{R}o\gg1,\,{R}e\gg1}\sim{R}o^{-1}\,{R}e^{-1/2}. \end{equation}

Let us now consider the low-${R}o$ limit in which significant axial variations of the flow field exist only within the Ekman layer. The dominant balance for $\omega _z$ then reads $2\varOmega \,\partial _zU_z\approx -\nu \,\nabla ^2\omega _z$. Near the equatorial plane, axial variations at radial positions $\sigma \approx a$ in the Ekman layer take place over distances of the order of the sphere radius, so that $\partial _z\approx a^{-1}$. The axial velocity being of the order of $U_\infty$ at the outer edge of that layer, one has $\partial _zU_z\sim U_\infty /a$. Since $U_z$ is almost symmetric with respect to the sphere's equator, the dominant contribution to the source term in the $\omega _z$ balance is an odd function of $z$, and so is the leading contribution to $\omega _z$, say $\omega _z^oU_\infty /a$. The scaling of $\omega _z^o$ results from the balance $2\varOmega U_\infty /a\sim \nu (U_\infty /a)\omega _z^o/(a\delta _B)^2$, which yields $\omega _z^o\sim 1$. If the Rossby number is small but finite, then advection past the sphere brings a small correction to the $\omega _z$ distribution through the term $\{U_\sigma \,\partial _\sigma \omega _z^o +(U_z+\sigma \,\partial _\sigma U_z)\,\partial _z\omega _z^o-\omega _z^o\,\partial _zU_z\}U_\infty /a$, which is almost an even function of $z$. To balance this term, an even correction to $\omega _z$ is required, say $\omega _z^eU_\infty /a$, and is provided by the corresponding viscous term, $\nu ( U_\infty /a)\,\nabla _{\sigma,z}^2\omega _z^e$. Still in the vicinity of the equatorial plane, $\partial _\sigma \sim (a\delta _B)^{-1}$ in the Ekman layer, and $U_\sigma \sim \delta _B U_\infty$ at its outer edge. Therefore, the above inertial source term is dominated by the contribution $U_\infty (\sigma /a)\,\partial _\sigma U_z\,\partial _z\omega _z^o\sim (U_\infty /a)U_\infty /(a\delta _B)\omega _z^o$, and the above balance implies $U_\infty /(a\delta _B)\omega _z^o\sim \nu \omega _z^e/(a\delta _B)^2$, i.e. $\omega _z^e\sim {R}e\,\delta _B\,\omega _z^o$. According to the scalings obeyed by $\delta _B$ and $\omega _z^o$, this yields $\omega _z^e\sim {R}e\,{T}a^{-1/2}$. Hence $\omega _z^e/\delta _B\sim {R}e$ and

(4.5)\begin{equation} \left. C_T\right|_{{R}o\ll1,\,{T}a\gg1}\sim1, \end{equation}

indicating that the torque coefficient is now independent of the control parameters. Therefore, (4.4) and (4.5) predict that $C_T$ exhibits two different scaling laws, according to the magnitude of the Rossby and Reynolds numbers. In rotation-dominated regimes, where advective effects provide only a small correction to the dominant axial vorticity balance, $C_T$ is constant, whereas it decays with both ${R}o$ and ${R}e$ in advection-dominated regimes.

Figure 15 shows how numerical results compare with the above predictions. As long as ${R}o\,{R}e^{1/2}$ is less than $\approx 5$, the torque coefficient agrees with the prediction (4.5), with $C_T= 0.56 \pm 0.06$. Beyond ${R}o\,{R}e^{1/2}\approx 5$, i.e. when inertial effects are moderate to large, and effects of rotation are moderate or weak, $C_T$ is approximated closely by the decay law $C_T=2.8\,{R}o^{-1}\,{R}e^{-1/2}$, in agreement with (4.4). As the colours of the symbols in figure 15 indicate, the first regime coincides with that in which figure 13 revealed that the drag coefficient agrees well with the prediction (1.3). Conversely, the second regime is that in which $C_D$ departs from this prediction.

Figure 15. Torque coefficient $C_{T}$ as a function of ${R}o\,{R}e^{1/2}$. Symbols $\blacklozenge$ and $\bigstar$ correspond to results obtained on domains with $\mathcal {L}=180$ and $\mathcal {L}\approx 10^3$, respectively. Blue (resp. orange) symbols are associated with configurations in which $C_D<6$ (resp. $C_D\geq 6$).

The numerical results for the torque coefficient may be used to estimate the differential rotation of the sphere, say $\varOmega _s$, required to satisfy the torque-free condition. Indeed, this rotation induces an azimuthal velocity $\varOmega _sa\sin \theta$ at the sphere surface, so that the dimensionless velocity gradient involved in the definition (4.3) of $C_T$ becomes $(\partial _{r/a}U_\phi ^e)|_{r/a=1}\sim (\omega _z^e-({a\varOmega _s}/{U_\infty })\sin \theta )/\delta _B$. The approximate change in the torque coefficient is then $4({R}e\,\delta _B)^{-1}({a\varOmega _s}/{U_\infty })\int _0^{\rm \pi} \sin ^3\theta \,{\rm d}\theta =\frac {16}{3}({R}e\,\delta _B)^{-1}{a\varOmega _s}/{U_\infty }$. Inspection of figure 7$(i)$ allows the distance to the sphere surface at which the swirl velocity $U_\phi /\sigma$ reaches its extremum in the equatorial plane to be determined. This leads to the approximate estimates $\delta _B\approx 2.0\,{T}a^{-1/2}$ and $\delta _B\approx 2.0\,{R}e^{-1/2}$ in the low- and high-${R}o$ regimes, respectively. In the former regime, numerical results showed that $C_T\approx 0.56$, so that the torque-free condition is achieved with ${a\varOmega _s}/{U_\infty }\approx -\frac {3}{16}\times 0.56\,{R}e\,\delta _B\approx -0.2\,{R}e\,{T}a^{-1/2}$. Similarly, in the high-${R}o$ regime, we found $C_T\approx 2.8\,{R}o^{-1}\,{R}e^{-1/2}$, so that ${a\varOmega _s}/{U_\infty }\approx -1.0\,{R}o^{-1}$. Normalized with respect to the imposed rigid-body rotation rate, these estimates become

(4.6a,b)\begin{equation} \left.\frac{\varOmega_s}{\varOmega}\right|_{{R}o\ll1,\,{T}a\gg1}\approx-0.2\,{R}o^{3/2}\, {R}e^{1/2}\quad\mbox{and}\quad \left.\frac{\varOmega_s}{\varOmega}\right|_{{R}o\gg1,\,{R}e\gg1}\approx-1.0\,. \end{equation}

The differential rotation is predicted to be very small in the low-${R}o$ regime, with for instance $\varOmega _s/\varOmega =-1.7\times 10^{-3}$ in the configuration of figure 5(g). In contrast, the sphere is predicted to have a negligible rotation with respect to the laboratory frame when the Rossby and Reynolds numbers are both large (figures 5b,c,f). Of course, these are only rough estimates, since we used crude approximations to evaluate the velocity gradient in (4.3), and the sphere rotation is expected to induce slight modifications in the distribution of the azimuthal vorticity in the sphere vicinity. Let us also mention that, in the limit ${R}e\ll {T}a^{1/2}\ll 1$, Childress (Reference Childress1964) predicted $\varOmega _s/\varOmega \approx -\frac {1}{5}\,{R}e^{2}\,{T}a^{-1/2}$. This prediction differs from those obtained here, in particular in the low-${R}o$ regime, since the first estimate in (4.6a,b) may be rewritten in the form $\varOmega _s/\varOmega \approx -0.2\,{R}e^{2}\,{T}a^{-3/2}$, which corresponds to the same dependence with respect to ${R}e$ but a faster decay with ${T}a$. This is no surprise, as we assume the Taylor number to be large, which makes all processes governing the body rotation controlled by the Ekman layer when ${R}o$ is low. In contrast, Childress’s analysis assumes ${R}e\ll {T}a^{1/2}\ll 1$, so that inertial effects responsible for this differential rotation manifest themselves only at large $O({T}a^{-1/2})$ dimensionless distances from the body.

To check the above prediction for $\varOmega _s$ and assess the influence of this rotation on the drag, we carried out additional simulations corresponding to the torque-free condition. This condition was enforced iteratively, and convergence was considered to be reached when the final torque was less than $1\,\%$ of its stationary value in the case $\varOmega _s=0$. We ran these simulations for the nine cases for which the flow structure is displayed in figure 5. These configurations span the range of conditions considered in this work, especially the two regimes exhibited in figure 15. Results of these simulations are summarized in table 1. It is seen that the sphere rotation has a negligible influence on the drag (i.e. the two values of $C_D$ differ by less than $2\,\%$) in all cases with ${R}o\lesssim 0.5$. This influence is larger with moderate to low rotation levels, as could be expected on the basis of the $O(1)$ values of $\varOmega _s/\varOmega$ predicted by (4.6a,b) in this regime. Nevertheless, the relative difference between the two $C_D$ values never exceeds $10\,\%$. Hence, replacing results of figures 11 and 13 obtained with $\varOmega _s=0$ with those corresponding to the torque-free condition makes no visual difference, even in the most inertial regimes. The estimates (4.6a,b) predict $\varOmega _s/\varOmega =-0.0017$ in case $(g)$ and $\varOmega _s/\varOmega =-1.0$ in case $(c)$, which compares well with the values reported in the top and bottom rows of table 1.

Table 1. Influence of the sphere rotation on the drag. Results are sorted by increasing values of ${R}o\,{R}e^{1/2}$. The labels in the first column refer to the plots in figure 5. For each set of conditions, the drag coefficient in the first row ($C_D^{0}$) was obtained with $\varOmega _s=0$, while that in the second row ($C_D^{tf}$) corresponds to the torque-free condition; $\Delta C_D=(C_D^0/C_D^{tf}-1)\times 100$ is the percentage difference between the two drag coefficients. Cases (g–i) were computed using the extended domain with $\mathcal {L}\approx 10^3$.