Triadic resonances in nonlinear simulations of a fluid flow in a precessing cylinder

Instabilities and waves in rotating fluids are important in numerous technical applications. In most cases it is important to avoid these phenomena in order to ensure the stability of fast spinning liquid filled bodies like gyroscopes (Pfeiffer 1974), spacecrafts (Bao and Pascal 1997) or projectiles (Miller 1982). Waves in rotating fluids are also of general interest because of their fundamental character in oceanographic and atmospheric flows (Rossby 1939) like, e.g., zonal flows in the atmospheres of Saturn or Jupiter that may develop via a modulation instability of meridional Rossby waves (Connaughton et al 2010). Moreover, the identification of inertial waves in the Earth's fluid outer core (Aldridge and Lumb 1987) allows conclusions on the core dynamics and the associated dynamo process that is responsible for the generation of the Earth's magnetic field. Nowadays, it is believed that in general planetary dynamos are driven by rotating thermal or compositional convection (King et al 2010, D Breuer et al 2010, M Breuer et al 2010). On the other hand, other types of flows cannot be ruled out, and, for example, a superposition of random inertial waves in a rotating conducting fluid is capable of transferring energy in a magnetic field as well (Moffatt 1970). However, in the model of Moffatt (1970), field and flow must decay in the long term because of the lack of energy sources required for a steady driving. This problem can be overcome by more complex models in which, for example, the equatorial heat flux in the Earth's outer core provides a persistent excitation mechanism for inertial waves with sufficient helicity for the generation of planetary magnetic fields (Davidson 2014). Inertial waves can be excited by precessional forcing as well, and it has long been discussed whether precession of the Earth can provide the necessary power to drive the geodynamo (Stewartson and Roberts 1963, Busse 1968, Malkus 1968).

In a next-generation dynamo experiment currently under development at Helmholtz-Zentrum Dresden-Rossendorf (HZDR) a precession driven flow of liquid sodium will be used to validate its suitability to drive a dynamo (Léorat et al 2006, Léorat 2006, Stefani et al 2012, 2015). Since the precessional forcing provides a natural driving mechanism, it represents a qualitatively new approach compared to previous dynamo experiments, where an optimized flow was driven either by impellers (Gailitis et al 2000, Monchaux et al 2007) or by a system of electromagnetic pumps (Stieglitz and Müller 2001).

A small-scale water experiment is currently running at HZDR in order to examine the hydrodynamic properties of a precession driven flow in a cylindrical container with strong forcing and a large precession angle. Measurements of pressure fluctuations and of axial velocity profiles show three distinct flow regimes in dependence of the precession ratio: a laminar state (with only the forced mode) followed by a nonlinear regime with various modes superimposed on the fundamental mode, eventually leading to a permanent chaotic state when a critical precession ratio is exceeded (Herault et al 2015, see also Goto et al 2014 for the spherical case). Presumably, the ability of the flow to drive a dynamo will be different in the three regimes and although the precessional dynamo experiment at HZDR will allow magnetic Reynolds numbers of the order of Rm ≈ 700 (based on the rotation velocity and the radius of the cylindrical container) it is a priori not obvious whether a magnetic field will be self-generated at all. There are promising indications for dynamo action driven by precession both from liquid sodium experiments by Gans (1971) which—despite of the rather small precession ratio—achieved an amplification of an applied field by a factor of three, as well as from simulations in spheres (Tilgner 2005), spheroids (Wu and Roberts 2009), ellipsoids (Ernst-Hullermann et al 2011), cubes (Krauze 2010) and cylinders (Nore et al 2011). However, kinematic dynamo simulations demonstrate that the primary forced mode alone, which consists of a single, non-axisymmetric mode with azimuthal wave number m = 1 and axial wave number k = 1, cannot drive a dynamo at a reasonable Rm (Giesecke et al 2015). This coincides with Herreman and Lesaffre (2011) who showed that a linear wave packet of inertial waves cannot drive a mean-field dynamo. However, this does not generally preclude the possibility of quasi-laminar (in the sense of Dudley and James 1989) or small scale dynamo action driven by inertial waves or another precession induced instability with sufficient intricate topology.

Experimental studies in a weakly precessing cylinder filled with a non-conducting fluid showed a large diversity of flow phenomena with complex structures frequently resulting in a vigorous turbulent state (Gans 1970, Manasseh 1992, 1994, Kobine 1995, Manasseh 1996). Other experiments revealed large scale structures like a system of intermittent, cyclonic vortices that may provide a strong source of helicity and hence be beneficial for a dynamo (Mouhali et al 2012), or free Kelvin modes with an azimuthal wave number m = 5 and m = 6 propagating around the symmetry axis of a weakly precessing cylinder (Lagrange et al 2008, 2011).

Free Kelvin modes are the natural eigenmodes in a rotating cylinder (Thomson 1880, Greenspan 1968, McEwan 1970) and are promising candidates for driving quasi-laminar dynamo action in case of precessional forcing because at least a subclass of these modes has a structure similar to the columnar convection cells that are responsible for dynamo action in convection driven models of the geodynamo (Lorenzani and Tilgner 2003). Free Kelvin modes may emerge from nonlinear interactions of a forced Kelvin mode with itself or as a parametric instability that involves the forced mode and two free Kelvin modes (Kerswell 1999, Lagrange et al 2008, 2011). The latter case requires appropriate combinations of wave numbers, frequencies and a geometry such that all three modes become resonant simultaneously which usually is named a triadic resonance. Triadic resonances have been observed in simulations of precession driven flow in a spheroid (Lorenzani and Tilgner 2003) and, experimentally, in weakly forced precession in a cylinder (Lagrange et al 2008) and a cylindrical annulus (Lin et al 2014). They may be involved in various phenomena in rotating fluids, e.g., the generation of mean flows or the transition to turbulence by providing a possibility for a secondary instability that may finally end up in small-scale turbulent behavior like it is observed in the resonant collapse (Kerswell 1999).

In the present study we conduct numerical simulations of a precessing flow far below the transition to the chaotic state observed in water experiments, but with sufficient forcing to excite free Kelvin modes with azimuthal wave numbers m > 1. The essential aim of our study is to examine the characteristics of the precession driven flow. This concerns mainly the energy that can be injected into the flow by the precessional driving (and thus would be at least partly available for the generation of a magnetic field) and the geometric structure in terms of large scale non-axisymmetric flow patterns. We believe that such flow patterns are beneficial for dynamo action and we expect that the excitation of triadic resonances would be a promising scenario to efficiently trigger such flow components. The computational costs for performing the simulations presented in our study are very high, and a similar work would no longer be feasible when performing full magnetohydrodynamic simulations. Therefore we have restricted ourselves to the hydrodynamic case in order to find guidance in which parameter regime a dynamo is most likely to occur according to the above-mentioned criteria. We focus on the impact of the aspect ratio which—at least in the linear approximation—determines whether an inertial wave becomes resonant. We identify different free Kelvin modes by means of their spatial structure, and we compute their frequencies from the azimuthal drift motion in order to conclude whether the free Kelvin modes constitute a triadic resonance. The results are used to further constrain the modes that may be observed in the water experiment at HZDR³ and will be applied in future simulations of the magnetic induction equation in order to proof whether free Kelvin modes are suitable to drive a dynamo.

The flow of a fluid with kinematic viscosity ν in the frame of a closed precessing cylinder is described by the Navier–Stokes equation including terms for the Coriolis force and the Poincaré force (see e.g. Tilgner 1998):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}+2({{\boldsymbol{\Omega }}}_{{\rm{p}}}+{{\boldsymbol{\Omega }}}_{{\rm{c}}})\times {\boldsymbol{u}}=-{\rm{\nabla }}P+\nu {{\rm{\nabla }}}^{2}{\boldsymbol{u}}-({{\boldsymbol{\Omega }}}_{{\rm{p}}}\times {{\boldsymbol{\Omega }}}_{{\rm{c}}})\times {\boldsymbol{r}}.\end{eqnarray} \tag{ 1 }$$

Here, ${\boldsymbol{u}}$ denotes the velocity field, which additionally obeys the incompressibility condition ${\rm{\nabla }}\cdot {\boldsymbol{u}}=0,\;$ ${\boldsymbol{r}}$ is the position vector, P is the reduced pressure (including the centrifugal terms), ${{\boldsymbol{\Omega }}}_{{\rm{c}}}$ is the angular frequency of the cylindrical container and ${{\boldsymbol{\Omega }}}_{{\rm{p}}}$ denotes the precession (see figure 1). The precession axis is time dependent in the cylinder frame, and the unit vector ${\widehat{{\boldsymbol{\Omega }}}}_{{\rm{p}}}$ that denotes the orientation of the precession axis is given by

$$\begin{eqnarray}&&{\widehat{{\boldsymbol{\Omega }}}}_{{\rm{p}}}=\mathrm{sin}\alpha \mathrm{cos}{{\rm{\Omega }}}_{{\rm{c}}}t\widehat{{\boldsymbol{x}}}-\mathrm{sin}\alpha \mathrm{sin}{{\rm{\Omega }}}_{{\rm{c}}}t\widehat{{\boldsymbol{y}}}+\mathrm{cos}\alpha \widehat{{\boldsymbol{z}}}\end{eqnarray} \tag{ 2 }$$

with α the angle between rotation axis and precession axis. We consider a cylinder with aspect ratio Γ = H/R, where H and R are height and radius so that $r\in [0;R]$ and $z\in [0;{\rm{\Gamma }}R].$ In the following we non-dimensionalize all quantities using R as length scale and ${{\rm{\Omega }}}_{{\rm{c}}}^{-1}$ as time scale.

Zoom In Zoom Out Reset image size

In the inviscid limit a linear solution that fulfills the boundary conditions $\hat{{\boldsymbol{z}}}\cdot {\boldsymbol{u}}=0$ at the top and the bottom is given by a Kelvin mode (Thomson 1880, Greenspan 1968):

$$\begin{eqnarray}&&{{\boldsymbol{U}}}_{m,k,n}(r,z,\varphi ,t)={{\boldsymbol{u}}}_{m,k,n}(r,z){{\rm{e}}}^{{\rm{i}}({{\rm{\Omega }}}_{m,k,n}t+m\varphi )}+{\rm{c.c.}}\end{eqnarray} \tag{ 3 }$$

with a frequency Ω_{m, k, n} that depends on the azimuthal wave number m, the axial wave number k and a third number n that counts the roots of the dispersion relation

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{j}{\lambda }_{j}{J}_{m-1}({\lambda }_{j})+m\left(2-{{\rm{\Omega }}}_{j}\right){J}_{m}({\lambda }_{j})=0{\rm{}}\;{\rm{with}}\;{\rm{}}{{\rm{\Omega }}}_{j}=\pm 2\sqrt{{\left(1+{\left(\displaystyle \frac{{\lambda }_{j}}{{\rm{\Gamma }}k\pi }\right)}^{2}\right)}^{-1}},\end{eqnarray} \tag{ 4 }$$

in which J_m denotes the Bessel function of order m and the triple index (m, k, n) is replaced by j. Positive (negative) frequencies correspond to retrograde (prograde) propagation. Both signs correspond to different radial wave numbers (different solutions of the dispersion relation) and hence to a different spatial structure of the corresponding Kelvin mode.

The dispersion relation ensures the fulfillment of the radial boundary conditions $\hat{{\boldsymbol{r}}}\cdot {\boldsymbol{u}}=0$ at the sidewalls, and although λ_j is not an integer it plays a role similar to a radial wave number with its position in the sequence of zeros corresponding to the number of half-cycles in the radial direction. The spatial structure of a Kelvin mode in a meridional plane ${{\boldsymbol{u}}}_{j}(r,z)$ is given by

$$\begin{eqnarray}{u}_{j}^{r}(r,z) & = & \displaystyle \frac{-{\rm{i}}}{(4-{{\rm{\Omega }}}_{j}^{2})}\left[{{\rm{\Omega }}}_{j}{\lambda }_{j}{J}_{m-1}({\lambda }_{j}r)+\displaystyle \frac{m(2-{{\rm{\Omega }}}_{j})}{r}{J}_{m}({\lambda }_{j}r)\right]\mathrm{cos}\left(\displaystyle \frac{\pi {kz}}{{\rm{\Gamma }}}\right),\\ {u}_{j}^{\varphi }(r,z) & = & \displaystyle \frac{1}{(4-{{\rm{\Omega }}}_{j}^{2})}\left[2{\lambda }_{j}{J}_{m-1}({\lambda }_{j}r)-\displaystyle \frac{m(2-{{\rm{\Omega }}}_{j})}{r}{J}_{m}({\lambda }_{j}r)\right]\mathrm{cos}\left(\displaystyle \frac{\pi {kz}}{{\rm{\Gamma }}}\right),\\ {u}_{j}^{z}(r,z) & = & -{\rm{i}}\displaystyle \frac{k\pi }{{{\rm{\Omega }}}_{j}}{J}_{m}({\lambda }_{j}r)\mathrm{sin}\left(\displaystyle \frac{\pi {kz}}{{\rm{\Gamma }}}\right),\end{eqnarray} \tag{ 5 }$$

with λ_j and Ω_j taken from the solutions of the dispersion relation (4).

In the inviscid linear approximation the total flow of the forced mode with m = 1 is given by

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{m=1}(r,z,\varphi )=\displaystyle \sum _{k}\displaystyle \sum _{n}{A}_{1{kn}}{{\boldsymbol{u}}}_{1{kn}}(r,z){{\rm{e}}}^{{\rm{i}}\varphi +{{\rm{\Omega }}}_{1{kn}}t},\end{eqnarray} \tag{ 6 }$$

where we replaced the index j by the individual contributions and ${A}_{1{kn}}$ represents the amplitude of the forced mode which can be computed by a projection of the applied forcing onto a particular mode (Manasseh 1994). For a precession ratio (or Poincaré number) Po = Ω_p/Ω_c and a precession angle α the corresponding calculation yields (Liao and Zhang 2012)

$$\begin{eqnarray}&&{A}_{j}=\displaystyle \frac{-{\rm{Po}}\mathrm{sin}\alpha {{\rm{\Gamma }}}^{3}(1-{(-1)}^{k}){(2+{{\rm{\Omega }}}_{j})}^{2}(2-{{\rm{\Omega }}}_{j})}{{\lambda }_{j}^{2}[{(k\pi /{\rm{\Gamma }})}^{2}+\frac{1}{2}(2-{{\rm{\Omega }}}_{j})](1-{{\rm{\Omega }}}_{j}){J}_{1}({\lambda }_{j})}.\end{eqnarray} \tag{ 7 }$$

If the wave length that corresponds to an axial wave number of a mode matches exactly the height of the cylinder, this mode becomes resonant with an eigenfrequency Ω_j = 1 and the expression for the inviscid amplitude (7) diverges. For each Γ there exists, in principle, an infinite number of resonant modes. From the dispersion relation (4) we find that the primary forced mode, i.e. the mode with m = 1, the axial wave number k = 1 and a radial wave number corresponding to the first root, becomes resonant at Γ = 1.9898174 which is rather close to the so called spherical cylinder with height equals diameter (Γ = 2). Further resonances for increasing radial wave number occur at Γ = 0.9560735 (for n = 2) and at Γ = 0.6206981 (for n = 3). The computation of the amplitude at resonance requires a consideration of viscous effects in the bulk and in terms of boundary layers with associated Ekman layer suction. Then the amplitude can be calculated by matching the Ekman layer suction to the precessional force (Gans 1970). A more general approach for the computation of the amplitudes of the response created by the precessional forcing is given by Liao and Zhang (2012) who derived an asymptotic solution including viscous effects that is valid at and away from resonance.

The precessional forcing only excites modes with m = 1 and k odd. Higher azimuthal modes or modes with even axial wave number must be triggered by nonlinear interactions, e.g., in terms of triadic resonances that involve this forced Kelvin mode and two free Kelvin modes (m_a, k_a) and (m_b, k_b) or by the nonlinear interaction of the forced Kelvin mode with itself (Meunier et al 2008). In the following we assume that the resonant case is most promising in view of the dynamo problem due to the expected larger amplitudes of the velocity field.

The free Kelvin modes are solutions of the linearized Navier–Stokes equation with a Coriolis term as a restoring force but without a precessional driving on the right-hand side:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+2({{\boldsymbol{\Omega }}}_{{\rm{p}}}+{{\boldsymbol{\Omega }}}_{{\rm{c}}})\times {\boldsymbol{u}}=-{\rm{\nabla }}P.\end{eqnarray} \tag{ 8 }$$

The principle of the interaction of free and forced modes becomes evident when we consider an exact triad with a forced mode ${{\boldsymbol{U}}}_{{\rm{f}}},$ with m_f = 1 and k_f = 1 and two free Kelvin modes ${{\boldsymbol{U}}}_{a}$ with m_a and k_a and ${{\boldsymbol{U}}}_{b}$ with m_b and k_b. We ignore viscosity and further assume a constant amplitude A_f of the forced mode so that we can write the total flow as

$$\begin{eqnarray}&&{\boldsymbol{U}}={A}_{{\rm{f}}}{{\boldsymbol{U}}}_{{\rm{f}}}+a(t){{\boldsymbol{U}}}_{a}+b(t){{\boldsymbol{U}}}_{b},\end{eqnarray} \tag{ 9 }$$

where a(t) and b(t) denote the amplitude modulations of the free mode (without the harmonic oscillations Ω_a and Ω_b) which are given by Kerswell (1999)

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}a}{{\rm{d}}t} & = & \langle {{\boldsymbol{U}}}_{a},{{\boldsymbol{U}}}_{b}\times ({\rm{\nabla }}\times {{\boldsymbol{U}}}_{{\rm{f}}}^{*})+{{\boldsymbol{U}}}_{{\rm{f}}}^{*}\times ({\rm{\nabla }}\times {{\boldsymbol{U}}}_{b})\rangle b(t),\\ \displaystyle \frac{{\rm{d}}b}{{\rm{d}}t} & = & \langle {{\boldsymbol{U}}}_{b},{{\boldsymbol{U}}}_{a}\times ({\rm{\nabla }}\times {{\boldsymbol{U}}}_{{\rm{f}}}+{{\boldsymbol{U}}}_{{\rm{f}}}\times ({\rm{\nabla }}\times {{\boldsymbol{U}}}_{a})\rangle a(t).\end{eqnarray} \tag{ 10 }$$

The square brackets $\left\langle\cdot ,\cdot \right\rangle$ denote a scalar product defined as $\left\langle{{\boldsymbol{U}}}_{a},{{\boldsymbol{U}}}_{b}\right\rangle=\int {{\boldsymbol{U}}}_{a}^{*}\cdot \;{{\boldsymbol{U}}}_{b}{\rm{d}}V,$ and a normalization of the free modes is assumed such that $\left\langle{{\boldsymbol{U}}}_{i},{{\boldsymbol{U}}}_{j}\right\rangle={\delta }_{{ij}}.$ In order to achieve non-trivial solutions of (10), and to constitute a triadic resonance, the wave numbers and frequencies of the free modes a and b must fulfill the conditions

$$\begin{eqnarray}\delta m & = & | {m}_{b}-{m}_{a}| & = & {m}_{{\rm{f}}} & = & 1,\\ \delta k & = & | {k}_{b}-{k}_{a}| & = & {k}_{{\rm{f}}} & = & 1,\\ \delta {\rm{\Omega }} & = & | {{\rm{\Omega }}}_{b}-{{\rm{\Omega }}}_{a}| & = & {{\rm{\Omega }}}_{{\rm{c}}} & = & 1,\end{eqnarray} \tag{ 11 }$$

where m_f = 1, k_f = 1 and Ω_c = 1 denote the properties of the forced mode. In contrast to sine and cosine-functions, there are no corresponding addition theorems for Bessel functions so that no further restrictions are imposed on the interaction of different radial wave numbers.

In the following, we solve the dispersion relation (4) for k = 1, 2 and $m=1,...,7.$ We assume that higher radial modes are damped and only consider the first root for each combination of m and k. Figure 2 shows the frequency difference δ Ω of two free Kelvin modes with δ m = 1 and δ k = 1 versus the aspect ratio in the vicinity of the critical value δΩ = 1. The crossings of the curves with the horizontal line at δΩ = 1 mark aspect ratios for which an exact triadic resonance may occur. These aspect ratios at which δΩ = 1 are listed in table 1 together with the corresponding frequencies of the resonant free Kelvin modes.

Zoom In Zoom Out Reset image size

A triad is always formed by free Kelvin modes consisting of one prograde and one retrograde mode (Waleffe 1993). Note that the frequencies of free Kelvin modes that constitute different triads often are close (if we restrict to triads with k_a = 1 and k_b = 2, see table 1) making it difficult to reliably identify individual Kelvin modes in realistic setups with potential frequency shifts due to viscous and/or nonlinear effects.

We conduct numerical simulations in a cylindrical geometry with the code SEMTEX that applies a spectral element Fourier approach for the numerical solution of the Navier–Stokes equation and has been thoroughly validated against analytical solutions (Blackburn and Sherwin 2004). Coriolis forcing was added later, and was shown to converge exponentially against a constructed, 3d, time-dependent exact solution in both cylindrical and Cartesian coordinates. Further, numerical results are in excellent agreement with experimental data (Albrecht et al 2015b, figure 5) and theory (Albrecht et al 2015a). In order to simplify the implementation of the precessional forcing we switch to the precessional frame in which the precession is stationary and the cylinder rotates with a frequency Ω_c. Thus only a term for the Coriolis force appears in the Navier–Stokes equation which now reads

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}+2{{\boldsymbol{\Omega }}}_{{\rm{p}}}\times {\boldsymbol{u}}=-{\rm{\nabla }}P+\nu {{\rm{\nabla }}}^{2}{\boldsymbol{u}}.\end{eqnarray} \tag{ 12 }$$

Due to the rotation of the cylinder in the precessing frame the boundary conditions for the azimuthal flow change to u_φ = rΩ_c (at endcaps and sidewall) whereas no-slip conditions are still applied for u_z and u_r.

The problem is described by four parameters, the Reynolds number defined with the angular velocity of the container ${\rm{Re}}={{\rm{\Omega }}}_{{\rm{c}}}{R}^{2}/\nu ,$ the precession ratio (or Poincaré number) Po = Ω_p/Ω_c, the aspect ratio Γ = H/R and the precession angle α with $\mathrm{cos}\alpha ={{\boldsymbol{\Omega }}}_{{\rm{p}}}\cdot {{\boldsymbol{\Omega }}}_{{\rm{c}}}/| {{\boldsymbol{\Omega }}}_{{\rm{p}}}| | {{\boldsymbol{\Omega }}}_{{\rm{c}}}| .$ In the present study, we keep the precession ratio fixed at Po = 0.0141 and the precession angle is set to α = π/4. All simulations are performed at Re = 6500 except for one run at Re = 10 000 to briefly examine the impact of increasing Re. The aspect ratio is varied in the range ${\rm{\Gamma }}\in [1.7;2.4]$ in order to find the container geometry for which individual modes or triads become resonant.

4.1. Pattern of the total flow

The simulations are started from an initial state with pure solid-body rotation ${\boldsymbol{u}}=r{\widehat{{\boldsymbol{e}}}}_{\varphi }.$ After switching on the Coriolis force at t = 0 the fluid shows a direct response in form of a large scale mode with m = 1. In dependence of the aspect ratio this is only a transient state until after a finite time (which may take up to 1000 rotation periods or more) a quasi-steady state is reached. The typical flow pattern for Γ = 2 and Re = 6500 is presented in figure 3 where the nested isosurfaces show a snapshot of the axial velocity u_z at 30% (60%, 90%) of its maximum value. The flow is clearly dominated by a velocity mode with m = 1 and k = 1. This is the primary Kelvin mode driven by the precessional forcing. Superimposed on the m = 1 mode we see some non-axisymmetric (time-dependent) disturbances that will be analyzed below.

Zoom In Zoom Out Reset image size

4.2. Kinetic energy

4.2.1. Forced mode

In the following we discuss the behavior of the flow by means of azimuthal Fourier modes ${\tilde{{\boldsymbol{u}}}}_{m}(r,z)$ with the total flow given by ${\boldsymbol{u}}(r,z,\varphi )=2\pi /({N}_{\varphi }+1){\displaystyle \sum }_{m=0}^{{N}_{\varphi }}{\tilde{{\boldsymbol{u}}}}_{m}(r,z){{\rm{e}}}^{{\rm{i}}m\varphi }.$ Figure 4 shows the history of the kinetic energy of the m = 1 mode for different aspect ratios. A final state arises either after a short transitional period (Γ = 2.000) or after one or more bifurcations (see e.g. the blue curve for Γ = 2.150). The final state can be chaotic (Γ = 2.000, black curve), oscillatory (Γ = 2.15, blue), quasi-periodic with collapses (Γ = 2.2, red), quasi-periodic with bursts (Γ = 2.225, green) or stationary (Γ = 1.8, yellow). We also find transient periods with constant energy (e.g. Γ = 2.15 for t = 300...600), and it may take up to 1200 revolutions of the cylinder (Γ = 2.15, blue curve) till a final quasi-stationary state is reached. This saturation time also depends on the Reynolds number and decreases with increasing Re.

Zoom In Zoom Out Reset image size

The bifurcations shown by the m = 1 mode go along with the emergence of non-axisymmetric modes with m ≥ 2 (figure 5). These modes become unstable only after quite a long time, e.g., for Γ = 1.825 it takes up to 600 rotation periods of the cylinder till modes with m ≥ 2 start to grow, and it may take up to twice that time till a quasi-steady state is reached. The non-axisymmetric modes with m ≥ 2 essentially proceed similar to the behavior of the m = 1 mode, i.e., we see a chaotic behavior when the m = 1 mode behaves chaotic (Γ = 2.0, figure 5, central panel) and we see periodic growth and decay when the m = 1 mode is periodic (Γ = 1.825, top panel in figure 5 and Γ = 2.200, bottom panel in figure 5). In the periodic regime the amplitude of the energy oscillations of the modes with m ≥ 2 is rather large so that the kinetic energy may vary in time by more than one order of magnitude.

Zoom In Zoom Out Reset image size

In the following we compare analytic expressions for the energy of the forced mode with the time-averaged energy taken from the quasi-stationary period of the final state in our simulations. This comparison is only partly justified because the analytical expressions are based on a number of serious simplifications, such as the linearization of the Navier–Stokes equation, the limit ${\rm{Po}}\to 0,$ and the neglect of the time dependent part of the Coriolis force which is justified only for small precession angles α. Furthermore, the final state in our simulations is instable and additional caution is advised when comparing simulations with a linear time independent solution. Nevertheless, such an analysis is helpful to identify the regimes in which, e.g., a maximum response of the flow can be expected, or to localize the aspect ratios at which triadic resonances are possible.

Analytically, the energy of the forced Kelvin mode is given by (Liao and Zhang 2012)

$$\begin{eqnarray}&&{E}_{m=1}={{\rm{\Gamma }}}^{2}\displaystyle \sum _{j}\;\displaystyle \frac{{(k\pi /{\rm{\Gamma }})}^{2}+\frac{1}{2}(2-{{\rm{\Omega }}}_{j})}{{{\rm{\Omega }}}_{j}^{2}(4-{{\rm{\Omega }}}_{j}^{2})}{| {A}_{j}| }^{2}{J}_{1}^{2}({\lambda }_{j}),\end{eqnarray} \tag{ 13 }$$

where the index j now denotes a combination of the axial wave number k and the radial wave number n, and A_j is the inviscid linear amplitude given by (7) with the azimuthal wave number fixed at m = 1. A_j can also be calculated by the method of Liao and Zhang (2012) which includes viscosity and is valid at and off the resonance. However, in that case equation (13) is accurate only up to an error of $O\left({{\rm{Re}}}^{-\frac{1}{2}}\right)$ since equation (13) does not consider the energy from the flow in the boundary layers. Figure 6 shows the time-averaged energy taken from the simulations (black curve) compared with the analytic solutions (green curve, inviscid solution and blue curve, viscous solution). For the calculation of the analytical solutions equation (13) is truncated at n = 6 and k = 6 which is sufficient to reach convergence (off-resonance in the inviscid case) and avoids the occurrence of too many resonances in the inviscid case (green curve) which anyway are suppressed in the more realistic viscous or nonlinear computations. Comparing the numerical solutions (black curve) with the linear solutions we see significant deviations around the primary resonance (k = 1, n = 1, Γ_res ≈ 1.9898) while outside of this regime the agreement of all curves remains good (say for Γ ≲ 1.85 and for Γ ≳ 2.05). Two key features characterize the behavior of the mean energy in our simulations. On the one hand we see a clear shift of the maximum to smaller aspect ratios and the maximum energy arises rather close to the aspect ratio expected for the resonance with k = 3 and n = 3. However, this resonance is suppressed by viscous effects and we do not see strong contributions with k = 3 and/or n = 3 in our simulations (see also section 4.5.1). Therefore the correlation is either coincidental or we have to assume an indirect impact of the (k = 3, n = 3) resonance. The second scenario is supported by the abrupt transition from the absolute maximum to a rather low energy state right below Γ ≈ 1.871. This jump is also connected to a modification of the character of the m = 1 mode which changes its behavior below Γ ≈ 1.86 from the chaotic regime to the periodic regime. Abrupt transitions of the energy of higher azimuthal modes can also be seen at the other (theoretical) resonances of the forced mode (k = 5, n = 5 and, a little further away, k = 5, n = 4, see figure 6) which will be revisited in the next paragraph. Beside the shift of the maximum to smaller aspect ratios the energy obtained in the simulations remains considerably smaller compared to the linear viscous solution. The maximum of the energy based on the amplitude of Liao and Zhang (2012) is ${E}_{{\rm{kin}}}^{{\rm{max}}}\approx 0.021$ (at Γ ≈ 1.9814) which is roughly two times larger than the maximum kinetic energy obtained in our simulations (${E}_{{\rm{kin}}}^{{\rm{max}}}\approx 0.0093$ at Γ = 1.871). The differences are not surprising because the prerequisites for the computation of the amplitudes A_j are not well met in our simulations. In particular, we observe a significant impact of nonlinear interactions in the simulations in terms of higher non-axisymmetric modes and the forming of an azimuthal shear flow that both draw energy from the forced m = 1 mode.

Zoom In Zoom Out Reset image size

4.2.2. Higher azimuthal wave numbers

We only find higher non-axisymmetric modes with a noteworthy amount of energy between Γ ≈ 1.81 and Γ ≈ 2.24. For the sake of brevity we will call this regime the nonlinear regime, assuming that the higher modes are triggered by nonlinear effects. The kinetic energy of the modes with m ≥ 2 is roughly 2 orders of magnitude smaller than the energy of the forced mode but essential characteristics such as the maximum at Γ = 1.871 and the abrupt drop of kinetic energy below this threshold are also reflected in the behavior of the higher non-axisymmetric modes (figure 7). The sudden drop of the energy for Γ ≤ 1.871 pertains to all modes which indicates that this breakdown might be a global flow property that is not connected to an individual mode.

Zoom In Zoom Out Reset image size

A similar steep decrease of energy can also be found on the two outer edges of the nonlinear regime. It is striking that all the abrupt transitions are close to theoretical resonances of the forced mode ($k=3,n=3\Rightarrow {{\rm{\Gamma }}}_{{\rm{res}}}\approx 1.8621,k=5,n=5\Rightarrow {{\rm{\Gamma }}}_{{\rm{res}}}\approx 1.8142,$ and with some greater distance $k=5,n=4\Rightarrow {{\rm{\Gamma }}}_{{\rm{res}}}\approx 2.2911,$ see figure 6). These resonances are, however, suppressed and cannot be directly observed in the simulations. So far, we cannot conclude whether the correlation of the energy jumps and the occurrence of (theoretical) resonances is a coincidence, or whether there is some systematic relationship.

Away from the primary resonance of the forced mode we find two further regimes with local maxima of the higher modes. For 1.80 ≲ Γ ≲ 1.85 we see a narrow window with the maximum around Γ ≈ 1.825 in which the m = 6 and m = 7 mode dominate (blue and yellow curve in figure 7). Likewise, but more explicit, we find the m = 4 and m = 5 mode becoming dominant for 2.05 ≲ Γ ≲ 2.25 with a maximum around Γ ≈ 2.2 (green and red curves in figure 7).

The local maxima around Γ ≈ 1.825 and Γ ≈ 2.2 for Fourier modes with δm = 1 indicate the presence of corresponding triadic resonances and roughly agree with predictions from the dispersion relation (see table 1 and figure 2). From the linear predictions we would also expect triads with m = 5 and m = 6 in the intermediate regime (for 1.85 < Γ < 2.05) but obviously there are other flow contributions with at least comparable energy which may disguise the energetic signature of the involved free Kelvin modes (e.g., we see relatively high energies of the m = 4 and the m = 5 modes, whereas in comparison the m = 6 mode is clearly suppressed).

4.3. Spatial structure of the Fourier modes

A qualitative impression of the flow structure is provided by isosurfaces of the individual Fourier modes of the axial velocity shown in figure 8 for Γ = 2 (for an impression of the temporal behavior see the snapshots of the m = 5 mode in figure 9 and the movie at stacks.iop.org/NJP/17/113044/mmedia with the visualization done in the precession frame of reference). The m = 1 mode is obviously dominated by a k = 1 contribution and shows comparatively small variations in time. Regarding the higher non-axisymmetric modes we see a complex pattern for the m = 2 and m = 3 mode with irregular temporal and spatial fluctuations and a more regular behavior of the modes m = 4, 5 and 6 with typical signatures of k = 1 and k = 2. It is striking that the m = 6 mode, which should become resonant around Γ ≈ 2, remains weak compared to the m = 4 or the m = 5 mode.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The amplitude of the individual contributions with m = 4, 5, 6 varies in time and the resulting flow perpetually changes its axial structure. The azimuthal mode resulting from the superposition of k = 1 and k = 2 contributions has a distinct and time dependent asymmetry with respect to the equatorial plane of the cylinder, so that the associated flow is concentrated alternately in both halves of the cylindrical container (see time series of the m = 5 mode in figure 9). The symmetry breaking has been shown to be essential for the functioning of a dynamo in simulations of Tilgner (2005) in a precessing sphere whereas it seemed less important in the direct numerical simulations of Nore et al (2011). The typical timescale for this process is of the order of the rotation period 2π/Ω_c but we can not derive a reliable periodicity from our simulations.

Away from the primary resonance of the forced mode the structure of higher azimuthal modes is much more regular and—except for periodic variations of the amplitude—essentially remains constant over time. Figure 10 shows the dominant azimuthal modes (beyond m = 1) at Γ = 1.825 (m = 6 and m = 7) and at Γ = 2.200 (m = 4 and m = 5). In both cases, the axial behavior of the modes shows clear indications of axial wave numbers k = 1 and/or k = 2 and we will show below that these modes indeed fulfill all conditions of equation (11) and thus constitute triadic resonances. In contrast to the behavior at Γ = 2, the geometric structure of the dominant modes exhibits only minor temporal variations except strong oscillations of the amplitude (see animation at stacks.iop.org/NJP/17/113044/mmedia with the visualization done in the precession frame of reference).

Zoom In Zoom Out Reset image size

4.4. Radial dependence

In order to quantify the visual observations made in the previous section and to estimate amplitude and frequency of individual Kelvin modes we use a discrete sine-transformation applied to each Fourier mode

$$\begin{eqnarray}&&{\widehat{{\boldsymbol{u}}}}_{m}^{k}(r)=\displaystyle \sum _{q=1}^{{N}_{z}}\;{\tilde{{\boldsymbol{u}}}}_{m}(r,{z}_{n})\mathrm{sin}\left(\displaystyle \frac{{qk}\pi }{{N}_{z}+1}\right)\end{eqnarray} \tag{ 14 }$$

with N_z the number of points in axial direction. Figure 11 shows radial profiles for the dominant contributions to the axial velocity at Γ = 1.825 (a) and at Γ = 2.200 (b) where we only present the three leading modes for each case. The solid thick curves present the time average and the bright colors in the background represent the variation of the radial profiles in time. For comparison, the thick dashed curves show the analytical behavior $\propto {J}_{m}({\lambda }_{j}r)$ predicted by the linear in-viscid approximation (5) with λ_j the first root of the dispersion relation (4). The differences between the numerical data and the analytical curves are surprisingly small with deviations mainly in the boundary layers which can be explained by the different boundary conditions for the simulations (no slip) and the in-viscid solution (free slip). Except for (considerable) changes in the amplitude the pattern of the radial profiles exhibits only minor temporal changes and shows no significant contributions from higher radial modes (see color coded contour plots in figure 11 that show the variations of the radial profiles in time).

Zoom In Zoom Out Reset image size

The behavior becomes more complex at Γ = 2 at which we see many modes with different m and k with nearly the same amplitude and clear changes in the radial structure. Figure 12 shows characteristic radial profiles for the mode (m = 5, k = 1) and their temporal behavior as a typical example. The amplitude exhibits considerable fluctuations and, unlike in the previous cases, we see additional changes in the radial structure. When the amplitude is around a (local) maximum, the radial velocity profile is in accordance with a wavenumber that corresponds to the first root of the dispersion relation and the radial profiles can nicely be described by a function $\propto {J}_{5}({\lambda }_{1}r)$ (see red curves in figure 12(a)). However, when the amplitude is weak we see clear signatures of higher radial wave numbers (see figures 12(b) and (d)). Even though the higher radial modes remain weak, they might represent a channel for the transfer of energy from the (periodically occurring) maxima of a mode into smaller scales.

Zoom In Zoom Out Reset image size

4.5. Amplitudes and frequencies

4.5.1. Height equal diameter (Γ = 2)

In the following, the frequencies of individual non-axisymmetric modes are specified in the mantle system. Figure 13 shows the temporal behavior of $| {\widehat{u}}_{m,z}^{k}|$ at r = 0.67 for m = 1 to m = 7 (from top to bottom) and for k = 1 to k = 4 (blue, red, green, yellow curve). The peak values of the amplitudes are additionally listed in table 2 (Γ = 2) and in table 3 (Γ = 1.825 and Γ = 2.2) together with the frequencies obtained from the time derivative of the azimuthal phase of each mode. Note that we only list modes with a regular behavior of the frequencies, i.e., contributions which exhibit a steady and unique azimuthal drift that allows a conclusive computation of the time derivative of the azimuthal phase. The results quantitatively complement the observations made in the previous paragraph. For Γ = 2 (central column and table 2) we see a dominant forced mode with m = 1 and k = 1 and the contributions with k > 1 are negligible (smaller by a factor of 20). The higher azimuthal modes ($m=2...7$) reach approximately 10% of the amplitude of the m = 1 mode with k = 1 and k = 2 slightly prevailing over k = 3 and k = 4. A striking property of the flow at Γ = 2 is the common orientation of the azimuthal drift motion of all modes (see table 2). We only find modes with a retrograde azimuthal drift so no combination is possible that fulfills the requirements for a triadic resonance (in all cases δΩ $\ne$ 1, see table 2). The presence of distinct frequency signals that fit to higher radial modes (which however are hardly evident in the radial profiles for the (m = 5, k = 1) mode, see figure 12) shows that there must be further contributions with chaotic behavior which are probably dominant and having no regular azimuthal drift.

Zoom In Zoom Out Reset image size

Table 2.  Amplitudes and frequencies from the simulations at Γ = 2. Only modes with well-defined drift behavior are listed. The columns denote (from left to right): Reynolds number, azimuthal wave number, axial wave number, maximum of the amplitude, azimuthal drift frequency from the numerical data, and theoretical frequency from the dispersion relation (4) for larger radial wave numbers n that fit best to the observed frequencies.

Re	m	k	$\mathrm{max}| {u}_{m}^{k}|$	Ωnum	Ωtheo
6500	1	1	0.2673	−1.0000	−1.0000 (forced)
	3	1	0.0412	+0.2911	+0.3342 (for n = 2)
	3	2	0.0266	+0.3154	+0.3333 (for n = 5)
	4	1	0.0536	+0.0357	no corresponding mode
	4	2	0.0525	+0.3593	+0.3647 (for n = 4)
	5	1	0.0491	+0.2278	+0.2030 (for n = 3)
	5	2	0.0578	+0.4692	+0.5157 (for n = 2)
	6	1	0.0392	+0.2006	+0.1868 (for n = 3)
	6	2	0.0391	+0.6239	+0.6309 (for n = 1)
10000	1	1	0.2775	−1.0000	−1.0000 (forced)
	3	1	0.0329	+0.3090	+0.3342 (for n = 2)
	3	2	0.0302	+0.3013	+0.3333 (for n = 5)
	4	1	0.0562	+0.1687	+0.1532 (for n = 5)
	4	2	0.0420	+0.3511	+0.3647 (for n = 4)
	5	1	0.0570	+0.2001	+0.2030 (for n = 3)
	5	2	0.0574	+0.5168	+0.5157 (for n = 2)
	6	1	0.0417	+0.3452	+0.3204 (for n = 1)
	6	2	0.0437	+0.6410	+0.6309 (for n = 1)

Table 3.  Amplitudes and frequencies from the simulations for the dominant modes at Γ = 1.825 and Γ = 2.200. The columns denote (from left to right): aspect ratio, azimuthal wave number, axial wave number, maximum of the amplitude, azimuthal drift frequency from numerical data, and theoretical frequency from the dispersion relation (4) at the aspect ratio at which the linear inviscid approximation predicts a triadic resonance (always supposing a radial wave number corresponding to the first root of the dispersion relation which is confirmed in figure 11).

Γ	m	k	$\mathrm{max}| {u}_{m}^{k}|$	${{\rm{\Omega }}}_{{\rm{num}}}$	Ωtheo
1.825	1	1	0.1113	−1.0000	−1.0000 (forced mode)
	6	1	0.0272	−0.2979	−0.3045 (at Γ = 1.80074)
	7	2	0.0382	+0.7010	+0.6956 (at Γ = 1.80074)
	5	1	0.0094	−0.2596	−0.3511 (at Γ = 1.94082)
	6	2	0.0144	+0.7405	+0.6489 (at Γ = 1.94082)
2.200	1	1	0.1095	−1.0000	−1.0000 (forced mode)
	4	2	0.0596	−0.6368	−0.6650 (at Γ = 2.17832)
	5	1	0.0423	+0.3632	+0.3350 (at Γ = 2.17832)

We also performed simulations at a slightly larger Reynolds number Re = 10 000 (at Γ = 2) and found only minor changes with respect to the run at Re = 6500 except for the modes (m = 4, k = 1) and (m = 6, k = 1) which have lower frequencies at Re = 6500. However, so far our data is not sufficient to explain the impact of Re on the drift frequencies or to establish a conclusive scaling towards more realistic Re that will be reached in the experiment. Hence, we refrain from any further discussion of the properties of the flow at larger Re.

4.5.2. Triadic resonances at Γ = 1.825 and Γ = 2.200

The behavior of the amplitude confirms the modified characteristics of the flow when the forced mode is off-resonance (Γ = 1.825, left column in figure 13, or Γ = 2.2, right column in figure 13; see also table 3). For Γ = 2.200 we find a well-defined single triadic resonance with unique wave numbers (m = 4, k = 2 and m = 5, k = 1) and frequencies (Ω_a = −0.6368 and Ω_b = 0.3632). These values fulfill the conditions for a triadic resonance and are quite close to the theoretical values obtained from the dispersion relation (4) (see table 3). The emergence of a triadic resonance is less explicit at Γ = 1.825. However, at this aspect ratio we find even two triadic resonances (see figure 13) with (m, k) = (6, 1) and (7, 2) and a second resonance with (m, k) = (5, 1) and (6, 2) (with much weaker amplitude, see table 3). In all cases the amplitude of the free modes that constitute a triadic resonance strongly oscillates with a maximum of up to 50% of the forced mode (Γ = 2.2, m = 4, k = 2) which corresponds to roughly 6% of the angular velocity of the cylindrical container.

4.6. Azimuthal shear flow

The growth of the triads goes along with a growth of an axisymmetric azimuthal shear flow that is mostly geostrophic (see left panel in figure 14). The axisymmetric mode emerges with a slight delay with respect to the free Kelvin modes which indicates a saturation process by a detuning of the resonance frequencies (see right panel in figure 14). The induced axisymmetric flow component is negative on average and, hence, causes a breaking of the solid body rotation which increases with increasing precession ratio. In the extreme case the breaking of the induced axisymmetric flow entirely cancels the solid body rotation, giving the impression that in the laboratory frame the cylindrical container is rotating around a standing fluid. Kerswell (1999) showed that a geostrophic instability even if it is subdominant to the triadic instability should finally outweigh the triad if Re is large enough because it does not suffer from detuning effects. This may be important for the observation of large mean ciculations in experiments but cannot explain alone the abrupt transition to a chaotic flow which is determined by the precession ratio and hardly depends on Re (Kong et al 2014, Herault et al 2015). This does not mean that the geostrophic mode does not matter for this phanomenon, but in order to be able to judge this more simulations are required at larger Γ.

Zoom In Zoom Out Reset image size

We performed numerical simulations of a precession driven flow in cylindrical geometry. The simulations confirm that the energy that can be injected via precessional forcing is very sensitive to the aspect ratio of the cylindrical container. Significant contributions of higher non-axisymmetric components appear only in a limited range of aspect ratios (for ${\rm{\Gamma }}\in [1.81;2.24]$) and it turns out that in that regime the behavior of the forced mode (with m = 1) cannot be described by a linear theory. Nevertheless, in all cases the forced Kelvin mode (m = 1) dominates and free Kelvin modes with higher azimuthal wave number appear as distortions of the forced mode. The maximum response of the flow takes place at Γ ≈ 1.871 which is remarkable far away from the resonance predicted by the linear in-viscid approximation (at Γ = 1.98982). At the resonance maximum our results yield an amplitude for the z-component of the flow of $| {u}_{z}^{{\rm{max}}}| \approx 0.27,$ which—assuming isotropy—roughly agrees with the estimations of Léorat (2006). The amplitude of the forced mode decreases with increasing distance from the resonance maximum, and in the regimes with triadic resonances, we observe values of $| {u}_{z}^{{\rm{max}}}| \approx 0.10$ (in terms of the azimuthal velocity of the container). In real units (assuming ${{\rm{\Omega }}}_{{\rm{c}}}=2\pi \cdot 10\;{\rm{Hz}}$ and R = 1 m) this would correspond to fluid velocities of 15 m s⁻¹ (at resonance) and 6 m s⁻¹ (off resonance) which is of the same order as the typical flow speed in the Riga dynamo experiment where a helical flow is produced by an optimized system of a propeller and pre- and post-propeller vanes. Regarding the higher Kelvin modes, the maxima of the saturated amplitudes of the strongest free Kelvin modes are mostly independent of the forced mode and achieve values of about 0.05–0.06 in terms of the angular velocity of the container corresponding to ∼3.5 m s⁻¹ in the experiment.

Although at present we do not know whether the flow fields that emerge in our simulations will be able to drive a dynamo, we may try to estimate typical magnetic field strengths that can be expected in the planned dynamo experiment. We assume a saturation of the magnetic energy at roughly 5% of the kinetic energy of the hydrodynamic flow which is substantiated by the saturation behavior of the Riga dynamo (Gailitis et al 2008). We further presuppose that the internal velocity caused by the precessional forcing linearly scales with the rotation of the container (and hence with Re) and assume an angular velocity of Ω_c = 2π · 10 Hz, a radius R = 1 m, a height H = 2 m and a density of liquid sodium of $\rho \approx 930\;{\text{kg m}}^{-3}$ (at ∼400 K).

Using this values the kinetic energy obtained in our simulations at Γ = 2 (E_kin ≈ 0.006) can be scaled to real units which approximately yields a magnetic field strength of $B\approx \sqrt{0.05{\mu }_{0}\rho {E}_{{\rm{kin}}}{{\rm{\Omega }}}_{{\rm{c}}}^{2}}\approx 40\;{\rm{mT}}.$ This is in the range of the Riga dynamo which is not particularly surprising since the typical velocities obtained in our simulations indeed match the velocity arising in the Riga dynamo. However, given that the corresponding m = 1 mode has produced no dynamo in the kinematic simulations (Giesecke et al 2015) it might be more realistic to refer to the free Kelvin modes with higher m which have less kinetic energy. In that case the typical magnetic field strength amounts to only B ≈ 5 mT which, however, is still in the range of the values obtained at the Von-Kármán-sodium dynamo (Monchaux et al 2007). We expect a stronger response and hence a larger saturation field strength for increasing Po, at least as long as we remain below the critical Poincare number for the sudden transition to a chaotic state (Po^crit ≈ 0.0725 at Re = 5.65 × 10⁵, Herault et al 2015). However, until we know whether or not the precession driven flow will have the right structure and sufficient magnitude to excite dynamo action at all, these estimations are only useful to demonstrate that the reference values assumed for the planned dynamo experiment will allow magnetic fields with a reasonable field strength.

Regarding the structure of the flow, we find various manifestations of non-axisymmetric contributions in terms of free Kelvin modes with higher azimuthal wave numbers m. We find free Kelvin modes in resonance with the forced mode only when the forced mode is sufficiently far away from its primary resonance, i.e., when the amplitude of the forced mode is not too strong. Triadic resonances can be clearly identified using a combined Fourier-discrete sine transformation in the azimuth and along the axis with structure and frequencies close to predictions from linear theory. The modes, which we can identify, fulfil the requirements of a triadic resonance consisting of a forced mode with m = 1, k = 1 and ω = 1 and two free Kelvin modes that obey Δ m = 1, Δk = 1 and Δ ω = 1 and dominate the non-axisymmetric flow contributions with m ≥ 2. The amplitude of the higher azimuthal modes (with m ≥ 2) may become quite large in terms of the container velocity. However, regarding their kinetic energy, the results of our study are somewhat disappointing because the higher m modes are mostly concentrated within a rather narrow radial region so that the vast majority of the kinetic energy remains in the forced mode (with m = 1) which is supposed not to drive a dynamo. The modes with m ≥ 2 provide, at best, only a few percent of the total energy, and it seems to be difficult to branch off more energy from the fundamental mode via precession. This behavior may change when increasing the precession ratio, but in that case the computational costs increase drastically making a detailed investigation as performed within our present study at least very long lasting if not almost impossible.

The triadic instability requires a rather long time to emerge and saturates to a final state with slow periodic growth and decay. This period is likely related to the precession time scale but a conclusive correlation requires further simulations at different Po. The periodic growth followed by a fast decay reminds of the so called resonant collapse (McEwan 1970, Manasseh 1992). However, in our simulations the energy is essentially distributed among very few, large scale modes (the fundamental forced mode, two free Kelvin modes that establish the triad and an axisymmetric, essentially geostrophic mode). We do not spot any turbulence-like behavior during the decay of the free Kelvin modes, i.e. there is no transition from the dominance of the large scale mode into small scale turbulent flow which would be characteristic for the collapse phenomenon. This behavior might change when the Reynolds number approaches larger (i.e. more realistic) values. This is substantiated by the simulations at Γ = 2. Here we see no triadic resonance, but rather a quasi-periodic occurrence of free Kelvin modes that exhibit a trend towards larger radial wave numbers during their decay which may represent (part of) a cascade to small scale structures.

Actually, we were also expecting triads around Γ ≈ 2 with m = 5 and m = 6. Indeed, we do find patterns of free Kelvin modes in that regime, but no combination of them satisfies the triadic resonance conditions. Instead, we see free Kelvin modes solely with retrograde drift motion which on the one hand agrees with experimental observations of inertial waves in a spherical Couette experiment with liquid sodium (Kelley et al 2007). On the other hand, they seem to contradict with the results from Albrecht et al (2015a) who show that at Re = 2482 (close to the threshold of the instability) the growth rate is very small resulting in a very long time required for the triadic resonance to emerge (which explains why Kong et al (2014) did not find the triadic resonances at the same parameters). However, the reason for the disappearance of the triads at larger Reynolds number is probably an instability of the m = 1 mode, which occurs close to the fundamental resonance due to the increased efficiency of the forcing, while around Γ = 1.825 or Γ = 2.2 the m = 1 mode may saturate by viscosity before becoming unstable.

The authors acknowledge support from the Helmholtz-Allianz LIMTECH. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915. AG acknowledges the experiences he could make in participating in the program on Wave-Flow Interaction in Geophysics, Climate, Astrophysics, and Plasmas at the Kavli Institute for Theoretical Physics, UC Santa Barbara. The authors kindly acknowledge the discussions with Andreas Tilgner and Oliver Göpfert from the Institute of Geophysics at the University of Göttingen.