Experimental study of the flows in a non-axisymmetric ellipsoid under precession

Abstract

Precession driven flows are of great interest for both, industrial and geophysical applications. While cylindrical, spherical and spheroidal geometries have been investigated in great detail, the numerically and theoretically more challenging case of a non-axisymmetric cavity has received less attention. We report experimental results on the flows in a precessing triaxial ellipsoid, with a focus on the base flow of uniform vorticity, which we show to be in good agreement with existing theoretical models. As predicted, the uniform vorticity component exhibits two branches of solutions leading to a hysteresis cycle as a function of the Poincaré number. The first branch is observed at low forcing and characterized by large amplitude of the total fluid rotation and a moderate tilt angle of the fluid rotation axis. In contrast, the second branch displays only a moderate fluid rotation and a large tilt angle of the fluid rotation axis, which tends to align with the precession axis. In addition, we observe the occurrence of parametric instabilities early in the first branch, which saturate in the second branch, where we observe the same order of the kinetic energy in the base flow and instabilities.

1. Introduction

The term precession denotes the slow, gyroscopic motion of a spinning object, resulting from a torque that tends to tilt the object's rotation axis. Precessing, fluid-filled containers are found in industrial contexts, for example as soft mixers in bio-engineering (Meunier 2020) or in spacecraft with a liquid payload, for which minimization of mechanical energy dissipation is crucial for stability (Vanyo & Likins 1971). In the context of geo- and astrophysical fluid dynamics, it has been suggested that precession driven flows could act as a stirring mechanism for planetary dynamos (Bullard 1949; Malkus 1968) or could be responsible for power dissipation in liquid planetary cores and subsurface oceans (Yoder & Hutchison 1981; Williams et al. 2001; Lin, Marti & Noir 2015; Cébron et al. 2019).

The problem of flows inside a precessing spheroidal container has been studied theoretically since the end of the 19th century, starting with the seminal works of Hough (1895), Sloudsky (1895) and Poincaré (1910). Under the assumption that the vorticity is uniform and steady in the frame of precession, they derived an inviscid solution in the form of a tilted solid body rotation, which is complemented by a gradient flow. This solution is often referred to as Poincaré flow and is closely related to the so-called spin-over mode, a solid body rotation around an equatorial axis. Later, viscosity was reintroduced through the viscous torque in the boundary layer (Busse 1968; Noir et al. 2003; Kida 2020) and the main characteristics of the viscous solutions have been confirmed experimentally and numerically for spherical and spheroidal cavities (e.g. Vanyo et al. 1995; Noir et al. 2001a; Noir, Jault & Cardin 2001b; Tilgner & Busse 2001; Noir et al. 2003). While the viscous solution in the sphere is always unique, in spheroids two stable solutions may be found over a finite range of precession rates, leading to a hysteresis cycle (Noir et al. 2003; Cébron 2015; Nobili et al. 2020).

In non-axisymmetric, i.e. triaxial, ellipsoids, the flow of uniform vorticity has been investigated theoretically and numerically by Noir & Cébron (2013) and experimentally and numerically by Cebron, Le Bars & Meunier (2010). The former treat the true geophysical problem of a solid container with a fixed shape in the rotating frame, while the latter consider a deformable container with a shape fixed in the frame of precession, yet both approaches share similar dynamics. As for the spheroidal geometry, multiple solutions are expected over a range of precession rates, but have not yet been observed either experimentally or numerically. In contrast to the solutions in an axisymmetric spheroid, the uniform vorticity flow in a triaxial ellipsoid is not steady, which renders the analytical treatment more cumbersome.

The base flow of uniform vorticity is prone to instabilities, in both the boundary layer and the bulk of the fluid. Instabilities of the Ekman layer are observed when the local Reynolds number $Re_{bl}$, defined on the boundary layer thickness and a proxy for the free-stream velocity, exceeds a critical value, typically between 50 and 100 (Lorenzani & Tilgner 2001). Also, Sous, Sommeria & Boyer (2013) derived and validated experimentally the onset criteria for the primary instability ($Re_{bl}\gtrsim 55$) and the subsequent transition to turbulence ($Re_{bl}\gtrsim 150$) for a steady Ekman layer. However, more recently Buffett (2021) argued that a boundary layer driven by precession should be treated as a time dependent Ekman layer, for which their numerical simulations yield turbulence only above $Re_{BL}>500$. Concerning bulk instabilities, the irrotational component of the base flow and its viscous correction associated with the critical latitude in the boundary layer can drive bulk parametric instabilities, eventually leading to space filling turbulence (Malkus 1968; Kerswell 1993; Goto et al. 2007; Lin et al. 2015; Kida 2020). However, instabilities due to the elliptical distortion of the Poincaré solution (uniform vorticity base flow) are also observed in the purely inviscid case, as reported by Roberts & Wu (2011). Finally, Giesecke et al. (2019) recently proposed that transition to turbulence in a precessing cylinder may be triggered by a centrifugal instability of the background zonal flow profile. While there is now a considerable understanding of the destabilizing mechanisms, the saturation of those instabilities and the associated kinetic energy remain poorly understood, in particular the fundamental nature of the turbulence, quasi-geostrophic or three-dimensional wave turbulence, remains unpredictable for a given range of parameters (Le Reun, Favier & Le Bars 2019).

Here, we experimentally investigate the flows driven in a precessing triaxial ellipsoid for which we present the theoretical foundations in § 2. Using the experimental device introduced in § 3, we first investigate the uniform vorticity component of the flow (§ 4). In the second part we focus on instabilities (§ 5) paying particular attention to the distribution of the kinetic energy between base flow and instabilities. Finally, we draw some conclusion on the applicability of our findings for planetary core dynamics and industrial applications in § 6.

2. Problem formulation

2.1. Governing equations and uniform vorticity base flow

Let us consider a precessing triaxial ellipsoid with semi-major axes $a\neq b \neq c$, filled with an incompressible fluid of uniform density $\rho$ and kinematic viscosity $\nu$. The ellipsoidal cavity is spinning around the $c$-axis as $\boldsymbol {\varOmega }_{\boldsymbol {s}}=\hat {\boldsymbol {k}}_{\boldsymbol {s}}\varOmega _s$, which itself precess as $\boldsymbol {\varOmega }_{\boldsymbol {p}}=\hat {\boldsymbol {k}}_{\boldsymbol {p}}\varOmega _p$ fixed in the inertial frame of reference (see figure 1a). Here, $\hat {\boldsymbol {k}}_{\boldsymbol {s}}$ denotes the rotation vector and $\hat {\boldsymbol {k}}_{\boldsymbol {p}}$ the precession vector. For the description of the problem we choose a right-handed Cartesian coordinate system in which the $z$-axis is always aligned with the spin axis along $c$. Using $\varOmega _s^{-1}$ as the time scale and $R=(abc)^{1/3}$ as the length scale, the dimensionless Navier–Stokes equations in the frame of the container are

(2.1)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + 2(\hat{\boldsymbol{k}}_{\boldsymbol{s}} + Po \hat{\boldsymbol{k}}_{\boldsymbol{p}}) \times \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-} \boldsymbol{\nabla} p - Po(\hat{\boldsymbol{k}}_{\boldsymbol{p}} \times \hat{\boldsymbol{k}}_{\boldsymbol{s}}) \times \boldsymbol{r} + E \boldsymbol{\Delta} \boldsymbol{u}, \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

where $\boldsymbol {u}$ is the flow velocity and $p$ is the reduced pressure including centrifugal forces. The two non-dimensional numbers in (2.1), are the Poincaré number $Po = \varOmega _p/\varOmega _s$, that measures the strength of the precession forcing and the classical Ekman number $E = \nu /\varOmega _s R^{2}$, with $\nu$ the kinematic viscosity.

Figure 1. (a) Sketch of the experimental device. (b) Close-up photography of the container and definition of the coordinate system. The $x$-axis is pointing towards the reader.

Following Hough (1895), Sloudsky (1895) and Poincaré (1910), let us search for a solution in the form of a uniform vorticity component written in the frame of the container as

(2.3)\begin{equation} \boldsymbol{U} = \boldsymbol{\omega} \times \boldsymbol{r} + \boldsymbol{\nabla} \psi, \end{equation}

where $\boldsymbol {\omega } \times \boldsymbol {r}$ represents a rigid body rotation, with r the position vector $(x, y, z)$ and ω the fluid rotation vector ($\boldsymbol{\omega}_x$, $\boldsymbol{\omega}_y$, $\boldsymbol{\omega}_z$) and $\boldsymbol {\nabla } \psi$ the irrotational stretching term necessary to enforce the non-penetrating boundary condition $\boldsymbol {u}\boldsymbol {\cdot }\hat {\boldsymbol {n}} = 0$ at the wall of the container. Using a set of coordinates fixed in the frame attached to the container with $\hat {\boldsymbol {x}}$ along $a$, $\hat {\boldsymbol {y}}$ along $b$ and $\hat {\boldsymbol {z}}$ along $c$, the velocity field (2.3) is expressed as

(2.4)$$\begin{gather} U_x = \omega_y \frac{2a^2}{a^2+c^2}z - \omega_z \frac{2a^2}{a^2+b^2}y, \end{gather}$$

(2.5)$$\begin{gather}U_y = \omega_z \frac{2b^2}{a^2+b^2}x - \omega_x \frac{2b^2}{c^2+b^2}z, \end{gather}$$

(2.6)$$\begin{gather}U_z = \omega_x \frac{2c^2}{b^2+c^2}y - \omega_y \frac{2c^2}{a^2+c^2}x. \end{gather}$$

Substituting this ansatz in (2.1), Noir & Cébron (2013) derived the equations governing the evolution of the three components of the rotation vector $\boldsymbol {\omega }$ in the frame of reference attached to the container, which we present here for a triaxial ellipsoid precessing at an angle of $90^\circ$:

(2.7)\begin{align} \frac{\partial \omega_x}{\partial t}&= \left[\frac{2a^{2}}{a^{2}+c^{2}}-\frac{2a^{2}}{a^{2}+b^{2}}\right]\omega_z\omega_y + Po \sin(t) \frac{2a^{2}}{a^{2}+b^{2}}\omega_z\nonumber\\ &\quad + \frac{2a^{2}}{a^{2}+c^{2}}\omega_y+ Po \sin(t)+\mathcal{L}\varGamma_x, \end{align}

(2.8)\begin{align} \frac{\partial \omega_y}{\partial t} &= \left[\frac{2b^{2}}{a^{2}+b^{2}}-\frac{2b^{2}}{b^{2}+c^{2}}\right]\omega_x\omega_z + Po \cos(t) \frac{2b^{2}}{a^{2}+b^{2}}\omega_z\nonumber\\ &\quad - \frac{2b^{2}}{b^{2}+c^{2}}\omega_x+ Po \cos(t)+\mathcal{L}\varGamma_y, \end{align}

(2.9)\begin{align} \frac{\partial \omega_z}{\partial t} &= \left[\frac{2c^{2}}{b^{2}+c^{2}}-\frac{2c^{2}}{a^{2}+c^{2}}\right]\omega_x\omega_y - Po \cos(t) \frac{2c^{2}}{a^{2}+c^{2}}\omega_y \nonumber\\ &\quad - Po \sin (t)\frac{2c^{2}}{b^{2}+c^{2}}\omega_x +\mathcal{L}\varGamma_z, \end{align}

where $\boldsymbol {\varGamma }$ is the viscous torque and $\mathcal {L}$ is a matrix accounting for the geometry of the cavity

(2.10)\begin{equation} \mathcal{L}=\frac{15}{16\pi} \begin{bmatrix} \dfrac{b^2+c^2}{b^2c^2} & 0 & 0 \\ 0 & \dfrac{a^2+c^2}{a^2c^2} & 0 \\ 0 & 0 & \dfrac{b^2+a^2}{b^2a^2} \end{bmatrix} . \end{equation}

Assuming a laminar Ekman boundary layer, Busse (1968), Noir et al. (2003) and Noir & Cébron (2013) derived expressions for the viscous torque in the asymptotic limit of small Ekman and Poincaré numbers. The viscous torque arises from the differential rotation between the fluid and the container, $\boldsymbol {\delta \omega }$, which we decompose into an axial ($\boldsymbol {\delta \omega }_{ax}$) and equatorial component ($\boldsymbol {\delta \omega }_{eq}$) with respect to the fluid rotation axis. The two components are given by

(2.11)\begin{align} \boldsymbol{\delta \omega}_{ax}=\left(\frac{\boldsymbol{\varOmega}-\hat{\boldsymbol{k}}_s }{\varOmega^2} \boldsymbol{\cdot} \boldsymbol{\varOmega} \right)\boldsymbol{\varOmega}, \end{align}

(2.12)\begin{align} \boldsymbol{\delta \omega}_{eq}=\boldsymbol{\varOmega}-\hat{\boldsymbol{k}}_s - \boldsymbol{\delta \omega}_{ax}, \end{align}

where we have introduced

(2.13)\begin{equation} \boldsymbol{\varOmega}=\boldsymbol{\omega}+\hat{\boldsymbol{k}}_s, \end{equation}

the fluid rotation vector as seen from the frame of precession, i.e. the turntable.

Assuming that the viscous torque acting on the fluid prevents the growth of a spin-over mode from the equatorial differential rotation and the spin-up of the fluid from the axial differential rotation, Noir & Cébron (2013) derived the following expression for the viscous term:

(2.14)\begin{equation} \mathcal{L}\boldsymbol{\varGamma}=\sqrt{E\varOmega}\left( \lambda_r \begin{bmatrix}\varOmega_x\\\varOmega_y\\\varOmega_z-1\end{bmatrix} + \frac{\lambda_i}{\varOmega}\begin{bmatrix}\varOmega_y\\-\varOmega_x\\0\end{bmatrix} +\varepsilon ( \lambda_{sup} -\lambda_{r} ) \begin{bmatrix}\varOmega_x\\\varOmega_y\\\varOmega_z\end{bmatrix} \right), \end{equation}

where $\lambda _r$ and $\lambda _i$ are the real and imaginary parts of the viscous correction to the spin-over mode and $\lambda _{sup}$ represents a spatially integrated decay rate of the axial rotation. Furthermore, we have introduced

(2.15)\begin{equation} \varepsilon=\frac{\varOmega^2-\varOmega_z}{\varOmega^2}, \end{equation}

a measure of the no spin-up condition introduced by Noir et al. (2003). We refer to expression (2.14) as the full damping.

In the limit $\varepsilon <\omega /\varOmega$, when the axial torque is negligible in comparison with the equatorial one, we may neglect the last term in (2.14), to obtain

(2.16)\begin{equation} \mathcal{L}\boldsymbol{\varGamma}=\sqrt{E\varOmega}\left( \lambda_r \begin{bmatrix}\varOmega_x\\\varOmega_y\\\varOmega_z-1\end{bmatrix} + \frac{\lambda_i}{\varOmega}\begin{bmatrix}\varOmega_y\\-\varOmega_x\\0\end{bmatrix} \right), \end{equation}

which is equivalent to the expression of the viscous torque in a spheroid from Busse (1968) and Noir et al. (2003). We refer to (2.16) as the reduced damping.

While the inviscid left-hand side of (2.7)–(2.9) is exact, all the assumptions are contained in the parameterization of the viscous effects, and it is thus important to make a few remarks regarding their range of validity: first, the approximation of the torque by the decay of spin-over and axial spin-up is only valid for a laminar boundary layer, and we shall later see that this assumption is satisfied in our experiments. Additionally, the parameters for the viscous correction of the spin-over mode $\lambda _r$ and $\lambda _i$ can be calculated only when the fluid rotation axis is aligned with one of the principal axes of the ellipsoid (Vantieghem 2014). Finally, to the best of our knowledge, there is no analytical expression of $\lambda _{sup}$ in a non-axisymmetric ellipsoid and one has to rely on experimental estimates for this parameter. In conclusion, we note that, while the model will prove to be robust for small to moderate values of $Po$, we should keep in mind these limitations when discussing experiments with a large differential rotation between the fluid and the container or with a large tilt of the fluid rotation axis. In the remainder of this article, we will use $\lambda _r=-2.55$ and $\lambda _i=0.79$ as calculated from Vantieghem (2014) for the geometry of our ellipsoid and for $\lambda _{sup}$ we use an experimentally determined value of $-2.50$, which is explained in further detail in Appendix A.

2.2. Numerical integration

In contrast to the spheroidal geometry, the governing equations for a non-axisymmetric ellipsoid are not in the form of an implicit equation and we must solve the system by time stepping the system of ordinary differential equations (ODEs) (2.7)–(2.9). To capture potential multiple solutions of the system, we randomly select $300$ initial rotation vectors $\boldsymbol {\omega }$ by uniformly sampling the three variables $(\omega,\theta,\phi )$ in the range $([0,1],[0,2\pi ],[0,2\pi ])$ from which we construct the initial rotation vectors as $\omega _x = \omega \sin (\theta )\cos (\phi ),\omega _y = \omega \sin (\theta )\sin (\phi )$ and $\omega _z = \omega ^2 \cos (\theta ))$. Starting from the respective initial condition, we numerically integrate the (2.7)–(2.9) for 30 viscous spin-up times. Figure 2 represents the time averaged amplitude of the fluid rotation vector $\langle | \boldsymbol {\omega } | \rangle _t$ at the same conditions as our experiment. A more detailed discussion of the time evolution of the uniform vorticity in our numerical model can be found in Appendix B. In a range from $Po_c^1 = 0.075\pm 0.005$ to $Po_c^2 = 0.23\pm 0.005$, the system admits two solutions at each $Po$, as previously reported for axisymmetric spheroids (Cébron 2015; Nobili et al. 2020).

Figure 2. Time averaged amplitude of the fluid rotation vector from the numerical integration of (2.7)–(2.9): two distinct branches of solution are identified and bi-stability (grey) is observed for $0.075\pm 0.005< Po<0.23\pm 0.005$. The colours represent the number of randomly selected initial conditions converging to the respective solution. The values of $a$, $b$ and $c$ are the exact same ones as in the experiment.

3. Experimental set-up and flow measurements

3.1. Experimental set-up

The experimental set-up is sketched in figure 1(a). An ellipsoidal acrylic container of semi-major axes $a = 0.125\ \text {m}$, $b= 0.104\ \text {m}$ and $c = 0.078\ \text {m}$ is rotating around $c$ using a servomotor (type: SGMGV-09DDA6F, produced by YASKAWA) and a belt driving system. The rotation rate $\varOmega _s$ is set to 1.57 rad s$^{-1}$. To create a precessional forcing, the rotating tank and the drive system are mounted on a turntable driven by another belt system connected to a second servomotor (type: SGMGH-30DCA61, produced by YASKAWA). The precession rate $\varOmega _p$ can be varied between 0.018 and 3.674 rad s$^{-1}$ in increments of 0.018 rad s$^{-1}$. The angle between $\boldsymbol {\varOmega }_{\boldsymbol {p}}$ and $\boldsymbol {\varOmega }_{\boldsymbol {s}}$, i.e. the precession angle, is fixed to $90^\circ$. The container is completely filled with water at room temperature with a kinematic viscosity $\nu =1\times 10^{-6}\ \text {m}^2\ \text {s}^{-1}$. The experimental values of the non-dimensional parameters introduced above are $6.3\times 10^{-5}$ for the Ekman number and the Poincaré number is varied between $0.035$ and $0.33$.

3.2. Flow measurements

We measure the flow inside the ellipsoid with ultrasonic Doppler velocimetry (UDV) from which we recover velocity profiles along selected chords in the fluid. To that end, the fluid is seeded with a mixture of 2AP1 Particles produced by Griltex, of sizes 50 $\mathrm {\mu }$m (60 % by weight) and 80 $\mathrm {\mu }$m (40 % by weight). The density of the tracer particles is $1.02$ g cm$^{-3}$. We use a DOP4010 system produced by Signal Processing SA,1073 Savigny, Switzerland, which allows us to simultaneously measure three velocity profiles at a sampling rate up to 12 Hz on each probe. Each profile consists of $\sim$500 individual sampling points along the chord with a spatial resolution of 0.5 mm.

The probe location has been carefully chosen to allow a direct measurement of the uniform vorticity component (2.4)–(2.6). The precise locations of the transducers are presented in table 1. To illustrate the concept of our determination of the fluid rotation $\omega _x,\omega _y$ and $\omega _z$, let us consider the probe measuring $u_x$, i.e. probe 1 in figure 1(b), located at $y=0$ and $z=-0.47$ (in dimensionless units) and measuring in direction of $x$. In the direction of measurement, the uniform vorticity component $U_x$ is constant along the chord and takes the value $U_x=\omega _y (2a^2)/(a^2+c^2)(-0.47)$. Hence, we can use the velocity average along the chord $\langle u_x\rangle _x$, where $\langle \rangle _x$ denotes the spatial average along $x$, as a proxy for $\omega _y$. Placing two additional probes in the same fashion, we can reconstruct all three components of $\boldsymbol {\omega }$ as a function of time. In the rest of the paper, unless explicitly stated, all velocities are in dimensionless units, i.e. re-scaled by $\varOmega _s R$.

Table 1. Position of the UDV probes (centre of the front lens), measured velocity components ($\boldsymbol {u}_{\boldsymbol {i}}$) and associated component of the fluid rotation vector ($\boldsymbol {\omega }_{\boldsymbol {i}}$). All coordinates are non-dimensional.

Figure 3(a) illustrates a typical measurement of the velocity components $u_x, u_y$ and $u_z$ and the corresponding time series of the deduced uniform vorticity at $Po=0.05$ (figure 3b). At leading order, the velocity is constant along the profile, in agreement with the uniform vorticity assumption. The uniform vorticity is dominated by the equatorial components $\omega _x$ and $\omega _y$, which are oscillating at $\varOmega _s$, while the axial component $\omega _z$ is characterized by a small, steady retrograde rotation and a oscillation at $2\varOmega _s$.

Figure 3. Typical velocity measurements in the experiment: (a) velocity profiles at $Po = 0.05$. (b) Computing the spatial average along each profile allows us to reconstruct the three components of the uniform vorticity $\boldsymbol {\omega }$ as a function of time. The error bars display the standard deviation of the spatial average. The Ekman number is $6.30\times 10^{-5}$ and $x$, $y$ and $z$ in (a) are made dimensionless by $R$ and the time in (b) is made dimensionless by $\varOmega _s$.

3.3. Experimental protocol

We apply the following experimental protocol: the container is set into rotation at $\varOmega _s$ and we wait until the fluid co-rotates with the container, which is indicated by a vanishing velocity on all UDV channels. Subsequently, we start the motion of the turntable at $\varOmega _p$ and wait until a statistically steady state is reached. We then start our measurements collecting 3000 samples on each channel at a sampling frequency of 12 Hz. Finally, we change the precession rate to the next value of $\varOmega _p$ and wait again until the fluid motion at the increased precession rate is in a steady state. We perform two sets of experiments: in the first one we subsequently increase $Po$ from $Po=0.035$ to $Po=0.33$, in the second set of experiments we start from $Po=0.33$ and decrease the precession rate down to $Po=0.035$.

4. The uniform vorticity

4.1. Overview of the results

We start with a discussion of the uniform vorticity component and define four important quantities that characterize the rotation of the fluid: first, the time averaged total rotation of the fluid $\langle \varOmega \rangle _t$, second, the time averaged differential rotation between the fluid and the container, $\langle \delta \omega \rangle _t$, third, the time average rotation of the fluid along the spin axis $\langle \varOmega _z\rangle _t$ and, finally, the time averaged tilt angle of the fluid rotation axis with respect to the container axis, $\langle \theta \rangle _t$. These quantities are defined in the frame of precession, i.e. viewed from the turntable, and can be obtained from our measurements of $\omega _x$, $\omega _y$ and $\omega _z$ as follows:

(4.1)$$\begin{gather} \langle \varOmega\rangle _t=\langle \sqrt{\omega_x^2+\omega_y^2+(\omega_z+1)^2}\rangle _t , \end{gather}$$

(4.2)$$\begin{gather}\langle \delta \omega \rangle _t=\langle \sqrt{\omega_x^2+\omega_y^2+\omega_z^2}\rangle _t , \end{gather}$$

(4.3)$$\begin{gather}\langle \varOmega_z\rangle _t=\langle \omega_z+1\rangle _t, \end{gather}$$

(4.4)$$\begin{gather}\langle \theta\rangle _t= \left\langle \arccos \frac{\varOmega_z}{\varOmega}\right\rangle _t, \end{gather}$$

where $\langle \rangle _t$ denotes the average in time. The experimental data for the four quantities as a function of $Po$ are presented in figure 4, together with the prediction from the model (2.7)–(2.9), using the full damping (blue curve) and the reduced damping (red curve). As we have seen above, the model predicts two distinct branches of solutions in the range $Po_c^1=0.075\pm 0.005< P_o< Po_c^2=0.23\pm 0.005$, indicated by the grey area. We call branch 1 the solutions characterized by a large total fluid rotation, a small differential rotation, a large axial rotation and a moderate tilt angle. In contrast, the solutions of branch 2 have a small total fluid rotation, a large differential rotation of ${O}(1)$, a vanishing axial rotation and a large tilt angle.

Figure 4. Fluid rotation viewed from the precession frame: (a) total fluid rotation. (b) Differential rotation $\langle \delta \omega \rangle _t$ between the fluid and the container. (c) Axial component of rotation of the fluid. (d) Angle between the fluid and container axis of rotation. The points represent the experimental data, circles when increasing $Po$, triangles when deceasing $Po$. The solid and dashed lines represent the model branches 1 and 2, respectively, for the full damping (blue) and the reduced damping (red). The shaded area represents the range of $Po$ for which the model accepts two solutions. The two bounds of this region are $Po_c^1=0.075\pm 0.005$ and $Po_c^2=0.23\pm 0.005$. The arrows materialize the experimental transition from branch 1 to branch 2 and vice versa. All experiments are performed at $E=6.3\times 10^{-5}$ and the parameters for the damping in the model are $\lambda _r=-2.55$, $\lambda _i=0.79$ and $\lambda _{sup}=-2.50$. The error bars are representative of the standard deviation of the displayed quantities over the entire time series.

Experimentally, as we start from the smallest values in $Po\approx 0.03$ and subsequently increase $Po$ between the individual experiments, we follow the solutions of branch 1 until a critical value $Po_{exp}^2=0.18\pm 0.006$ is reached. At this point, the uniform vorticity component in the experiment transits to branch 2, where it remains for all larger $Po$ values. In contrast, when decreasing $Po$ after starting from large values $Po\approx 0.33$, the experimental data follow branch 2 down to $Po_{exp}^1=0.11\pm 0.006$, from where they transit to branch 1. Outside of the hysteresis region ranging from $Po_{exp}^1$ to $Po_{exp}^2$, the measured quantities for increasing and decreasing $Po$ overlap almost perfectly. Comparing the model with the experimental data, we observe a very good agreement in branch 1 for all four quantities. However, in branch 2 the model systematically underestimates the differential rotation and the tilt of the fluid rotation axis. Moreover, the model predicts a small positive axial rotation whereas we observe negative values in the experiment. Nevertheless, it is fair to say that the model still captures the fundamental properties of this second branch correctly.

We notice that $Po_{exp}^{1}$ and $Po_{exp}^{2}$ in the experiment do not match exactly the boundaries of the hysteresis domain of the model. For both, increasing and decreasing $Po$, the experimental data transit to the other branch before the end of the hysteresis range in the model. As we shall see later, both transitions occur in a range of $Po$ where fully developed instabilities are observed, with amplitudes comparable to that of the underlying uniform vorticity flow. In such conditions it is likely that finite amplitude perturbations may trigger early transitions, in comparison with the uniform vorticity model, which is free of any perturbations or instabilities by construction. Alternatively, since the fluid rotation is no longer close to the container axis, it may be argued that the asymptotic values of the damping coefficients $\lambda _r,\lambda _i$ and $\lambda _{sup}$ are not representative of the boundary layer dynamics in this range of $Po$. In Appendix C we consider finite variations of the three damping coefficients, demonstrating that no set of values ($\lambda _r,\lambda _i,\lambda _{sup}$) can consistently explain our observations. Using different values of $\lambda _r,\lambda _i$ and $\lambda _{sup}$ for each $Po$ may seem somewhat more physical and might even yield a better agreement between experimental data and model, however, there is a high risk of over-fitting the data, especially in the absence of an underlying physical model. We believe it is more appropriate to use the model with the asymptotic values and keep in mind the limitations of its validity range.

4.2. Spin-over and spin-up damping contribution

Integrating the system of ODEs with the full and the reduced damping yields almost identical solutions, as indicated by the similarity between the blue and red curves in figure 4. Both damping models fit equally well the data at the smallest $Po$ in branch 1 while showing a less good agreement at larger $Po$, which suggest that the viscous torque is dominated by the spin-over component over the entire range of explored $Po$. From the expression of the viscous term (2.14), we calculate the time averaged amplitude of the spin-over contributions $\langle | \mathcal {L}\varGamma _{so} |\rangle _t$ and spin-up contribution $\langle | \mathcal {L}\varGamma _{sup} |\rangle _t$ as

(4.5)$$\begin{gather} \langle|\mathcal{L}\boldsymbol{\varGamma}_{sup} |\rangle _t= \langle |\varepsilon\sqrt{E\varOmega} \lambda_{sup} \boldsymbol{\varOmega} |\rangle _t, \end{gather}$$

(4.6)$$\begin{gather}\langle |\mathcal{L}\boldsymbol{\varGamma}_{so} |\rangle _t =\langle | \mathcal{L}\boldsymbol{\varGamma} - \mathcal{L}\boldsymbol{\varGamma}_{sup}|\rangle _t. \end{gather}$$

The results for branch 1 and branch 2 solutions are depicted in figure 5. Over the entire range of investigated $Po$ values, both the model and the experiment are dominated by the spin-over viscous term, which explains the close similarity between the full and reduced model. In branch 1, the damping from the model is in excellent agreement with the experimental data, with a spin-over contribution one order of magnitude larger than the spin-up. In branch 2, the model still captures the dominance of the spin-over contribution correctly but the actual amplitude of both torques is larger in the experiment, which may be a signature of nonlinear contributions which are not accounted for in the model.

Figure 5. Comparison of the spin-over (red) and spin-up (green) damping term in the experiment and in the model. (a) Branch 1 solutions (solid line). (b) Branch 2 solutions (dashed line).

These findings suggest that the so-called no spin-up condition in a spheroid as derived by Busse (1968) and Noir et al. (2003), i.e. $(\varOmega ^2-\varOmega _z)/\varOmega ^2=0$, may also hold at small $Po$ in an non-axisymmetric ellipsoid. Figure 6 shows $\langle |\varepsilon |\rangle _t= \langle |(\varOmega ^2-\varOmega _z)/\varOmega ^2 |\rangle _t$ as a function of $Po$ for the solutions of branches 1 and 2 from the model and the experiment, confirming that the no spin-up condition is a reasonable assumption at small $Po$ in branch 1. In branch 2, this condition is clearly violated, as the experimental data and the model display values $>1$.

Figure 6. Measure of the no spin-up condition $\langle |\varepsilon |\rangle _t= \langle |(\varOmega ^2-\varOmega _z)/\varOmega ^2 |\rangle _t$ in the experiment. The error bars are representative of the standard deviation of the displayed quantities over the entire time series.

5. Non-uniform vorticity flow

5.1. Kinetic energy in the uniform and non-uniform vorticity flow

By construction, the semi-analytical model only accounts for the uniform vorticity component of the flow, which may not be dominant at large $Po$. Considering the flow in the frame of precession, the velocity can be decomposed into its uniform vorticity component, hereafter referred to as the base flow U_b and the non-uniform vorticity components $\boldsymbol{u}_{non}$:

(5.1)\begin{equation} \boldsymbol{u}=\boldsymbol{U}_b +\boldsymbol{u}_{non}, \end{equation}

with

(5.2)\begin{equation} \boldsymbol{U}_b=\boldsymbol{\varOmega}\times \boldsymbol{r}+\boldsymbol{\nabla}\psi. \end{equation}

We define the associated mean kinetic energy as

(5.3)$$\begin{gather} E_b=\frac{1}{V}\int \frac{\boldsymbol{U}_b^2}{2} \,\textrm{d}v, \end{gather}$$

(5.4)$$\begin{gather}E_{non}=\frac{1}{V}\int \frac{(\boldsymbol{u}-\boldsymbol{U}_b)^2}{2} \,\textrm{d}v. \end{gather}$$

In figure 7(a) we present the time averaged ratio $\langle E_{non}/E_{b}\rangle _t$ when increasing and decreasing $Po$ in branches 1 and 2. The dynamics evolves from a largely dominant uniform vorticity flow in branch 1 to a state of almost equal kinetic energy in the uniform and non-uniform vorticity flow. Comparing the individual energies $\langle E_{b}\rangle _t$ and $\langle E_{non}\rangle _t$ we observe two distinct behaviours: the base flow amplitude hardly evolves over the entire branch 1, but abruptly drops by a factor of 10 across the transition to branch 2 at $Po_{exp}^2$ (figure 7b). In contrast, $\langle E_{non}\rangle _t$ increases continuously by almost two orders of magnitude in branch 1, reaching a saturation value that also holds in branch 2 with a much smaller jump at the branch transition (figure 7c). The increasing amplitude of the fluctuating velocity component $u_{non}$ at larger forcing amplitude is also reflected in the more complex flow regime, as can be anticipated from movies of qualitative flow visualization using rheoscopic fluid (Borrero-Echeverry, Crowley & Riddick 2018), which we provide as supplementary movies are available at https://doi.org/10.1017/jfm.2021.932.

Figure 7. (a) Time averaged ratio of the non-uniform vorticity flow and the base flow mean kinetic energy $E_{non}/E_{b}$ as a function of $Po$; (b) uniform vorticity base flow mean kinetic energy $E_{b}$ as a function of $Po$; (c) non-uniform vorticity flow mean kinetic energy $E_{non}$ as a function of $Po$. The colour scheme characterizes the branch of solution, the two symbols represent the increasing/decreasing $Po$ experiments. The error bars are representative of the standard deviation of the displayed quantities over the entire time series. The unit of energy is $(\varOmega _sR)^2$.

5.2. Instabilities

After characterizing its amplitude, we now aim at shedding light on the dynamical nature of the non-uniform vorticity flow. Three instability mechanisms have been identified in precessing fluid cavities: parametric resonances (Kerswell 1993; Goto et al. 2007; Lin et al. 2015), viscous boundary layer instabilities (Tilgner & Busse 2001; Cébron et al. 2019; Nobili et al. 2020; Buffett 2021) and, more recently, experiments and numerical simulations suggest that centrifugal instabilities may occur in precessing cylinders leading to space filling turbulence (Giesecke et al. 2018, 2019). Our spatially limited UDV measurements may not allow us to completely disentangle these mechanisms, for example we cannot get access to the radial profile of the angular momentum that plays a crucial role in the identification of the centrifugal instability reported by Giesecke et al. (2018). Nevertheless, we can still gain insights into the evolution of instabilities from the spectral content of $\boldsymbol {u}_{non}$. To that end, we calculate the discrete Fourier transform (DFT) in time of $\boldsymbol {u}_{non}$, at all positions along the three velocity profiles. We then take the spatial average along each chord, to obtain the space averaged spectral content of $u_{non,x}, u_{non,y}$ and $u_{non,z}$. Finally, we stack all three DFT components and the resulting spectra are presented in figure 8a) for increasing $Po$ and figure 8b) for decreasing $Po$. At very small $Po\leq 0.11$, we observe only frequencies corresponding to the forcing at dimensionless frequency $\varpi /\varOmega _s=1$ and $\varpi /\varOmega _s=2$. Following branch 1 for $Po\geq 0.12$ up to the transition at $Po=0.18$ (figure 8a) we distinguish individual peaks which satisfy the typical parametric resonant conditions $\varpi _1 \pm \varpi _2=\varpi$, with $\varpi = \varOmega _s$ or $\varpi =2\varOmega _s$ near onset (Kerswell 1993; Lin et al. 2015), where the instability saturates via viscous and/or detuning effects. After the transition to branch 2, the spectra are again characterized by dominating peaks $\varpi /\varOmega _s=1$ and $\varpi /\varOmega _s=2$, but this time significant energy is widely spread over all frequencies. This observation is characteristic for all experiments of branch 2, also for the experiments at decreasing $Po$ depicted in figure 8b)

Figure 8. Stacked DFT in time of $u_{non}$. Each spectrum is normalized by the respective Po. (a) DFT for increasing Po; (b) DFT decreasing Po. Blue and orange colours mark branch 1 and branch 2, respectively.

We observe typical growth and collapses of the instability as illustrated in figure 9 for the smallest $Po$ of branch 2, i.e. at $Po = 0.13$, where we observe a decay of the base flow during the growth of the instability until a state is reached, where base flow and instability have almost the same kinetic energy. At that point, the instability collapses and allows the base flow to recover its initial amplitude. These cycles of growth and collapse are typical of weakly super critical parametric instabilities (Malkus 1989; Eloy, Le Gal & Le Dizes 2003; Herreman et al. 2010) and somewhat reminiscent of the intermittent solution as predicted by the amplitude equations for the weakly nonlinear instability in a precessing cylinder reported by Lagrange et al. (2011). At larger $Po$, the temporal dynamics of the flow becomes more and more complex, with shorter period of growth and collapse and eventually indistinguishable cycles.

Figure 9. Time series of $E_{non}$ and $E_b$ for the lowest $Po=0.129$ of branch 2. The background grey curves correspond to the unfiltered data, the coloured curves represent a running averaging of the energies with a window of 1.5 rotations.

Concerning boundary layer instabilities, it is suggested that instabilities of a steady Ekman boundary layer should occur at $Re_{bl}>55$ and a transition to turbulence at $Re_{bl}>150$ (Caldwell & Van Atta 1970; Sous et al. 2013), where $Re_{bl}$ is the Reynolds number based on the thickness of the Ekman layer. However, the recent numerical study of Buffett (2021) suggests that, for oscillating flows, such as the ones driven by precession, the Ekman boundary layer would only become turbulent at $Re_{bl}>500$. We define the boundary layer Reynolds number based on the amplitude of the differential rotation between the wall and the fluid $\delta \omega$ as $Re_{bl}=\delta \omega E_f^{-1/2}$, where $E_f$ is the effective Ekman number based on the fluid rotation $\varOmega$. With our notation it follows that

(5.5)\begin{equation} Re_{bl}=\langle \delta \omega \rangle _t (E/\langle \varOmega\rangle _t)^{{-}1/2}. \end{equation}

As shown in figure 10, Ekman boundary layer instabilities may occur in branch 2, yet never reach the onset condition for turbulence. Consequently, the instabilities reported above starting at $Po\sim 0.117$, could not be underlined by a boundary layer mechanism.

Figure 10. Boundary layer Reynolds number $Re_{bl}=\delta \omega (E\varOmega )^{-1/2}$ as a function of $Po$ in branches 1 and 2. The dashed black line represents the onset range $Re_{bl}>55$ for the initial boundary layer instability.

To conclude this section, we shall briefly discuss the evolution of the Rossby number based on the fluctuating velocity component $\boldsymbol {u_{non}}$ and the fluid rotation

(5.6)\begin{equation} Ro=\frac{\langle (2E_{non})^{1/2}\rangle _t}{\langle \varOmega\rangle _t R}. \end{equation}

In figure 11 we display $Ro$ as a function of $Po$ for our experimental data. In branch 1, $Ro$ increases from $\sim$0.01 to $\sim$0.06, in agreement with the increasingly nonlinear dynamics observed from the frequency spectra, yet remaining less than unity. In contrast, branch 2 is characterized by a Rossby number greater than one, giving rise to a strongly nonlinear dynamics. Given the reported values of $Ro>0.5$ it is difficult to apprehend this regime, as it may no longer follow the classical rotating nor non-rotating dynamics. Further investigations will be necessary to fully comprehend this regime, in particular the emergence of a regime where $E_{b}$ and $E_{non}$ are of the same order is of great interest for dynamo experiments and industrial applications.

Figure 11. Effective Rossby number based on the fluctuating velocity component $u_{non}$ and the total rotation of the fluid.

6. Conclusions

We have experimentally investigated the flows in a triaxial ellipsoid subject to precession at an angle of $90^\circ$. We first focused on the base flow of uniform vorticity and, by exploiting symmetry properties of this flow, we measure the three components of the uniform vorticity flow using UDV. In agreement with a semi-analytical model for the evolution of the uniform vorticity proposed by Noir & Cébron (2013), we observe two distinct branches of solutions: branch 1 is typically observed at small to moderate forcing ($Po<0.18$) and is characterized by a large amplitude of the total fluid rotation ($\sim$1 in dimensionless units), as well as a small angle between the fluid rotation and the rotation axis of the container. In contrast, branch 2 is characterized by a much smaller total rotation of the fluid and a large tilt of the fluid rotation axis, which tends to align with the axis of precession. Furthermore, we find a hysteresis cycle of the uniform vorticity component as a function of $Po$, from both the model and the experimental data. A remarkable agreement between the uniform vorticity model and experimental data is observed in branch 1, especially at small values of $Po$. This might justify the application of the model in the context of planetary core and subsurface ocean dynamics, which are typically characterized by a small amplitude of the forcing.

Although the model still captures the fundamental features of branch 2, the agreement with experimental data is less good. This is to be expected as the semi-analytical model is only valid at small $Po$, yet the predictions remain in good enough qualitative agreement with the observations to capture the first-order dynamics of the system.

It is of great interest for both planetary science and industrial applications to characterize different instability mechanisms in precessing flows, which could lead to space filling turbulence. Our results suggest that a parametric instability mechanism occurs early in branch 1 and evolves towards a complex saturated state in branch 2 with a kinetic energy comparable to that of the underlying base flow. It should be emphasized that our observation of a hysteresis in the amplitude of the uniform vorticity base flow is not necessarily related to that of turbulence (as reported by e.g. Malkus 1968; Herault et al. 2015; Horimoto et al. 2018; Komoda & Goto 2019), while the opposite is likely true. In fact, one could imagine a scenario where branch 2 of the base flow is not necessarily unstable. Nevertheless the observation of the fluctuating velocities being of the same order as the base flow in branch 2, as well as the qualitative flow visualization (see supplementary movies), point towards the first scenario in our experiment. Of course it would be quite valuable to draw a conclusive picture of the transition to turbulence in various container shapes, (i.e. in cylinder, spheres, spheroids and ellipsoids) but this is beyond the focus of the present article.

The observation of a state with same order of kinetic energy in the base flow and fluctuating velocity component provides a reasonable estimate for the typical $u_{non}$ in industrial applications operating in similar range of parameters, i.e. $E\gtrsim 10^{-6}$ and $Po\sim 1$. In planetary settings, i.e. $E\sim 10^{-14}$, it has been suggested that the saturation mechanism of the instabilities might be fundamentally different, i.e. inertial wave turbulence instead of geostrophic turbulence (Le Reun et al. 2019). Whether a state of quasi-equal kinetic energy holds in planetary cores or subsurface oceans remains an open question. Further investigations at much lower Ekman numbers will be necessary, such as the one accessible in the upcoming Dresdyn experiment (Stefani et al. 2015).