2. Configuration

We investigate numerically the streamwise impingement of a canonical Batchelor vortex onto an elongated axisymmetric obstacle. A cylindrical coordinate system $(R,\theta ,Z)$ is introduced, whose axial component is aligned along the obstacle centreline. In this coordinate system, the fluid motion, associated with the state vector $\boldsymbol {q}=(U_{R},U_{\theta },U_{Z},P)$, is governed by the Navier–Stokes equations and characterized by a homogeneous and incompressible fluid of kinematic viscosity $\nu$. The variables are scaled by the physical length and velocities defined as the vortex core size $R_v$ and a uniform axial velocity $U_{\infty }$ at the inlet, respectively. These reference scales refer to the upstream vortex defined as a Batchelor vortex, which writes in dimensionless variables

(2.1)\begin{equation} U_R(R)=0,\quad U_{\theta}(R)=S\frac{\left(1-\textrm{e}^{-R^{2}}\right)}{R},\quad U_Z(R)=1,\quad \text{on}\ \varGamma_{in}, \end{equation}

where $S$ is the swirl number defined as $S = R_v \varOmega _0 / U_{\infty }$ with $\varOmega _0$ the vorticity and $R_v$ the radius of the vortex. Using this scaling, the dimensionless Navier–Stokes equations writes in vector formulation as

(2.2)\begin{equation} \left.\begin{gathered} \frac{\partial \boldsymbol{U}}{\partial t} + \left( \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{ U} = -\boldsymbol{\nabla} P + Re^{-1} \nabla^{2}\boldsymbol{U}, \quad \text{in}\ \varOmega,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \quad \text{in}\ \varOmega, \end{gathered}\right\} \end{equation}

with the Reynolds number defined as $Re = R_v U_{\infty } / \nu$. We solve these equations in a large computational domain $R_{max}=10 R_s$ and $Z_{free} = 19 R_{s}$, where $R_s$ is the obstacle radius. These conditions ensure the vortex in front of the obstacle to be free of inlet and radial disturbances, by isolating the vortex from the free outflow condition imposed on the external boundary $\varGamma _{ext}$ (see Pasche et al. Reference Pasche, Gallaire and Avellan2018b). The Batchelor vortex, (2.1), is imposed at the inlet $Z=0$ and can be aligned or shifted by an offset $D$ regarding the obstacle axis. The obstacle extends up to the outlet boundary $\varGamma _{out}$ to avoid flow unsteadiness due to the wake of the obstacle. The Navier–Stokes equations are solved in a 2-D axisymmetric domain and in a 3-D domain. In the 2-D axisymmetric domain, a free outflow ($(-P\boldsymbol {I} + Re^{-1}\boldsymbol {\cdot }(\boldsymbol {\nabla } \boldsymbol {V}))\boldsymbol {\cdot }\boldsymbol {n}=0$) is imposed on the outlet boundary $\varGamma _{out}$, while in the 3-D domain, a convective boundary condition is imposed on $\varGamma _{out}$ ($\partial _t \boldsymbol {V} + \boldsymbol {U}_c\boldsymbol {\cdot } \partial _{\boldsymbol {n}} \boldsymbol {V} = 0$). We fix the convective velocity to be equal to the free-stream velocity $\boldsymbol {U}_c=(0,0,1)$. This change of the outlet boundary condition does not impact the flow dynamics as it is imposed far away from the obstacle nose $Z_{obs}=10R_s$.

In this configuration, the obstacle that the vortex must accommodate is defined by a swirl profile radius ratio, in short aspect ratio $h = R_v/R_s$. The different parameters investigated are summarized in table 1. We consider only obstacles that have a larger radius than the vortex core inducing a strong topology change of the vortex. The adverse pressure gradient modifies both the viscous core and the potential vortex region to force the deflection of the vortex circumventing the obstacle.

Table 1.

Parameters of the present numerical analysis, $h$ the aspect ratio between half the vortex core size and the obstacle radius, $S$ the swirl number, $Re$ the Reynolds number based on the vortex and $D$ the lateral displacement of the vortex.

We introduce a second non-dimensional length scale $\tilde {Z} = Z \times h$ and $\tilde {R}=R \times h$, which is more appropriate to portray the flow solution because it enables a generalized representation of the computational domain for any aspect ratio studied. Thus, in the figures of the present study, the obstacle radius becomes unity and the vortex core size adjusts consequently. Any variables specified with $\tilde {(\cdot )}$ use the newly introduced length scale, which is just a rescaling factor.

3. Methods of analysis

We consider both 2-D axisymmetric configuration and 3-D configuration of the computational domain. The 2-D axisymmetric analysis allows us to perform bifurcation, modal and weakly nonlinear (WNL) analysis. The 3-D DNS analysis validates the 2-D axisymmetric simulations and provides a deeper understanding of symmetry breaking of the flow when the vortex is laterally displaced.

3.1. Two-dimensional axisymmetric analysis

We look for laminar solutions of the flow by solving the 2-D axisymmetric steady Navier–Stokes equations for a wide range of Reynolds numbers and swirl numbers. These solutions are explored by performing a bifurcation analysis on the Reynolds number for several specified swirl numbers. An arc-length continuation method is used to follow the stable and unstable branches of the base flow. As an axisymmetric coordinate system is used, critical points in the space of the bifurcation parameters can be associated with different azimuthal wavenumbers $m$. The boundaries of these critical points are tracked using a fold tracking algorithm and a Hopf tracking algorithm (Salinger et al. Reference Salinger, Burroughs, Pawlowski, Phipps and Romero2005), for the azimuthal wavenumber $m=0$ and $m=1$, respectively. The higher azimuthal wavenumbers are not tracked in the present study because the azimuthal wavenumber $m=2$ becomes dominant for high swirl number, as already observed in Sarpkaya (Reference Sarpkaya1974) and Meliga et al. (Reference Meliga, Gallaire and Chomaz2012). The results of a convergence study are presented in appendix A. These modes are highlighted in the modal analysis too. Linear global stability analysis is performed (see Pasche et al. (Reference Pasche, Gallaire and Avellan2018b) for details and validation), which allows us to access the full spectrum of the flow while tracking algorithms follow the marginally unstable eigenmode of specific azimuthal wavenumbers.

In the linear global stability analysis, the flow is decomposed into a base flow $\boldsymbol {q}_0$ and a disturbance $\boldsymbol {q}'$. The disturbance $\boldsymbol {q}'$ is expanded in normal modes for different azimuthal wavenumbers $m\in \mathbb {Z}$,

(3.1)\begin{equation} \boldsymbol{q}'(R,\theta,Z,t)=\boldsymbol{q}_1(R,Z)\exp(-\textrm{i}\omega t + \textrm{i}m \theta)+\textrm{c.c.}, \end{equation}

where c.c. is the complex conjugate. The substitution of this definition in the governing equations (2.2), leads to an eigenvalue equation expressed in a compact form as

(3.2)\begin{equation} \left( (\omega_i-\textrm{i} \omega_r)\mathcal{N} + \mathcal{L}_m(\boldsymbol{U}) \right)\boldsymbol{q}_1=\boldsymbol{0}, \quad \text{in}\ \varOmega_{a}, \end{equation}

where $\mathcal {L}_m$ is the operator for the linearized Navier–Stokes equations of azimuthal wavenumber $m$. This eigenvalue problem results in eigenvalues and eigenmodes that will lead to self-sustained instability for unstable cases. This problem is solved using the finite element library FreeFem$++$ (Hecht Reference Hecht2012), see appendix B.

3.2. Weakly nonlinear analysis

A multiple time scale analysis is performed at the critical Reynolds value $Re_c$, of the marginally unstable eigenmode $\boldsymbol {q}_{1}$ of eigenvalue $-\textrm {i}\omega _c$, to investigate the WNL mechanism behind the flow bifurcations. An asymptotic expansion of the flow field $\boldsymbol {q} = (U_R,U_{\theta },U_Z,P)$ regarding the small parameter $\epsilon$ is considered

(3.3)\begin{equation} \boldsymbol{q} = \boldsymbol{q}_0+\epsilon \boldsymbol{q}_1+\epsilon^{2} \boldsymbol{q}_2+\epsilon^{3} \boldsymbol{q}_3+ \text{h.o.t}, \end{equation}

where h.o.t stands for higher-order terms, and where $0<\epsilon \ll 1$ and defines the distance from the thresholds (3.4). In addition, we assume that the Reynolds number departs from criticality at order $\epsilon ^{2}$, yielding

(3.4)\begin{equation} \frac{1}{Re} = \frac{1}{Re_c}-\epsilon^{2}. \end{equation}

The unsteady dynamics is decomposed into a fast time scale $t_0$ and slow time scale $T_{2i}$, the latter start acting at $\epsilon ^{2}$ too:

(3.5)\begin{equation} t = t_0 + \sum_{i \geqslant 1} \epsilon^{2i}T_{2i}. \end{equation}

This scaling has been used for axisymmetric flows (see Meliga et al. Reference Meliga, Gallaire and Chomaz2012; Citro et al. Reference Citro, Tchoufag, Fabre, Giannetti and Luchini2016). Substituting the asymptotic expansion, (3.3) in the governing equations (2.2) and using the chain rule for the derivative involving multiple time scales results in systems of equations at different epsilon orders that are solved. These equations are summarized in appendix B and the reader is invited to read Meliga et al. (Reference Meliga, Gallaire and Chomaz2012) and Sipp & Lebedev (Reference Sipp and Lebedev2007) for additional explanations.

The WNL analysis leads to the definition of an amplitude equation that describes slow modulation in time and space of patterns out of equilibrium. This approach is only valid for WNL conditions, which means close to the threshold of instability $Re\approx Re_c$, and provides only qualitative results in the stronger nonlinear regime. Nevertheless, the amplitude equation describes completely the pattern forming in its range of validity. In the present case, we will show that the flow dynamics follows a Landau amplitude equation, which writes, up to the fifth order (see Fujimura Reference Fujimura1991; Dusek, Gal & Fraunié Reference Dusek, Gal and Fraunié1994),

(3.6)\begin{equation} \frac{\partial a}{\partial t} - \epsilon^{2}(\lambda_1+\epsilon^{2}\lambda_2)a - (\mu_1+ \epsilon^{2} \mu_2) a|a|^{2} - \gamma a|a|^{4} =0, \end{equation}

with $\lambda$, $\mu$ and $\gamma$ complex coefficients and $a$ the complex amplitude.

The amplitude $a$, solution of (3.6), allows a direct comparison to DNS results when the energy norm of the eigenvector matches the energy of the DNS disturbance. The DNS solution is a real value velocity field and the WNL analysis is performed in the complex number set. Thus, only half the energy of the WNL analysis has to be used to compare with the DNS, i.e. the complex conjugate part of the product is dropped. The final solution of the WNL analysis writes

(3.7)\begin{align} \boldsymbol{q} &= \boldsymbol{q}_0 + a \boldsymbol{q}_{1a} \exp({\textrm{i}\theta-i\omega_c t}) + a|a|^{2} \boldsymbol{q}_{3a|a|^{2}} \exp({\textrm{i}\theta-\textrm{i}\omega_c t})\nonumber\\ &\quad + a|a|^{4}\boldsymbol{q}_{5a|a|^{4}} \exp({\textrm{i}\theta-\textrm{i}\omega_c t}) + \text{h.o.t} + \textrm{c.c.} \end{align}

Assuming the decomposition: $a= \rho (t)\,\textrm{e}^{\textrm {i}\phi (t)}$, the following equation holds for the disturbance amplitude:

(3.8)\begin{equation} \sqrt{(\boldsymbol{q}-\boldsymbol{q}_0)^{2}} = \sqrt{\rho \rho^{*} (\textrm{e}^{\textrm{i}\phi} \boldsymbol{q}_{1a} \ \textrm{e}^{\textrm{i}\theta-\textrm{i}\omega_c t}, \textrm{e}^{-\textrm{i}\phi} \boldsymbol{q}^{*}_{1a} \ \textrm{e}^{-\textrm{i}\theta+\textrm{i}\omega_c t})} + {\rm h.o.t}, \end{equation}

where the left-hand side is the disturbance amplitude of the DNS and the right-hand side is the disturbance amplitude of the WNL analysis. It results in

(3.9)\begin{equation} \sqrt{(\boldsymbol{q}-\boldsymbol{q}_0)^{2}} = \sqrt{\rho\rho^{*} + \rho\rho^{*} |\rho|^{2} + {\rm h.o.t} } \end{equation}

using the following normalization of the adjoint eigenmode:

(3.10)\begin{equation} (\boldsymbol{q}_{1a},\boldsymbol{q}^{*}_{1a}) = \int_\varOmega \boldsymbol{q}_{1a} \boldsymbol{\cdot} \boldsymbol{q}^{*}_{1a} = 1, \end{equation}

which approximates very well the nonlinear disturbance amplitude at threshold. However, it can differ further away from the threshold because of nonlinear interactions. Once again, the computation of the WNL expansion is performed using the FreeFem$++$ library, involving the resolution of eigenvalue, adjoint eigenvalue and resolvent problems, which are expressed in appendix B.

3.3. Three-dimensional direct numerical flow simulations

Direct numerical flow simulations of the configuration displayed in figure 1 have been performed using NEK5000, a spectral element solver developed by Fischer, Lottes & Kerkemeier (Reference Fischer, Lottes and Kerkemeier2008). An O-grid mesh is used to mesh the obstacle. We use a convective boundary condition at the outlet that avoids backward waves and helps convergence for this swirling flow. Similar simulations of the vortex breakdown have already been performed by the authors in Pasche et al. (Reference Pasche, Gallaire and Avellan2018b). A convergence analysis is presented in appendix A.

Figure 1.

Half of the computational domain.

4. Bifurcation analysis

Let us start by introducing an axisymmetric flow solution for $S = 1.05$ and $Re = 255.5$ with and without a downstream obstacle having an aspect ratio of $h=0.5$, see figure 2. The flows are illustrated by surface streamlines coloured by the axial velocity magnitude on the upward $R$ coordinate and by the azimuthal velocity magnitude on the downward $R$ coordinate. In figure 2(a), the flow is purely columnar and the vortex progressively dissipates by viscous diffusion as it propagates downstream. On the contrary, the formation of a recirculation bubble is observed when an obstacle is located on the vortex trajectory for the parameters considered. This obstacle generates an adverse pressure gradient which originates in this change of flow topology. The recirculation bubble is reminiscent of a bubble vortex breakdown, which is characterized by a negative velocity on the axis defining the starting location of the bubble (see Leibovich & Stewartson Reference Leibovich and Stewartson1983). Behind the first vortex breakdown bubble, a typical swirling wake flow develops, which attaches to the obstacle, generating a second recirculation region.

Figure 2.

Two-dimensional axisymmetric vortex flow solutions at $S = 1.05$, $Re = 255.5$ without impingement (a), and with impingement onto a downstream obstacle $h=0.5$ (b). On the upward $R$ coordinate, the figure shows surface streamlines coloured by the axial velocity magnitude, and on the downward $R$ coordinate the figure shows the azimuthal velocity component. Both examples are stable solutions. Panel (b) represents a typical solution observed in the present numerical study.

The adverse pressure gradient decelerates mainly the axial velocity on the centreline, while the azimuthal velocity remains almost constant before the vortex breaks down (see Delery Reference Delery1994; Pasche et al. Reference Pasche, Gallaire, Dreyer and Farhat2014). In an axisymmetric system, the azimuthal velocity is mainly affected by the radial momentum through the radial pressure gradient. The transport of the azimuthal momentum is ensured by nonlinear advection. In contrast, the axial pressure gradient modifies straightforwardly the axial momentum. Thus, the axial velocity is first decelerated, and then the azimuthal velocity is modified as the recirculation bubble is reached.

In figure 2(b), we define the recirculation length, which measures the length between the obstacle front and the upstream zero value of the axial velocity on the axis. This recirculation length is used, in the present study, as a bifurcation parameter to identify the different bifurcation branches of the flow. The associated bifurcation diagrams are presented for two configurations: an aspect ratio $h=0.5$ in figure 3(a) and an aspect ratio $h=0.25$ in figure 3(b). The black curves, which represent bifurcation branches, show the recirculation length of the bubble on the vertical axis as a function of the Reynolds number $Re$, for different values of the swirl number $S$. As the Reynolds number is increased, the recirculation bubble propagates upstream, following a cusp shape in the bifurcation diagram. Such a transition is similar to the saddle-node bifurcation observed by Lopez (Reference Lopez1994) resulting from the presence of a radial contraction in a numerical tube experiment. In figure 3, the forward branches represent stable solutions, while the backward branches represent unstable solutions. For $h=0.5$, the transition is smooth, while for $h=0.25$ several unstable branches appear, increasing the complexity of the bifurcation diagram. In both cases, as the swirl number is increased, the first bifurcation appears at smaller Reynolds numbers and the cusp shapes are damped, tending to the bubble vortex breakdown solution of the canonical vortex observed in tube experiments (see Meliga & Gallaire Reference Meliga and Gallaire2011). At lower Reynolds number, we observe in figure 3 that the bifurcation curves are flat, implying that the flow is attached to the obstacle. On the contrary, at the largest displayed Reynolds numbers, the recirculation lengths remain constant because the back flow approaches the inlet of the computational domain. A larger obstacle size shifts the saddle-node bifurcation to smaller Reynolds number as already reported by Delery (Reference Delery1994) and Pasche et al. (Reference Pasche, Gallaire, Dreyer and Farhat2014).

Figure 3.

Axisymmetric bifurcation map of the vortex impinging a downstream obstacle for two aspect ratios; $h=0.5$ (a) and $h=0.25$ (b). The red dots show the saddle-node bifurcation threshold, and the blue dots show the Hopf bifurcation threshold.

In figure 3, we also present the outcome of the turning point and the Hopf point tracking algorithms. The red dots are the neutrally stable limit of the flow for the azimuthal wavenumber $m=0$. These curves accurately fall on the stable–unstable limit of the bifurcation branches, which has already been observed for swirling flows (Meliga & Gallaire Reference Meliga and Gallaire2011; Pasche, Avellan & Gallaire Reference Pasche, Avellan and Gallaire2018a). The blue dots represent the neutrally stable limit of the flow concerning the azimuthal wavenumber $m=1$, i.e. limit cycle solutions.

The detailed bifurcation curves for a $S=1.0$ are shown in figure 4. The unstable branches are represented by discontinuous curves and the solid lines represent stable branches. The critical Reynolds number at which the Hopf bifurcation occurs is indicated with blue crosses. In both cases, the turning points giving birth to the second stable branches appear before the critical Reynolds number of the Hopf bifurcation.

Figure 4.

Bifurcation curves of the axisymmetric vortex flow impinging an obstacle of aspect ratio $h=R_v/R_s =0.5$ (a) and $h=0.25$ (b), for a swirl number of $S=1.0$. The dashed curves represent the unstable branch, while the solid curves represent the stable branches. The blue crosses exhibit the critical Reynolds number $Re_c$ of the Hopf bifurcation, and the green dots exhibit the selected bifurcation points represented in the subsequent figures 5 and 6.

The flow solutions of the first and second stable branches are illustrated in figure 5, and correspond to the bifurcation parameters of the green dots in figure 4. On the first stable branch, at low Reynolds number, the classical stagnation point is recovered. As the Reynolds number increases, a small recirculation region is observed in front of the obstacle as illustrated in figure 5(a) at $Re=350$, $S=1.0$ and $h=0.5$. This recirculation slowly propagates upstream until the turning point is reached, and remains attached to the obstacle. On the unstable branch, the recirculation zone increases without the formation of a separate bubble, to merge with the first axisymmetric solution on the second stable branch, see the turning point in figure 5(b) at $Re=263.2$ and $S=1.0$. Then, the recirculation zone separates from the obstacle and defines a detached bubble as illustrated in figure 5(c) at $Re=350$ and $S=1.0$ (note the relative analogy with the upper branch solution depicted in figure 2(b) for $S=1.05$ and $Re=255.5$). The transition from the attached to the detached recirculation bubble occurs smoothly as the Reynolds number increases. The second stable branch is, therefore, associated with the development of the bubble vortex breakdown. Regarding the smaller aspect ratio $h=0.25$, whose first stable branch solution is displayed in figure 5(d) at $Re=325.1$ and $S=1.0$, a similar scenario is observed up to the second stable branch: the recirculation bubble detaches from the obstacle and moves upstream letting a second recirculation region appear, which remains attached to the obstacle.

Figure 5.

Axisymmetric vortex flow solutions for $h=0.5$ and $S=1.0$, on the first stable branch $Re = 350.0$ (a), at the threshold of the second stable branch $Re=263.2$ (b) and on the second stable branch $Re= 350$ (c), as well as the solution for $h=0.25$ and $S=1.0$ on the first stable branch $Re=325.1$ (d). These solutions correspond to the green dots in figure 4.

Moreover, we present the eigenvalue spectra of the solutions corresponding to the green dots displayed in figures 4 and 5. On the first stable axisymmetric branch at $h=0.5$ and $Re=350$ (see figure 6a) the eigenvalue spectrum displays an unstable eigenvalue of azimuthal wavenumber $m=1$. This eigenvalue growth rate is positive in agreement with the critical Reynolds number at $Re_c=324.5$ of the Hopf tracking method. The eigenvalue spectrum at $h=0.5$ and $Re = 263.2$ for the second stable branch displays several unstable eigenvalues of azimuthal wavenumber $m=1$ and $m=2$ showing the hydrodynamically unstable flow field (figure 6b). As we move along the second stable branch, the number of unstable eigenmodes increases and their growth rates are strengthened (figure 6c). The eigenvalue spectrum at $h=0.25$ displays two unstable eigenvalues on the first stable branch: a first unstable eigenmode of azimuthal wavenumber $m=1$ and a second of $m=2$, see figure 6(d). The latter has a smaller growth rate and a lower frequency than the unstable mode $m=1$. We observe on the eigenvector plot (not shown here) that the second unstable mode develops in the recirculation region in front of the obstacle while the first unstable mode is the helical unstable mode of the spiral vortex breakdown.

Figure 6.

Eigenvalue spectra of the axisymmetric vortex solution $h=0.5$, $S=1.0$ for the first stable branch at $Re = 350.0$ (a), for the second stable branch at $Re = 263.2$ (b) and for the second stable branch at $Re=350$ (c), as well as the first stable branch for $h=0.25$, $S=1$, $Re=325.1$ (d). These solutions correspond to the green dots in figure 4.

We observe, finally, that in contrast to the saddle-node bifurcation of Lopez (Reference Lopez2006), the second stable axisymmetric branch of a vortex impinging an obstacle is strongly unstable and cannot exist in 3-D configurations. Only the first axisymmetric stable branch can exist. While the Reynolds or swirl numbers are increased, a recirculation zone ahead of the obstacle develops, which is submitted to a Hopf bifurcation of azimuthal wavenumber $m=1$ as a first instability.

5. Subcritical Hopf bifurcation

The flow, issued from the 3-D direct numerical flow simulations, is illustrated by the lambda 2 criterion in figure 7 for $h=0.25$ and $S=1.0$. The development of a single helical spiral in front of the obstacle is shown. This mode is also obtained in the 2-D axisymmetric solution represented by the eigenmode of the most unstable eigenvalue having the same winding direction and temporal rotation as the 3-D solution. In addition, the spiral suddenly forms from the centreline of the vortex without the development of a recirculation bubble. Such a spiral vortex breakdown has the same origin as the single spiral appearing in the wake of a recirculation bubble, (see Gallaire et al. Reference Gallaire, Ruith, Meiburg, Chomaz and Huerre2006). Both are hydrodynamic instabilities of a swirling wake flow.

Figure 7.

Axial vorticity of the 3-D solution for $S=1.0$, $Re=224.6$ and $h=0.25$.

The time series of the flow and the Fourier amplitude spectrum of this 3-D solution is portrayed in figure 8. A limit cycle solution is observed with a well-defined frequency. However, the periodic solution displayed in figure 8 at $Re=218.75$ appears for a smaller Reynolds number than the critical Reynolds number $Re=224.6$ identified by modal analysis. This shift is not a numerical inaccuracy. The transition from a steady to an unsteady dynamics of the vortex streamwise impingement onto an axisymmetric obstacle is caused by a subcritical Hopf bifurcation, in contrast to supercritical Hopf bifurcation of the spiral vortex breakdown of Gallaire et al. (Reference Gallaire, Ruith, Meiburg, Chomaz and Huerre2006).

Figure 8.

Time series (a) and amplitude Fourier spectrum (b) of the 3-D DNSs at $(\tilde {X},\tilde {Y},\tilde {Z})=(0.5,0.0,18.0)$ for the axial velocity and $S=1.0$ and $h=0.25$, and for $Re=218.75$ (a,b), respectively.

To fully characterize the nature of the Hopf bifurcation, we perform a WNL analysis. It provides, in addition, a powerful tool to obtain a model for this bifurcation with the resulting amplitude equation. Subcritical Hopf bifurcation requires to perform the expansion up to the fifth order. At the fifth order, the amplitude equation of a Hopf bifurcation is the Landau equation, which writes as

(5.1)\begin{equation} \cfrac{\partial a}{\partial t} - \epsilon^{2}(\lambda_1+\epsilon^{2}\lambda_2)a - (\mu_1+ \epsilon^{2} \mu_2) a|a|^{2} - \gamma a|a|^{4} =0 ,\end{equation}

with $\lambda$, $\mu$ and $\gamma$ complex coefficients and $a$ the complex amplitude. However, the radius of convergence regarding the bifurcation parameter decreases as we compute higher-order expansions. In both cases investigated, the WNL analysis fails to predict accurate coefficients at the fifth order, see figure 18 in appendix B. We, therefore, present the third-order coefficients which are combined with a numerical fit of the 3-D DNS to obtain the fifth-order coefficients, see table 2. This fit is based on the fluctuating kinematic energy.

Table 2.

Coefficients of the fifth-order Landau equations for $h=0.5$ and $h=0.25$. The third-order coefficients are computed by solving the WNL analysis. The fifth-order coefficients are fitted on the DNS data.

Regarding the coefficients of the Landau equation, the $\mu _1$ coefficients are always positive meaning that the quadratic nonlinearities destabilize the linear eigenmode in addition to the linear growth rate, represented by $\lambda _1$. Only $\gamma$, the fifth-order term, stabilizes the amplitude equation. This is typical of a subcritical bifurcation. Such a bifurcation type is sensitive to external perturbation and the location of the vortex breakdown in experimental cases can be, therefore, very sensitive.

The amplitude equations are displayed in figure 9(a) and the corresponding frequencies are displayed in figure 9(b). The third-order solution of the WNL solution is represented by the dashed blue curves. The fluctuating kinematic energy of the 3-D DNS are displayed by the red dots, and the fifth-order amplitude equation is displayed by the red curves. The solid curves represent the stable branches and the dashed curves are unstable branches. As the Hopf bifurcation is subcritical, the third-order amplitude equation is unstable and the fifth-order amplitude equation stabilizes at an earlier Reynolds number than the critical Reynolds number. The frequency remains almost constant for the Reynolds number investigated.

Figure 9.

Limit cycle amplitude (a,c), and frequency (b,d) predicted by the third- and fifth-order amplitude equations for $S=1.0$, $h=0.5$ (a,b) and $h=0.25$ (c,d), where $a = \rho \ \textrm {e}^{\textrm {i}\phi }$ as defined in (5.1).

The linear threshold characterized by the critical Reynolds number does not correspond to the only possible transition of the flow, the nonlinear threshold defining the hysteresis region appears for smaller Reynolds number and develops possibly first. The hysteresis range spans over ${\rm \Delta} Re = 7$ for $h=0.25$ and ${\rm \Delta} Re = 10$ for $h=0.5$. In both cases considered, the nonlinear thresholds appear for Reynolds larger than the turning point giving rise to the second, stable, branch in axisymmetric bifurcation diagram (see figure 5). Nevertheless, the 3-D solution shows only the spiral mode and no recirculation regions are observed. Due to the nature of the bifurcation, which is very sensitive to disturbance, the eventuality to see the recirculation bubble for such a swirling flow, submitted to an adverse pressure gradient, is unlikely for the considered 3-D configurations.

6. Lateral displacement of the vortex

The vortex dynamics induced by a lateral displacement of the obstacle is summarized in the bifurcation diagram portrayed in figure 10. This bifurcation diagram is produced by displaying the successive minima and the maxima of a time series recorded at the location $(\tilde {R},\tilde {\theta },\tilde {Z}) = (0.5,0.0,18.0)$. Therefore, limit cycle solutions appear as two points, while more complex orbits such as quasi-periodic states or chaotic states display a dispersed distribution of points.

Figure 10.

Bifurcation diagram of the min-max temporal series of the axial velocity of the 3-D DNS at $(\tilde {X},\tilde {Y},\tilde {Z})=(0.5,0.0,18.0)$ as a function of the obstacle lateral displacements $D$ at $S=1.0$, $Re=224.6$ and $h=0.25$.

Figure 10 shows that the vortex dynamics preserves its periodicity below the lateral displacement $D<0.1$. In this limit, which can be considered as small, the obstacle appears still as a flat wall for the vortex and does not disturb the sub-critical unstable mode of the vortex. The displacement required to fully stabilize the unstable spiralling mode is $D=0.2$, where the vortex is seen to manage to escape the obstacle. A single point is displayed in figure 10, meaning that a steady-state solution is retrieved.

Complex orbits are observed in this bifurcation diagram for a lateral displacement between $D=0.1$ and $D=0.16$. These quasi-periodic or chaotic states exhibit the nonlinear interactions induced by the misalignment of the vortex centre and the obstacle. The associated time series and amplitude Fourier spectrum is displayed in figure 11 for the representative case of $D=0.1$. The transient phase, which ends at $t=1600$, is removed from the Fourier analysis. The flow solution evolves toward a more complex orbit, where the main frequency of the spiralling motion is kept, and low-frequency modulations of the signal develop.

Figure 11.

Time series (a) and amplitude Fourier spectrum (b) of the 3-D DNS at $(\tilde {X},\tilde {Y},\tilde {Z})=(0.5,0.0,18.0)$ for $S=1.0$, $Re=224.6$, $h=0.25$, and $D=0.1$.

Visualizations of the flow solutions using the lambda 2 criterion are displayed in figure 12. Figure 12(a) represents the flow solution for a lateral displacement of $D=0.16$. The main spiral is observed, with smaller vortices inside the main recirculation region. In addition, the steady flow solution for lateral displacement of $D=0.2$ is displayed in figure 12(b). A horseshoe vortex develops in front of the obstacle with a stronger branch where most of the vorticity is redirected as the vortex circumvents the obstacle.

Figure 12.

Lambda 2 criterion of the 3-D DNS for $S=1.0$, $Re=224.6$, $h=0.25$ and $D=0.16$ (a), $D=0.2$ (b), and their associated time series at the probe location $(\tilde {X},\tilde {Y},\tilde {Z})=(0.5,0.0,18.0)$.

An amplitude Fourier spectrum is displayed in figure 13 for $D=0.16$. The low-frequency oscillations are axisymmetric pulsations of the flow with residual oscillations of a double helix that can be associated with a flow structure developing in the recirculation zone between the obstacle and the spiralling vortex. It should be noticed that this double helix, which rotates at low frequency, is a mode reported in the modal analysis for the axisymmetric solution in figure 6(d). The main frequency peak is mainly composed of the single helix spiral and partially of the previously mentioned axisymmetric pulsation. Such a type of scenario is typically due to the quadratic nonlinearity of the Navier–Stokes equations. The lateral displacement of the obstacle can be associated with a steady mode of azimuthal wavenumber $m=1$, which interacts with the main single spiral mode generating axisymmetric oscillations. This scenario has already been reported in Francis turbines with the generation of a synchronous pressure (Nishi et al. Reference Nishi, Matsunaga, Kubota and Senoo1984; Pasche et al. Reference Pasche, Gallaire and Avellan2019) and in the spiral vortex breakdown where several single helical modes generate axisymmetric pulsations (Pasche et al. Reference Pasche, Avellan and Gallaire2018a).

Figure 13.

Amplitude Fourier spectrum decomposed in azimuthal wavenumber of the 3-D DNS at $(\tilde {X},\tilde {Y},\tilde {Z}) =(0.24,0,18.75)$ for azimuthal velocity and $S=1.0$, $Re=224.6$, $h=0.25$ and $D=0.16$.