[ Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil
KTH Royal Institute of Technology, Linné FLOW Centre, SE-10044, Stockholm, Sweden
?; revised ?; accepted ?. - To be entered by editorial office

This work deals with the closed-loop control of streaky structures induced by free-stream turbulence (FST) in a zero-pressure gradient, transitional boundary layer, by means of localized sensors and actuators. A linear quadratic gaussian regulator is considered along with a system identification technique to build reduced-order models for control. Three actuators are developed with different spatial supports, corresponding to a baseline shape with only vertical forcing, and to two other shapes obtained by different optimization procedures. A computationally efficient method is derived to obtain an actuator which aims to induce the exact structures which are inside the boundary layer, given in terms of their first spectral proper orthogonal decomposition (SPOD) mode, and an actuator that maximizes the energy of induced downstream structures. Two free-stream turbulence levels were evaluated, corresponding to 3.0% and 3.5%, and closed-loop control is applied in large-eddy simulations of transitional boundary layers. All three actuators lead to significant delays in the transition to turbulence and were shown to be robust to mild variations in the free-stream turbulence levels. Differences are understood in terms of the SPOD of actuation and FST-induced fields along with the causality of the control scheme when a cancellation of disturbances is considered. The actuator optimized to generate the leading downstream SPOD mode, representing the streaks in the open-loop flow, leads to the highest transition delay, which can be understood due to its capability of closely cancelling structures in the boundary layer. However, it is shown that even with the actuator located downstream of the input measurement it may become impossible to cancel incoming disturbances in a causal way, depending on the wall-normal position of the output and on the actuator considered, which limits sensor and actuator placement capable of good closed-loop performance.

On the role of actuation for the control of streaky structures in boundary layers] On the role of actuation for the control of streaky structures in boundary layers K. Sasaki, P. Morra, A. V. G. Cavalieri, A. Hanifi, D. Henningson]Kenzo Sasaki,\nsPierluigi MorraAndré V. G. Cavalieri,\nsArdeshir HanifiDan S. Henningson \volume??? \pagerange???–???

1 Introduction

The control of transitional and turbulent boundary layers in wall-bounded flows could potentially lead to high benefits in terms of energy saving, owing to the fact that up to 50% of drag and therefore fuel consumptions of modern aircraft is to overcome the skin friction due to turbulent boundary layers (schrauf2005status), a ratio which is expected to increase due to the growing wing aspect ratio of new generation of aircraft and corresponding reduction in wave and induced drag components. Any reduction on the skin friction will result in significant savings for the operational cost of such aircraft.

1.1 Transition to turbulence

In the classical route to turbulence, which occurs in a low-perturbation scenario, a laminar boundary-layer solution becomes unstable to infinitesimal perturbations, which will grow exponentially in the form of two-dimensional Tollmien-Schlichting (TS) waves. When a critical amplitude for such fluctuations is reached, nonlinear interactions start to occur, which will eventually lead to three-dimensionality and breakdown to turbulence, a process which is thoroughly described in the review of kachanov1994physical.

However, if the zero-pressure gradient laminar boundary layer is subject to levels of free-stream turbulence higher than the transition to turbulence will occur via a different mechanism, which “bypasses” the classical Tollmien-Schlichting case (matsubara2001disturbance). Such behaviour is explained by the non-normality of the Orr-Sommerfeld/Squire operator which describes the flow dynamics, which is associated to non-orthogonal eigenmodes (reddy1993energy; schmid2012stability). Such non-orthogonality may lead to strong transient amplifications, which may occur even if the flow is stable. In the case of boundary layers, the upstream perturbations which undergo the highest transient amplifications take the form of streamwise elongated structures with comparably narrow spanwise scales. Such streaky structures are sometimes referred to as the Klebanoff mode, referencing to the experiments of klebanoff1971effect; more recent experiments have also identified such structures for different levels of free-stream turbulence (westin1994experiments; matsubara2001disturbance). In the works of andersson1999optimal and luchini2000reynolds it has been shown that the free-stream turbulence generates streaky structures matching those generated by the optimal perturbation, calculated from a transient growth analysis.

The physical origin of these streaks may be explained by the lift-up effect (ellingsen1975stability; landahl1980note), where wall-normal velocity disturbances cause the movement of fluid across the boundary layer; low-speed fluid is pushed away from the wall and high-speed fluid is pushed towards it. This movement creates quasi-periodic low and high-speed streaks which will grow linearly in the streamwise direction. With the growing intensity of the streaks, they become susceptible to higher-frequency secondary instabilities (brandt2002transition; brandt2002weakly), which will develop in turbulent spots, localized regions of chaotic motion (brandt2004transition; schlatter2008streak). When these spots merge, they will lead to a fully developed turbulent boundary layer.

1.2 Control

The high-dimensionality and inherent nonlinearity of the Navier-Stokes equations cause the computational requirements both of the simulated system and online actuation calculation to rapidly become intractable with the size of the calculation domain. The usual strategy to overcome these difficulties consists in the “reduce-then-design” approach (semeraro2013riccati), where the control-laws are designed off-line in a reduced-order model and tested a posteriori in the full nonlinear system, either a simulation or experiment (bagheri2009input; belson2013feedback; semeraro2013transition; semeraro2011feedback) .

Once the reduced-order model is available, a common strategy for control of boundary layer transition is to place the actuation in a region where the amplitude of the perturbations is small and to account for the convective nature of the flow via a feedforward scheme, where the actuator is placed downstream of the input and upstream of the objective position. The control action is then decided by means of measuring the input and acting to minimize a given quantity at the objective position. This can be accomplished using static compensators, such as the Linear Quadratic Gaussian (LQG) regulator (barbagallo_sipp_schmid_2009; barbagallo2011input; juillet2014experimental; schmid2016linear; fabbiane2017energy).

The previously cited works deal with the control of the transition induced by modal instabilities, such as TS waves. The control of non-modal structures is more rare and applications may be found in the works of jacobson1998active; hanson2010transient; papadakis2016closed; bade2016reactive, all of which deal with isolated streaks. This implies that the streaks are generated inside the boundary layer, either via roughness elements or via the inclusion of pairs of modes in a numerical simulation.

In less artificial studies, lundell2007reactive and monokrousos2008dns used blowing and suction at the wall and wall-shear stress measurements combined with feedforward control to delay the transition induced by free-stream turbulence, which inherently considers a much greater number of upstream modes. However, Lundell tuned the control effort for one specific configuration and Monokrousos et al. used spatially extended actuators with many degrees of freedom which would be prohibitive for practical implementations. lundell2009feedback demonstrated the drawbacks of currently available actuators and suggested they pose a considerable limitation for the control of streaky structures in flow control applications (see the review of cattafesta2011actuators for an overview of actuators for flow control applications).

The difficulty in the control of bypass transtion is that, differently from the Tollmien-Schlichting case, which corresponds to a definite modal instability, a family of streaks may be generated inside the boundary layer. Even though the resulting structure will correspond to the one generated by the optimal perturbation, as shown in luchini2000reynolds, its precise shape will be different depending on where it is generated. This poses a challenge to the actuator which in practice has to be located in an specific position.

1.3 Contribution of the present work

The present study tackles the mitigation of unsteady streaks, generated by means of free-stream turbulence, which penetrates the boundary layer via the receptivity mechanism (brandt2004transition). We assess the role of actuation by considering different strategies for the design of the resulting forcing, which gives insight into the physics behind the active control of streaks. Such strategies are useful for the design and evaluation of actuators for the active control of streaky structures.

The paper is organized as follows. Section 2 introduces the flow configuration control and estimation methods. Spectral proper orthogonal decomposition (SPOD) is applied to the open-loop data in Section 3, the result being compared to the optimal perturbation. The methods for the design of actuators are given in Section 4 with the results and discussion in Sections 5 and 6, respectively; finally, conclusions are drawn in Section 7. The appendix presents the specifics of the SPOD calculation and a detailed description of the adjoint optimization methods considered in the design of the forcings.

2 Flow configuration, control methods and estimation tools

2.1 Flow configuration

The incompressible Navier-Stokes equations model the flow,


where and are the velocity and pressure, respectively, at each time step and position , taken in the cartesian coordinates.

The same flow configuration of the parallel investigation of morra2019_inpress will be considered here. A plate of semi-infinite length lies in the plane, where no-slip conditions are enforced at . The control action is analysed via large-eddy Simulations (LES) with the pseudo-spectral code SIMSON (chevalier2007simson), which gives a high numerical accuracy. The flow is periodic along the spanwise direction and a fringe forcing, given as , is introduced in the last 20% of the domain to ensure periodicity also along the streamwise direction. Spatial coordinates and velocities are non-dimensionalized using the displacement thickness in the entrance of the domain and the free-stream velocity , respectively. The resulting Reynolds number, defined as , where is the kinematic viscosity, is 300. The computational domain for the 3D simulation is of in the , and directions, with and Fourier modes discretizing the plane and Chebyshev polynomials in the vertical direction.

A volume forcing is used to perform the control action, and its spatial shape will be obtained by three different methods, which will be introduced in section 4.

At the fringe region a number of modes from the continuous branch of the Orr-Sommerfeld Squire operators (which will be referred to as OSS modes) is forced. The prescribed perturbation takes the form of


where , the prime indicates a fluctuation and represents the eigensolution of the Orr-Sommerfeld Squire eigenvalue problem for the velocity fluctuations for a parallel flow, , and are the streamwise and spanwise wavenumber and the angular frequency, respectively. For further details concerning the method, the reader is referred to the work of brandt2004transition. A number of 200 modes, with an integral length scale of and turbulent intensity of 3.0% or 3.5% will be considered in this work. The characteristic spectrum of the free-stream turbulence seeks to represent the von Karman spectrum and is the same as in brandt2004transition and also used in morra2019_inpress to produce homogeneous isotropic turbulence. For further details the reader is referred to the previously cited works.

A localized measurement of the streamwise skin friction is used to define the inputs given by sensors , and downstream objective, . Three rows of 36 equispaced independent objects are placed with a transverse separation of , which is adequate to resolve the spanwise wavenumber content of the fluctuations considered here. Measurements are taken at and , defining input and objective, respectively. Actuation is performed at . Alternatively, streamwise positions will sometimes be referred to by the local Reynolds number based on , . For sensor, actuation and objective positions, is equal to 105000, 127000 and 150000, respectively. Figure 1 presents a scheme for the current simulation and coordinates considered in this paper.

Figure 1: Scheme for the 3D simulation of the flat plate considered. The blue and red circles represent the input sensors and actuators, respectively.

2.2 Estimation tools and control methods

For the development of the control law, the same approach as per morra2019_inpress will be followed with the construction of a linear quadratic gaussian (LQG) regulator (bagheri2009input2; fabbiane2015role; sasaki2018wave), using the eigensystem-realization-algorithm (ERA) (juang1985eigensystem; ma2011reduced) to supply a state-space representation of tractable dimension for the design of the LQG controller.

The actuation is computed from a convolution of the measurements with a kernel . The spanwise direction is discretized considering the position of the localized sensors and actuators, such that each actuator will behave as


where the index refers to each actuator and sensor, such that all the sensor measurements are considered in the computation of the actuation signal of each actuator.

The design of the LQG regulator involves the solution of two Riccati equations, one for the Kalman gain and one for the actual control kernel. Such calculation requires a state-space description of the problem, which is given in terms of a matrix , describing the system dynamics, matrices and which describe the effect of actuation of the disturbance and matrices and determining the actual measurements,


white noise is also assumed to be present in the measurement sensor. The numerous degrees of freedom of typical fluid mechanics problems require the usage of a reduced-order model for the description of equation (5). As in previous works by this group (sasaki2018wave), the eigensystem realization algorithm (ERA) (juang1985eigensystem) was chosen for this task. ERA involves the singular value decomposition of a Hankel matrix, formed by the impulse responses of all the inputs of the system which, for this case, correspond to the disturbances and actuation . For the details concerning the construction of the Hankel matrix, the reader is referred to sasaki2018wave.

The difficulty here is that the considered disturbance is formed by a great number of OSS modes, which implies that the computation of each individual impulse response is not feasible computationally. Furthermore, such impulse responses are not available for the case of an experimental implementation. We therefore proceed by a somewhat different strategy, defining a new set of “dummy” measurements which is placed upstream of the and measurements. Empirical transfer functions are then calculated between and or , following the procedure introduced in morra2019_inpress,


where and or are the auto and cross-spectra of the dummy measurement and the measurements, and , and are calculated via ensemble averaging. The discrete spanwise wavenumbers are defined by considering each localized actuator as a discrete measurement at a given position , , where .

Inverse Fourier transforming the quantities defined in equation (6),


where is the number of discrete transverse wavenumbers considered, provides empirically identified impulse responses which may be directly applied in the ERA method for the construction of ROMs for designing LQG. This procedure based only on the measured signal, which permits the definition of LQG controllers even for experimental implementations, was firstly introduced in morra2019_inpress and the reader is directed for such work for further details.

Application of ERA for this problem results in a system with modes, which accurately reproduces the empirically identified transfer functions.

Such empirically derived transfer functions may also be used to predict the time and spanwise behaviour of the measurement, at the objective position, from the measurement, when the actuator is not active in the system. The empirical transfer function is then


with resulting from the double inverse Fourier transform, as per equation (7). The prediction is then taken as the double convolution of with the measurements . This procedure may be applied to any streamwise separated measurements. Figure 2 presents a sample of the prediction of from the measurements for the set-up considered in this paper, and validates the procedure. For more details concerning the application of the proposed methodology for the time-domain prediction of streaky structures induced by free-stream turbulence, the reader is referred to the work of morra2019_inpress, which firstly introduced the method for this type of application.

Figure 2: Comparison between (a) LES field at and (b) its prediction from the empirical transfer function using wall-shear stress measurements at .

3 Spectral proper orthogonal decomposition applied to transitional streaks

In what follows, spectral proper orthogonal decomposition (SPOD) is applied to fluctuation data at the cross-stream plane at the streamwise objective position, , without any control action taking place. SPOD has been used in previous studies (picard2000pressure; cavalieri2013wavepackets; semeraro2016stochastic; towne2018spectral) with the objective of extracting the most energetic, and probable, structures in the flow, for each combination. Here the SPOD modes will be used to extract the dominant structure in the flow, to determine the best suited actuator for this application, and, finally, to evaluate how the closed-loop actuation is attenuating the streaks in the flow.

SPOD is applied to the velocity fluctuations such that they are optimal modes to represent the turbulent kinetic energy, where the modes are defined from the solution of the following integral equation,


where will correspond to an eigenfunction (SPOD mode) with corresponding , eigenvalue, and is the two-point cross spectral density, which is defined from the Fourier transform of the correlation tensor,


The correlation tensor is obtained by:


with , the three velocity components, and the expectation operator, representing the expected value of a given realization of the flow.

Equation (9) may be replaced by an eigenvalue problem (towne2018spectral) which reads:


where the elements of are calculated via an ensemble averaging,


where . For a detailed description of SPOD, the reader is referred to the work of towne2018spectral, in Appendix A, a brief description of the application of SPOD to data is provided.

The elements in equation 13 are determined by means of the Welch method, as outlined in Appendix A, with a triangular window and 80% overlap of the segments. Each segment had 100 points with a time discretization of . The total number of segments in the averaging was 90. These choices were seen to be adequate for the current application to accurately resolve the frequencies and wavenumbers of the structures in the flow, exemplified in figure 3.

The SPOD modes are compared to the flow response to the optimal upstream perturbation, which is calculated by using direct-adjoint power iterations via the boundary layer equations, as in levin2003exponential. The objective of such comparison is to determine whether the free-stream turbulence modes are inducing the optimally growing structures, which correspond to streaks for this application. The optimal perturbation is made for a given , and the comparison made to the most amplified case. The calculation is performed for different streamwise positions and the perturbation which is most amplified with respect to its initial position is chosen for comparison. We have also obtained the flow response to the optimal forcing, adapting the formalism in levin2003exponential for resolvent analysis, as shown in Appendix B. The resulting fluctuation at the final integration position is approximately the same for the optimal upstream perturbation and optimal forcing, given that they are both generated at the same streamwise position.

Figure 3 presents the comparison of the leading SPOD mode with the result of the optimal perturbation, which is found to be generated at . The behaviour of the first SPOD mode for the streamwise velocity fluctuation in the plane is also shown and highlights the characteristic streaky behaviour of the flow. The calculation was made for , as this corresponds to the most amplified frequency/wavenumber pair, as shown in figure 4.

Figure 3: Comparison of the SPOD modes in the wall-normal direction induced by the FST with the result of the optimal perturbation (a) and characteristic streaky behaviour observed in the streamwise velocity fluctuations (b), the pseudocolours vary between -1 and 1 and the arrows correspond to the velocity field induced by the wall-normal and spanwise velocity fluctuations. Modes were normalized to present unitary amplitude. Comparison for

As highlighted in figure 3, there is a good correspondence between the first SPOD mode of the velocity fluctuations induced by FST and the optimal perturbation. A similar feature had already been observed in other works (luchini2000reynolds), where the optimal perturbation is seen to be approximately independent of Reynolds number and to match the structures induced by free-stream turbulence modes.

Finally, figure 4 presents the behaviour of the first SPOD eigenvalue as a function of the the frequency and transverse wavenumber. Such analysis is necessary for the definition of the pair which will be considered in the optimization of the actuator in sections 4.2 and 4.3. Although not shown here, the first eigenvalue dominates the dynamics of this flow, being approximately one order of magnitude higher than the subsequent modes. It is clear that the dominating structures are present for , which will therefore be targeted by the optimization techniques presented herein.

Figure 4: Behaviour of the first SPOD eigenvalue as a function of the frequency and spanwise wavenumber, for the Tu 3.0 % case. Similar characteristics for the most amplified frequency/wavenumber pair are also observed for Tu 3.5%.

4 Actuators

A total number of 36 elements is considered in the row of actuation, . Each element adds a body force to the flow with a given spatial support which is modulated by a time signal ,


and the role of the control law is then to determine the time modulation, for each element. Three different actuators will be evaluated which vary in terms of their spatial support.

4.1 Vertical force only - actuator

The first actuator corresponds to a vertical body force only and it seeks to mimic the effect of ring plasma actuators, such as in the works of kim2016report and shahriari2018control who deal with a plasma actuator with a similar spatial support, acting on the wall-normal direction. The effectiveness of such actuator is related to the lift-up effect, which is a known trigger of streaks (brandt2014lift). For such actuator, we define the following spatial support, leading to a force only in the wall-normal direction,


with the other components of the forcing equal to zero, corresponding to the position of actuation, , and . The resulting spatial support along the wall-normal direction will be shown in figure 8 in comparison with the two other cases evaluated here.

The impulse response measured at the objective location and its corresponding frequency content are shown in figure 5. It should be noted that the frequency content of such actuator is concentrated close to , with a preferable spanwise wavelength , which corresponds to the most amplified streaks generated by the free-stream turbulence. The delay observed in the time-domain is in accordance with the group velocity of such structures and the streamwise separation between the row of actuators and sensors.

Figure 5: Impulse response of the blowing actuator, considering wall-shear stress as the measured quantity, in the time (a) and frequency domains (b), respectively.

4.2 Optimal forcing actuator

The second actuator to be considered corresponds to a spatial support given in terms of the optimal forcing, which is calculated at the position of actuation, . It should be noted that such forcing is different from the one considered in section 3 where the streamwise dependence on the generation of the forcing is also considered. The method to calculate the optimal forcing is outlined in the appendix and corresponds to a modification of the procedure described in levin2003exponential, using adjoint methods for constrained optimization; the goal is to obtain the forcing that leads to the highest energy gain at the position of objective, at , for the most amplified spanwise wavenumber, .

The actuation is restricted to spatially localized upstream areas by inclusion of a Gaussian mask in the optimization procedure. This avoids a spatially extended forcing which would be impractical in experimental applications, for example. As shown in Appendix B, the spanwise spatial support is imposed by a Gaussian function given by , with . The resulting spatial support along the wall-normal direction is shown in figure 8 for the span and wall-normal components; the contribution of the streamwise forcing is irrelevant, as it is comparably inefficient for the generation of streaks.

The corresponding impulse response in the time and frequency domains are shown in figure 6. As before, the most amplified streaks are located at , in accordance with the performed optimization.

Figure 6: Time (a) and frequency-domain (b) behaviour of the impulse response of the optimal forcing actuator considering wall-shear stress as the measured output quantity.

4.3 Identified actuator

Finally, the third actuator to be considered is calculated targeting the specific shape of the structures present at the objective position , given in terms of their first SPOD mode for the most amplified pair, as described in section 3. The actuator will be referred to as “identified” as it targets the structures which were previously identified at the position of objective. This procedure is also outlined in Appendix B and it is inspired in the work of tissot2017sensitivity. This actuator is expected to be the most efficient, as it targets the specific structures present in the flow and should therefore lead to their best cancellation, in accordance with the physical mechanisms behind active flow control. The resulting impulse response is shown in figure 7. The optimization was performed exclusively with the wall-normal and spanwise direction forcings, which should dominate the generation of streaks.

Figure 7: Time (a) and frequency-domain (b) behaviour of the impulse response of the identified actuator, wall-shear stress as considered as the measured output quantity.

As before, the maximum of the frequency content is consistent with the targeted streaks; the consideration of the spanwise component of the forcing causes the behaviour to be non-symmetric along the direction.

4.4 Comparison of the different forcings

The main difference on the spatial support of the forcings is on their wall-normal behaviour, as the same Gaussian mask was considered in the span and streamwise directions. The three different cases are shown in figure 8 for the wall-normal and spanwise components. The spatial support was normalized such that the energy content of the different forcings is the same.

Figure 8: Spatial support of the three forcings considered along the wall-normal direction for the streamwise (a) and wall-normal (b) components, respectively.

The two optimization techniques lead to the typical behaviour for the optimal forcing shape as in monokrousos2010global. The main difference between the optimal forcing and SPOD-based optimization is that the latter presents a peak at higher wall-normal positions, a feature which is seen to be related to the streaks the actuator generates.

Finally, the energy of the fluctuation,


resulting from the application of these forcings is shown in figure 9. The result is compared to the calculation of the optimal forcing at different streamwise locations. There is a strong dependence of the fluctuation energy on the position where the optimization is performed, as previously observed in other works (levin2003exponential). The optimal position for the generation of the optimal forcing is at and the fluctuation resulting from such forcing approximately matches the FST induced streak at the objective position for control (), as previously seen in figure 3.

Figure 9: Comparison of the energy fluctuation as a function of the streamwise position for the different forcings considered; optimal forcing calculated for different streamwise positions (black solid), optimal forcing at the most amplified position (red crosses), optimal forcing calculated at the position of actuation (blue dashed), vertical forcing (pink dotted) and SPOD-based identification (green dash-dotted). A zoom of the area of interest () is shown in the inset.

The optimal forcing calculated at the , where the actuation is actually performed for control, leads to a much higher energy at the objective position when compared to the other two approaches. However, as it will be shown later, it leads to a thinner streak when compared to the actual structures inside the boundary layer.

5 Results

The attention will be focused on two turbulence intensity cases, and , both of which present challenging scenarios on which transition to turbulence will be induced by the free-stream disturbances and where some nonlinearity is already present in the sensor/actuator region, posing a limitation in the accuracy of the considered reduced-order models. The same kernels and actuators, designed for were used for both cases which also allows an evaluation of the robustness of the controllers and optimization methods considered. Different kernels were calculated for each actuator, a necessary step since the actuators present significant differences in their impulse responses, illustrated in figures 5 to 7.

Figure 10 presents the performance indices, which consider the mean square values of the output, for the two cases evaluated considering the three different actuators. It should be noted that the performance index takes into account exclusively the wall-shear stress and therefore does not represent a metric for the disturbances throughout the boundary layer. On the other hand this parameter serves as a good evaluation of the effectiveness of the controllers themselves, as their single objective is the minimization of the wall-shear fluctuations at the objective position.

(a) Performance index for
(b) Performance index for
Figure 10: Performance indices for the two turbulence intensity cases evaluated.

All three methods present an adequate attenuation of the objective, with the identified actuator outperforming the other two. Table 1 summarizes the average reduction for the different cases. The better performance observed for the case is related to a better accuracy of the reduced-order models, where less non-linear effects are present in comparison to .

- Only
Optimal Forcing
Table 1: Summary of the closed-loop cases evaluated.

In order to evaluate the energy spent in actuation, the following metric, which is related to the energy budget for actuation, is defined;


This metric considers both the amplitude modulation , which is calculated by the control law, and the spatial support of the forcing . The behaviour of for the different cases is shown in figure 11. The energy budget of the optimal forcing actuator is about one order of magnitude lower than the one corresponding to the other two cases, a fact which is related to the streaks induced by it presenting the highest possible growth rates for the specific position where they are generated, as illustrated in figure 9.

(a) Energy budget for
(b) Energy budget for
Figure 11: Energy budget for the different turbulence intensity and actuators evaluated.

We now evaluate the effectiveness of closed-loop control with the different actuators in delaying transition. Figure 12 shows the friction coefficient and maximum root-mean square (RMS) values for the three evaluated actuators for the turbulence intensity levels of 3.0% and 3.5 %. The corresponding behaviour of the RMS values of the streamwise velocity fluctuation in the plane are shown in figure 13, for the case and in figure 14 four streamwise positions are shown in order to better highlight the effect of the different actuation schemes; similar results are also observed for the higher turbulence intensity value.

(a) Friction coefficient,
(b) Maximum RMS values along the wall-normal,
(c) Friction coefficient,
(d) Maximum RMS values along the wall-normal direction,
Figure 12: Friction coefficient and maximum rms values for the streamwise velocity fluctuation and the different actuation schemes considered in this work. The dashed line in the friction coefficient plots gives reference values for the laminar and turbulent cases respectively.
(a) , uncontrolled case
(b) , controlled with forcing
(c) , controlled with the optimal forcing actuator
(d) , controlled with the identified actuator
Figure 13: Behaviour of the RMS of the streamwise velocity fluctuation for the uncontrolled and different controlled scenarios evaluated.
Figure 14: RMS value of the streamwise velocity fluctuations at four streamwise positions as a function of the wall-normal direction for the uncontrolled and different actuators evaluated in this work. corresponds to the objective position.

The identified actuator considerably outperforms the other two on what concerns the delay in transition, a feature observed from the friction coefficient and maximum RMS values in figure 12, which take much longer to increase to values typical of turbulent boundary layers. The RMS values in figure 13 indicate that the effect of the identified actuator is more extended along the wall-normal direction, since, at the position of actuation, there is a significant decrease of RMS levels throughout the boundary layer. Furthermore, the actuation energy of the identified actuator is similar to that of the -only actuator and it is also more robust to the evaluated changes in turbulence intensity, leading also to significant delays in transition for the two evaluated cases.

As for the optimal forcing actuator, although it is capable of delaying the transition in the two FST intensity cases, it presents the lowest performance in terms of the transition delay. The RMS values indicate that its effect is more localized, which leads to an imperfect cancellation of the incoming streak. These characteristics will be further explored in section 6 by means of an evaluation of the SPOD of the actuation effect. In spite of these characteristics, the optimal forcing actuator leads to the lowest actuation energy for the evaluated cases, about one order of magnitude lower than the other two. This is related to the fact that it excites streaks with the highest energy growths and it is therefore capable of leading to a cancellation, specifically at the objective position, with less energy spent.

Finally, the -only actuator presents an intermediary behaviour, leading also to a significant delay in transition. However, it leads to the highest actuation energy.

6 Discussion

Two approaches will be used to better understand the results of the previous section: the SPOD of the open and closed-loop cases; and the possibility of causal streak cancellation, which is related to the different delays of the actuators in exciting a response in the output. Both of these analysis are related to the fact that the streamwise velocity fluctuations present their peak value above the wall, as it will be explored next. All results in this section are shown for the case. Trends for are similar and will not be displayed for brevity.

6.1 Correlations along the wall-normal direction

We first consider the RMS value of the streamwise velocity fluctuations along the wall-normal direction, given in figure 15. It is noticeable that the maximum value of the RMS is above the wall, located at . Therefore, in order to obtain a more global effect on the field the actuator should lead to changes over such higher wall-normal positions rather than in the near-wall region only. Furthermore, the fluctuations at such positions should be predictable given the considered input sensor, which corresponds to measurements of wall shear stress.

Figure 15: (a) Correlation between estimated and LES streamwise velocity fields - position of input and output were fixed at 250 and 400, respectively, and the wall-normal position was varied. Dashed line indicates input/output positions at the same wall-normal positions. (b) Root-mean square values of the streamwise velocity fluctuation as a function of the wal-normal direction

In order to evaluate if the streamwise velocity fluctuations at may be predicted from wall-measurements, the correlation between the estimated and LES fields was calculated as a function of the wall-normal positions of input and output measurements. The streamwise positions were kept at and 400, for input and output respectively, in accordance with the and measurements. The predicted field was obtained by means of the empirical transfer function approach, as outlined in section 2. The result is shown in figure 15 and demonstrates that there is a strong correlation between wall shear stress measurements () and the streamwise velocity fluctuations until , indicating that the currently considered sensors are adequate for this type of application. Differences are therefore only accountable for the considered actuators, as the method to determine the control law was the same for all of them. In what follows, all the analysis is made with the input sensor corresponding to wall-shear stress.

6.2 SPOD in open-loop

Spectral POD is calculated at the objective position, , in the plane defined by the wall-normal and spanwise coordinates. The decomposition is made per pair and the results shown here will consider , which corresponds to approximately the most amplified case used in the actuator optimization techniques. The eigenvalues resulting from such decomposition, for the open and closed-loop cases, considering the different actuation strategies evaluated in this paper, are shown in figure 16.

Figure 16: Eigenvalues of the SPOD modes calculated at position , for , for the open (circles) and closed-loop cases (squares, diamonds and crosses, for the blowing, identified and optimal forcing actuator, respectively).

The first mode is approximately two orders of magnitude higher than the second one, a fact that indicates its dominance in the flow, at the evaluated streamwise position. The three evaluated actuators lead to an attenuation of such mode in closed-loop, along with the subsequent modes higher than two. The identified actuator presents the best performance, leading to a reduction of about five times in the magnitude of the first mode eigenvalue.

In order to better distil the effect of each actuator in the flow, simulations were performed without free-stream turbulence and with the actuators with a white-noise time modulation. SPOD was then applied to the resulting data at the position of objective at the pair considered here; the objective is to extract the exact structures that are being excited by the actuators and compare them to the one corresponding to the free-stream-turbulence-induced streaks. Figure 17 presents the first SPOD mode for the different cases as a function of the wall-normal direction, for the streamwise velocity fluctuation, which is the velocity component that dominates the fluctuation.

Figure 17: First SPOD mode for the uncontrolled (solid line) and different actuation strategies evaluated in this work (dash, dash-dot and dots for the -only, identified and optimal forcing actuator).

It is clear that whereas the identified actuator leads to fluctuations in a compelling agreement with the fluctuations inside the boundary layer, the other two excite “thinner” structures with a peak value which occurs closer to the wall than the peak of the streaks inside the boundary layer. This behaviour leads to an imperfect cancellation of the incoming streak. Since a destructive interference is the physical mechanism behind flow control of convectively unstable flows (sasaki2018wave; morra2019_inpress), the -only and optimal forcing actuators lead to a lower performance in terms of transition delay. It should also be noted that the identified actuator presents a peak in its spatial support which is at a higher wall-normal location than the other two (as shown in figure 8); this is probably related to the higher peak of the streak it is identifying inside the flow.

6.3 The role of causality in streak cancellation

Figure 18 shows the impulse responses for the estimation (taken between , using wall-shear stress as the measurement, and or , with streamwise velocity as the output measurement, where is the first point of the grid above the wall, using the empirical transfer function) and actuation (taken between and or , considering streamwise velocity as the output measurement). The impulse response for estimation is representative of open-loop disturbances, whereas the one for actuation highlights properties of streaks induced by the control law.

A few characteristics can be observed in figure 18; First, the estimation impulse responses, taken for an input further upstream, present a time delay comparable to the actuation cases, indicating that the group velocity of open-loop streaks is higher than the one corresponding to the actuator-induced disturbances. The delay changes between wall and measurements, indicating a tilting of the structures, which reach the higher wall-normal position at before reaching the wall at the same streamwise position. Finally, from the time delays of impulse responses it is inferred that the three actuators induce streaks with different group velocity, with the identified case presenting the highest value and optimal forcing the lowest. The relevance of such time delays for closed-loop control are related to the possibility of a causal cancellation of incoming disturbances, as discussed by sasaki2018_inpress, and can be summarised, in simplified manner, as follows. Once a given structure is detected by the upstream sensors , it is estimated, through the transfer functions in figure 18, that it will reach the downstream objective after a time delay . The actuator should cancel this disturbance, but this cannot occur instantly, since the actuation-induced structures take a time delay to reach the objective location. If such cancellation is feasible, whereas in the opposite case it becomes impossible to cancel the incoming streaks; in the latter situation, once an upstream streak is detected in , it is already too late to attempt to cancel it in .

Considering the time delays for which the impulse reponse peaks, all actuators are able to generate disturbances that reach the objective position located at the wall before the one related to the estimated field, which can be seen by the characteristic time delays , lower than in figure 18 for all cases. The same is not true when we consider the output at , particularly for the optimal forcing actuator. This characteristic will prevent the optimal forcing case of acting where the highest energy of the streak is present, at the objective position.

(a) Estimated field.
(b) Estimated field.
(c) Optimal forcing actuator.
(d) Optimal forcing actuator.
(e) -only actuator.
(f) -only actuator.
(g) Identified actuator.
(h) Identified actuator.
Figure 18: Behaviour of the impulse responses of the estimated field and different actuators considered here. The measurement was located at the wall (left column) and at (right column), and the streamwise position was located at , corresponding to the objective of the control law. The colorbar was adjusted for each plot in order to better visualize the behaviour. The dashed line indicates the time delay for the maximum of the impulse response.

This characteristics are summarized in figure 19, where the streamwise position is given as a function of the time-delay for it to be reached for the estimation and each actuator. The delay was considered as the time when the peak value is reached, for each impulse response at the considered position. The slope of lines in the plot can thus be related to the group velocity of disturbances in the open-loop case, given by the estimation transfer function, and of the ones resulting from the three actuators. The values at the wall (considered as ) and at are shown in figure 19. The input measurements were considered as , for estimation, and , for actuation. For all actuators and both wall-normal positions it is noticeable that the impulse response of the estimation has a higher velocity. Among the three actuators, the identified one leads to structures with higher group velocity, which is particularly clear for the case, the most important one in terms of the energy content of the fluctuations. This higher velocity can be related to generation of streaks at higher wall-normal positions (figure 17), and may further justify the effectiveness of the identified actuator.

It should be noted that when the curve corresponding to the estimated field surpasses those of the impulses, these positions correspond to uncontrollable cases, as the control-induced streaks will reach a downstream position after the incoming structures one wishes to attenuate. The considered objective position of is highlighted and the impulse responses of all actuators reach it before the estimated field, which therefore leads to a causal kernel which is capable of the attenuations reported here, which result in a good attenuation of the objective quantity.

However, when one considers the case, where most of the fluctuation energy is contained, the streamwise objective position of does not correspond to a causal behaviour when the optimal forcing actuator is considered, which explains why it presents the worst performance among the different spatial supports considered here. The higher group velocity of streaks for larger prevents a perfect cancellation particularly for the optimal forcing and vertical forcing-only actuators, which present considerably lower values of group velocity.

One could try to compensate for the lower group velocity of the actuator-induced streaks by moving the objective position upstream, closer to the actuator. This would cause all schemes to be causal, both at wall and at . However, since the free-stream turbulence is continuously forcing the boundary layer, it is expected that the lower group velocity of the actuator-induced structures, particularly for the vertical-force and optimal cases, will eventually play a role in downstream areas of the flow.

(a) Measurement at .
(b) Measurement at .
Figure 19: Streamwise position for the different actuators and estimated field as a function of time. Input positions () were kept fixed at , with wall-shear stress, and , using streamwise velocity, for estimation and actuation cases, respectively. The transverse position of the measured impulse corresponds to (a) and (b) .

7 Conclusions

Three methodologies have been considered for the design of localized actuators for the control of streaky structures, induced by free-stream turbulence. A set-up close to practical applications was considered, where a large number of OSS modes was used to model free-stream turbulence. This makes it infeasible to compute individual impulse responses of each disturbance, and an empirical transfer function was derived to obtain a reduced-order model, for which application of LQG led to control laws.

Two of the methods of actuator design corresponded to optimization procedures, where either the energy of the actuator induced disturbances at the position of objective, or the difference with respect to a previously measured structure between actuator induced and open-loop disturbances was considered as the cost function, which was maximized in the first case (leading to an optimal forcing) and minimized in the second (with an identified, tailored actuator that optimally target open-loop streaks). The resulting direct/adjoint iteration algorithm was computationally efficient and led to desired results when applied in the nonlinear simulation. A third actuator, corresponding to a vertical forcing-only served as a baseline case, in line with the recent results of shahriari2018control where a ring of plasma actuators was used to excite a vertical body force in a boundary layer.

Closed-loop control with all the actuators led to significant delay in transition, and this was shown to be robust to mild changes in the free-stream turbulence level, a desired characteristic for real-life applications. Differences between the three cases were understood in terms of the SPOD of estimation and actuation fields, in open-loop, which highlighted the dissimilarities between the structures induced by the actuators and the one actually present in a boundary layer. Here, an important difference between control of TS waves and streaks appears. Whereas in the former case any actuator leads to exactly the same TS waves at downstream positions, as these are the only structures in spatial growth, for streaks a whole family of disturbances can be generated by actuators. It thus becomes important to target precisely the streaks that are actually expected in a given transitional boundary layer, and thus the identified actuator obtains a closed-loop performance superior to the other ones, since it cancels more accurately the open-loop streaks.

The distinct velocities of structures induced by the three actuators and the streaks induced by the free-stream turbulence, along with the tilting of the structures along the wall-normal direction, also plays a role in this type of application, where it may become impossible to obtain a causal cancellation of incoming disturbances, even if the actuator is downstream of the input measurement. Causality will also depend on the wall-normal position under consideration, a feature which had not yet been observed or quantitatively computed.

Finally, an evaluation of the correlations along the wall-normal direction indicated that wall measurements were adequate for the prediction of the output signals. Such technique allowed similar conclusions to observability tools without the need to perform adjoint simulations. These analysis were possible by means of the empirically calculated impulse responses, which permitted an exploration of the parameters of the problem reducing the number of required non-linear simulations.

On what concerns the design of actuators for experimental applications, it was shown that a vertical forcing only, which is currently possible to implement, should be adequate for the control of streaky structures. Better results may be obtained both in terms of the transition delay and energy budget of actuation when access to the open-loop data is available prior to the design of the actuator, where the methods outlined here for the evaluation of the forcing and optimizations should aid in the design of new actuators for flow control.

The authors would like to acknowledge the VINNOVA Projects PreLaFlowDes and SWE-DEMO and the Swedish-Brazilian Research and Innovation Centre CISB for funding. Moreover, part of this work was performed during an exchange programme at KTH, for which Kenzo Sasaki received a scholarship from Capes, project number 88881.132008/2016-01. Kenzo Sasaki is also funded by a scholarship from FAPESP, grant number 2016/25187-4. André V. G. Cavalieri was supported by a CNPq grant 310523/2017-6. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC, HPC2N and PDC.

Appendix A Calculation of SPOD modes from data

We briefly outline the approach to compute SPOD modes from snapshots taken from a simulation or experimental data. As in the main body of the paper, , which is then Fourier transformed from to ,


The spanwise direction is discretized and the continuous Fourier transform becomes a discrete Fourier transform (DFT), which is evaluated at the discrete wavenumbers . Given the periodicity of this coordinate, the DFT is regarded as the discretization of the coefficients of the corresponding Fourier series.

Consider now the discretization in time, the quantity represents the instantaneous state of . If a total number of snapshots is used, the signal may be regarded as:


where is , representing the number of spatial grid points times the number of physical quantities considered (on this case, the three velocity components and pressure). Application of the DFT directly into the lines of matrix should not be performed as the result will not converge with the number of snapshots (bendat2011random), and the order of magnitude of the error could be as high as the corresponding magnitude of the spectrum. Therefore, in order to obtain converged values of the spectral density, for calculation of the spectral density tensor, it is necessary to average the spectra over multiple realizations of the flow. This may be accomplished by application of Welch’s method (welch1967use).

Start by partitioning the full signal into blocks, each with elements, the nth block is then given as:


such that each block can be regarded as a realization of the flow. Overlapping the blocks with adjacent elements is possible and allows a higher number of blocks for the same length of the original signal, permitting a faster convergence of the statistics. The th entry in the th block is then given as , where is the number of overlapping snapshots. The DFT is then calculated for at each block,


where the th element of the block is then given as:


with and . Equation (22) represents the discrete Fourier transform of the block, with the addition of the weights which allow the application of a window function, used to reduce spectral leakage due to the non-periodicity of the block. The normalization factor ensures the transform is unitary for a square window. is the kth element of the DFT of the nth block, with a corresponding frequency ,




Finally, the cross-spectral density tensor can be estimated at a frequency and spanwise wavenumber by averaging the blocks,


and . Each Fourier coefficient at a frequency for each block, at a given , can be arranged in a matrix form,


where and is . The cross-spectral density tensor is then written compactly as,


The calculation of the cross-spectral density tensor then converges as the number of blocks and snapshots at each block is increased together (bendat2011random).

Defining the positive-define Hermitian matrix , , to account for the weight and the numerical quadrature for performing an integral on a discrete grid, the SPOD eigenvalue problem reduces to an matrix eigenvalue problem, at each frequency and transverse wavenumber,


The SPOD modes are then given in the columns of , ranked accordingly to their corresponding eigenvalues, which are in the diagonal matrix .

Appendix B Derivation of the Optimization Schemes

In this section, the two schemes for the calculation of the actuators considered in subsections 4.2 and 4.3 will be outlined. We follow the works of andersson1999optimal and levin2003exponential in a scheme appropriate for algebraically growing disturbances along the streamwise direction, for a slowly divergent mean flow. The pressure and velocity fluctuations follow the boundary layer equations which, written in matrix form, are given as:


where is a forcing applied in the three directions and , and the hat indicates quantities given in the domain. Equation (29) results from the application of the Ansatz into the linearized Navier–Stokes equations. The operators are then given as:


where and are the mean velocity components in the streamwise and wall-normal direction, respectively and is the identity matrix. The same non-dimensionalizations as in the remaining of the paper are considered here. It should also be noted that the streamwise wavenumber, normally referred to as , is not present in this equations as there will be no exponential dependence of the fluctuations and all streamwise variation is to be absorbed in , which implies that Tollmien-Schlichting waves will not be considered here. This equations are appropriate for algebraically growing disturbances, as the streaky structures induced by the free-stream turbulence.

Equation 29 can be written in compact form as


where the spatial, frequency and wavenumber dependences have been absorbed into as


where the derivative operators and represent discretized derivative operations in the wall-normal and streamwise directions, respectively.

Equation (29) is integrated in the streamwise direction using a first or second order explicit Euler method. We solve the problem subject to an initial condition at , and, since the equation is parabolic, downstream spatial marching can be performed. The discretization over the wall-normal direction is made by means of Chebyshev polynomials considering 300 points. Dirichlet boundary conditions are applied to the velocity components at the wall and at .

The objective is to minimize a cost function which considers the difference between the calculated fluctuation, at the objective position, and the SPOD of the field in open-loop. A similar approach has been used by tissot2017sensitivity to identify forcing terms in a turbulent jet. The considered cost-function is then:


it should be noted that by considering and performing a maximization, rather than a minimization, the usual procedure to obtain the optimal forcing is recovered. This will be treated as a particular case of this procedure.

We then follow the approach of pralits2000sensitivity by defining an extended Lagragian functional which includes the cost function and the constraint, which is given in terms of equation 34, with the addition of a Lagrange multiplier, which plays the role of the adjoint variable, .


The brackets represent the inner product, which is defined for two arbitrary functions as:


We then take the variation of equation 37, which leads to:


the sensitivity of the Lagrangian functional to infinitesimal changes , and . When the variation becomes zero, the Lagrangian is minimized or maximized. In this case, the third term is equal to zero, such that the state equation in (34) is satisfied, as desired. Zeroing the second term will result in the adjoint problem and corresponding boundary and initial conditions, as follows.

The operators are moved to the right side of the inner product by considering the following property of the adjoint,


where the star , refers to a conjugate transpose, when applied to a matrix operator. The derivatives are moved from the direct to the adjoint problem by integrations by parts. We then obtain:


where the subscripts and indicate that the corresponding operator has been derived with respect to or . Setting this to zero for arbitrary leads to the adjoint boundary layer equations, given by


where the subscripts and represent derivatives along the corresponding directions.

The term corresponds to four integrals. Once set to zero, it will supply the boundary conditions for the adjoint boundary layer equations. Writing explicitly the first three, we have:


The boundary conditions are then given as,




Finally, zeroing the fourth boundary term, which comes from the streamwise derivative, supplies the initial condition for the adjoint variables, set in the final point of the domain, which corresponds to the objective position:


such that, at , we have:


and .

After setting all these terms to zero, the variation of remains with a single term;


Since no source terms are being considered on the wall, the variation of may be written in terms of the gradient of the objective with respect to the forcing term,