2. Simulation framework

2.1. Governing equations

The flow is governed by the non-dimensional continuity and incompressible Navier–Stokes equation

(2.1) \begin{align} \nabla \cdot \textbf {u} &=0, \end{align}

(2.2) \begin{align} \frac {\partial \textbf {u}}{\partial t} +\textbf {u}\cdot \nabla \textbf {u} &= -\nabla p + \frac {1}{{\textit {Re}}} \nabla ^2 \textbf {u}, \end{align}

where $\textbf {u}$ is the velocity vector and $p$ is the pressure. In the present study, the inner sphere is rotating at constant angular velocity ( $\omega _i$ ) while the outer sphere is held stationary. The natural choice of length scale and velocity scale are $R_i$ and $R_i\omega _i$ , respectively, where $R_i$ is the radius of the inner sphere. The parameter $\nu$ is the kinematic viscosity of the fluid. The Reynolds number ( $Re$ ) is defined as $Re=R_i^2\omega _i/\nu$ . The governing equations in spherical coordinates, where the velocity vector $\textbf {u}=(u_\phi$ , $u_\theta$ , $u_r)$ , are solved using Dedalus, a Python-based framework for solving partial differential equations using spectral methods.

The boundary conditions are:

(i) at the inner sphere, $r=R_i$

(2.3) \begin{align} u_{\theta }=u_r=0 \quad and \quad u_{\phi }=1 ;\end{align}

(ii) at the outer sphere, $r=R_o/R_i$

(2.4) \begin{align} u_{\theta }=u_r=u_{\phi }=0 .\end{align}

2.3. Grid dependence study and validation

A grid independence study is performed for the simulations of a spherical annulus in spherical coordinates with the inner sphere rotating with different Reynolds numbers. The solutions for $ Re = 2000$ and $7000$ are obtained for three different grids, G1, G2 and G3, with varying divisions in the radial, polar and azimuthal directions. The maximum and minimum values of the azimuthal and radial components of the velocity field in the domain are compared for these three grids. Table 1 shows the grid dependency for both Reynolds numbers. The results are obtained for a variable time step between $10^{-2}$ and $10^{-4}$ . The simulations were run with an adaptive time-stepping algorithm to ensure numerical stability and accuracy. The exact time-step size varied depending on the Reynolds number and flow condition, with a smaller time-step size used for higher Reynolds numbers or flow states with higher vorticity. For example at $Re =6300$ the time step is $5.58\times 10^{-3}$ . The differences in the maximum and minimum values for the G2 and G3 grids are within 5% for both the Reynolds numbers. Therefore, G2 is chosen as the optimal grid for all the simulations presented in this study.

Table 1.

Grid dependence study for $\beta$ = 0.24 for $Re = 2000$ and $7000$ . Maximum and minimum contour values in the domain have been compared for three different grids for the azimuthal and radial components of the velocity.

Figure 2.

Visualised flow pattern of five spiral waves from (a) Wulf et al. (Reference Wulf, Egbers and Rath1999) and (b) the azimuthal component of the vorticity calculated using the present numerical model for a gap ratio of $\beta =0.5$ and $Re=1320$ .

The results from the present computations are validated qualitatively as well as quantitatively with the existing literature for a range of Reynolds numbers, gap ratios and rotation rates. Figure 2 shows the azimuthal component of the vorticity from the present study for a gap ratio of $0.5$ and $Re=1320$ . The results are compared with the flow visualisation of experimental results from Wulf et al. (Reference Wulf, Egbers and Rath1999), who observed a periodic supercritical flow of wavenumber $k=5$ . While Wulf et al. (Reference Wulf, Egbers and Rath1999) provides experimental evidence, Abbas et al. (Reference Abbas, Yuan and Shah2018) has numerically observed $k=5$ spiral waves for $Re=1280$ using an artificial compressibility method. The present study reproduces the same flow pattern and vorticity distribution as Wulf et al. (Reference Wulf, Egbers and Rath1999), confirming the validity of the proposed numerical method.

Figure 3.

Comparison of the present azimuthal and radial components of velocity with the results of Hollerbach et al. (Reference Hollerbach, Junk and Egbers2006) at (a) $Re=11270$ and $\beta =0.1$ and (b) $Re=489$ and $\beta =1$ .

Figure 3(a) shows the azimuthal and radial components of velocity at the equatorial plane along the radial direction for a gap ratio of $0.1$ and $Re=11270$ . The present study results are in perfect agreement with the results of Hollerbach et al. (Reference Hollerbach, Junk and Egbers2006), who used a finite difference method. The azimuthal velocity is higher than the radial velocity, indicating a strong shear layer near the inner sphere. Figure 3(b) shows the same components of velocity for a gap ratio of 1 and $Re=489$ . Again, the present study results match well with the results of Hollerbach et al. (Reference Hollerbach, Junk and Egbers2006). The azimuthal velocity is higher than the radial velocity, but the difference is smaller than in figure 3(a), indicating a weaker shear layer.

Figure 4.

Comparison of flow structures in the $r-\theta$ plane for $\beta =0.14$ and $Re=1000$ ; (a) azimuthal component of vorticity contour and (b) streamlines. For both plots, the left half is from Abbas et al. (2021) and the right half from the present simulations.

Further, figure 4 compares the streamlines and azimuthal component of the vorticity line contour for a steady two-vortex flow observed for a gap ratio of $0.14$ and $Re=1400$ . The present study results are consistent with the previous study by Abbas et al. (2021), who reported a similar flow structure and vorticity pattern using different numerical schemes. These comparisons demonstrate that the proposed numerical method can accurately simulate the spherical Couette flow with inner sphere rotation for different flow regimes and parameters. The results also agree well with the existing theoretical and experimental studies on the topic, validating the physical relevance and significance of the present study.

3. Details of simulations

Figure 5.

Variation of radial velocity component, $\langle u_r\rangle _{\phi }$ , for various $Re$ . Axisymmetric represents the branch on which the flow remains axisymmetric up to the largest Reynolds number studied in the present study. TWI represents the branch on which the flow is dominated by the travelling wave instability. EI represents the branch dominated by equatorial instability. EI_Backward is the backward extension of the EI branch by reducing the Reynolds number of the flow on the EI branch, showing the hysteresis. The topology of the flow at different branches for various $Re$ is shown using contours of $\theta$ component of vorticity.

The present work analyses the flow for a range of Reynolds numbers $(400{-}6300)$ for a fixed gap ratio $\beta =0.24$ . Figure 5 shows a radial component of velocity, $u_r$ , averaged in the $\phi$ direction, plotted for a range of Reynolds numbers. This figure shows the bifurcations in spherical Couette flow for $\beta =0.24$ with $Re$ as a parameter. We have performed mainly three sets of simulations to obtain various branches of the bifurcation curve of the flow shown in figure 5. In the first set, all the simulations are initialised with the velocity field of the previous lower $Re$ simulation. Results of these simulations form branch I (upper branch) of figure 5 (black circles). Hereafter, this branch is referred to as the axisymmetric branch.

For the second set of simulations, all the simulations start with zero initial field and zero initial perturbation. The results of these simulations form branch II (middle branch) of figure 5 (blue crosses). This branch is referred to as the travelling wave instability (TWI) branch in the following discussions. For the third set, each simulation is started with the following perturbation (Yuan Reference Yuan2012):

(3.1) \begin{align} \hat {v}_r&=-4 \varepsilon _1 \frac {\left (r-R_i\right )\left (r-R_o\right )}{\left (R_o-R_i\right )^2} \cos \left [\pi \left \{1-\frac {R_m\left (\frac {\pi }{2}-\theta \right )}{\left (R_i-R_o\right )}-\varepsilon _2 \sin \left (m_a \phi \right )\right \}\right ], \end{align}

(3.2) \begin{align} \hat {v}_\theta &=\varepsilon _1 \sin \left [\frac {\pi }{2} \frac {\left (r-R_i\right )\left (r-R_m\right )\left (r-R_o\right )}{\left (R_o-R_i\right )^3}\right ] \nonumber \\ &\quad\times \sin \left [\pi \left \{1-\frac {R_m\left (\frac {\pi }{2}-\theta \right )}{\left (R_i-R_o\right )}-\varepsilon _2 \sin \left (m_a \phi \right )\right \} \frac {1}{\beta }\right ]. \end{align}

Here, $(\hat {\cdot })$ denotes the perturbation, $R_m = 0.5(R_i+R_o)$ is the mean radius, $\varepsilon _1$ is the amplitude of the perturbation velocities and $\varepsilon _2$ is the amplitude of the azimuthal perturbation. The values of the perturbations are chosen as $\varepsilon _1=10^{-4}$ , $\varepsilon _2 = 0.4$ and $m_a=6$ in accordance with the literature (Yuan Reference Yuan2012; Abbas et al. Reference Abbas, Yuan and Shah2018). This perturbation satisfies the no-slip boundary conditions at the inner and the outer spheres. A similar form of perturbation is used in Taylor–Couette flow and mimics the Taylor vortices (Schröder & Keller Reference Schröder and Keller1990). This branch is referred to as the equatorial instability (EI) branch.

All the three branches of the curve separate from each other approximately at $Re=4500$ . Once the first point of a new branch is obtained, the subsequent points on that particular branch are obtained using the axisymmetric branch’s strategy. In addition to these three sets of simulations, we have carried out simulations to extend the EI branch backward. The Reynolds number is reduced to check the backward extent (hysteresis) of the EI branch, until the flow jumps back to the axisymmetric branch. In the following sections, we describe the dynamics and topology of the flow on each branch as the Reynolds number is varied.

5. Comment on equatorial instability and polar instability

Figure 20.

Time series data of azimuthal vorticity contours at $r=(R_i+R_o)/2$ ; $\phi = 0$ is plotted along $\theta$ from 0 to $\pi$ for Reynolds numbers (a) 4500, (b) 4700, (c) 5340 and (d) 6300 on the TWI branch.

Both the TWI branch and the EI branch show equatorial instability with wavenumber $k=7$ . However, the nature of this instability is different in both branches. Figure 20 shows the contours of the azimuthal vorticity along the $\theta$ direction, at the mid-radius and $\phi =0$ with time, for various Reynolds numbers on the TWI branch. As mentioned in § 4.2, a rotating equatorial instability appears at $Re=4500$ along with a spiral instability of wavenumber $k=7$ that extends to the polar region of the sphere, as shown in figure 20(a). The rotating spiral instability starts to change its propagation direction at $Re=4700$ . As mentioned earlier, this transition is completed at $Re=5340$ . During this transition, the equatorial instability also exhibits multiple wavenumbers. These wavenumbers are presented in table 2 for various Reynolds numbers. The table shows only the dominant wavenumbers in the equatorial instability and the bold numbers show the most dominant wavenumbers. As the reversal of the propagation direction of the travelling spiral wave instability is completed at $Re=5340$ (cf. § 4.2), we observe that the equatorial instability exhibits multiple wavenumbers with $k=8$ being the most dominant, as can be also seen in figure 20(b) for $Re=5340$ . Table 2 shows that, at larger Reynolds numbers, the equatorial instability supports a range of wavenumbers with multiple dominant modes. This multi-mode equatorial instability can be seen in figure 20(c) for $Re=6300$ , which is the largest Reynolds number analysed on the TWI branch in the present work. We can also observe a multi-mode rotating spiral instability at larger $Re$ , as shown in figures 20(b) and 20(c), which results in the chaotic state of the flow, discussed in the following section.

Table 2.

Dominant wavenumbers present in the equatorial instability for various Reynolds numbers on the TWI branch. The bold number represents the most dominant mode/s, and the dot-dot indicates many other modes.

Figure 21 shows the contours of the azimuthal vorticity along the $\theta$ direction, at the mid-radius and $\phi =0$ with time, for various Reynolds numbers on the EI branch. As mentioned previously (refer to § 4.3), the EI branch exhibits a two-vortex structure at the equator at $Re = 4500$ , resulting from the radial jet hitting the outer sphere at the equator. This results in twin jets above and below the equator (cf. § 4.3), as shown in figure 17. At $Re = 4500$ , the flow is steady and axisymmetric, as shown in figure 21(a). At $Re = 4750$ , the equatorial jet and the twin jets become unstable due to the oscillatory instability in the radial jet. This instability is clearly evident in figure 21(b). At the Reynolds number is increased to $Re = 5750$ , we also observe a weak rotating azimuthal instability appearing in the vicinity of the equator, which can be seen in figure 21(c). This structure is formed due to twin jets above and below the equator (cf. § 4.3), as shown in figure 17. This instability is different from the one observed on the TWI branch, as the flow exhibits a two-vortex structure for all the cases on the EI branch. Above $Re=4750$ , the two-vortex region is analogous to an unsteady double-gyre flow, as can be seen in figures 17(f), 17(g) and 17(h). No such structure is observed on the TWI branch. At larger Reynolds numbers ( $Re=5750$ and $6250$ ), we also observe a weak spiral instability in the flow away from the equator, as shown in figures 21(c) and 21(d). Further studies are required to check whether this spiral instability grows and results in a chaotic flow, something that is not observed for the range of Reynolds numbers investigated in the present work. In the next section, we briefly discuss the phase space dynamics on the three branches of the flow.

Figure 21.

Time series data of azimuthal vorticity contours at $r=(R_i+R_o)/2$ ; $\phi = 0$ is plotted along $\theta$ from 0 to $\pi$ for Reynolds numbers (a) $4500$ , (b) $4750$ , (c) $5750$ and (d) $6250$ on the EI branch.

6. Phase space dynamics

Figure 22.

Phase space trajectories showing a limit-cycle behaviour for the periodic flow at $Re=6250$ on the axisymmetric branch, where time delay $\tau = 350$ .

We have used the numerical probe data of various cases to reconstruct the phase space of the flow as the $Re$ value of the flow is varied. The phase space for each case is reconstructed using Cao’s time delay method (Cao Reference Cao1997). More details of the method can be found in the literature (Kennel et al. Reference Kennel, Brown and Abarbanel1992; Cao Reference Cao1997). In the present work, we have plotted three-dimensional phase space diagrams from the velocity magnitude timeseries. As discussed in the previous section, a weakly periodic, axisymmetric solution emerges when the Reynolds number is increased to $Re=6250$ , on the axisymmetric branch. Figure 22 shows the phase space diagram corresponding to $Re=6250$ on the axisymmetric branch. One axis of the phase space diagram is the velocity magnitude $u(t)$ while the other two are $u(t+\tau )$ and $u(t+2\tau )$ , where $\tau$ is the delay time, $\tau =350$ . We see that the phase space shows a limit-cycle behaviour corresponding to the axisymmetric periodic solution. This is the largest Reynolds number studied on the axisymmetric branch, and we have not seen any other bifurcation on this branch.

On the TWI branch, we have observed an azimuthal spiral instability of $k=7$ at $Re=4500$ . The flow is periodic and non-axisymmetric due to this mode. The timeseries data from a numerical probe has shown the presence of a fundamental frequency and its first harmonic (cf. figure 11), indicating it is a period-doubling bifurcation. The phase space corresponding to the $Re=4500$ exhibits closed orbits, shown in figure 23(a). During the direction-reversing bifurcation that starts at $Re\approx 4700$ , the orbits in the phase space diagram shown in figure 23(b) are spread out as other harmonics of the fundamental mode also appear in the flow (cf. figure 12). After the completion of the direction-reversing bifurcation, the flow exhibits multiple spiral instability modes as it reaches a chaotic state. The diagram for $Re=6300$ , shown in figure 23(c) exhibits the orbits filling the entire phase space indicating the chaotic nature of the flow.

Figure 23.

Phase space trajectories for (a) $Re=4500$ , (b) $Re=4700$ and (c) $Re=6300$ on the TWI branch.

Figure 24.

Phase space trajectories for (a) $Re=5000$ , (b) $Re=5750$ and (c) $Re=6250$ on the EI branch.

The EI branch exhibits a steady and axisymmetric flow at $Re=4500$ . The first bifurcation, which is a symmetry-breaking bifurcation, occurs at $Re=4750$ , where the flow exhibits an equatorial instability of wavenumber $k=7$ . The timeseries of the velocity magnitude shows a fundamental frequency along with its first harmonic (cf. Figure 18), indicating a period-doubling bifurcation. Figure 24(a) shows the phase space diagram for $Re=5000$ , exhibiting period-2 orbits corresponding to this period-doubling bifurcation. Figures 24(b) and 24(c) show the phase space plots respectively for the Reynolds numbers $5750$ and $6250$ . The phase spaces still exhibit a period-2 nature. As discussed in the previous section (cf. § 4.3), the equatorial instability of $k=7$ persists up to the largest $Re=6250$ studied in the present work. The flow at these large $Re$ exhibits weak spiral instability in the polar region but does not show a chaotic nature, as evident from the phase space diagrams of the flow at large $Re$ .

7. Summary and conclusion

This work examines the spherical Couette flow for a narrow-gap ratio of $\beta =0.24$ for a range of Reynolds numbers. We find that the flow is sensitive to initial conditions, and various combinations of the initial conditions are used to obtain different branches of the bifurcation diagram. Such strategies have been used previously in the literature for medium-gap ratios of Spherical Couette flows (Wimmer Reference Wimmer1976). We identify three branches in the bifurcation diagram, plotted in the $\langle u_r \rangle _{\phi }-{\textit {Re}}$ plane, namely: the axisymmetric branch, the TWI branch and the EI branch. Backward marching in Re along the EI branch shows hysteresis in the flow.

As the name suggests, the flow on the axisymmetric branch remains axisymmetric and steady up to a large range of Reynolds numbers $ ({\textit {Re}}\lt 6250 )$ . The flow undergoes a Hopf bifurcation at $Re=6250$ , where an axisymmetric periodic solution emerges. The flow at this branch shows a dominant equatorial jet for all the Reynolds numbers analysed on this branch in the present work. At $Re=4500$ , a branch, called TWI, separates from this axisymmetric branch, where the flow supports travelling spiral waves and an equatorial wave.

The flow on the TWI branch is dominated by an azimuthal spiral instability of wavenumber $k=7$ , along with an equatorial instability of the same wavenumber. As the $Re$ value of the flow is increased, the direction of propagation of the travelling wave reverses. This flow transition starts at $Re\approx 4700$ and completes at $Re\approx 5340$ . During this transition, the equatorial wave first exbhibits a dominant wavenumber of $k=8$ and later, at larger Reynolds numbers, supports multiple wavenumbers. Further analysis is required to understand the relation between the instability of the travelling wave and the polar instability on the EI branch. At large Reynolds numbers, the flow also exhibits a chaotic behaviour, which results from the appearance of the multiple wavenumbers in both spiral and equatorial instabilities in the flow.

The EI branch, which gets separated from the axisymmetric branch at $Re\approx 4500$ , exhibits a two-vortex structure at the equator. This two-vortex structure is caused by two jet streams above and below the equator. As the Reynolds number is increased, the flow exhibits period doubling, the twin jet streams become unstable and the flow becomes non-axisymmetric. Similar to the TWI branch, we observe a weak azimuthal wave of wavenumber $k=7$ at the equator. Unlike the TWI branch, we do not observe a chaotic state for the range of Reynolds numbers investigated in the present work. The flow also exhibits hysteresis on the EI branch. Reducing $Re$ on this branch does not result in the flow jumping back immediately on the axisymmetric branch for $Re\lt 4500$ . Instead, the EI branch gets extended backward to a very low $Re=400$ before the flow jumps back to the axisymmetric branch. A similar analysis needs to be done for other gap ratios to generalise the inferences made here. Larger gaps will be more susceptible to instability, and we expect the transition to turbulence to happen early. It will be interesting to see the equatorial jet instability for these larger annular gaps. Since the initial conditions significantly influence the flow topology, a similar investigation is required with alternative initial conditions to broaden the understanding of bifurcation behaviour and the potential emergence of additional branches.