Open-loop and closed-loop flow control based on Van der Pol modeling

A real-time modeling and control architecture for coupled shear layer von Kármán instabilities on a wing section wake is developed. Previous wind tunnel experiments highlighted that these wake instabilities give rise to limit cycles on both local and resulting airloads. The two-state Van der Pol system is employed to model the self-sustained oscillations exhibited by the actual flow. Schemes featuring either one or two coupled Van der Pol oscillators in parallel are taken under consideration. Open-loop and closed-loop control laws to manipulate and reduce the self-sustained oscillations are extensively studied. An adaptive control scheme is developed to optimize the performance of the controller, according to the flow conditions. For the open-loop control scheme, the consistency with experimental tests performed on a wind tunnel model is shown.

of attack-is observed. Le Pourhiet et al. [11] developed a VdP oscillator model to reproduce and control closed-loop self-sustained oscillations on physical systems. The feasibility of this approach was proved by means of applications to mathematical models of buffet phenomena on a bidimensional wing and to surge effects on an axial compressor.
Formulations based on differential equations with expressions similar to the VdP system have been also used for hemodynamic problems. Petrov and Edissonov [14] used a first-order nonlinear differential equation, resembling the first derivative of the VdP model, to study the influence of the bulk elasticity, hydraulic, and peripheral resistances on the flow rate and pressure in the arterial system. In a successive work Petrov [15] used a nonlinear second-order differential equation like that of the VdP oscillator to describe the left ventricle pulsations in terms of contractile length and internal ventricular pressure. It was shown that, under defined conditions, the model behaves as a relaxant system. Bifurcation analysis of the model transitions from relaxant behavior to other ones is used to propose a way to predict the genesis for arterial hypertension and hypotension.
This modeling approach has been extensively employed for structural problems as well, see, for instance, Refs. [16,17]. In Ref. [8], the forced Mathieu equation was used to assess the limit cycles featured by the blades of a wind turbine due to wind shear, tower shadowing, and cyclic gravitational axial loading, all occurring at the rotational characteristic frequency. Additionally, Dotson et. al [18] employed a nonlinear second-order differential equation very similar to that of the VdP oscillator to model the aeroelastic response to nonlinear aerodynamic forces typical of transonic flow conditions. The behavior of arising limit cycle oscillations was assessed, and parametric studies of the modeling were carried out.
The work presented here aims to build up a low-order dynamic system capable to model the effects of wake flow instabilities experimentally observed by wind tunnel flow visualizations discussed in Ref. [19] and the pressure measurements presented in Ref. [20]. The two-state VdP system is herein used to model this specific limit cycle behavior. Open-and closed-loop control law architectures-to manipulate and reduce the self-sustained oscillations of wake turbulent structures-are then developed. The formulation and the results obtained with an adaptive control scheme are finally given.
The work is organized as follows. Section 2 discusses the development of the Van der Pol model for the turbulent flow past the wing section model. The block scheme and the formulation in the state space are provided. The unforced response of a single VdP system, as well as that of two oscillators coupled in parallel, is presented in detail.
Section 3 deals with the open-loop control of the single and of the coupled VdP oscillators. The responses to sinusoidal inputs in the time and in the frequency domain are presented and discussed. Additionally, it is proved that the forced response of the model is consistent with experimental tests on a wing section controlled in open-loop. Section 4 illustrates how the closed-loop feedback control is built up and presents the results of its application to the single and to the coupled VdP oscillators.
Section 5 provides details of the development of the adaptive control law applied to the single VdP system. The adaptation is performed either on the feedback gain, or on the feedback phase delay. Section 6 draws conclusions to the present works and proposes future recommendations within the framework of the overall research project.

Van der Pol modeling of the wake flow
The development of a proper control scheme for a complex system-like the turbulent wake past a wing section-makes worthwhile building an appropriate analytical model of the physical phenomena.
Experimental tests were carried out on a wing section model within the framework of the RTRA-STAE research project [22]. The behavior of the wake was investigated at Reynolds numbers ranging from ∼ 10 5 to ∼ 10 6 . Flow visualizations by particle image velocimetry surveys [19] and pressure measurements [20] were performed for this purpose. Figure 1 is an example of the power spectral density (PSD) of the local velocity measured downstream the trailing edge. Coupled shear layer/von Kármán instabilities are highlighted by the two peaks clearly visible in the PSD of the velocity. Namely, f SL indicates the shear layer instability peak. On the other hand, f I arises  [20]; Re ∼ 10 6 ; angle of attack: 10 • . Pressure is achieved with a Kulite transducer located at 90% of the chord in the mid-span on the model of Ref. [20] from a linear combination of the von Kármán and shear layer instability peaks. Such wake instabilities are found to give rise to limit cycles in the local, as well as in the resulting aerodynamic loads on the wing [19,25]. Figure 2 shows the PSD of pressure measured on the wing section model of Ref. [20]. The experimental campaign was performed within the framework of the RTRA-STAE DYNAMORPH research project as well [22]. Pressure is achieved with a Kulite XCQ063 transducer, located at approximately 90% of the chord on the suction side at the mid-span of the model sketched in Fig. 1a of Ref. [20]. Details of the experimental rig and of the measurements campaign are not the purpose of the present paper. More information is reported in Ref. [20] and briefly presented in Sect. 3. Two peaks at approximately 125 and 175 Hz-to be modeled with VdP oscillators-are clearly visible. Notice that these peaks are close to those of Fig. 1, achieved with PIV measurements taken during a previous experimental campaign [19,22,25].
A VdP oscillator-based model is employed here in order to reproduce the self-sustained oscillations, featured by the flow under consideration. The corresponding block scheme is reported in Fig. 3.
The block L(s) in Fig. 3 is the transfer function of an unstable linear dynamic system and reads: where ω 0 is the natural circular frequency of the system and ξ is the nondimensional damping coefficient. The system is stabilized, for large y, by a static nonlinear feedback, v = ϕ(y) = K s y 3 , with gain K s . As a result, a limit cycle appears in the system, being this latter unstable for small values of y and stable for large ones. The variable u is the external forcing, possibly present, and G is the corresponding gain.
The state space representation for the flow model, with its control actuation, of Fig. 3, is the following: where x 1 , x 2 are the states of the system, y is the measured output-here the time history of the local pressure measured on the wing-and u is the control input, i.e., the signal for the piezo-actuators on the wing trailing edge of [23]. An overview of the physical operation of these actuators is provided in the following. More details can be found in Ref. [23]. Specific information on the analytical derivation and the tuning of the model parameters is provided in Ref. [12]. Notice that the control system developed in this work is conceived to be operating in real time within the test bench at the IMFT (Institute de Mécanique des Fluides de Toulouse) wind tunnel. Therefore, in order to reproduce the digital signal processing, the VdP model is simulated in discrete time as well. Classical Runge-Kutta-like integration methods are not appropriate for this purpose. Details on the used discrete time integration scheme are provided likewise in Ref. [12]. Figure 4 shows the time history and the PSD of the unforced response, achieved with the present VdP modeling. The natural frequency of the model f 0 = ω 0 /2π is set to 100 Hz, whereas the nondimensional damping ξ is fixed to 0.07. Recalling Fig. 3, the forcing u is here set to zero. The simulated time is 30 s; however, the time histories plots are reported blowed up in the range 0.8-1.5 s, to highlight the details of the signals. The time step for the simulation is 10 −4 s.
Notice that the unforced system reaches its limit cycle oscillations regime after approximatively 1.2 s. The dynamic of the response is mostly determined by the damping coefficient of the model. The permanent regime  In order to try to reproduce self-sustained oscillations featuring two peaks-like those exhibited by the flow and plotted in Fig. 1-two VdP oscillators are coupled in parallel. The corresponding block scheme is reported in Fig. 5. The parameters of the two systems are the same, except of the natural frequency and the quantities depending on it. The natural frequencies for the first and the second VdP oscillators are set to 100 and 200 Hz, respectively. The outputs of the two systems are equally weighted, with gains of 0.5 each. Figure 6 shows the PSD of the output for the unforced coupled VdP systems (y in Fig. 5). The simulated time is 30 s, and the time step for the simulation is 10 −4 s. The two peaks of equal amplitude corresponding to the natural frequencies of the oscillators are clearly visible. This result shows how the VdP model is appropriate for modeling the wake instabilities of Fig. 1. Of course the larger the number of oscillators in parallel employed for the modeling, the more accurate is the representation of the multi-harmonic turbulent flow under consideration. Within the limit of infinite oscillators in parallel, the model tends to match exactly the spectrum of the flow.

Open-loop control scheme
Open-loop controls are first applied to the single and to the coupled VdP oscillators. Namely, the response of the system is investigated in the time and in the frequency domain, for sinusoidal external forcing. Figure 7 shows an example of time history and PSD of the VdP oscillator, when a sinusoidal forcing is applied. The external forcing has unitary amplitude and frequency of 50 Hz, see Fig. 7a, on the top. The simulated time is 30 s, but the time histories plots are reported as blowed up in the range 0-0.5 s, to highlight the details of the signals. The time step for the simulation is again 10 −4 s. The time history of the output shows that the system responds very rapidly to the input. Namely, Fig. 7a, bottom, shows that, after a first oscillation with frequency 100 Hz-from 0 to 0.01 s-the period of the response slides to that of the forcing input. The permanent regime response features oscillations of unitary amplitude and frequency of 50 Hz. The PSDs of the input and of the output are reported in Fig. 7b. Notice that the PSD of the output is computed on the signal in the permanent regime. The PSD of the system response resembles the counterpart of the forcing input, in terms of location-50 Hz-and amplitude-∼ 0.4-of the peak. Consistently with the formulation of Refs. [11,12], only the frequency of the forcing is expected to appear in the system response.
Notice that, when the forcing frequency is equal to the natural frequency of the system, one has ω 1 = ω 0 . Indeed, the formulation developed in Ref. [12] yields G = 0. That is, regardless of the amplitude of the forcing signal, the output of the system will correspond to the unforced response, see Fig. 3.
The response of the two coupled VdP oscillators to a sinusoidal input is investigated as well. Figure 8 shows the time histories and the PSDs of the input and the resulting output for the scheme of Fig. 5. Namely, a sinusoidal input with unitary amplitude and frequency 150 Hz is applied. It is worth recalling that the natural frequencies of the two VdP systems are again 100 and 200 Hz, respectively.
Similarly to what was seen for the single oscillators, the forced system rapidly switches to the harmonic of the forcing. Indeed, the PSD of the output clearly shows a single peak at 150 Hz. Therefore, properly forcing a multi-harmonic system-a low-order model of the turbulent flow past a wing section-allows to withdraw the peaks of the unforced response. Refs. [11,12] show that, according to the choice of the control and the model parameters, the response to a sinusoidal input could feature: (i) only the harmonic of the external forcing, which is the case discussed above; (ii) both the natural and the forcing frequencies; and (iii) only the natural

Experimental tests on a wind tunnel wing section model
Experimental tests have been carried out on a wind tunnel wing section model equipped with piezo-actuators on the trailing edge. The operation of these actuators, alongside a detailed description of the experimental rig, is beyond the purposes of the present paper and was described in Ref. [20]. For convenience, it appears useful to mention that the experimental model features a span-wise uniform airfoil corresponding to an A320 wing section. The chord length is 700 mm, and the span length is 590 mm. Piezoelectric actuators are located at the trailing edge, both on the pressure and on the suction side of the model. They are span-wise uniformly distributed. Macrofiber composite patches-glued on a metallic support-are activated alternatively on pressure and suction side to provide upward/downward trailing edge displacements. The amplitude of the maximum trailing edge displacement corresponds to ± 0.5 mm. The actuation is conceived to yield combined vortex breakdown and eddy blocking effects on the wake. That is, trailing edge vibration causes alternate shedding of counter-rotating vortical structures. These latter interact with the existing vortices related to the upper (suction side) and lower (pressure side) shear layers. When the trailing edge actuation is tuned appropriately in terms of frequency, amplitude, and phase, this interaction provides a beneficial vortex breakdown effect, characterized by a loss of coherence of the most energetic, large-scale flow modes. As a result, the energy of the largest  [20]. Open-loop harmonic forcing at 150 Hz; Re ∼ 10 6 ; angle of attack: 10 • . Pressure is achieved with a Kulite XCQ063 transducer located at 90% of the chord on the suction side in the mid-span on the model of [20] coherent structures in the wake is reduced. At the same time, an eddy blocking effect is achieved. Namely, trailing edge actuation introduces energy into the wake, specifically at the interface between the upper and the lower shear layers. As a consequence, smaller-scale structures are activated in the inertial range between the two shear layers. An upscale energy transfer toward larger-scale coherent eddies occurs. The resulting shear-sheltering mechanism leads to a constriction of the turbulent-nonturbulent interfaces, and the thickness of the wake is in turn reduced [26]. The effects in terms of aerodynamic integral loads consist of a reduction in drag and an increase in lift, relative to the nonactuated counterpart, see [20]. Figure 9 shows the PSD of the pressure fluctuations measured with a Kulite XCQ063 transducer located at ∼ 90% of the chord on the suction side at the mid-span. The solid black line corresponds to the nonactuated configuration, depicted also in Fig. 2. The white dot-dash line represents the actuated configuration. Piezoactuation is provided in the form of a sinusoidal law, featuring maximum displacement of the trailing edge-i.e., ± 0.5-and frequency of 150 Hz. Consistently with the effects described above, excitation at one frequency causes a redistribution of the flow modes. The two peaks related to von Kármán and Kelvin Helmholtz instabilities are remarkably reduced. As expected, the peak at 150 Hz is increased, relative to the nonactuated counterpart. It is worth remarking that the effects achieved on the experimental model are consistent with those depicted in Fig. 8 for the two coupled oscillators controlled in open loop. This proves that the proposed low-order modeling is appropriate not only to represent the dynamics of the unforced flow/loads, but also to simulate the effects of an active control on the turbulent wake. Indeed, when forcing the two coupled VdPs at a certain frequency, the response of the system switches to the harmonic content of the control input. The eigenfrequencies of the unforced system are almost completely withdrawn. This corresponds indeed to what was observed in the experiment. That is, the von Kármán/Kelvin Helmholtz instability peaks are remarkably attenuated, whereas the peak at the forcing frequency is-as expected-increased. This is a very promising result in view of the transposition of the closed-loop control architectures-presented in Sects. 4 and 5-to the experimental model.

Closed-loop feedback control scheme
The general expression of the closed-loop control input which is considered in this work is: being u the control input, y the measured signal, K the feedback gain, and τ the time delay of the command law. For the application to the actual flow, the control input is meant to be the motion law for the trailing edge piezoelectric actuators mentioned in Sects. 1 and 3. The output signal is expected to be the pressure measured locally on the wing section model of [20]. Figure 10 depicts the block diagram employed for the computation of the feedback control input. The CONTROL LAW block takes as input the measured signal and provides the real-time control law. The VAN DER POL SYSTEM block simulates the flow over the wing section. Either one VdP at natural frequency Targeted tests are carried out with time varying input gain K = K (t) or delay τ (t), to find out a preliminary optimum-i.e., y minimizing-parameters for the control law in the chosen flow conditions. The optimum values achieved for the feedback gain and delay are K = −0.2520, with τ = 0 s, and K = 0.2520, with τ = 2 · 10 −3 s, respectively. Figure 11 shows the time history of the control input u and output measure y obtained with the present approach. The plots on the right-hand side are a blow-up of the left-hand side counterpart in the range [2-2.4] s. It is clearly visible that the oscillations in the outputs are effectively withdrawn in less than 0.2 s past the application of the control input.
A proportional feedback control is then applied to the coupled VdP oscillators. Notice that, in this case, the time step of the control law is set to T c = 2 · 10 −4 s and that of the VdP system to T s = 1 · 10 −4 s. Having larger time resolution allows for better assessing the response, expected to be more complicated than for the single VdP scheme.
The simulated time is kept 30 s. For the first 5 s of the simulation, the system is kept unforced, to allow for the complete development of the self-sustained oscillations. Subsequently, a proportional feedback control with gain of K = 0.25 and no delay is applied. Figure 12 shows the time history of the input and of the output for the coupled VdP oscillators. The natural frequencies of the two systems are again set to 100 and 200 Hz, respectively.
Notice that the resulting output is almost completely withdrawn within less than 1 s since the application of the feedback.

Adaptive control algorithm
Because the flow conditions are likely far from being constant during the flight, an adaptation of the parameters of the control law is also studied. The adaptive control is presented hereinafter for the single VdP systems, but a similar scheme is employed for the coupled oscillators. Figure 13 shows the block diagram of the adaptive closed-loop control scheme developed in this work.  The real-time pressure measurements, i.e., the outputs of the VdP flow model, are used as inputs to the CRITERION block, where the variance J of the signal is computed. Specifically, the mean μ y and the variance σ 2 yy of the measure y are computed as a first-order and a second-order low-pass filter. To this aim, a dynamic system is built up, with input y and output J = σ 2 yy . The two states μ y and σ 2 yy are described by the following differential equations: where 1/T is the filtering frequency, such that T >> T c . The corresponding discrete time representation is given by: Within the ADAPTATION block, a gradient-based algorithm is employed, to calculate the optimum parameter θ which minimizes the variance J = σ 2 yy . The adaptation parameter is a generic quantity of the control law and is computed at each adaptation step as follows: where α, such that 0 ≤ α < 1, defines the criterion rate of decrease. The optimum criterion is initialized as The initialization of the criterion gradient g = d J dθ is performed by imposing two values of θ , θ 1 , and θ 2 . The corresponding values of the criterion J 1 and J 2 are computed as J 1 = J | θ 1 and J 2 = J | θ 2 . The resulting gradient reads: To avoid achieving unfeasible values for the adaptation parameters, constraints are imposed to θ such that θ min ≤ θ ≤ θ max . At every step of the adaptation, the current J current is compared to J opt . If J current < J opt , then θ opt = θ current , J opt = J current , and g opt = g current . Otherwise, θ opt , J opt , and g opt are left unchanged.
Currently, tests are performed by defining θ as either the gain, or the time delay for the input signal. Specifically, referring to Eq. (1), the following two test cases are considered: The CONTROL block takes as input the current measure y and the adaptation parameter θ and computes the control input expressed in Eq. 1. Figure 14 shows an example of the results obtained by the adaptation algorithm, with u(t) = K (t), y(t) = θ(t) y(t), i.e., with no delay in the control law. The first plot from the top represents the time history of the adaptation parameter θ(t) = K (t), computed during a 220-s-lasting simulation. The variable θ is initialized to 0 and changes during the simulation, finally converging to the value which minimizes the variance of the measure. The second plot from the top shows the variance of the output computed during the simulation and changes accordingly with y. Recalling that the adaptation algorithm aims to minimize J with respect to θ , it can be observed that the variance approaches zero, after the transitory consequent to the initialization. The third plot from the top represents the time history of the control input. The amplitude is maximum at the very beginning of the control application. Then, when a proportional feedback control is applied, u converges to zero as the output suppression is achieved. The plot in the bottom of Fig. 11 shows the time history of the output measure. Its behavior proves ultimately the effectiveness of the adaptation algorithm here developed. As the control input is applied to the dynamic system, the output amplitude decreases, approaching zero within about 100 s. It is worth noting that the adaptive control is effective since the very first instants past the application of the input u. Less than 10 s is enough to achieve a 50% reduction in the amplitude of the unforced system.
The adaptation is performed with a period of 3 s, to allow for the setting of the system to the current control input. The algorithm exploits the first-order expansion of the criterion with respect to the variable θ , and the gradient is computed as the incremental ratio between two consecutive values of J and θ . Therefore, two values of θ are required to initialize the gradient. Those reported in Fig. 14 are specifically selected to depart from the expected optimum. This aims to demonstrate the algorithm robustness with respect to the model parameters and to their initialization. Indeed, from Fig. 14, it is clearly visible that the adaptation changes the control gain K progressively, to ultimately minimize the measure y.

Concluding remarks
A low-order nonlinear model is developed, based on the Van der Pol oscillator. The model is meant to represent the flow instabilities, previously measured past a wing section model at Reynolds number ∼ 10 5 [19] and ∼ 10 6 [20]. Specifically, the model is capable to reproduce the limit cycle oscillations exhibited by the wake velocity and in turn by the wing pressure. The development of the model, together with a detailed discussion of the unforced response, is given at a first instance. In particular, it is shown how the coupling of two Van der Pol oscillators in parallel allows for reproducing the measured velocity oscillations, related to coupled Von Kármán shear layer instabilities. Open-loop and closed-loop control laws are then applied on the single as well as on the coupled Van der Pol oscillators. Applying a sinusoidal forcing to the system allows for effectively reducing the self-sustained oscillations. It is shown that the open-loop control is effective if its harmonic content differs from the natural frequency of the system. Experimental results obtained on an open-loop actively controlled wing section model are also provided. It is shown that the response of the physical system to a harmonic control input is consistent with that achieved on the forced VdP-based model. This means that the proposed low-order model is suitable not only to characterize the main physics of the unforced flow, but also to reproduce the forced response of the physical system. Therefore, an effective control architecture for the flow can be designed on this 2-degrees-of-freedom (DOFs) model, without requiring time-and money-consuming CFD computations or wind tunnel tests-at least for the first phases of the design. Feedback control laws are then presented. Applications to the single and to the coupled Van der Pol systems are provided, to demonstrate the effectiveness of the proposed control scheme. Finally, the development of an adaptive control scheme, capable of computing in real time the optimum feedback gain, is proposed. This control architecture is particularly appropriate for on-board utilizations on actual wings, as the flow conditions are expected to change significantly during flight. The novelty of this work stems from two principally intertwined elements: (i) Van der Pol (VdP) systems are widely used as models of turbulent flows; yet this is one of the first attempts to develop an open-and closed-loop real-time wake control architecture based on this formulation. It is highlighted that not only the VdP oscillator is suitable for modeling the unforced response of a turbulent flow, but also an effective control architecture for the flow dynamics can be built up on a very simple scheme, as the present two DOFs system. This result is regarded as quite remarkable, because it is shown that the VdP oscillator is a feasible basis on which to design and to assess control architectures. This avoids-at least at preliminary design stages-expensive and time-consuming CFD computations and wind tunnel tests. (ii) The closed-loop control scheme developed here exploits a nonmodel-based approach, capable of operating in real time. On the contrary, the largest part of closed-loop flow control applications is based on complex physical models-built up, for example, on the basis of modal decompositions, linear stochastic estimations, Galerkin models, and black box approacheswhich are issued, for instance, from high-order 3D computations or PIV measurements performed a priori, see, for example, Ref. [26]. As a consequence, the resulting closed-loop control schemes are actually based on scheduled data. This reduces significantly the effectiveness of the control, due to the uncertainties/disturbances of the system, and to unexpected responses-not included in the scheduled feedbacks-which can always occur. Here a single stochastic variable, representative of the physical phenomenon, is used to compute the control feedback, whereas many DOFs are usually required for CFD-or PIV-based closed-loop strategies.
Tests of the closed-loop feedback control scheme, as well as the of adaptation algorithm, are planned for the near future.