4.1 Computation of base flows near the bifurcations

Within a time-stepper approach, the 3D steady base-flow solution can be obtained by simply integrating the time-dependent Navier–Stokes equations (2.1) over a long time interval when the flow is in subcritical conditions, i.e. before the occurrence of a bifurcation. However, if the Reynolds number exceeds the instability threshold, we cannot use a simple time integration to compute the unstable steady base flow that is required for the global stability analysis. It is, therefore, necessary to use a stabilizing numerical procedure able to retrieve the unstable reference state. In the case of low-dimensional systems, on one hand, there exist several continuation and bifurcation packages like AUTO (Doedel Reference Doedel1986) or CONTENT (Kuznetsov & Levitin Reference Kuznetsov and Levitin1996) based on direct solvers for the linear algebraic systems and eigenvalue problems involved in the computations. On the other hand, since the discretization of a full three-dimensional problem often leads to a very large discrete system, the extension to high-dimensional problems is not straightforward.

It should be noted that computation of the steady solution using a Newton algorithm, which is the most commonly used method in global stability studies of 2D flows, is not suitable here, as it necessitates a matrix inversion which would require too much memory and time in a 3D case. An alternative, however, would be to adopt a Krylov-subspace method like the classical GMRES (Saad Reference Saad2003). Other approaches can be found in the literature. In Åkervik et al. (Reference Åkervik, Brandt, Henningson, Hoepffner, Marxen and Schlatter2006), for instance, the authors managed to reach the steady state by damping the unstable frequencies using a dissipative relaxation term proportional to the high-frequency content of the velocity fluctuation field. Shroff & Keller (Reference Shroff and Keller1993), instead, were able to compute unstable states by using a projection onto the small unstable subspace coupled with a Newton procedure.

As discussed in the previous section, the base flow is obtained in two different ways. For the Cartesian multigrid method, details have been discussed above. For the Nek5000 code, we adopt a novel efficient algorithm (Citro et al. Reference Citro, Luchini, Giannetti and Auteri2016), inspired by the Krylov-subspace methods, to compute unstable steady states of the inherent dynamical system. This method, like GMRES, is based on the minimization of the residual norm at each integration step. It allows us to obtain the unstable steady field by correcting the new iteration of the numerical procedure using a linear combination of previous fields. The key steps of this procedure can be found in appendix A of Citro et al. (Reference Citro, Giannetti, Luchini and Auteri2015). However, we tested our algorithm in computing the flow past a circular cylinder. In this case, the maximum difference between the base flow computed using our algorithm and the classical Newton–Raphson method is less than  $10^{-10}$ .

4.2 Convergence analysis and validation

The computational domain is chosen to be a large rectangular box surrounding the bluff body (see figure 1). In our case, it extends from $x_{inlet}$ to $x_{outlet}$ , from $y_{down}$ to $y_{up}$ and is symmetric in the $z$ direction with a width of $z_{lat}$ .

Table 1.

The different meshes used to validate our numerical simulations. Nek5000: $n_{elem.}$ is the total number of spectral elements inside the computational domain; $N_{pol.ord.}$ is the polynomial order adopted for the spatial discretization. $IBM$ code: $N_{x}$ , $N_{y}$ , $N_{z}$ are respectively the total number of points in the $x$ , $y$ and $z$ directions.

Figure 2.

Truncation error ( $L_{2}$ norm) as a function of the number of points for the finite-difference second-order immersed-boundary multigrid code. The solution on a grid of $1024\times 512\times 512$ is assumed as the reference solution. In particular, we report here the number of points in the $x$ direction, $N_{x}$ ; the other two directions are characterized by $N_{y}=N_{z}=0.5N_{x}$ . The dashed line represents the theoretical second-order scaling of the error.

Table 1 shows the details of the several meshes used for the analysis (for both the IBM and the SEM code). The adopted grids present a clustering of points near the sphere surface and are symmetrically located with respect to the plane ( $x$ – $y$ ). For the Cartesian IBM code, Roberts stretching transformations (Tannehill, Anderson & Pletcher Reference Tannehill, Anderson and Pletcher1997, pp. 336–337) are used to smoothly cluster the points close to the sphere centre. A similar approach has been successfully used in Giannetti & Luchini (Reference Giannetti and Luchini2007).

In the case of Nek5000, the domain is discretized into hexahedral elements, which are deformed by isoparametric mappings. Within each element, the unknown is represented in local Cartesian coordinates by tensor-product Lagrange polynomials of degree $N_{pol.ord.}$ (and $N_{pol.ord.}-2$ for the pressure field). The upper part of table 1 shows the characteristics of these meshes. The grids adopted in the immersed-boundary calculations are also reported in the same table. In order to validate our numerical methods, we performed several tests to study the effects of the resolution and of the size of the computational domain. A grid convergence analysis was performed for the new immersed-boundary code. Figure 2 depicts the evolution of the error as a function of the number of points adopted in the discretization (in the $x$ direction). We found perfect agreement between the second-order theoretical scaling of the error and the evolution of the error computed by our multigrid code. Furthermore, a validation analysis has been carried out in order to select the best grid in terms of speed and accuracy. Table 2 shows the comparison between the results obtained by the present numerical approaches and the dedicated literature. We found small changes in the lift coefficient $C_{L}$ and in the critical Reynolds number $Re^{SS}$ for the first bifurcation. Thus, for efficiency reasons, we chose the meshes $m_{2}$ and $M^{2}$ to compute all of the results documented in the present work.

Table 2.

Validation of the present numerical approaches. In particular, we compare the lift coefficient $C_{L}$ and the critical Reynolds number for the first bifurcation $\mathit{Re}^{SS}$ with values reported in the literature. (FD $=$ finite-difference method, SEM $=$ spectral element method, FV $=$ finite-volume method, FEM $=$ finite-element method).

Figure 3.

Flow past a sphere immersed in a uniform stream: contour plot of the pressure fields for several flow conditions. The white lines represent the streamlines of the flow. The fixed sphere case ( $\unicode[STIX]{x1D6FA}=0.0$ ) is depicted using two cases: (a) axisymmetric state at $\mathit{Re}=200$ and (b) asymmetric flow at $\mathit{Re}=275$ (onset of mode I without rotation). Furthermore, we select four cases to show the spatial structure of the flow for mode I: (c) $\unicode[STIX]{x1D6FA}=0.01$ , $\mathit{Re}=200$ , (d) $\unicode[STIX]{x1D6FA}=0.01$ , $\mathit{Re}=270$ , (e) $\unicode[STIX]{x1D6FA}=0.2$ , $\mathit{Re}=200$ , (f) $\unicode[STIX]{x1D6FA}=0.2$ , $\mathit{Re}=240$ ; and two cases for mode II: (g) $\unicode[STIX]{x1D6FA}=0.9$ , $\mathit{Re}=200$ , (h) $\unicode[STIX]{x1D6FA}=0.9$ , $\mathit{Re}=270$ .

6.1 Weakly nonlinear approach

Following Meliga et al. (Reference Meliga, Chomaz and Sipp2009a ), we pose $\unicode[STIX]{x1D716}^{2}=1/Re^{SS}-1/Re$ and assume this parameter to be small. The flow is thus expanded in terms of this parameter as follows:

(6.1) $$\begin{eqnarray}\boldsymbol{q}=[\boldsymbol{u},p]=\boldsymbol{q}_{0}+\unicode[STIX]{x1D716}\boldsymbol{q}_{1}+\unicode[STIX]{x1D716}^{2}\boldsymbol{q}_{2}+\unicode[STIX]{x1D716}^{3}\boldsymbol{q}_{3}+\cdots \,.\end{eqnarray}$$

In order to introduce the effect of rotation in the method, we also have to assume that the rotation is small. Since in the non-rotating case the amplitude equation derived in Meliga et al. (Reference Meliga, Chomaz and Sipp2009a ) arises from the resolution of the problem at order $\unicode[STIX]{x1D716}^{3}$ , we anticipate that $\unicode[STIX]{x1D74E}$ has to be of the same order of magnitude, so that we set

(6.2) $$\begin{eqnarray}\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D74E}_{0},\end{eqnarray}$$

where $\unicode[STIX]{x1D74E}_{0}$ is the rescaled rotation rate of order one.

Following the standard weakly nonlinear expansion for such flows (Sipp & Lebedev Reference Sipp and Lebedev2007; Meliga et al. Reference Meliga, Chomaz and Sipp2009a ), we rely on a multi-time scale technique and assume $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}^{2}t$ , where $\unicode[STIX]{x1D70F}$ and $t$ are respectively the slow and fast time scales. Injection of this decomposition into the governing equations (2.1) along with the boundary conditions (2.2) leads to various problems to be solved for $Re=Re^{SS}$ at different orders $\unicode[STIX]{x1D716}^{j}$ , with $j=0,1,2,3,\ldots .$

1. (i) At order $\unicode[STIX]{x1D716}^{0}$ , we obtain the nonlinear problem governing $\boldsymbol{q}_{0}$ , the axisymmetric base flow with no rotation.

2. (ii) At order $\unicode[STIX]{x1D716}$ , we obtain a linear problem which can be written in a symbolic way as ${\mathcal{L}}\boldsymbol{q}_{1}=0$ . Since we are at the threshold $Re^{SS}$ , this equation has non-trivial solution with azimuthal wavenumbers $\pm 1$ (the marginal eigenmodes). The solution at this order is thus a linear combination of the global modes at $Re=Re^{SS}$ , and is taken as

   (6.3) $$\begin{eqnarray}\boldsymbol{q}_{1}=\hat{A}(\unicode[STIX]{x1D70F})\hat{\boldsymbol{q}}(r,x)\text{e}^{\text{i}\unicode[STIX]{x1D719}}+\text{c.c.},\end{eqnarray}$$
   where $\hat{\boldsymbol{q}}$ is the eigenvector of the azimuthal wavenumber $(m=+1)$ and $\hat{A}$ is a complex amplitude that depends on the slow time scale. It should be noted that the eigenvector $\hat{\boldsymbol{q}}$ can be normalized so that the real and imaginary parts of the amplitude $A$ correspond to the $y$ and $z$ components of the lift (see Fabre et al. Reference Fabre, Auguste and Magnaudet2008).

3. (iii) At order $\unicode[STIX]{x1D716}^{2}$ , the solution contains base-flow viscous corrections due the departure from $Re^{SS}$ as well as quadratic terms resulting from the self-interactions of the global mode. Using the notations of Meliga et al. (Reference Meliga, Chomaz and Sipp2009a ), this reads as

   (6.4) $$\begin{eqnarray}\boldsymbol{q}_{2}=\boldsymbol{q}_{\unicode[STIX]{x1D6FF}}+\boldsymbol{q}_{AA\ast }|A|^{2}+(\boldsymbol{q}_{AA}A^{2}\text{e}^{2\text{i}\unicode[STIX]{x1D719}}+\text{c.c.}).\end{eqnarray}$$

4. (iv) At order $\unicode[STIX]{x1D716}^{3}$ , one obtains a linear problem which can be written in a symbolic way as follows:

   (6.5) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{q}_{3}+(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}A\hat{\boldsymbol{q}}\text{e}^{\text{i}\unicode[STIX]{x1D719}}+\text{c.c.})={\mathcal{N}}(\boldsymbol{q}_{2},\boldsymbol{q}_{1})+{\mathcal{L}}\boldsymbol{q}_{3},\end{eqnarray}$$
   where ${\mathcal{N}}(\boldsymbol{q}_{2},\boldsymbol{q}_{1})$ represents the nonlinear interaction between the solutions at the two previous orders, and ${\mathcal{L}}$ is the same linear operator as for the problem at order one. This problem is again similar to the case of Meliga et al. (Reference Meliga, Chomaz and Sipp2009a ), except for the fact that $\boldsymbol{q}_{3}$ has to verify the boundary condition corresponding to a steady rotation, namely
   (6.6) $$\begin{eqnarray}\boldsymbol{u}_{3}=\unicode[STIX]{x1D74E}_{0}\boldsymbol{z}\times \boldsymbol{r}.\end{eqnarray}$$
   As the linear operator governing the system (6.5) is singular, the usual procedure to remove the secular terms consists of taking a scalar product of the forcing terms with the adjoint of the problem. This yields a single differential equation for the amplitude $A$ , with the form
   (6.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D716}^{2}A-\unicode[STIX]{x1D707}\left|A\right|^{2}A+a\unicode[STIX]{x1D74E},\end{eqnarray}$$
   where $\unicode[STIX]{x1D706}$ , $\unicode[STIX]{x1D707}$ and $a$ are all real. Here, the coefficients $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}$ contain the effect of the forcing terms corresponding to the viscous correction ( $\boldsymbol{q}_{\unicode[STIX]{x1D6FF}}$ ) and the quadratic terms constituting $\boldsymbol{q}_{2}$ . These terms are computed exactly as in Meliga et al. (Reference Meliga, Chomaz and Sipp2009a ), so their detailed expression need not be given. On the other hand, the coefficient $a$ comes from the non-homogeneous boundary condition.

   Let $\langle \boldsymbol{v}_{1},\boldsymbol{v}_{2}\rangle _{\unicode[STIX]{x2202}\mathscr{D}_{sph}}=\int _{\unicode[STIX]{x2202}\mathscr{D}_{sph}}\boldsymbol{v}_{1}^{\ast }\boldsymbol{\cdot }\boldsymbol{v}_{2}$ , i.e. the Hermitian scalar product between two complex vectors $\boldsymbol{v}_{1}$ and $\boldsymbol{v}_{2}$ on the surface $\unicode[STIX]{x2202}\mathscr{D}_{sph}$ . Then, the coefficient $a$ is computed as follows:

   (6.8) $$\begin{eqnarray}a=\left<\frac{1}{Re}\unicode[STIX]{x1D735}\boldsymbol{u}^{\dagger }\boldsymbol{\cdot }\boldsymbol{n},\frac{1}{\unicode[STIX]{x1D74E}_{0}}\hat{\boldsymbol{u}}_{3}\right>_{\unicode[STIX]{x2202}\mathscr{D}_{sph}},\end{eqnarray}$$
   where $\hat{\boldsymbol{u}}_{3}=\unicode[STIX]{x1D74E}_{0}(x/2\boldsymbol{e}_{r}+\text{i}x/2\boldsymbol{e}_{\unicode[STIX]{x1D719}}-r/2\boldsymbol{e}_{3})$ (in cylindrical coordinates) because of the no-slip condition on the sphere at order $\unicode[STIX]{x1D716}^{3}$ .

All of the calculations presented in this section are carried out using a finite-element method based on the $\text{Freefem}++$ software, along the lines of Meliga et al. (Reference Meliga, Chomaz and Sipp2009a ). The methodology was adapted and validated in our team for various problems including spheres (Fabre et al. Reference Fabre, Tchoufag, Citro, Giannetti and Luchini2016), bubbles (Tchoufag, Fabre & Magnaudet Reference Tchoufag, Fabre and Magnaudet2013) and disks (Tchoufag et al. Reference Tchoufag, Fabre and Magnaudet2014b ). The reader is referred to these latter papers for more details, especially the supplementary material related to Tchoufag et al. (Reference Tchoufag, Fabre and Magnaudet2014b ) where the weakly nonlinear approach was set up for the more general case of a freely moving object, and the associated code was made available to potential users on an open-source basis. For the present case, convergence of the results was checked by varying the mesh details, including the grid density in the close vicinity of the sphere and the dimensions of the domain, along the lines of what was done for ellipsoidal bubbles (table II of Tchoufag et al. (Reference Tchoufag, Fabre and Magnaudet2013)) or disks (table II of Tchoufag, Fabre & Magnaudet (Reference Tchoufag, Fabre and Magnaudet2014a )). The retained mesh is the same as used in Fabre et al. (Reference Fabre, Tchoufag and Magnaudet2012, Reference Fabre, Tchoufag, Citro, Giannetti and Luchini2016). It provides the critical Reynolds for the onset of instability with four digits of accuracy, i.e. $Re^{SS}=212.6$ , which is identical to the value given in Meliga et al. (Reference Meliga, Chomaz and Sipp2009b ). Numerical values of the coefficients entering the amplitude equations are obtained with three-digit accuracy, namely $\unicode[STIX]{x1D706}=126$ , $\unicode[STIX]{x1D707}=94.3$ , $a=-0.0247$ .

To study the dynamics predicted by the amplitude equation (6.7), we first look for equilibrium solutions, then consider their asymptotic trends as $\unicode[STIX]{x1D74E}\rightarrow 0$ and finally investigate their stability by examining the behaviour of small-amplitude perturbations with respect to them. For the latter purpose, one must ensure that the amplitude $A$ is defined as a complex number, the real and imaginary parts contributing to the lift in the $y$ and $z$ directions, so that the stability with respect to both in-plane and out-of-plane perturbation has to be distinguished. The details are presented in appendix A. In summary, it is shown that the system possesses three branches of equilibrium solutions which verify the following properties.

1. (i) The first one, denoted $A_{01}$ , is real and positive, and exists for all values of $\unicode[STIX]{x1D716}^{2}$ . For very small rotation, it displays the following asymptotic behaviours:

   (6.9a,b ) $$\begin{eqnarray}A_{01}\approx -\frac{a\unicode[STIX]{x1D74E}}{\unicode[STIX]{x1D706}\unicode[STIX]{x1D716}^{2}}\quad \text{as }\unicode[STIX]{x1D74E}\rightarrow 0\text{ with }\unicode[STIX]{x1D716}^{2}<0;\quad A_{01}\approx \sqrt{\frac{\unicode[STIX]{x1D706}\unicode[STIX]{x1D716}^{2}}{\unicode[STIX]{x1D707}}}\quad \text{as }\unicode[STIX]{x1D74E}\rightarrow 0\text{ with }\unicode[STIX]{x1D716}^{2}>0.\end{eqnarray}$$
   This solution is stable to both in-plane and out-of-plane perturbations.

2. (ii) The second branch, denoted $A_{02}$ , is real and negative, and emerges from a saddle–node for $\unicode[STIX]{x1D716}^{2}>3\cdot 2^{-2/3}a^{2/3}\unicode[STIX]{x1D707}^{2/3}\unicode[STIX]{x1D74E}^{2/3}/\unicode[STIX]{x1D706}$ . For very small rotations, it displays the following asymptotic behaviour:

   (6.10) $$\begin{eqnarray}A_{02}\approx -\frac{a\unicode[STIX]{x1D74E}}{\unicode[STIX]{x1D706}\unicode[STIX]{x1D716}^{2}}\quad \text{as }\unicode[STIX]{x1D74E}\rightarrow 0\text{ with }\unicode[STIX]{x1D716}^{2}>0.\end{eqnarray}$$
   This second branch is unstable to both in-plane and out-of-plane perturbations.

3. (iii) The third branch, denoted $A_{03}$ , originates from the same saddle–node as the previous one, and for very small rotations, it has the asymptotic behaviour

   (6.11) $$\begin{eqnarray}A_{03}\approx -\sqrt{\frac{\unicode[STIX]{x1D706}\unicode[STIX]{x1D716}^{2}}{\unicode[STIX]{x1D707}}}\quad \text{as }\unicode[STIX]{x1D74E}\rightarrow 0\text{ with }\unicode[STIX]{x1D716}^{2}>0.\end{eqnarray}$$
   This last branch is stable with respect to in-plane perturbations but unstable with respect to out-of-plane perturbations.

It should be noted that the fact that the three solutions are all real means that the corresponding structures are symmetrical with respect to the transverse ( $x$ – $y$ ) plane. Accordingly, these states are characterized by a constant lift in the $y$ direction, given by their amplitude $A$ in the original scaling.

Figure 9.

The lift force as a function of the Reynolds number for different rotation rates. The symbols correspond to DNS and the lines correspond to the normal form (6.7) of the imperfect pitchfork bifurcation. The solid lines represent the stable $A_{01}$ solution, while the dashed lines represent the unstable $A_{02}$ and $A_{03}$ solutions.

6.2 Comparison with 3D numerical simulations

We recall that due to the normalization of the direct modes, the lift force is directly given by $A$ , allowing one to compare, directly, the solutions of the amplitude equations with the numerical results of figure 4. The comparison is presented in figure 9, which shows that the base flow discussed in § 5.1 is in perfect agreement with the predictions of branch $A_{01}$ of the weakly nonlinear model. In addition, the weakly nonlinear approach allows the unveiling of the existence of other disconnected branches of steady solutions for a rotating sphere, corresponding to branches $A_{02}$ and $A_{03}$ . As shown in figure 9, the global approach also allows the existence of these states to be confirmed, at least for $\unicode[STIX]{x1D74E}=0.001$ and $\unicode[STIX]{x1D74E}=0.0001$ . On the other hand, for larger values of the rotation rate, these branches become difficult to track. Since they are always unstable, we disregarded these branches in the stability approach described in § 5.2 and concentrated only on the primary branch $A_{01}$ .

As expected, the discrepancies between the theory and the DNS data in figure 9 increase with the rotation rate. As the sphere spins faster, the departure from the axisymmetric base flow is greater even for Reynolds numbers far smaller than $Re^{SS}$ , and the steady global mode bifurcating from the fixed sphere configuration progressively loses its relevance. It should be noted also that, the weakly nonlinear expansion having been derived in the vicinity of the steady bifurcation, its comparison with DNS data is expected to be less satisfactory as the departure from criticality increases. Figure 9 confirms such behaviour, for all three solution branches and for all rotation rate values.