2.1 Governing equations

We consider the dimensionless incompressible Navier–Stokes equations in a semi-infinite domain $\unicode[STIX]{x1D6FA}$ ,

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{V}}{\unicode[STIX]{x2202}t}+(\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{V}=-\unicode[STIX]{x1D735}P+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{V}\quad \text{in}~\unicode[STIX]{x1D6FA}\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}=0\quad \text{in}~\unicode[STIX]{x1D6FA},\end{array}\right\}\end{eqnarray}$$

where $(\boldsymbol{V},P)$ are the state variables and $Re$ the Reynolds number, built on the fluid’s kinematic viscosity and on the vortex core radius and incoming centreline streamwise velocity prevailing at inlet. This system is written in compact form as

(2.2) $$\begin{eqnarray}{\mathcal{N}}\unicode[STIX]{x2202}_{t}\boldsymbol{Q}+{\mathcal{M}}(\boldsymbol{Q})=\mathbf{0}\quad \text{in}~\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $\boldsymbol{Q}=(\boldsymbol{V},P)$ , ${\mathcal{M}}$ is the nonlinear evolution operator and ${\mathcal{N}}$ is the singular operator premultiplying the time derivative. We consider the dynamics governed by (2.1) and resulting from an inlet (Grabowski & Berger Reference Grabowski and Berger1976) vortex profile

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}U_{R}(R)=0,\quad U_{\unicode[STIX]{x1D703}}(R)=S\unicode[STIX]{x1D6F9},\quad U_{Z}(R)=1,\quad \text{on}~\unicode[STIX]{x1D6E4}_{in}\\ \unicode[STIX]{x1D6F9}(R\leqslant 1)=R(2-R^{2}),\quad \unicode[STIX]{x1D6F9}(R>1)=1/R.\end{array}\right\}\end{eqnarray}$$

This vortex is composed of a unitary uniform dimensionless axial velocity component $U_{Z}=1$ , a null radial velocity component and a tangential velocity component with maximum intensity assigned by the swirl number $S$ and composed of a vortex core of dimensionless radius $1$ in solid body rotation and a potential decay outside the vortex core. In the present study, we fix the swirl number to $S=1.095$ but let vary the Reynolds number up to $Re=500$ , thereby restricting our investigation to the spiral vortex breakdown case, in contrast to pure bubble vortex breakdown or multiple helix vortex breakdown observed respectively at lower and larger swirl numbers (Ruith et al. Reference Ruith, Chen, Meiburg and Maxworthy2003; Meliga et al. Reference Meliga, Gallaire and Chomaz2012a ).

2.2 Direct numerical flow simulation

Direct numerical flow simulations have been performed on a 3-D cylindrical domain based on a Cartesian coordinate system $(X,Y,Z)$ . The main axis of the cylinder is oriented along the $Z$ -component and its origin is located at the centre of the inlet, see figure 1. A second, cylindrical, coordinate system is introduced $(R,\unicode[STIX]{x1D703},Z)$ , which is more convenient to study vortex flows, its axial axis $Z$ and origin remaining identical. The transformation matrix between the two reference frames is defined as follows: the $X$ -axis is oriented along $R$ for $\unicode[STIX]{x1D703}=0$ , while the $Y$ -axis completes the reference frame to obtain a right-handed coordinate system and corresponds therefore to $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ . The radial extension of the domain equals $R_{max}=50$ while its length is $Z_{max}=40$ . Such dimensions were necessary to guarantee an unconfined radial boundary condition (see § A.1) because the vortex breakdown resulting from these flow conditions is known to be very sensitive to confinement effects, as experimented by Ruith et al. (Reference Ruith, Chen, Meiburg and Maxworthy2003). The computational domain is bounded by the Grabowsky and Berger vortex (2.3) on the inlet boundary $\unicode[STIX]{x1D6E4}_{in}$ , a free-outflow condition $(-P\boldsymbol{I}+Re^{-1}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{V}))\boldsymbol{\cdot }\boldsymbol{n}=0$ on the external boundary $\unicode[STIX]{x1D6E4}_{ext}$ and a convective condition $\unicode[STIX]{x2202}_{t}\boldsymbol{V}+\boldsymbol{U}_{c}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{n}}\boldsymbol{V}=0$ on the outlet boundary $\unicode[STIX]{x1D6E4}_{out}$ . We have fixed the convective velocity to be equal to the free-stream velocity $\boldsymbol{U}_{c}=\boldsymbol{e}_{Z}$ . The governing equations are solved numerically with NEK5000 a spectral element solver developed by Fischer, Lottes & Kerkemeier (Reference Fischer, Lottes and Kerkemeier2008). The temporal discretization of the nonlinear terms is treated explicitly by a third-order backward-differentiation scheme (BDF3) combined with a third-order extrapolation scheme (EXT3). The linear terms are treated implicitly in time and a pressure–velocity decoupling method is used for the spatial discretization. The velocity and pressure space are represented by a tensor-product array of Gauss–Lobatto–Legendre (GLL) and Gauss–Legendre (GL) points of polynomial orders $N$ and $N-2$ . In the present study, $P_{10}-P_{8}$ polynomial orders for velocity-pressure and 11 040 hexahedral elements are required to compute the flow field. This represents a computation of 8 million of degree of freedoms. The validation of the 3-D DNS is presented in § A.1.

Figure 1.

Three-dimensional schematic of the flow configuration.

2.3 Axisymmetric flow

The steady axisymmetric flow solution developing from the inlet Grabowsky and Berger vortex (2.3) is also investigated, with governing equations written in compact form as

(2.4) $$\begin{eqnarray}\tilde{{\mathcal{M}}}(\boldsymbol{Q}_{B})=\mathbf{0}\quad \text{in}~\unicode[STIX]{x1D6FA}_{a},\end{eqnarray}$$

with $\boldsymbol{Q}_{B}=(\boldsymbol{V}_{B},P_{B})$ the base flow state variables. An axisymmetric domain $\unicode[STIX]{x1D6FA}_{a}$ , with the same dimensions as those of the 3-D domain is used $R_{max}=50$ and $Z_{max}=40$ , except that a sponge region is added along the $Z$ -component at the end of the domain. The length of this sponge region is equal to $L_{sponge}=60$ , using the same smoothing function as in Meliga et al. (Reference Meliga, Gallaire and Chomaz2012a ). The inlet and external boundary conditions are equivalent to those of the 3-D case but a symmetry condition is applied on the axis $\unicode[STIX]{x1D6E4}_{axis}$ that equals to $V_{BR}=V_{B\unicode[STIX]{x1D703}}=0,\unicode[STIX]{x2202}_{R}V_{BZ}=0$ , while the outlet boundary has a free-outflow condition.

The flow solution is carried out by the means of the finite element library Freefem++ (Hecht Reference Hecht2012). The steady axisymmetric incompressible Navier–Stokes equations are first premultiplied by $R$ to avoid the axis singularity and then solved numerically via a Newton–Raphson iterative method. The related linear system is carried out by the unsymmetric multifrontal sparse LU factorization package (UMFPACK). The computational domain is meshed by 408 866 triangular Taylor–Hood elements, $P_{2}-P_{1}$ polynomial-order elements for the velocity–pressure unknowns. A convergence criterion of $10^{-8}$ on the $H_{1}$ -norm is reached in the Newton–Raphson iterative method, which leads to the base flow solution.

2.4 Global stability analysis

The small disturbance dynamics, $\boldsymbol{q}^{\prime }=(\boldsymbol{v}^{\prime },p^{\prime })$ , of the Navier–Stokes equations is investigated by global linear stability analyses around the base flow $\boldsymbol{Q}_{B}$ and the mean flow $\bar{\boldsymbol{Q}}=(\bar{\boldsymbol{V}},\bar{P})$ in the axisymmetric domain $\unicode[STIX]{x1D6FA}_{a}$ . The disturbances $\boldsymbol{q}^{\prime }$ are expanded in normal modes for different azimuthal wavenumbers $m\in \mathbb{Z}$ ,

(2.5) $$\begin{eqnarray}\boldsymbol{q}^{\prime }(R,\unicode[STIX]{x1D703},Z,t)=\boldsymbol{q}(R,Z)\exp (\text{i}(m\unicode[STIX]{x1D703}-\unicode[STIX]{x1D714}t))+\text{c.c.},\end{eqnarray}$$

where c.c. is the complex conjugate. We assume mean flow solutions unperturbed at first order by the Reynolds stresses (Barkley Reference Barkley2006, Meliga et al. Reference Meliga, Pujals and Serre2012b , Mantič-Lugo et al. Reference Mantič-Lugo, Arratia and Gallaire2014 for discussion of this hypothesis). Therefore the eigenvalue problems related to the base and mean flow state become similar and may be written in a compact form using a generic flow solution $\boldsymbol{Q}_{0}=(\boldsymbol{V}_{0},P_{0})$ as

(2.6) $$\begin{eqnarray}[-\text{i}(\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}){\mathcal{N}}+{\mathcal{L}}_{m}(\boldsymbol{V}_{0})]\boldsymbol{q}=\mathbf{0},\quad \text{in}~\unicode[STIX]{x1D6FA}_{a},\end{eqnarray}$$

where ${\mathcal{L}}_{m}$ is the operator for the linearized Navier–Stokes equations of azimuthal wavenumber $m$ and $\boldsymbol{Q}_{0}$ the linearization point that is equal to $\boldsymbol{Q}_{B}$ for the base flow eigenvalue problem and $\bar{\boldsymbol{Q}}$ for the mean flow eigenvalue problem. The time-averaged flow distribution $\bar{\boldsymbol{Q}}$ is computed on the fly in NEK5000 with a minimum of 160 periods of the characteristic oscillation, after the transient phase has died out. In addition, an interpolation step is needed to transfer the mean flow from the 3-D mesh to the axisymmetric mesh. This is carried out by a spectral interpolation of $10^{-13}$ accuracy integrated in NEK5000. At this level the mean flow is assumed purely axisymmetric and only a slice at $X=0$ and $Y>0$ is interpolated and exported to the eigenvalue solver.

In this context we use the definition that a positive azimuthal wavenumber $m>0$ corresponds to a spiral coiling in space in the anti-trigonometric direction and therefore, since $S>0$ , in the opposite direction as that of the ambient flow while $m<0$ designates a spiral coiling in the same direction as the flow. The boundary conditions applied to the disturbances are homogeneous Dirichlet conditions when Dirichlet conditions are imposed on the base flow and remain unchanged in case of a Neumann condition. Specific boundary conditions, which are recalled in table 1, are applied on the symmetry axis due to the azimuthal expansion. They are obtained by symmetry considerations of the perturbations (Khorrami Reference Khorrami1991).

Table 1.

Boundary conditions on the axisymmetric axis applied to the disturbances for different azimuthal wavenumbers.

The global stability analysis is carried out by finite element method implemented in the Freefem++ software. The eigenvalue problem is first premultiplied by $R$ and then the discretized equations are solved by the implicit restarted Arnoldi method of the ARPACK library. The eigenvalues are obtained with a tolerance of $10^{-6}$ of the ARPACK solver. The axisymmetric mesh described in § 2.3 is used to compute the eigenvalue problem of the base and mean flow, except that the sponge region is removed for the latter. A validation of the eigenvalue calculations of the base and mean flows is presented in §§ A.2 and A.3 respectively.

3.1 Three dimensional flow

The 3-D instantaneous flow computations are illustrated for two specific Reynolds values, $Re=180$ and $Re=300$ , by axial vorticity contours in figures 2(a) and 2(b), respectively. These contours highlight the recirculation bubble and the spiralling motions of the flow. In the case of $Re=180$ , the recirculation bubble is accompanied by a single helical vortex structure coiling in its wake, figure 2(a). The spiral spins temporally in the same direction as the Grabowsky and Berger vortex but coils spatially in the opposite direction, which is consistent with the previous results of Ruith et al. (Reference Ruith, Chen, Meiburg and Maxworthy2003). The temporal discrete Fourier amplitude spectrum of the radial velocity at the location $(R,\unicode[STIX]{x1D703},Z)=(0.1,0.0,5.0)$ , see figure 3(a), exhibits a quasi-sinusoidal signal at the frequency $f=0.189$ . The system has reached a limit cycle that arises from a self-sustained instability. In contrast, at $Re=300$ , the wake of the bubble is followed by two helical spirals, which have the same spatial coiling and temporal spinning as for $Re=180$ . The axial vorticity contours in figure 2(b) exhibit a second, slaved, spiral in the vicinity of first, master, spiral. They are differentiated by their axial vorticity magnitudes, the slaved structure exhibiting a modulated axial vorticity iso-contour when compared to the master structure, characterized by a dense iso-contour, at a value of 2.4. The temporal discrete Fourier amplitude spectrum is richer than that prevailing at $Re=180$ , see figure 3(b). This amplitude spectrum exhibits four peaks: a main one with a frequency of $f=0.204$ corresponding to the master helical vortex, a smaller peak close to the main peak with a frequency of $f=0.167$ , a low frequency peak at a frequency of $f=0.003$ as well as the main peak’s second harmonic at $f=0.408$ . The low frequency component and the broad base of the main peak, which latter is not due to a lack of fast Fourier transform (FFT) accuracy, are markers of the nonlinear interactions between the two spirals.

Figure 2.

Iso-surfaces of axial vorticity contour of the 3-D DNS flow solution for (a) $Re=180$ and (b) $Re=300$ , for $S=1.095$ .

Figure 3.

Discrete Fourier amplitude spectrum of the time signal of the 3-D DNS of the radial velocity component located at $(R,\unicode[STIX]{x1D703},Z)=(0.1,0.0,5.0)$ for $Re=180$ and $Re=300$ , for $S=1.095$ .

The axisymmetric time-averaged flow fields, to be investigated by global stability analysis, are displayed in figure 4. These figures are composed by the tangential velocity component on the upper part, $R>0$ , and streamsurfaces coloured by the magnitude of the axial velocity component on the lower part, $R<0$ . The time-averaged streamsurfaces exhibit a recirculation bubble, characteristic of vortex breakdown. At $Re=180$ , the bubble is elongated and lifted up from the centreline in its downstream part while the bubble has a more round shape for $Re=300$ , and additionally a second decelerated region is observed at $Z=5$ . We have observed that the extension of the time-averaged bubble decreases with the Reynolds number while the intensity of the undulating wake in the lee of the bubble increases with the Reynolds number.

3.2 Axisymmetric flow

The axisymmetric base flow, governed by (2.4), is displayed in figures 5(a) and 5(b) for $Re=180$ and $Re=300$ , respectively, for $S=1.095$ . At a Reynolds value of $Re=180$ , the difference between the mean flow (figure 4 a) and the base flow (figure 5 a) is almost invisible. The mean flow correction resulting from the saturation of the self-sustained instability is therefore best observed at $Re=300$ , see figures 5(b) and 4(b), where the mean flow strongly deviates from the base flow, supporting the observations already reported by Ruith et al. (Reference Ruith, Chen, Meiburg and Maxworthy2003).

Figure 4.

Time-averaged solution of the 3-D Navier–Stokes equations, $\bar{\boldsymbol{V}}$ , for $S=1.095$ and $Re=180$ (a) and $Re=300$ (b). The upper part of the graph shows the tangential velocity component and the lower part, the streamsurfaces coloured by the magnitude of the axial velocity component.

Figure 5.

Solution of the steady axisymmetric Navier–Stokes equations, $\boldsymbol{V}_{B}$ , for $S=1.095$ and $Re=180$ (a) and $Re=300$ (b). The upper part of the graph shows the tangential velocity component and the lower part, the streamsurfaces coloured by the magnitude of the axial velocity component.

3.3 Global stability analysis

The global stability of the base flow and the mean flow have been investigated, with a focus on the first helical mode $m=1$ . A representative example is presented in figure 6 at a Reynolds value of $Re=300$ . Up to three unstable eigenmodes, labelled $B_{n}$ for $n=1,2,3$ can be identified in the eigenvalue spectrum of the base flow, while only two unstable modes $M_{1}$ and $M_{2}$ are detected in the mean flow eigenvalue spectrum. The associated growth rates are significantly larger when the linearization is performed around the base flow than around the mean flow, in agreement with mean flow modification based saturation mechanisms reported in different flow geometries (Barkley Reference Barkley2006; Mantič-Lugo et al. Reference Mantič-Lugo, Arratia and Gallaire2014).

Figure 6.

Eigenvalue spectra of the global linear stability of mode $m=1$ about the base (a) and mean (b) flow at $Re=300$ and $S=1.095$ . Unstable eigenvalues are labelled respectively $B_{1}$ , $B_{2}$ , $B_{3}$ and $M_{1}$ and $M_{2}$ .

We have then systematically compared the frequency prediction of the base and mean flow stability analyses for different Reynolds values $Re=100{-}500$ by keeping constant the swirl number $S=1.095$ , focusing only on unstable modes and on the azimuthal wavenumber $m=1$ , since other azimuthal wavenumbers were all found to be linearly damped. The instability threshold was found at $Re=143.5$ for the base flow and approximately at $Re=143.7$ for the mean flow, where the dying out of the transient in the DNS become prohibitively long close to the threshold. In the vicinity of the threshold, the linear frequency predictions (respectively in red and blue for the base and mean flow) are in agreement with the limit cycle frequency extracted from discrete Fourier transform of the 3-D DNS signal (in green), see figure 7(a). As the Reynolds number increases, the growth rate of the dominant unstable mode linearized about the base flow ( $B_{1}$ ) increases rapidly, followed by a second unstable mode ( $B_{2}$ ) bifurcating at $Re=200$ and subsequently emerging unstable eigenvalues $B_{3}$ and $B_{4}$ . At $Re=350$ , where the base flow stability analysis has been stopped, the base flow is found to be unstable to no less than four unstable modes, as seen in figure 7(b). The frequencies of the four eigenvalues considered are all around $\unicode[STIX]{x1D714}_{r}\cong 1.08$ at their respective instability thresholds, beyond which they increase linearly with the Reynolds number. The base flow frequency prediction of the dominant unstable mode remains acceptable before it crosses the non-monotonic limit cycle frequency evolution (in green) around $Re=235$ . This appears as an a posteriori justification of the validity range of the linear stability analysis about the base flow of Meliga et al. (Reference Meliga, Gallaire and Chomaz2012a ) and Qadri et al. (Reference Qadri, Mistry and Juniper2013).

Figure 7.

Growth rate and angular frequency of the global linear stability of the mean, $M_{1}$ , $M_{2}$ and base flow, $B_{1}$ , $B_{2}$ , $B_{3}$ , $B_{4}$ , for azimuthal wavenumber $m=1$ compared to the 3-D angular frequency from the DNS, computed by discrete Fourier transform at location $(R,\unicode[STIX]{x1D703},Z)=(0.1,0.0,5.0)$ .

The limit cycle frequency (in green) indeed first strongly increases away from the threshold until it drops for $Re>230$ and finally slowly increases with the Reynolds number for $Re>240$ . The frequency of the leading eigenvalue ( $M_{1}$ ) of the mean flow (in blue) closely follows the limit cycle frequency (in green) along the whole Reynolds number range, demonstrating a correct frequency prediction by the global stability analysis around the mean flow, as found in the canonical cylinder wake flow (Barkley Reference Barkley2006). Similarly the maximum growth rate of the mean flow remains much smaller that the growth rate of the dominant mode of the base flow that stays close to zero. Above $Re>300$ , a second unstable eigenvalue is found with a frequency that corresponds to the second peak observed in the temporal discrete Fourier amplitude spectrum (the angular frequencies $\unicode[STIX]{x1D714}_{r}$ should be divided by $2\unicode[STIX]{x03C0}$ to compare to the frequencies obtained in figure 3, for example at $Re=300$ , $f_{M_{2}}=0.156$ , and $f_{DNS}=0.167$ ). The most unstable eigenmode is likely to correspond to the master helix. The second unstable eigenmode would be associated with the slaved helix. These frequency predictions, together with the global stability of the base and mean flows with respect to $m=2$ (not shown for brevity) are a confirmation that the snapshot of figure 2(b) at $Re=300$ does not correspond to an $m=2$ instability but rather to the superimposition of two $m=1$ spirals, each corresponding to a separate frequency peak and without any phase correlation.

These results on the frequency prediction of the stability analysis about the mean flow support the recent extension of this technique to high Reynolds number industrial flows, such as the frequency prediction of the precessing vortex core in a combustion chamber (Oberleithner et al. Reference Oberleithner, Sthr, Im, Arndt and Steinberg2015; Paredes et al. Reference Paredes, Terhaar, Oberleithner, Theofilis and Passchereit2015; Tammisola & Juniper Reference Tammisola and Juniper2016) or of the part load vortex rope in Francis hydraulic turbines (Pasche et al. Reference Pasche, Avellan and Gallaire2017). Moreover, we identify up to two unstable eigenmodes, which may offer a relevant test case for the need to take into account mode interaction in generalizing the self-consistent model for saturation dynamics developed by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014), in the spirit of the generalized quasi-linear models developed by Marston, Chini & Tobias (Reference Marston, Chini and Tobias2016).

4.1 Optimal control theory

Although it can be beneficial to improve mixing, the self-sustained instability in the spiral vortex breakdown is a source of perturbations that is undesirable in several applications (Rheingans Reference Rheingans1940; Syred Reference Syred2006; Gursul et al. Reference Gursul, Wang and Vardaki2007; Nishi & Liu Reference Nishi and Liu2013). The finite amplitude nonlinear predictive control of such perturbations is investigated for two specific Reynolds values, $Re=180$ and $Re=300$ , with the aim of stabilizing the periodic fluctuations by applying a constant volume force. In the present study, we have chosen to minimize two different objectives to achieve this goal, a ‘flow stabilization control’ which targets a zero value of the growth rate of the most unstable eigenvalue declined into two versions controlling the base and mean flow eigenvalue, and a ‘flow manipulation control’ which targets a minimal radial velocity intensity. The rationale behind this latter choice is to minimize the radial flow associated with the bubble formation and thereby to maintain the flow as columnar as possible. This strategy is based on the physical mechanism underpinning spiral vortex breakdown: the wake induced by the bubble triggers an absolute instability region that leads to the spiral vortex oscillation (Gallaire et al. Reference Gallaire, Ruith, Meiburg, Chomaz and Huerre2006; Qadri et al. Reference Qadri, Mistry and Juniper2013). The mathematical framework required in these minimization problems is the optimal control theory. The objective functions for the two control strategies, to be defined in (4.1) and (4.2), are derived for the axisymmetric domain where the base flows and the relevant base and mean flow stability analyses were computed, bypassing larger and more expensive 3-D computations. This dimensionality reduction from the 3-D domain to a 2-D axisymmetric domain is performed thanks to the azimuthal Fourier decomposition naturally appearing in the eigenvalue analysis involved in the ‘flow stabilization control’. In the context of the ‘flow manipulation control’, only the axisymmetric base flow is considered, which is two-dimensional by definition. To be more specific the following objective functions are defined,

Flow manipulation control:

(4.1) $$\begin{eqnarray}\displaystyle \text{Min}\,{\mathcal{J}}_{Manip}(\boldsymbol{F})=\frac{1}{2}\int \!\!\!\int _{\unicode[STIX]{x1D6FA}_{a}}V_{R}(\boldsymbol{F})^{2}\,\text{d}\unicode[STIX]{x1D6FA}_{a}+\frac{\unicode[STIX]{x1D6FC}}{2}\int \!\!\!\int _{\unicode[STIX]{x1D6FA}_{a}}\Vert \boldsymbol{F}\Vert ^{2}\,\text{d}\unicode[STIX]{x1D6FA}_{a}, & & \displaystyle\end{eqnarray}$$

Flow stabilization control:

(4.2) $$\begin{eqnarray}\displaystyle \text{Min}\,{\mathcal{J}}_{Stab}(\boldsymbol{F})=|\text{Im}(\unicode[STIX]{x1D714})|+\frac{\unicode[STIX]{x1D6FC}}{2}\int \!\!\!\int _{\unicode[STIX]{x1D6FA}_{a}}\Vert \boldsymbol{F}\Vert ^{2}\,\text{d}\unicode[STIX]{x1D6FA}_{a}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{J}}$ is the objective function, $\boldsymbol{F}$ the control variable (here a volume force), $(\boldsymbol{V},P)$ are the state variables, $\unicode[STIX]{x1D714}$ is the eigenvalue of the most unstable eigenmode linearized around the base $\unicode[STIX]{x1D714}_{B_{1}}$ or mean flow $\unicode[STIX]{x1D714}_{M_{1}}$ and $\unicode[STIX]{x1D6FC}$ is the scalar weight penalizing of the control cost (here the $L^{2}$ -norm of the volume force). We control the flow on the entire axisymmetric domain $\unicode[STIX]{x1D6FA}_{a}$ through the radial and axial volume force components only. We made this choice to avoid trivial control solutions, which would simply decrease the tangential velocity magnitude, and to consider simple control distributions that are easier to visually translate in terms of potential appendage control devices. At convergence, the distributed control investigated will end to an optimal volume force $(F_{R}^{\ast },F_{Z}^{\ast })$ and optimal state variables $(\boldsymbol{V}^{\ast },P^{\ast })$ .

The gradient direction of these objective functions allows us to look for the minimum of ${\mathcal{J}}$ through an iterative algorithm such as a conjugate-gradient method. This gradient is computed by the use of the Frechet derivative applied to a specific Lagrangian functional related to each optimal control problems. These Lagrangian functionals, which are expressed in § B.1, use the adjoint variables as Lagrange multipliers, $(\boldsymbol{V}^{\dagger },P^{\dagger })$ for the Navier–Stokes equations and $(\boldsymbol{v}^{\dagger },p^{\dagger })$ for the eigenvalue problem. The gradient direction is defined by (4.3) for the distributed control

(4.3) $$\begin{eqnarray}\unicode[STIX]{x1D735}{\mathcal{J}}=\boldsymbol{V}^{\dagger }+\unicode[STIX]{x1D6FC}\boldsymbol{F}\quad \text{in}~\unicode[STIX]{x1D6FA}_{a}.\end{eqnarray}$$

We should mention that the computation of the ‘flow stabilization control’, where the leading eigenvalue directly intervenes, uses the sensitivity to a steady force $\unicode[STIX]{x1D735}_{F}\unicode[STIX]{x1D714}_{i}$ as the gradient direction, see Giannetti & Luchini (Reference Giannetti and Luchini2007) and Marquet et al. (Reference Marquet, Sipp and Jacquin2008). Since however the steady force is no more expected to remain infinitesimal, the nonlinearity of the controlled Navier–Stokes equations is included in our optimization strategy, which can therefore be seen as a ‘nonlinear sensitivity’ to a steady force. Since in the case of the canonical cylinder wake flow (Strykowski & Sreenivasan Reference Strykowski and Sreenivasan1990; Marquet et al. Reference Marquet, Sipp and Jacquin2008), the linear sensitivity to a steady force, mimicking a small control cylinder, was sufficient to stabilize the flow close to the instability threshold, we have also investigated the efficiency of this linear sensitivity approach in the case of the spiral vortex breakdown. In this ‘linear sensitivity-based finite amplitude control’, the first gradient direction $\unicode[STIX]{x1D735}_{F}\unicode[STIX]{x1D714}_{i}=\boldsymbol{F}^{(o)}$ is multiplied by a finite amplitude $A$ determined by a bisection method until the growth rate of the controlled leading eigenvalue is equal to zero. Next, we compute the ‘base flow manipulation control’, the ‘base flow stabilization control’, the ‘mean flow stabilization control’ and the ‘linear sensitivity-based finite amplitude control’ of the base and mean flow at a Reynolds value of $Re=180$ , before turning to a Reynolds value of $Re=300$ , further away from instability threshold, where only the ‘mean flow stabilization method’ and the ‘linear sensitivity-based finite amplitude control’ of the mean flow are determined.

4.2 Minimization and numerical methods

The optimal control problems are solved by a conjugate-gradient method where the direction is updated by a Polak & Ribiere (Reference Polak and Ribiere1969) method and the step length is defined by an Armijo line search method, see Polak (Reference Polak1997). The initial step $\unicode[STIX]{x1D70F}$ is dependent on the flow case and has been chosen to ensure convergence (see § B.2) and the control weight is set to $\unicode[STIX]{x1D6FC}=10^{-5}$ . Since the ‘base flow stabilization’ is a non-convex problem, the solution invariance was checked with respect to the parameter $\unicode[STIX]{x1D6FC}$ , the initial step lengths and the optimization methods. The results are presented in § B.2, where the solution has been found to be quasi-invariant under these choices in the methods.

The optimization algorithms for the ‘flow manipulation control’ and the ‘base flow stabilization control’ have been solved numerically with Freefem++ software in the same spirit as the base flow and eigenvalue problem computations. The same mesh and discretization technique are applied. In the case of the ‘mean flow stabilization control’, the time-averaged flow field is computed by DNS in NEK5000 before being transferred on the axisymmetric mesh. The other minimization steps of the conjugate-gradient algorithm are computed within Freefem++ until the next updated control (updated volume force) is obtained and returned to the 3-D DNS. The optimization process is stopped when the objective function ${\mathcal{J}}$ reaches a minimum, see § B.2.

5.1 Base flow manipulation control at $Re=180$

The predictive control of the axisymmetric bubble base flow minimizing the quadratic objective function equation (4.1) at $Re=180$ is displayed through the associated volume forcing in figure 8(a) and the flow solution $\boldsymbol{V}^{\ast }$ in figure 8(b). The radial velocity component exhibited by the streamsurfaces is clearly reduced by the minimization algorithm, starting with a maximum value of $V_{R}^{max}=0.3$ and ending with a value of $V_{R}^{max}=0.037$ at convergence, leading to a columnar solution contrasting with the uncontrolled bubble breakdown, see figure 5(a). Although the control target is not the recirculation bubble itself (in contrast to what was proposed in Boujo & Gallaire (Reference Boujo and Gallaire2014) for plane flows) or the recirculation intensity, we note that the use of our physical mechanism inspired surrogate objective entirely suppresses the recirculation: the minimum centreline axial velocity evolves from an uncontrolled value of $V_{Z}^{min}=-0.08$ to a controlled value of $V_{Z}^{min}=0.8$ . The volume force needed to restore a breakdown-free columnar solution is a predominating axial suction upstream of the nominal breakdown bubble combined with a subdominant radial correction, see figure 8(a), both resulting in a control amplitude at convergence of $\Vert \boldsymbol{F}^{\ast }||=0.084$ , see figure 18 in § B.2.

Figure 8.

Navier–Stokes solutions for the base flow manipulation control at $Re=180$ (b) and the associated volume force (a) at the minimum of the functional ${\mathcal{J}}_{Manip}$ . The upper part of the graph (b) shows the tangential velocity component and the lower part, the streamsurfaces coloured by the magnitude of the axial velocity component.

We have verified that the suppression of the recirculation bubble results in the stabilization of the spiral mode, as anticipated by the understanding of the underlying physical mechanism. Figure 9(a) (blue solid curve) demonstrates that the stabilization is achieved for a control amplitude of $\Vert \boldsymbol{F}\Vert =0.028$ , far beyond the minimizing control amplitude $\Vert \boldsymbol{F}^{\ast }\Vert =0.084$ . Note that this stabilization is accompanied by a 10 % frequency shift (see the blue solid curve in figure 9 b). While physically well motivated the ‘flow manipulation control’ will appear as suboptimal when compared to the ‘flow stabilization control’ presented in § 5.2. As we will see, the comparison of figures 8(a) and 10(a) will show that a lot of control authority is spent in reducing the radial velocity in regions of low receptivity of the instability. We therefore now turn our efforts to base flow stabilization control, trading the conservative but quadratic cost functional (4.1) for a more precise but non-convex cost functional (4.2).

5.2 Base flow stabilization control at $Re=180$

Despite the absence of mathematical proof, we have observed that the minimization algorithm targeting the eigenvalue growth rate of the spiral vortex breakdown mode developing on the base flow, herein referred to as ‘base flow stabilization control’, succeeds in stabilizing the flow. As seen by the red solid curve in figure 9(a), the normalized eigenvalue growth rate progressively decreases until stabilization for $\Vert \boldsymbol{F}^{\ast }\Vert =0.0058$ , i.e. 5 times less than the ‘flow manipulation control’ approach. Interestingly, a finite amplitude recirculation bubble is still present in the stable, optimally controlled flow reported in figure 10(b). The comparison with the uncontrolled base flow (see figure 5 a) shows that minimal changes in the recirculation bubble are sufficient to restore the stability of the flow and this stabilization goes together with a small frequency shift, reported in figure 9(b).

Figure 9.

Base flow dominant growth rate normalized by the most unstable uncontrolled eigenvalue at $Re=180$ as a function of the norm of the control (a) and of the angular frequency of the mode (b) for the ‘flow manipulation control’ (blue solid curve) and ‘base flow stabilization control’ (red solid curve), and the ‘linear sensitivity-based finite amplitude control’ (red dashed curve).

Figure 10.

Optimal ‘base flow stabilization control’ force (a) and associated flow solution (b) at $Re=180$ and $S=1.095$ . The upper part of the graph (b) shows the tangential velocity component and the lower part, the streamsurfaces coloured by the magnitude of the axial velocity component.

The determination of this optimal, nonlinear, predictive control requires a computationally expensive iterative minimization procedure. A linear sensitivity-based finite amplitude control therefore provides a promising alternative. The latter is included in the red dashed curve in figure 9(a), showing that, although linearity is preserved, flow stabilization is achieved for approximately two times more control amplitude. A comparison of both volume force distributions in figures 10(a) and 11(a) shows in addition that the first gradient direction may deviate quite substantially from the nonlinear optimal control. While the latter suggests to decelerate the flow in the inlet region, thereby pulling on the recirculation bubble, the sensitivity-based centreline axial force accelerates the flow upstream of the bubble, thereby pushing the bubble downstream. This results in slightly different flow topologies (see figures 10 b and 11 b), noting that in both cases, the bubble is slightly smaller than in the nominal, unstable, base flow reported in figure 5(a). In contrast, the volume forces acting in bubble wake around $Z=5$ are very similar and act to accelerate the flow in the wake region. This is also clearly visible from the resulting controlled flows (figures 10 b and 11 b) when compared to the base flow (figure 5 a). All these results illustrate that, even for moderate growth rates, the finite amplitude linear sensitivity to a steady force can remain significantly suboptimal. We now turn our attention to the mean flow stabilization control, first at $Re=180$ .

Figure 11.

Suboptimal ‘linear sensitivity-based finite amplitude control’ (a) and associated flow solution (b) at $Re=180$ and $S=1.095$ . The upper part of the graph (b) shows the tangential velocity component and the lower part, the streamsurfaces coloured by the magnitude of the axial velocity component.

Figure 12.

Mean flow dominant growth rate, normalized by the most unstable uncontrolled mean flow eigenvalue at $Re=180$ (red curves) and $Re=300$ (blue curves) as a function of the norm of the control (a) and of the angular frequency of the mode (b) for the nonlinear mean flow stabilization (solid curves) and the ‘linear sensitivity-based finite amplitude control’ (dashed curves).

5.3 Mean flow stabilization at $Re=180$

The control magnitude and eigenvalue path during the minimization process are exhibited in solid red lines in figure 12(a,b) for the ‘mean flow stabilization control’ at $Re=180$ , together with the ‘linear sensitivity-based finite amplitude control’ of the mean flow (red dashed lines). The optimal mean flow stabilization control is reached for a norm of $\Vert \boldsymbol{F}^{\ast }\Vert =0.0057$ while the linear sensitivity-based finite amplitude control stabilizes the flow for $\Vert \boldsymbol{F}^{\ast }\Vert =0.014$ . These values compare very well with those obtained for the base flow stabilization control. The volume force and controlled solution, see figure 13, are also very similar to those of the base flow control, see figure 10. This indicates the consistency of both approaches targeting the dominant growth rate of the mean and base flows respectively, for parameters close to the bifurcation threshold, even if the least unstable eigenvalues of the base and mean flow differ by a factor 10. The 3-D vorticity contour of the controlled flow solution is plotted in figure 16(a), demonstrating that the axisymmetry of the flow is restored and the inspection of the time-dependent traces remains stationary.

Figure 13.

Time-averaged 3-D Navier–Stokes solution of the mean flow stabilization control (b) and the related force applied (a) for $Re=180$ at convergence. The upper part of the graph (b) shows the tangential velocity component and the lower part, the streamsurfaces coloured by the magnitude of the axial velocity component.

The linear sensitivity-based finite amplitude control about the mean flow is presented in figure 14(a) and represents the initial sensitivity to a steady force for $Re=180$ . The receptivity distribution, while suboptimal, is very similar to the nonlinear sensitivity map, and enables a correct identification of the maximum and minimum locations.

Figure 14.

Suboptimal ‘linear sensitivity-based finite amplitude control’ about the mean flow at $Re=180$ (a) and $Re=300$ (b) for a swirl number of $S=1.095$ .

5.4 Mean flow stabilization at $Re=300$

At $Re=300$ , the minimization path is less straight, and the optimal solution (solid blue line in figure 12) seems to oscillate around the trajectory of the ‘linear sensitivity-based finite amplitude control’ (dashed blue line in figure 12) before eventually sensibly reducing the eigenvalue growth rate at the end of the algorithm. The control magnitude is $\Vert \boldsymbol{F}^{\ast }\Vert =0.019$ at convergence, to be compared with the value of $\Vert \boldsymbol{F}^{\ast }\Vert =0.027$ required for the mean flow ‘linear sensitivity-based finite amplitude control’. The optimal volume force as well as the resulting controlled flow solution are exhibited in figure 15. In contrast to the Reynolds number $Re=180$ case, most of the control effort is located in the recirculating bubble region while little forcing acts in the wake region. The upstream forcing blows the bubble downstream, as opposed to the optimal forcing at $Re=180$ , where the bubble was sucked. The controlled flow solution has no recirculation, but adopts instead a wavy form and demonstrates therefore that a large modification of the mean flow is required to restore its stability, while a small modification is required at a Reynolds value of $Re=180$ . The 3-D solution is displayed in figure 16(b), demonstrating the axisymmetrization of the controlled flow. Although several waves are observed on the vorticity contour, we have checked that the controlled flow remains steady. In addition, the ‘linear sensitivity-based finite amplitude control’ about the mean flow is presented in figure 14(b) and represents the initial sensitivity to a steady force for $Re=300$ . Similarly to the ‘mean flow stabilization control’ at a Reynolds value of $Re=180$ , the linear receptivity well identifies the maximum and minimum of the nonlinear receptivity map (see figure 15 a), although it remains suboptimal (see figure 12).

Figure 15.

Time-averaged 3-D Navier–Stokes solution of the mean flow stabilization control (b) and the related force applied (a) for $Re=300$ at convergence. The upper part of the graph (b) shows the tangential velocity component and the lower part, the streamsurfaces coloured by the magnitude of the axial velocity component.

Figure 16.

3-D Navier–Stokes vorticity of the optimal mean flow stabilization control solution for $Re=180$ (a) and $Re=300$ (b). These 3-D flows are steady solutions and they correspond to the mean flows (and therefore base) reported in figures 13 and 15.