Flow past a sphere translating along the axis of a rotating fluid: revisiting numerically Maxworthy's experiments

Abstract

We compute the flow induced by the steady translation of a rigid sphere along the axis of a large cylindrical container filled with a low-viscosity fluid set in rigid-body rotation, the sphere being constrained to spin at the same rate as the undisturbed fluid. The parameter range covered by the simulations is similar to that explored experimentally by Maxworthy (J. Fluid Mech., vol. 40, 1970, pp. 453–479). We describe the salient features of the flow, especially the internal characteristics of the Taylor columns that form ahead of and behind the body and the inertial wave pattern, and determine the drag and torque acting on the sphere. Torque variations are found to obey two markedly different laws under rapid- and slow-rotation conditions. The corresponding scaling laws are predicted by examining the dominant balances governing the axial vorticity distribution in the body vicinity. Results for the drag agree well with the semi-empirical law proposed for inertialess regimes by Tanzosh & Stone (J. Fluid Mech., vol. 275, 1994, pp. 225–256). This law is found to apply even in regimes where inertial effects are large, provided that rotation effects are also large enough. Influence of axial confinement is shown to increase dramatically the drag in rapidly rotating configurations, and the container length has to be approximately a thousand times larger than the sphere for this influence to become negligibly small. The reported simulations establish that this confinement effect is at the origin of the long-standing discrepancy existing between Maxworthy's results and theoretical predictions.

1. Introduction

The spectacular and subtle characteristics of the flow field generated by a rigid or deformable body translating in a rapidly rotating fluid have received much attention for more than a century, starting with the landmark investigations of Proudman (1916) and Taylor (1917). This configuration, which shares similarities with flows in stratified or magnetized fluids, is of practical relevance in problems where particles, drops or bubbles settle or rise in locally rotating flows, such as in the dynamics of rapidly rotating suspensions or in centrifugal separation techniques employed in two-phase flows (Ungarish 1993; Bush, Stone & Tanzosh 1994). It is also relevant in ocean and atmosphere dynamics (Loper 2001) and, combined with thermal or compositional convection, in astrophysics to understand the dynamics of liquid cores in terrestrial and rapidly rotating planets (Bush, Stone & Bloxham 1992; Cheng et al. 2015).

The flow disturbance generated by a rigid axisymmetric body with equatorial radius $a$ moving at speed $U_{\infty }$ in a Newtonian fluid of kinematic viscosity $\nu$ rotating as a whole with angular velocity $\varOmega$ depends on the Taylor number ${T}a\equiv {\varOmega a^2}/{\nu }$ and the Rossby number ${R}o \equiv {U_{\infty }}/{\varOmega a}$ (or equivalently the Reynolds number ${R}e \equiv {U_{\infty }a}/{\nu }= {R}o\,{T}a$). Pioneering experiments with a cylinder or a sphere translating in a viscous fluid set in rigid-body rotation were performed by Taylor, with the body translating either parallel to the rotation axis (Taylor 1922) or perpendicular to it (Taylor 1923). These experiments revealed the existence of slender recirculating fluid regions, later referred to as Taylor columns, confined within a cylinder circumscribing the body and having their generators parallel to the rotation axis. Later, Maxworthy repeated Taylor's 1922 experiments with a sphere translating along the rotation axis over a broad range of ${T}a$ at both low Reynolds number (${R}e \lesssim 0.5$; Maxworthy 1965) and moderate-to-large Reynolds numbers: $5\lesssim {R}e \lesssim 100$ (Maxworthy 1968), and $3\lesssim {R}e \lesssim 300$ (Maxworthy 1970). He confirmed Taylor's observations regarding the typical features of the flow structure, and found that the drag force on the sphere is generally increased by the fluid rotation, this increase scaling linearly with the Taylor number once the drag force has been normalized by the Stokes drag.

A sketch of the corresponding flow at a relatively large Taylor number (${T}a=193$) is depicted in figure 1. No fore–aft symmetry with respect to the sphere equator exists in this case, as advective effects are large (${R}e=93$). Two prominent recirculation regions standing upstream and downstream of the sphere may be observed. Existence of such recirculation regions in the present case is in line with the predictions of Tanzosh & Stone (1994) and Vedensky & Ungarish (1994), the latter for a disc, which indicate that these structures take place when ${T}a\gtrsim 50$, and their axial extent (normalized by the body radius) grows approximately as $0.052\, {T}a$. The second noticeable feature is the nearly geostrophic region in which the Taylor–Proudman theorem applies approximately (Moore & Saffman 1968). In this smaller region, the non-dimensional length of which is approximately $0.006\, {T}a$ (Tanzosh & Stone 1994), the fluid almost achieves a rigid-body rotation, the rotation rate being faster (resp. slower) than $\varOmega$ downstream (resp. upstream) of the body. This nearly uniform swirling motion is accompanied by a weak plug-like axial flow thanks to which a tiny flux is transmitted from one nearly geostrophic region to the other via the Ekman boundary layer surrounding the body. The last salient flow feature is the Stewartson layer that connects the outer flow to the Taylor column (which is the body of fluid made of the recirculation and nearly geostrophic regions, and the above Ekman boundary layer). In the Stewartson layer, which has a complex internal ‘sandwich’ structure made of three concentric sublayers, the thicknesses of which obey different scaling laws, an intense axial motion takes place while the swirl velocity varies rapidly in the radial direction (Baker 1967; Moore & Saffman 1969). This layer is the main region through which the fore and aft Taylor columns exchange fluid when the container is long enough for the end walls not to interact dynamically with these columns.

Figure 1. Qualitative flow structure past a sphere translating along the axis of a large fluid container set in rigid-body rotation. The flow is observed in the reference frame translating with the sphere. The thin white lines are streamlines obtained from a simulation at ${R}e=93$ and ${T}a=193$, i.e. ${R}o=0.48$.

Numerous studies have attempted to characterize the influence of the rigid-body rotation on the drag experienced by the sphere, in both finite-length and infinitely long containers. Stewartson (1952) considered the asymptotic limit of an impulsive but slow motion in an inviscid flow and an infinitely long container. Using a Laplace transform technique, he predicted that the drag force $F_D$ is

(1.1)\begin{equation} \frac{F_D}{F_{St}}=\frac{8}{9{\rm \pi}}\,{T}a \quad \Leftrightarrow \quad C_D=\frac{32}{3{\rm \pi}}\,{R}o^{-1}\approx 3.4\,{R}o^{-1} \quad ({T}a= \infty,\ {R}o\rightarrow0), \end{equation}

where $F_{St}=6{\rm \pi} \rho \nu a U_{\infty }$ stands for the Stokes drag (with $\rho$ the fluid density), and the drag coefficient $C_D$ is defined through the relation $F_D=\frac {1}{2}C_D{\rm \pi} a^2\rho U_{\infty }^2$. The above result was later confirmed by Moore & Saffman (1969) assuming small-but-finite viscous effects. Conversely, Childress (1964) considered the viscous regime and assumed ${R}e\ll {T}a^{1/2}\ll 1$. Making use of the matching asymptotic expansion technique, he obtained

(1.2)\begin{equation} \frac{F_D}{F_{St}}=1+\frac{4}{7}\,{T}a^{1/2} \quad \Leftrightarrow \quad C_D=\frac{12}{{R}e}+\frac{48}{7}\,({R}e\,{R}o)^{-1/2} \quad ({R}e\ll1, {T}a\ll1,\ {T}a/{R}e^2\gg1). \end{equation}

Interestingly, Childress's theory also predicts that the drag is smaller than that in a non-rotating fluid when ${T}a/{R}e^2\lesssim 0.2$, the largest reduction being ${\approx }5\,\%$ for ${T}a/{R}e^2\approx 0.09$. Later, Weisenborn (1985) and Tanzosh & Stone (1994) predicted the drag for arbitrary Taylor numbers, still assuming the Reynolds number to be negligibly small. While both used distinct approaches (the so-called ‘induced-force’ method and a boundary integral technique, respectively), the two sets of results are in agreement within 0.5 % up to ${T}a=10^4$, and both agree within 5 % with the semi-empirical law proposed by Tanzosh & Stone (1994), namely

(1.3)\begin{equation} \frac{F_D}{F_{St}}=1+\frac{4}{7}\,{T}a^{1/2}+\frac{8}{9{\rm \pi}}\,{T}a \quad \Leftrightarrow \quad C_D=\frac{12}{{R}e}+\frac{48}{7}\,({R}e\,{R}o)^{-1/2}+\frac{32}{3{\rm \pi}}\,{R}o^{-1} \quad ({R}e \ll1). \end{equation}

The prediction (1.3) is nothing but the linear combination of (1.1) and (1.2). Independently, Vedensky & Ungarish (1994) used a system of dual integral equations to predict the drag on a disc under similar conditions. Influence of the finite length of the container was considered by Moore & Saffman (1968), assuming small-but-finite viscous effects and neglecting inertial effects. Considering a container with rigid ends and a half-length $H$ such that $1\ll \mathcal {L}\equiv H/a\ll {T}a^{1/2}$, they showed that

(1.4) \begin{equation} \frac{F_D}{F_{St}}=\frac{43}{630}\,{T}a^{3/2} \quad \Leftrightarrow \quad C_D=\frac{86}{105}\,{R}o^{-1}\,{T}a^{1/2} \quad ({T}a \to \infty,\ {R}o\rightarrow 0). \end{equation}

Recently, Kozlov et al. (2023) performed experiments with a sphere rising in a rapidly rotating short container ($\mathcal {L}=9.4$) in the range ${T}a \in [250, 2.5\times 10^{4}]$, ${R}o \in [10^{-4}, 10^{-2}]$, and confirmed the ${T}a^{3/2}$ dependence predicted by (1.4). In this ‘short-container’ limit, the Ekman layers that develop along the two end walls interact directly with the Taylor columns and ensure a good part of the fluid transport between the fore and aft columns, making the drag coefficient depend on viscosity (through the Taylor number), in contrast to the ‘long-container’ limit. In the latter, characterized by container aspect ratios such that $\mathcal {L}\gg {T}a^{1/2}$, the radial flow in these two Ekman layers is very weak and plays no role. However, the end walls may still influence the internal structure of the Taylor columns through a purely kinematic ‘blocking’ effect, thereby modifying the drag. For this reason, Hocking, Moore & Walton (1979) considered finite values of the ratio $\delta =\mathcal {L}/{T}a$ (still in the limit on negligibly small Rossby numbers) and concluded that the drag increases monotonically as $\delta$ is reduced. For instance, when normalized by the prediction (1.1) corresponding to $\delta \rightarrow \infty$, they found that the drag on a sphere standing midway between the end walls increases by approximately $9\,\%$ for $\delta =1$, and $30\,\%$ for $\delta =1/4$.

The low-Reynolds-number drag measurements (${R}e \lesssim 0.5$, ${T}a \in [0.05, 0.7]$) carried out by Maxworthy (1965) agree within a few percent with (1.2). It is worth noting that these data also support Childress’s prediction that, at low enough ${T}a/{R}e^2$, the drag is smaller than that in a non-rotating fluid. Conversely, at large enough Reynolds and Taylor numbers (${R}e \in [3, 300]$, ${T}a \in [10, 450]$), the data reported later by the same author (Maxworthy 1970) follow the scaling (1.1), albeit with a significantly larger pre-factor. Based on the comparison between (1.1) and (1.4), Maxworthy suspected that the origin of the discrepancy may stand in the finite length of his container, which was such that $\mathcal {L}\approx 80$ or $\mathcal {L}\approx 120$, depending on the size of the particles used. Hence he corrected his results from end-wall effects using supplementary data, some of which, reported in Maxworthy (1968), were obtained in a much shorter container ($5\lesssim \mathcal {L}\lesssim 10$). Based on this correction, he concluded that his data may be extrapolated to an infinitely long container in the form

(1.5)\begin{equation} C_D = (5.2\pm 0.1) \times {R}o^{-1 \pm 0.01}. \end{equation}

However, the pre-factor involved in (1.5) is still nearly $50\,\%$ larger than that in (1.1). This discrepancy motivated the aforementioned extension of (1.1) to finite-length containers. However, the corresponding correction was found to reduce the discrepancy only slightly, making Hocking et al. (1979) conjecture that finite-${R}o$ effects not accounted for in their theory cannot be ignored.

The very first simulations of the problem under consideration based on the full Navier–Stokes equations, hence incorporating finite-${R}o$ effects, were carried out by Dennis, Ingham & Singh (1982). Computational limitations at that time restricted the explored parameter range to ${R}e\leq 0.5$ and ${T}a\leq 0.5$. Nevertheless, these simulations were able to confirm quantitatively the experimental findings of Maxworthy (1965) regarding the increase in drag with ${T}a$ in the range $0\leq {T}a\leq 0.5$. Rao & Sekhar (1995) explored a much broader range of Reynolds number (up to ${R}e=500$) but considered only Rossby numbers larger than $2$. They could observe the changes in the flow structure in the presence of moderate rotation effects, especially the shrinking and disappearance of the standing eddy at the back of the sphere when ${R}e\gtrsim 20$ and ${R}o$ is decreased from $O(10)$ to $O(1)$ values. They found that in this moderate-${R}o$, moderate-to-large-${R}e$ regime, rotation effects reduce the drag, a finding also noticed by Maxworthy (1970) and later reconfirmed numerically by Sahoo et al. (2021). Minkov, Ungarish & Israeli (2000, 2002) considered the case of a circular disc rising under low-${R}o$ conditions in short and long containers, respectively. They confirmed that the relative height of the container deeply affects the drag force. They also investigated the influence of the advective terms, i.e. finite-${R}o$ corrections, by exploring (in the long-container case) the range ${R}o\leq 0.25$ with ${T}a \in [100, 200]$, i.e. ${R}e \leq 50$. They concluded that these effects actually reduce the drag, thus further increasing the discrepancy with Maxworthy's data. Wang, Lu & Zhuang (2004) performed three-dimensional simulations of the same configuration for a sphere with or without a differential spin for ${R}e=100$ and $250$, and ${T}a \in [50, 6.25\times 10^3]$. They confirmed the characteristic features of the flow structure sketched in figure 1 at low Rossby number, and examined the influence of the control parameters on the inertial waves pattern. However, they did not report any drag value. Therefore, full Navier–Stokes simulations have not helped so far to reconcile the experimental findings of Maxworthy (1970) in the low-${R}o$ regime with theoretical predictions (1.1) or (1.3) for the drag. This is why the conclusion of Minkov et al. (2002) that ‘in any case, the major discrepancy between theory and experiments concerning the value of the drag force remains unresolved, and becomes even more puzzling in view of the present results’ still holds.

This intriguing and unexplained discrepancy was the main initial motivation for the present work. We use fully resolved simulations to get new insight into this issue, and more generally into the influence of rigid-body rotation, and viscous and advective effects on the organization of the flow past the body. The sphere is assumed to rotate at the same rate as the undisturbed flow, and we determine the corresponding drag force and torque, assessing the possible influence of axial confinement effects on the flow structure and the loads on the body. We consider Taylor numbers ${T}a \in [20, 450]$ and Reynolds numbers ${R}e \in [2, 300]$, yielding Rossby numbers in the range ${R}o \in [5\times 10^{-3}, 10]$, which corresponds to the parameter range covered in Maxworthy's 1970 experiments. The mathematical problem, the numerical set-up and a preliminary comparison with zero-${R}o$ results are presented in § 2. Characteristic features of the flow structure are discussed and compared with previous findings in § 3. Then the variations of the drag and torque with the control parameters are analysed in § 4. The main outcomes of the study and some avenues for future work are presented in § 5.

2. Problem formulation and numerical set-up

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 1997).

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 (1995) 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 1998) 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 (2019) 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. (1995). The detailed form of the governing equations expressed in this general coordinate system is also provided in Magnaudet et al. (1995). 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 (2007) and Auguste & Magnaudet (2018). Here, following Magnaudet et al. (1995), 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. (1995), Legendre & Magnaudet (1998) and Legendre, Magnaudet & Mougin (2003), 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 (1994), 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 (1994). 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 (1994) 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 (1994) for ${T}a=50$ and ${T}a=500$, respectively.

3. Flow field

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 (1970). 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 1993; Johnson & Patel 1999). In the presence of moderate rotation effects (${R}o=2$), Wang et al. (2004) 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. (2018) 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 2000). 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 2000; Poon et al. 2014) – 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 (1994) 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 2020; Torres et al. 2000). 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 (1994) 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 (1994). 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 (1994), 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 (1970), 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 1994) 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 (1994) in the zero-${R}o$ limit. We attributed this discrepancy to axial confinement effects, a track already suggested by Ungarish & Vedensky (1995). 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 (1970); $\triangleleft$ corresponds to boundary-integral simulations at ${R}o=0$ (Tanzosh & Stone 1994). The solid black line indicates the zero-${R}o$ limit $\ell _{s}=0.052\, {T}a$ (Tanzosh & Stone 1994). 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 1968). Taylor (1922) 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 1969). 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. (2015, 2018) using particle image velocimetry. We are not aware of similar experimental observations in the configuration considered here. The numerical investigation of Wang et al. (2004) 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 1967; Whitham 1974)

(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 (1922). 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 1982; Johnson 1982).

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. Loads on the body

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 1933), 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 (1970). 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 (1970) 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 (1995) and Sahoo et al. (2021). 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 (1970) 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 (1994). 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 (1970). 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. 2002). 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 1968). 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 (1995) 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 (1964) 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$.

5. Summary and concluding remarks

With the aid of numerical simulations, we revisited the classical problem of a rigid sphere translating steadily along the axis of a rotating container filled with a slightly viscous fluid. Assuming the flow to be axisymmetric and the sphere to rotate at the same rate as the container, we considered a large number of combinations in the range ${T}a \in [20, 450]$ and ${R}e \in [5, 300]$, covering the Rossby number range ${R}o \in [10^{-2}, 10]$. These conditions correspond to those explored experimentally by Maxworthy in his 1970 reference study (Maxworthy 1970).

Although the problem looks easy from a numerical point of view by today's standards, it is actually challenging regarding the computational domain and the discretization grid. The reason is that the flow has to be captured accurately both in the thin Ekman boundary layer surrounding the body and over very long distances upstream and downstream of it in the near-axis region corresponding to the Taylor column. This is presumably the reason why it took half a century to repeat Maxworthy's experiments on a computer. We dealt with this technical issue by making use of a boundary-fitted orthogonal curvilinear grid that combines the advantages of spherical coordinates in the sphere vicinity with those of cylindrical coordinates far from it.

Thanks to the design of this grid, the characteristics of the flow could be examined in detail throughout the desired parameter range, and several quantities were compared quantitatively with available predictions. In particular, we could observe the inertial wave pattern radiated by the sphere, and check that the associated wavelength agrees well with the inviscid theoretical prediction $\varLambda ={\rm \pi} \,{R}o$ (Taylor 1922). We also examined how the characteristics of the flow within the Taylor column vary with the control parameters. In particular, we found that, for ${T}a\gtrsim 100$, the length of the upstream recirculation region follows the law $\ell _s=0.052\,{T}a$ established by Tanzosh & Stone (1994) in the zero-${R}o$ limit. Using horizontal slices of the three velocity components at various altitudes, we could also clarify some interesting low-${R}o$ mechanisms, such as that leading gradually to a plug-like distribution of the angular swirl as one moves away axially from the body through the nearly geostrophic and recirculation regions. Slices in the equatorial plane also helped to highlight some consequences of inertial effects that break the symmetries inherent to the zero-Rossby-number limit. While these symmetries impose that the radial and azimuthal velocities are zero in that plane at ${R}o=0$, we found that these components develop large negative peaks within the Ekman layer, with respective minima of the order of $20\,\%$ and $60\,\%$ of the sphere speed for ${R}o\approx 0.5$. Conversely, these effects reduce drastically the magnitude of the large positive peak of the axial velocity encountered in that layer in the ${R}o=0$ limit, dividing it by a factor of two for ${R}o\approx 0.5$, which of course has direct consequences on the fluid exchange between the fore and aft Taylor columns.

We determined the drag experienced by the sphere for a large number of $({T}a, {R}e)$ sets, and performed a systematic comparison of the numerical results with Maxworthy's 1970 data and available predictions. Comparing low-${R}o$, large-${T}a$ results obtained on computational domains having $\mathcal {L}=O(10^2)$ with the zero-${R}o$ predictions of Tanzosh & Stone (1994) based on a boundary-integral approach (hence an infinite domain) made it clear that axial confinement effects dramatically enhance the drag in this regime, owing to the slight changes they induce in the structure of the Taylor column. To get rid of almost all of this undesired influence, we designed a grid with $\mathcal {L}=O(10^3)$ on which the drag coefficients were found to agree with the zero-${R}o$ prediction within a few percent. Hence this extreme sensitivity of the drag to axial confinement effects is the reason why Maxworthy's 1970 data (obtained in a container with $\mathcal {L}\approx 80$) stand systematically and significantly beyond theoretical predictions. Once these effects are eliminated, the drag coefficient agrees well with the semi-empirical law (1.3) that accounts for the combined effects of rotation, viscosity and weak inertia. Actually, the domain of validity of (1.3) was found to extend throughout the range of conditions under which $C_D\gtrsim 6$. Hence we could conclude that (1.3) is valid up to ${R}e= O(10^2)$, provided that rotation effects are large enough. In contrast, (1.3) overestimates the drag when inertial effects are ‘too’ dominant. Remarkably, in this high-${R}e$ and moderate-to-large-${R}o$ regime, the drag is also overestimated by the standard law corresponding to a sphere translating in a fluid at rest. The reason for this could be ascribed to the influence of (weak) rotation effects on the azimuthal vorticity in the sphere wake. Rotation contributing to increase this vorticity component in that region, it weakens the negative axial fluid velocity within the standing eddy, and therefore reduces the pressure drag, a scenario confirmed by the numerical velocity and pressure distributions at the back of the sphere.

Since the sphere was assumed to rotate at the same rate as the undisturbed fluid, we could determine the torque that it experiences. It turned out that this torque obeys two different scaling laws, depending on the flow regime. The torque coefficient is constant when rotation effects are dominant, more precisely as long as ${R}o\,{R}e^{1/2}\lesssim 5$. In contrast, this coefficient decays as ${R}o^{-1}\,{R}e^{-1/2}$ in inertia-dominated regimes. Interestingly, the conditions corresponding to the transition between the two scalings coincide with the threshold below which the drag law (1.3) ceases to be valid, i.e. $C_D\approx 6$. The two scaling laws were rationalized by examining the axial vorticity balance in the sphere vicinity and the associated symmetries with respect to the equatorial plane, from which the dominant scalings governing the symmetric vorticity component, which originates in advective effects, could be determined. Numerical results for the torque were used to infer the differential rotation of the sphere achieving the torque-free condition. It was found that the differential rotation rate, normalized by the rotation rate of the outer fluid, scales as ${R}o^{3/2}\,{R}e^{1/2}$ in the rotation-dominated regime, while it becomes constant in inertia-dominated regimes. Some simulations were carried out under the torque-free condition. They revealed virtually no influence of the sphere rotation on the drag as long as the Rossby number is less than unity, and a modest influence, with relative differences $\lesssim$10 %, in the most inertial regimes.

The present work calls for several extensions in at least three directions. First, three-dimensional effects were deliberately ignored here. Although their influence on the drag is presumably marginal in the parameter range that we explored, this secondary effect is worth quantifying. From a more fundamental point of view, determining the threshold ${R}e_c({R}o)$ beyond which the wake becomes three-dimensional, the nature of the corresponding bifurcation, and the spatial structure of the first three-dimensional mode, would be a significant addition to the current knowledge concerning high-${R}e$, low-to-moderate-${R}o$ flows past axisymmetric bodies. Numerical tools designed to perform global linear stability analysis in axisymmetric open flows are now mature and could be adapted easily to tackle this problem.

Second, although only results concerning the steady-state configuration were reported here, transient regimes are also worthy of investigation. In particular, examining how the flow structure and the drag change when the rotation rate is suddenly increased or decreased at a given Reynolds number is a relevant question to predict transient effects in rapidly rotating suspensions and centrifugation processes. Since such a variation induces a change in the drag, the settling or rise speed of the particle also varies. The force balance governing the velocity of particles moving under time-dependent conditions in a viscous fluid is usually split into several distinct contributions, although this splitting is justified rigorously only under creeping-flow conditions. Besides the net body weight and the steady (or quasi-steady) drag, one then finds an added-mass force that opposes the relative acceleration between particle and fluid, and a history force resulting from the unsteady transport of vorticity past the particle. The added-mass effect being due to the no-penetration of the fluid across the body surface, the corresponding force depends only on the instantaneous relative acceleration and on the body shape. Hence, for a given acceleration, it is unaffected by rotation effects, a conclusion that we could confirm numerically (Aurégan 2020). Things are different regarding the history contribution, the evolution of which depends on the past history of the relative acceleration weighted by a time-dependent kernel. This kernel expresses the way a change in the vorticity at the particle surface propagates in the flow under the combined influence of viscosity, inertia and possible non-conservative forces, here the Coriolis force. As such, this kernel is expected to depend on the Rossby number. In the aforementioned preliminary investigation, we could verify that this is indeed the case. Therefore, a systematic study of history effects in the presence of rigid-body rotation appears to be an important objective for future work. Such an investigation should presumably combine a theoretical approach in the zero-${R}o$ limit with numerical simulations to explore the influence of finite advective effects.

Last but not least, drops and bubbles offer challenging additional questions. The theoretical investigations of Bush et al. (1992) (in short containers) and Bush, Stone & Bloxham (1995) (in both short and long containers) performed in the zero-${R}o$ limit provide interesting insights into the effects of the drop-to-fluid viscosity ratio and the centrifugal-to-surface-tension force ratio (the so-called rotational Bond number). The drops were shown to take prolate shapes due to the centrifugal force; the larger the rotational Bond number, the more the drop elongates along the rotation axis. Remarkably, in long containers, the rise or settling speed was predicted to be nearly independent of the drop viscosity and detailed shape, and to depend essentially on its equatorial radius. The reason is that the drag results directly from the efficiency with which the fluid is transported from the fore to the aft Taylor column, and in long containers, this transport takes place mostly through the Stewartson layer rather than via the Ekman boundary layer. How these features are modified by advective effects is currently essentially unknown, but these effects are suspected to be in good part responsible for the significant overestimate of the drag predicted in the zero-${R}o$ approximation, compared to experimental data. Ungarish (1996) introduced a ‘quasi-geostrophic’ approximation incorporating some finite inertial corrections to remedy this problem, but this refinement reduced the disagreement only slightly. These are some of the reasons why the investigation carried out here should be repeated with drops and bubbles. Although this is technically challenging, a variety of numerical approaches now allow the efficient and accurate treatment of boundary conditions at a deformable interface with finite surface tension. Hence exploring how the zero-${R}o$ findings are altered by the presence of finite advective effects appears as an exciting and reachable continuation of the present work.