Optimal transient growth in an incompressible flow past a backward-slanted step

With the aim of providing a first step in the quest for a reduction of the aerodynamic drag on the rear-end of a car, we study the phenomena of separation and reattachment of an incompressible flow focusing on a specific aerodynamic geometry, namely a backward-slanted step at 25 degrees of inclination. The ensuing recirculation bubble provides the basis for an analytical and numerical investigation of streamwise-streak generation, lift-up effect, and turbulent-wake and Kelvin-Helmholtz instabilities. A linear stability analysis is performed, and an optimal control problem with a steady volumic forcing is tackled by means of variational formulation, adjoint method, penalization scheme and orthogonalization algorithm. Dealing with the transient growth of spanwise-periodic perturbations and inspired by the need of physically-realizable disturbances, we finally provide a procedure attaining a kinetic-energy maximal gain of the order of one million with respect to the power introduced by the external forcing.


Introduction
The research field of hydrodynamic stability has the objective of elucidating how structures of some specific temporal frequency and spatial scale are selected and emerge owing to the amplification of small-magnitude perturbations. The comprehension of these effects is of huge relevance, since many flows of practical interest are dominated by genuine instability mechanisms that can be either enhanced or alleviated to improve performances. Typical expected benefits consist in the reduction of the operational cost of vehicles by decreasing the skin friction or the aerodynamic drag, or the extension of the operating conditions of turbomachinery by increasing the surface heat flux. The investigation is based on structural sensitivity [1,2], a theoretical concept stemming from the framework of stability analysis in laminar flows. This allows one to identify beforehand which regions of a given flow are most sensitive to a prescribed actuation, without the needs of calculating the actual controlled flow and of resorting systematically to a trial-and-error procedure, which would represent an insurmountable bottleneck. Here we apply this concept to determine where and how to control efficiently the turbulent-flow separation occurring at the rearend of a ground vehicle. Such an approach can thus be used to obtain valuable information about the most sensitive regions for open-loop control based on the underlying physics.
Separated flows often arise in industrial applications, as results of an adverse pressure gradient stemming from either operating conditions or geometrical constraints (airfoil at high angles of attack, rear-end of a blunt body). They are usually associated with a loss of performance. For a ground vehicle, the flow separation taking place at its rear-end contributes to a huge increase in the drag force, and thus also in the fuel consumption and pollutant emissions. For instance, a drag increase by 10% is expected to augment the fuel consumption by 5% at highway speeds. Moreover, flow separation causes low-frequency instabilities, which can trigger the excitation of aeroacoustic noise (sunroof cavities, side mirrors). The implementation of efficient control strategies, aimed at preventing separation itself or -when this is prohibitively costly or inevitableat alleviating its detrimental consequences, is therefore a great environmental and economical issue. The dynamics is here triggered by complex interactions between small-scale structures inside the shear layer, huge flow separations, and trailing vortices expanding far in the wake.
Many complex phenomena are attacked by means of linear perturbation dynamics, aiming at describing the fate of infinitesimal disturbances superimposed on a steady basic flow, and providing a rigorous mathematical foundation to investigate the control of fluid systems. Various perspectives have emerged, depending on whether the disturbance growth is characterized at large or short times: on the one hand, the archetype of disturbance energy amplified over asymptotically large times is the occurrence of vortex-shedding in wake flows, a behaviour called a modal instability [1]; on the other hand, the transient amplification over finite times is typically observed in channel flows, and is referred to as a non-modal instability [3,4,5].
For boundary-layer-like flows exhibiting a marginal separation, as occurring at the rear-end of a vehicle with small slant angle, a non-modal theoretical analysis can identify flow regions where the transient amplification of streamwise streaks (by the lift-up effect) is most sensitive to steady spanwise periodic disturbances [6]. In the experiments, such disturbances can result from either steady jets or roughness elements positioned upstream of the separation location, reproducing a parietal or a volumic forcing respectively. Indeed, it is very-well-known that the global dynamics of complex flows can be modified by imposing local disturbances. Typical examples are the use of surface rugosities to delay the transition to turbulence in boundary layers, or the injection of fluid into the wake of a bluff body to alleviate unsteadiness [7]. We remind that the lift-up effect is related to the vertical mixing of large-speed fluid from higher layers to lower ones, and vice versa for small-speed fluid. A streamwise vorticity perturbation arises, and evolves into a set of streamwise stripes characterized by relevant variations of the streamwise velocity, possibly with a periodic structure in the spanwise direction: the streaks.
In practice, one can use the adjoint method to calculate the gradient of some objective function (energy gain over a specified time horizon, growth rate of unstable disturbances) with respect to each actuation parameter, thus making it possible to cover large parameter spaces with a limited number of computations. This capability is useful as an aid to guide the design of efficient, tractable control strategies. In the past, it has mainly been applied to related problems of vortex shedding in compressible or incompressible laminar wakes, and the agreement between the experimental results and the theoretical predictions is excellent near the threshold of instability. Its application to flow separation for ground vehicles of practical importance constitutes a major issue, since substantial developments are necessary in order to encompass the complexity of turbulent-flow regimes, where large intervals of temporal and spatial scales strongly interact.
The paper is organized as follows. In section 2 we describe our numerical approach and its validation. In section 3 we specify the geometry under consideration and the main equations into play. In section 4 we introduce the base flow we have adopted. In section 5 we perform the linear-stability analysis and we focus on the direct and adjoint perturbations. In section 6 we analyze the control mechanisms and the associated kinetic-energy gain. Conclusions and perspectives follow in section 7. The appendix A is devoted to showing some further details about boundary conditions and adjoint equations.

Description and validation of numerical tools
We have made use of the FreeFEM++ software [8] to build a Finite-Element Method code. Such a tool solves the continuity and Navier-Stokes equations in their variational formulation, with prescribed boundary conditions (Dirichlet, Neumann or mixed). We have implemented a P1 scheme for the pressure field, and a P1b scheme for the velocity field (for which we have also tested a P2 scheme, without any appreciable change in the results) [9].
We have performed two main validation tests. In both cases, the quantitative validation makes use of a software which, by scanning printed figures from scientific articles, gives the numerical values of plotted points or lines with great precision. First, we have considered the 2D open cavity from [10]. We have implemented the same exact geometry and mesh from figure 7 thereof, consisting of a long  [10]. Plotted in ordinate is the generalized displacement thickness (2) versus the streamwise coordinate (in abscissa). flat floor interrupted by a unit square excavated below it. We have performed a qualitative validation against their figures 8c&d and 10, for the direct & adjoint perturbation and the eigenvalues. More importantly, we have focused on their figure 9b reporting the generalized displacement thickness (to be defined more precisely in (2)), and we have found a good quantitative agreement as shown in our figure 1. (Notice that the prefactor 2.71 reported in the caption of their figure 9b in [10] is incorrect, as the correct coefficient from the theory and in the plot is actually 1.72 [11].) Second, we have considered the 3D backward-facing step (inclined at 90 • ) from [12]. We have implemented the same exact geometry from their figure 1, and (because of the different numerical scheme) an approximate mesh from their figure 2. We have performed a qualitative validation against their figures 3-4-7 for the base flow, the skin friction, and the eigenvalues. More importantly, we have focused on their figure 5 reporting the separation/reattachment points, and we have found a good quantitative agreement as shown in our figure 2.

Geometry and equations
We have focused on a geometry issued from a standard ERCOFTAC benchmark, namely a backward-slanted step with a slope of 25 • with respect to the horizon- Figure 2: Validation against figure 5 of reference [12]. As functions of the Reynolds number (in ordinates), plotted are the abscissae of the reattachment point on the lower wall ∀ Re (red plus signs and inverted triangles), and -only for Re ≥ 300 -of the separation (green times signs and standard triangles) and reattachment (blue stars and diamonds) points on the upper wall. tal surface. This configuration, plotted in figure 3, represents a simplification of the rear end of a car and of a portion of Ahmed's body [13]. The x and y axes correspond to the streamwise and wall-normal components, respectively, with the origin placed at the leftmost/uppermost point of the sloping zone. This picture is assumed as invariant in the spanwise z direction, which implies that the base flow is assumed as two-dimensional, while the perturbations can present a three-dimensional character. (Different geometries investigated through this scheme can be found e.g. in [14,15,16].) The linear density of meshing points for the automatic triangulation process has been assumed as 4 on segments UD, DC, CV; 14 on segments WU, UV, VZ; and 24 on segments EW, WX, XY, YZ, ZB, BA, AO, OI, IE -a finer grid is obviously required close to the lower physical boundary. The resulting number of triangular elements employed in the numerical simulations is about 5 × 10 5 . About the nondimensionalization, we have assumed as reference units the vertical projection of the step and the uniform inlet speed. As the two quantities have unitary values for us, the nondimensionalized kinetic viscosity ν equals the inverse of the Reynolds number based on the step height. We have imposed standard inlet and outlet conditions on the left and right boundaries, respectively, and a free-slip condition on the upper boundary. For what concerns the lower boundary (a physical wall), we have imposed the noslip condition, except for the very first portion EI on which a free-slip condition has been used, in order to allow for the evolution of a boundary-layer profile [17]. (See also appendix A.) Tests have also been made in which we have varied the streamwise length of the domain (both upstream and downstream), its normal height, and the length of the segment EI for the imposition of the boundary condition. The chosen reference geometry falls in a range where convergence has already taken place. The case in which the upper boundary is a physical wall has also been briefly investigated, both in the case of the standard segment DC, and in a modified domain where the latter has a curved S-shaped profile to simulate a streamline and study the influence of confinement [18,19]. The effect of the resolution has been tested as well, by implementing a discretization of up to almost 7 × 10 5 triangles, without appreciable changes.
The full incompressible flow u p (x, y, z, t), comprising both the velocity and the pressure fields, satisfies the Navier-Stokes and continuity equations, (For an interesting discussion on the role of compressibility -not considered here -see e.g. [20,21].) In what follows, we decompose it into a 2D steady solution plus a 3D small perturbation:

Base flow
We have assumed as our base flow, U P , a steady solution of the Navier-Stokes and continuity equations, satisfying the same boundary conditions as the full flow. We have obtained this flow numerically by means of Newton's iterative method [22]. (The relevance of small modifications in the base flow was studied e.g. in [23,24].) Notice that, because of mass conservation, this type of base flow presents a speed overshoot, i.e. for some range of y the horizontal velocity exceeds unity. The vertical profile is not monotonic as in a standard boundary layer, so that it is not appropriate to define a typical width as the height at which the velocity reaches a definite percentage of the far-field value. It is therefore more convenient to quantify the boundary layer by means of the so-called "generalized displacement thickness" [10]: where ω ≡ ∂ x V −∂ y U represents the vorticity of the base flow (a scalar quantity, i.e. the z component -the only non-zero -of the vector given by the curl of the base velocity). When this profile reaches the step, we find δ 1 (x = 0) ∈ [0.08, 0.18], depending on the Reynolds number. We have taken into consideration Reynolds numbers ranging from 500 to 3000, with increments of 500. The flow separates from the bottom boundary at the step and, with growing Reynolds number, a larger and larger recirculation bubble develops in the wake, until reattachment takes place. Figure 4 displays the dependence of the reattachment point (i.e. the abscissa, after the beginning of the step, at which the vertical derivative of the horizontal velocity at the lower wall turns from negative to positive) as a function of Re. A sketch of the base flow for Re = 1000 is presented in figure 5.

Linear stability analysis
After fixing our base flow, we perform a linear stability analysis (see e.g. [25,26]).
Owing to the steadiness of the base flow and to its invariance in the spanwise direction, we can expand the perturbation field from (1) onto an exponential basis in both t and z: The resulting linearized equation for the (direct) perturbation is: We then look for the spectrum of the eigenvalues, i.e. complex values of σ such that (3) has nontrivial solutions U P , which can accordingly be defined as (direct) eigenvectors. The baseline case, Re = 1000 and β = 0 (i.e. no spanwise  dependence), is plotted in figure 6, and is stable because all the eigenvalues have negative real part. Instability can be reached in two ways: either modifying the spanwise wavenumber (e.g. β = 1 in figure 7), or by augmenting the Reynolds number (e.g. Re = 3000 in figure 8). It is worth noticing that the former operation induces stationary instabilities -as the imaginary part of the  rightmost eigenvalue still vanishes -while the latter introduces unstationary ones ( (σ) = 0 for those points where (σ) > 0, and the picture is of course symmetric with respect to the horizontal axis). The critical Reynolds number, for which the flow develops its first linear instability at some value of β, is approximately 750. A sketch of the perturbation field (for the largest-real-part eigenvalue depicted in figure 7), at Re = 1000 and β = 1, is presented in figure 9. It is evident that this eigenvector is a physically meaningful one because it is concentrated in the recirculation bubble, which is the zone where instabilty develops.

Control and gain
In this section we follow [27,28] by introducing a forcing on the right-hand side of the Navier-Stokes equation, which we assume as steady and volumic: F (x, £ g t) = f (x, y)e iβz . A key point is that this allows us not only to leave the boundary conditions unchanged with respect to (3), but more importantly to confine ourselves to steady solutions U P . Indeed, we have already analyzed and found the temporal evolution of the general unforced solutions (eigenvalues and eigenvectors) in the previous section, and what we are looking for here is just a particular solution to a steady forcing, which can thus be assumed as time-independent. Therefore we focus on the equations: We define as gain the quantity and as optimal gain: The procedure to find the optimal forcing consists in an iterative algorithm making use of the adjoint variables coupled with suitable outlet boundary conditions (see e.g. [29]) specified in appendix A.

Standard (non-penalized) case
In our reference case (standard geometry with Re = 1000 and β = 1) the magnitudes of what we have obtained numerically as optimal response and corresponding optimal forcing are sketched in figure 10. A clear problem arises here: even if on the one hand the optimal response develops streamwise for the whole length, on the other hand the optimal forcing is localized in correspondence of the step and of the sloping portion of the wall. This is not what one would expect for physical realizability, as on the contrary a localization on the horizontal upstream part would be suitable. Indeed, our aim is to take advantage of the streak lift-up, i.e. the formation of counter-rotating longitudinal vortices (see e.g. [6,30,31]), and to let these latter interact with the recirculation bubble, in an interesting example of interaction between the Kelvin-Helmholtz and the wake instabilities. In [7,32,33] the formation of streaks was implemented experimentally by placing a series of small cylinders, Figure 10: Magnitude of the optimal response and of the corresponding optimal forcing fields, |u | and |f |, in the upper and lower panels respectively, at Re = 1000 and β = 1. The gain is maximized according to scheme (4).
acting as rugosity elements. 1 (Notice that the spanwise periodicity of this array of cylinders can be described effectively via our periodic expansion in z, i.e. by means of the wavenumber β which should equal 2π divided by the array spacing.) In principle, one could introduce step functions in the integrals defining the gain, and we briefly explored this option preliminarily. However, in the next subsection we are going to attack this problem by means of a penalization method (see e.g. [34,35]).

Penalized case
Let us momentarily introduce an effective-viscosity field, ν eff (x), defined as equal to ν upstream and until the beginning of the step, and to a value some orders of magnitude larger for abscissae downstream of it. We then focus on: In this way we obtain the optimal response and forcing sketched in figure 11, which should more precisely be defined as sub-optimal because of the penalization scheme. The shape of the forcing is now what we expect for physical realizability, due to its localization on the upstream portion of the wall, but the same cannot be said about the velocity response, due to its concentrated character very different in look from the envisaged streaks appearing in the previous 1 The modified flow that one would obtain after physically placing the roughness element is clearly not the same as the one in their absence. What is meant here is that this discrepancy must be small for our theory to work (which clearly poses severe restrictions on the applicable elements), so that the modified flow should be obtainable as the sum of the original basic flow plus the weak perturbations currently analyzed. Figure 11: Magnitude of the sub-optimal (penalized) response and forcing fields, |u | and |f |, in the upper and lower panels respectively, at Re = 1000 and β = 1. The gain is maximized following (6).
subsection. In figure 12 we plot the optimal gain G as a function of β at different Re. The maximum of the curve not only obviously grows with Re, but also shifts to the right. As on the contrary the boundary-layer thickness shrinks when increasing the Reynolds number, it is interesting to plot the product between the generalized displacement computed at the step, δ 1 | x=0 , and the optimal wavenumber  . The inset shows that, within a relative maximum error of less than 5%, this product is independent of Re. β opt . This is done in figure 13 and proves that the ratio between the thickness and the optimal spanwise wavelength is almost independent of Re.

Penalized control with non-penalized response
The way to circumvent the paradox presented in the previous subsection is very easy. One can indeed find the optimal control through the penalized scheme, but of course once this forcing has been found its real action on the physical velocity must be computed with the actual (space-independent) viscosity ν. We then implement what one could call "non-penalized response to penalized-optimal control": LOOP TO FIND F opt (7) In this way we obtain the response sketched in figure 14, together with the aforementioned penalized-optimal forcing field. The fact that streaks are actually generated is confirmed by figure 15, which represents vertical cuts of the domain y ∈ [0, 1] × z ∈ [0, 2π) at eight different streamwise locations. In figure 16 we plot a comparison for the optimal G as a function of β at Re = 500 according to the three aforementioned schemes. In this completely-stable situation, one can see that -moving from the initial non-penalized scheme (black) to the final scheme proposed in this subsection (blue) -the loss in the optimal Figure 14: Magnitude of the non-penalized response to penalized-optimal control, |u |, in the upper panel, according to the scheme (7), at Re = 1000 and β = 1. The gain is maximized according to (6) and the magnitude of the corresponding penalized-optimal forcing field |f | (the same as in figure 11) is reported again in the lower panel for the sake of simplicity.
gain consists in less than one order of magnitude, but with the advantage of delivering an optimal forcing definitely feasible in terms of physical realizability. The comparison between (6) and (7) is also plotted in figure 17 at Re = 3000. Notice that scheme (4) cannot be enforced here, because the Reynolds number is larger than the critical value. Also, since this configuration is unstable for some perturbations, we also plot (in black) the gain corresponding to a non-modal response field which is computed through an orthogonalization procedure, in order to exclude spurious peaks related to the modal amplification (not interesting here).

Conclusions and perspectives
We have studied the transient growth of perturbations in a separated boundary layer, namely in the wake of a backward-slanted step (at 25 • ). We have shown that, by means of a suitable penalization method, also unstable cases are tractable in our formalism. We have been able to find situations where the optimal control is spanwise-periodic and localized on the horizontal upstream portion of the wall (which experimentally can be reproduced by an array of rugosity elements), and the corresponding response is represented by streaks.
Among future perspectives, five of them come to our mind. First, one could change the nature of the forcing term, from volumic to parietal. This would represent an external blowing or suction on the lower wall before the step, which would be implemented by imposing on this piece of boundary a condition on the velocity, that should keep zero tangential component but . The horizontal axis is z ∈ [0, 2π) and the vertical one is y ∈ [0, 1]; notice that for the four latter plots the physical domain extends below the bottom border of the figure, namely for a depth −1 ≤ y ≤ 0 which exactly equals the height shown. The black horizontal lines represent the height of the generalized displacement thickness δ 1 (x) at each location; the line is not shown in the fifth panel because happening to be placed above the top border, i.e. δ 1 | x=5 > 1. Figure 16: Optimal gain (in ordinates) vs. spanwise wavenumber (in abscissa) at Re = 500, according to the three maximization schemes (4) (black), (6) (red) and (7) (blue). assume a certain nonzero normal component (a function of x). We have already made a preliminary test for this situation, but a more profound investigation is Figure 17: Optimal gain (in ordinates) vs. spanwise wavenumber (in abscissa) at Re = 3000, according to the two maximization schemes (6) (red) and (7) (blue). The black points are the result of a process of orthogonalization aimed at excluding spurious amplifications, as happening for the blue peak at β = 3. definitely required. Second, it would be interesting to relax the assumption of steady forcing and to investigate the problem also in the temporal domain. As a consequence, one should fix a finite time horizon for the optimization, and perform back-and-forth temporal loops until convergence [36]. This is due to the fact that, by keeping the time dependence in the equations, the evolution of the adjoint field corresponds to a well-known backward-in-time integration, with "final" conditions imposed on the final time horizon. Of course the evolution of the direct field is forward-in-time, and one has to perform an optimization on the initial conditions. Also, the incompressibility of the flow is a key ingredient for the results shown here. It might be worth investigating how they change if a compressible flow is considered instead. We expect the theoretical analysis to be much more difficult, in view of the necessity of introducing an equation of state as well. Moreover, the present is a linear study, which is rigorously valid only for infinitesimal perturbations. If the perturbations are small but finite, we expect our framework to be still in excellent agreement with the real picture. However, it is evident that this check can be done only numerically, by performing simulations of the full problem, in order to understand whether the nonlinear coupling in the Navier-Stokes advection term induces significant modifications [37]. This issue could e.g. be attacked by means of appropriate Large-Eddy Simulations and provide the basis for a more direct comparison with experiments [38]. Last, but definitely not least, a relevant question arises about the stability of the considered steady flows to perturbations involving large spatial scales. Further weakly-nonlinear analysis (as in [39] for hydrodynamic flows and in [40] for MHD flows) may reveal a complex dynamics of the large-scale perturbations affecting the performance of the vehicle.
The base flow U inherits the same exact conditions. On the contrary, the perturbation field must satisfy fully-homogeneous boundary conditions, so that all the formulae above hold also for the quantities with a prime except for the very first one, which becomes (at inlet ED) u = 0.
Equation (3) can be rewritten in terms of the field u p (x, y) as: The procedure to find the adjoint equations consists in performing the scalar product of equations (3) (considered as a hypervector) with the adjoint field U † P † , and in integrating by parts on the whole domain, exploiting the convenience of our boundary conditions. In particular, the ones for the adjoint variables are the same as for the direct counterpart, except for the outlet condition which reads (on segment BC): P † I + ν∇U † + U ⊗ U † = 0.
However, when a forcing term is also present in the direct equations (as in (4)), the procedure is more complex because it involves the whole formalism of Lagrange multipliers and functional derivatives [27,28]. We do not report it here, and we simply remind the final result (5).