Optimal flow control of vortex breakdown in a laminar swirling flow

In highly swirling flows, such as hydraulic turbines operating under part-load (PL) conditions, vortex breakdown occurs and performance is impaired. Consequently, it is imperative that mitigation measures are taken. In the present study, a laminar swirling flow with a vortex breakdown at a Reynolds number of 180 is investigated. At the inlet, a swirling velocity profile with a swirl number of 1.095 is set. A stability analysis is conducted to identify unstable modes based on the assumption that vortex breakdown is a global instability. The results indicate that spiral modes with wave number 1 are unstable. An optimal flow control method based on the Adjoint method is then utilized to mitigate vortex breakdown. In the present study, the control method targets vorticity using a minimization algorithm. Control variables include radial and axial body forces. According to the results, the method was effective in mitigating vortex breakdown. A stability analysis conducted during the control process revealed that as the vorticity decreased, the growth-rate of the eigenvalue decreased, indicating that the flow is stabilized.


Introduction
In engineering applications, swirling flows can be quantified by the Swirl number, which is defined as the ratio of axial flux of tangential momentum to axial flux of axial momentum [1].When the swirl number exceeds a critical value, there is a sudden change in the flow vortex core that leads to vortex breakdown.A vortex breakdown is defined by Leibovich [2] as the formation of a stagnation point on the axis of the vortex, followed by reversed flow in a region of limited axial extent.Sarpkaya [3] observed three forms of vortex breakdown: bubbles (axisymmetric vortex breakdown), spirals, and double helixes.Many industrial applications, such as delta-wings and hydraulic machines, are also subject to vortex breakdown, resulting in performance degradation, pressure fluctuations, and structural Numerous studies have been conducted on the formation of vortex breakdowns.A vortex breakdown was attributed by Squire [4] to the upstream propagation of standing waves and their accumulation upstream.Benjamin [5] considered the vortex breakdown phenomenon as a transition from supercritical to subcritical flow states, akin to the hydraulic jump in an open channel flow.Vortex rope formation is dependent on flow features such as stagnant points, flow separation, and recirculation zones in the draft tube of hydraulic turbines [6].According to Foroutan and Yavuzkurt, the rotating vortex rope is caused by Kelvin-Helmholtz instabilities.Susan-Resiga et al. [7] investigated the vortex rope inside hydraulic turbines.They demonstrated that at the critical swirl number, there would be a transition from a supercritical flow to a subcritical one.
Hydrodynamic instability is based on the theory that real fluid particles entering a flow domain are continuously subjected to infinitesimal random disturbances that can change the flow characteristics by the mechanism of exponential amplification.In order to examine the stability of the flow, smallamplitude perturbations are imposed to it, and their development is examined by determining whether they decay or grow over time.The behavior of the flow after introducing a disturbance determines whether it is convective or absolute instability.The downstream propagation of a disturbance indicates that the flow is convectively unstable, whereas the upstream and downstream propagation of the disturbance indicate absolute instability [8].Global stability analysis resembles absolute instability [8].Using a linear global stability analysis for an axisymmetric swirling base flow at Re = 200 and swirl number 1.095, Meliga et al. [9] determined that the flow is unstable to the spiral mode.Moreover, Pasche et al. [10] conducted a global stability analysis on the mean turbulent swirling flow inside the draft tube of a hydraulic turbine operating at part-load condition.An unstable mode was identified with a frequency close to the rotational frequency of the vortex rope and lower than the runner frequency.As a result, the instability was associated with the vortex rope.A linear global stability analysis was also performed by Seifi et al. [11] for the time-averaged turbulent flow field inside the straight draft tube of a swirl generator [12] in which a vortex rope is present for certain operating conditions.To investigate the formation of vortex rope in a straight draft tube, Seifi et al. [13] performed transient simulations from BEP to part-load conditions of an axial swirl generator.In addition, the researchers conducted modal stability analysis on the flow for different flow rates and determined the critical flow rate at which the vortex rope begins to form.It was found that the flow inside the draft tube is also sensitive to asymmetrical disturbances with a frequency close to the rotational component.Vortex breakdown instability can be mitigated with a variety of methods.
As part of a control strategy, a physical control target must be selected that will have a beneficial effect on mitigation.Krause [14] demonstrated that vortex breakdown will eventually occur if the vortex core grows in streamwise direction, i.e., if the radial velocity is always positive.In order to mitigate vortex breakdowns in a laminar swirling flow, Pasche et al. [15] employed both the concepts of minimizing the radial velocity and stabilizing hydrodynamic instabilities that contribute to vortex breakdowns.Despite the fact that both methods mitigated vortex breakdown, considering stabilization as a control object resulted in a smaller magnitude of body force.Pasche et al. [16] employed the same stabilization method to stabilize the unstable mode associated with the vortex rope in a hydraulic turbine operating at PL conditions.This method stabilized and mitigated the vortex rope, however, other unstable modes emerged during the control process which may cause additional challenges.In the present study, vorticity minimization is considered as physical control target.
It has been demonstrated in numerous studies that axial vorticity inside axisymmetric vortex breakdowns, whether in turbulent or laminar flows, decreases [17][18][19].Spall and Gatski [18] discovered a correlation between vortex breakdown and vorticity and concluded that the axial component of vorticity was redistributed to the radial and azimuthal components as the stagnation point approached.Furthermore, they observed a local intensification of total vorticity in the bubble simulation According to the studies, the axial vorticity after vortex breakdown increases considerably again in laminar flows, increasing the possibility of subsequent vortex breakdowns.Gartshore [20] believed that the vortex breakdown was caused by the diffusion of vorticity associated with the axial flow reversal rather than to the propagation of disturbances.He added that the flow downstream of this reversal might be unstable and proposed that the observed spiral flow is evidence of this instability.
In the present study the laminar swirling flow with a vortex breakdown is considered.The test case is similar to the test case of Pasche et al. [15] .Optimal flow control based on the adjoint method is used to mitigate the vortex breakdown.The control method is based on a minimization algorithm targeting the vorticity.The radial and axial force fields are considered as the control variables.The remainder of this paper is structured as follow.The computational domain is explained in Section 2. The governing equations of both stability analysis and optimal control are given in Section 3. Results are presented and discussed in Section 4. Finally, the main conclusions of the present study are summarized.

Computational test case
The axisymmetric computational domain is depicted in Figure 1.The non-dimensional radius of the domain is chosen to be 5, smaller than 50 used by Pasche et al. [15] .The computational domain consists of two regions.Due to the axisymmetric swirling flow, a sponge region with a slightly decreased Reynolds number based on Eq.1 [9] is required to prevent backflow.The non-dimensional length of the main domain is 40 and 30 for the sponge region.The mesh consists of triangular elements finer near the inlet and axis of the main domain and coarser in the sponge region, see Figure 2. The grid sensitivity analysis was conducted at Re=180.To this end, three different meshes of 50 000, 100 000 and 150 000 triangular elements have been considered.Figure 3 shows the axial velocity profile for each grid at 2 Z = where the bubble vortex breakdown appears.The velocity profile of the coarse mesh deviates significantly from the other two meshes.An increase from the medium sized mesh with 100 000 elements to the finest mesh with 150 000 elements has no tangible impact on the axial velocity profile.In the present work, the medium sized mesh has been considered for further simulations.

Base flow
Base flow is governed by the incompressible momentum and continuity equations in dimensionless form as follow: ( ) .0 where ( ) and p are the nondimensional velocity vector and pressure scalar, respectively.Re is the Reynolds number, based on the fluid's kinematic viscosity, the vortex core radius and incoming centerline streamwise velocity prevailing at the inlet.Since the flow is steady and axisymmetric, the derivatives with respect to time and azimuthal coordinate are set to zero.Base flow is solved using the finite element software Freefem++ [21] with Newton-Raphson method.A convergence criterion of 10-10 on the H1-norm is reached for the Newton-Raphson iterative method.Finite element space P2-P1 is applied for velocity and pressure, respectively.At the inlet, Grabowski & Berger [22] vortex profile based on Eq.3 is defined.This vortex consists of a unitary uniform dimensionless axial velocity component, a zero radial velocity component, and a tangential velocity component whose intensity is determined by the swirl number S. It is composed of a vortex core of dimensionless radius, 1 R = , in solid body rotation and a potential decay outside the vortex core.The other boundary conditions are shown in Table 1.The eigenvalue equation (Eq.6) is then obtained using the following assumptions.
• Flow is axisymmetric so that the circumferential derivatives are set to zero.
• The periodic terms are replaced with the Fourier transform, Eq.5, where m is an azimuthal wavenumber defining the spatial shape of instability in  direction and  is the complex frequency.
 ( )  ( ) ( ) ˆ, , , , , , exp .p z r t p z r im t c c .0 where the operator m  representing derivatives in the  direction is a multiple of the wavenumber (m).The real part of  represents the growth rate and determines the stability of the flow.A positive value shows that the flow is globally unstable, otherwise all perturbations decay over a long time.In addition, the imaginary part represents the frequency of the mode.Circumferential wavenumbers represent either an axisymmetric mode ( 0 m = ) or non-axisymmetric modes ( 0 m  ) which are categorized as single spiral ( 1 m = ) and double helix ( 2 m = ).Equations ( 4) are solved using the finite element method implemented in the software Freefem++ to obtain the eigenvectors and eigenvalues.The resulting eigenvalue equation is solved by ARPACK library and UMFPACK package for matrix inversion.Finite element space P2-P1 is applied for velocity and pressure, respectively.To avoid singularity at the axis of rotation the boundary conditions presented in Table 2 [11] are applied.

Wave number
Boundary condition

Optimal control methodology
The finite amplitude nonlinear predictive control of vortex breakdown is investigated for Re=180 by applying a constant volume force.In the present study, vorticity magnitude is chosen as the objective to be minimized.The objective function for the control strategy is defined as follow:  uF (7) where ( , ) is a control vector and  is the scalar weight penalizing the control cost.u is the norm of vorticity and Ω is the computational domain.To achieve optimality in the present study, the Lagrange multiplier method is utilized [23,24].Lagrange multiplier method is used to change a constrained optimization to an unconstrained one as follow: where NS is the abbreviation for Navier-Stokes equations. .To obtain the extremum, the gradient of the Lagrange multiplier with respect to the control variables, state and adjoint variables must be zero.Therefore, the optimality system is obtained as follow: ( ) .. Adjoint equation .0 The optimality system (Eqs.9-11) is solved iteratively, starting with (0,0) = F at the first iteration.At the end of each iteration the force is updated using the steepest descent method for the next iteration until it is converged.The step length is defined by Armijo line search method [25] and its initial value has been chosen 10 -4 .Moreover,

Base flow simulation
In Figure 4 the computed circumferential velocity contours and streamline for Re=180 are compared with those reported by Pasche et al. for the same flow condition.Both contours and streamlines are in good agreement.The axisymmetric vortex breakdown and the recirculation region inside it are observable.Figure 4 shows a bubble vortex breakdown and its associated recirculation region inside.A stagnation point which is the indication of the occurrence of a vortex breakdown is located upstream at 1 z = .At the tail of the bubble vortex breakdown, 3-D simulations of Pasche et al. [15] demonstrated the presence of a spiral vortex breakdown.

Global stability analysis
The global stability of the axisymmetric base flow has been investigated for 0,1 m = and 2. The perturbations can be equivalently represented by either ( )  − in which * stands for the complex conjugate.Therefore, it is possible to only consider 0 m  and skip 0 m  with no loss of generality [9].The results showed that the flow is only unstable to waves with 1 m = and stable to other wavenumbers similar to the results of Pasche et al. [15] .The unstable eigenvalue of the present study, 0.055 1.17i

 =+
, is in good agreement with the one obtained by Pasche et al. [15] which is 0.0549 1.173i  =+ .Pasche et al. showed that Re=143.7 is the critical Reynolds number since for the higher values, the eigenfrequency is far from the frequency extracted from the 3-D DNS simulation which equals 1.19 [15] .However, once again at Re=235, both eigenfrequency and numerical frequency have the same value.Therefore, they explained that the base flow frequency prediction of the dominant unstable mode remains acceptable for 235 Re  . Meliga et al. [9] defined it as a posteriori justification of the validity range of the linear stability analysis about the base flow.Being unstable to the spiral wave disturbance proves the hypothesis of Gartshore [20] who believed that downstream of flow reversal might be unstable due to the spiral vortex breakdown.In the next section the results of the control methodology are presented.During the control process, the stability analysis is conducted to quantify the effect of the vorticity reduction method on the unstable mode.

Optimal control of the axisymmetric base flow
The predictive control of the axisymmetric base flow minimizing the objective function of Eq.7 through radial and axial volume forces at Re=180 is conducted.The contours of the circumferential velocity and streamlines after the convergence of the control process are presented in Figure 5.The control methodology was successful in quenching the vortex breakdown.Figure 3a shows the volume force contours associated to the controlled flow.It appears that the flow control method is mostly based on the increase of axial momentum.The contours shown in figure 6 represent the magnitude and components of vorticity for uncontrolled and controlled flows.As a result of the control, the zero axial vorticity area occurring at the start of the bubble or stagnation point has been removed.Lopez et al. [26] explained that recirculation zones are caused by a negative azimuth component of vorticity, which produces an axial velocity opposite to the base flow.It has been observed by Spall and Gatski [27] that within the breakdown bubble, the level of axial vorticity has been reduced as a result of the vortex tilting, causing a redistribution of the vorticity to other components within the breakdown bubble.Furthermore, Leibovich [28] reported that vortex breakdown was accompanied by a change in azimuthal vorticity sign.Using these theories with the results in Figure 6 (left), negative azimuth and radial vorticities are formed outside of the recirculation region, resulting in the reversal of axial flow.Since negative azimuthal vorticity is a sign of bubble vortex breakdown, the optimality control system increased it to zero to eliminate the bubble (Figure 6f).A decrease in axial vorticity within the bubble recirculation region (12 z ) was transferred to the other components (Figure 6, left).Controlled flow involves increasing the axial vorticity inside the bubble while decreasing the other components.Vortex breakdown resulted in the opposite transfer of vorticity components.In Figure 7, the growth rate of the unstable global mode is shown as a function of the eigenfrequency through the control process.Growth rates are normalized by the initial maximum growth rate.In addition to mitigating vortex breakdown, the control strategy was able to stabilize the flow.The increase in eigenfrequency indicates that although the amplitude of the instability is decreasing, its frequency is increasing.Stabilizing the unstable spiral vortex breakdown requires much less force than suppressing the bubble.It can be seen in the results of Pasche et al. [15] where the force required to stabilize the flow was lower than the force to minimize radial velocity which means bubble suppression.The magnitudes of the force and vorticity at stabilization and at the end of control are shown as area averages in Table 3.To remove a bubble, a force 10 times greater is required.Additionally, it indicates a reduction in the magnitude of vorticity.Table 3. averaged values of force and vorticity at stabilization and at the end of mitigation.

Variables
Stabilizing values Controlled values

Conclusion
A predictive optimal control algorithm has been developed to control the vortex breakdown present in a swirling flow at Re=180 with S=1.095.By using radial and axial volume forces, the optimization method aims to minimize vorticity magnitude.The method was effective in mitigating vortex breakdown.Optimal forces were applied to the flow in order to convert it to a columnar flow without vortex breakdown.Furthermore, the growth rate of the unstable mode decreased through the optimization and reached negative values, indicating that the flow has become stable.Physically, volume force control may not be possible.As a result, it would be better to replace it with more realistic control variables that can be experimentally tested, such as boundary control variables.

Figure 2 .
Figure 2. Computational grid of the axisymmetric domain.

Table 1 .
Boundary conditions of base flowTo perform the linear stability formulation of a laminar flow, the velocity and pressure in the Navier-Stokes equations (Eq.2) are replaced by u + u and pp + , in which u and p are the base flow velocity and pressure fields extracted from numerical simulation.u and p represent the small-amplitude periodic velocity and pressure perturbations.Ignoring the nonlinear terms yield the following linear equation for u The optimality system of Eqs.9-11 in addition to the eigenvalue equation of Eq. 6 are solved using Freefem++ software.

Figure 5 .
Figure 5. Contours of (a) the radial (upper) and axial volume forces, (b) Contours of the circumferential velocity (upper) and surface streamlines (lower) colored by the magnitude of the axial velocity.

F 7 Figure 7 .
Figure 7. Growth rate and frequency variation through control.