Modeling and Control of Dynamic Stall Loads on a Smart Airfoil at Low Reynolds Number

: This article describes the development and testing of a modiﬁed, semi-empirical ONERA dynamic stall model for an airfoil with a trailing edge ﬂap—a “smart airfoil”—pitching at reduced frequencies up to 0.1. The Reynolds number is 10 5 . The model reconstructs the load ﬂuctuations associated with the shedding of multiple dynamic stall vortices (DSVs) in a time-marching solution, which makes it suitable for real-time control of a trailing edge ﬂap (TEF). No other model captures the effect of the DSVs on the aerodynamic loads on smart airfoils. The model was reﬁned and tuned for force measurements on a smart NACA 64 3 -618 airfoil model that was pitching with an inactive TEF and was validated against the measurements when the TEF was activated. A substantial laminar separation bubble can develop on this airfoil, which is challenging for modelers of the unsteady response. A closed-loop controller was designed ofﬂine in SIMULINK, and the output of the controller was applied to the TEF in a wind tunnel. The results indicated that the model has a comparable accuracy for predicting loads with the active TEF compared to inactive TEF loads. In the fully separated ﬂow regime, the controller performed worse when dealing with the development of the laminar separation bubble and DSVs.


Introduction
The interest in affordable wind energy production has led to technological innovations such as the smart rotor in which the time-dependent loads, induced by unsteady conditions, are controlled by active aerodynamic devices. For an airfoil in periodic motion, the unsteady flow is commonly quantified by the reduced frequency, k = ωc/2V ∞ , where ω is the angular frequency, c is the airfoil chord length, and V ∞ is the free-steam velocity [1]. There are many factors that cause unsteadiness, including the complex nature of wind, yaw, the aeroelastic response, particularly of the blades, and tower shadow effects [2,3]. As a consequence of all these factors, the blades are subjected to strong unsteady fluctuations in the angle of attack (α defined in Figure 1) and V ∞ , which can lead to dynamic stall. Dynamic stall delays flow separation to an angle that is much higher than the static stall angle (α ss ) and is followed by the development of a large dynamic stall vortex (DSV), which generates overshoots in the aerodynamic loads while convecting over the blade [4]. In addition, after the DSV leaves the blade, a secondary DSV might start to shed, inducing an additional overshoot of lift (C l ). These loads can excite the vibrational modes of the blades and increase the vibration. Although structural failure caused by such vibrations might not occur immediately, they will reduce the fatigue life of the blades and perhaps other components as well [5]. Utilizing an aerodynamic control device, such as a trailing edge flap (TEF), can alleviate fatigue loads through its high-frequency response and is additionally attractive because of its low power consumption. However, TEFs have yet to be included in wind turbine blades as more experiments are required to optimize the design parameters, to introduce an effective control strategy, and to provide modeling tools for the unsteady dynamics of the TEF [6]. This begins with numerically and experimentally examining smart airfoils, which is the generic name for any airfoil with a TEF or other control device, in a wind tunnel. The majority of unsteady aerodynamics models of smart airfoils start with Theodorsen's potential flow model of a flat plate with a TEF [1]. These models were introduced by Leishman [7] and Hariharan and Leishman [8] using an indicial function to model the TEF under unsteady conditions for helicopter airfoils. The unsteady aerodynamics of thick airfoils used in wind turbines are poorly approximated by these models. Therefore, the Leishman model was modified by Anderson [9] and Anderson et al. [10] for smart airfoils suitable for wind turbine blades. They validated the modified model against measurements of a smart Risø-B1-18 airfoil. Aside from the potential flow model, a measurement campaign was conducted by Medina et al. [11] on a NACA 0006 airfoil model to validate the accuracy of the Goman-Khrabrov state space model [12], which is based on modeling the location of the flow separation on the airfoil surface. Raiola et al. [13] used a reduced-order model to determine the steady and unsteady loads for a NACA 0018 airfoil with a 0.25c TEF. The overall loads were estimated from the loads for zero TEF angle, β, as defined in Figure 1), added to the contribution from deflecting the TEF. This model provided a fair correlation to the unsteady measurements despite neglecting the change in the effective angle of attack, in response to the TEF motion. In addition, it needs to be tested for cambered wind turbine airfoils over a wide range of operating conditions. Nevertheless, any nonlinearities in the steady curves of lift and drag, for example, due to the formation of a laminar separation bubble, might be poorly estimated [13]. A similar approach was provided by Loren [14] via collapsing the steady aerodynamic load curves at different TEF angles for the NACA 0012 and SC-1095 airfoils to be used as an input for the ONERA-EDLIN dynamic stall model [15]. Loren provided a single equation for the steady curves in terms of α and β. None of these models include the fluctuation in the aerodynamic forces generated by the DSV shedding for smart airfoils. Modeling the DSV effect on smart airfoils requires first predicting the onset of the dynamic stall resulting from the flow separation generated by the DSV, which causes rapid increases in the aerodynamic loads [16]. Modeling the onset of the dynamic stall for smart airfoils is more complex compared to clean airfoils due to the change in the effective angle of attack caused by the TEF motion.
Bak et al. [17] conducted wind tunnel experiments on a Risø-B1-18 airfoil model equipped with a morphed piezoelectric TEF. By introducing a phase shift between the sinusoidal motions of the TEF and the airfoil, a load reduction of 80% in C l was achieved when the airfoil was oscillating by α a = 0.7 • only. The impact of actuating the TEF on the DSV formation and detachment under deep dynamic stall conditions was investigated by Gerontakos and Lee [18]. They concluded that the TEF is not capable of controlling the development of the DSV, but the upward flap deflection at certain positive α results in a thinner and weaker DSV. Instead of using the typical approach of pitching the airfoil and utilizing the TEF to alleviate the unsteady loads, Frederick et al. [19] used a bluff body in front of a smart NACA 0012 airfoil model to generate a flow disturbance that simulates the unsteady loading of a wind turbine blade by atmospheric turbulence. In this experiment, the 0.04c long TEF was controlled using Proportional-Integral-Derivative (PID) and Linear Quadratic and Gaussian (LQG) controllers to mitigate the fluctuation in the unsteady loads. The error signal of the control model was the difference between the unsteady C l and the time-mean C l . The aerodynamic model used in this model was based on Leishman's modification [7]. Using a similar approach, Velte et al. [20] used a PID control scheme to mitigate the unsteady C l on a NACA 64418 airfoil model equipped with a 0.15c TEF. The inflow was oscillated harmonically by a collective pitching of two wings upstream to yield a change in the effective angle of attack of the airfoil between 0 • and 10 • . The unsteady thin airfoil potential flow theory model developed by Gaunaa and Andersen [21] was used in the experiment to model the unsteady C l on the airfoil. The PID-controlled TEF achieved a load reduction of 80% at a reduced frequency of k = 0.028. Nonetheless, the effectiveness of the TEF in alleviating the unsteady loads reduced with increasing k. Zhu et al. [22] duplicated this experimental configuration numerically in a finite volume incompressible solver. The maximum C l was reduced theoretically by 95% at a reduced frequency of k = 0.028, which is higher than the outcome of the experiments of Velte et al. [20] because the simulations did not involve any time delay in the response of the mechanical system and the controller. All of these studies are consistent in suggesting that the standard PID controller might be sufficient to alleviate the load changes in the attached flow regime. Nonetheless, the application of feedback-controlled TEF should be extended to include the separated flow regime, which can be achieved by utilizing a dynamic stall model that can predict the impact of the DSV shedding on the aerodynamic loads.
The present work provides one of the first applications of a dynamic stall model that can reproduce the load overshoots associated with the shedding of multiple DSVs for smart airfoil control in a time-marching solution. The model presented here is an extension of our previous work [23] to include the aerodynamic forces on smart airfoils. For the reasons given in [23], the model uses the chordwise and normal forces, shown in coefficient form, C a and C n , respectively, in Figure 1b instead of the commonly-used C l and C d . The model allows investigation of the control authority of the TEF for alleviating unsteady loads in the fully separated flow regime. Furthermore, the specific aim of the designed closedloop controller is to assess the model accuracy for the non-sinusoidal motion of the TEF. Therefore, direct force measurements in a wind tunnel were obtained on a smart NACA 64 3 -618 airfoil model, which will be shown to exhibit extreme aerodynamic behavior at low Re. Substantial laminar separation occurs at low α followed by an abrupt formation of a laminar separation bubble (LSB) at an Re-dependent α, which causes abrupt increases in both lift and drag. Previous studies of this airfoil were limited to the steady characteristics at low Re [24][25][26]; therefore, the airfoil selection allows investigation of a dynamic stall model that can handle very complex aerodynamics.
The following section describes the details of the experimental setup for directly measuring the steady and unsteady aerodynamic loads in the wind tunnel. In Section 3, the aerodynamic characteristics of a clean model of the airfoil are presented and discussed under steady and unsteady conditions. The proposed dynamic stall model and its application are detailed in Section 4. The model is tuned and validated in Section 5, while the integration of the model with a PID control scheme is presented in Section 6, followed by a conclusion section.

Experimental Setup
The wind tunnel experiment was conducted in a 1 × 1 m wind tunnel at the University of Calgary. Two models of a NACA 64 3 -618 airfoil were 3D-printed, with a chord length of c = 0.15 m, an aspect ratio of 4.8 and blockage ratio of 7.5% at α = 30 • . The first model represented the clean airfoil while the second was equipped with 0.2c plain full-span TEF, shown in Figure 1a, which can deflect up to an angle of β = ±20 • , defined in Figure 1b. The clean model and the main part of the other model consisted of five 3D-printed sections glued together on a 16 mm steel shaft. The final assembly was sanded and coated with WEST SYSTEM epoxy for a smooth surface finish. The contact between the main airfoil body and the TEF, shown in Figure 1b, was sealed to block any airflow from the pressure side to the suction side. The experiment was carried out at V ∞ = 10 m/s, corresponding to Re = 10 5 with a turbulence intensity of 0.3% and at 5 m/s for Re = 5 × 10 4 where the turbulence level was nearly the same. A timing-belt system, driven by a stepper motor, was used to pitch the airfoil model, as illustrated in Figure 1c. The stepper motor was integrated with an optical encoder that facilitated the acquisition of α with a resolution of 0.0625 • . A Trinamic stepper motor driver (TMCM-1110) received the step and direction commands from the Arduino controller board. The driver module has a resolution of 256 micro-steps, which corresponds to a motor angle of 0.00703 • . A Dynamixel RX-64 Smart Serial Servo, shown in the second photo of Figure 1c, with a continuous torque of 2.7 Nm at a speed of 27 rpm, was used to actuate the TEF. The smart motor has an internal potentiometer for position feedback, with a resolution of 0.3 • .
The models were mounted vertically at the center of the wind tunnel test section, between two ATI multi-axis force sensors, as illustrated in Figure 1c. The capacity of the lower force sensor is 32 N in the streamwise and cross-stream directions, with a force resolution of 0.00625 N corresponding to C n = 0.0007 at Re = 10 5 . The upper force transducer has a capacity of 20 N in the same directions, with a force resolution of 0.005 N corresponding to C n = 0.00056 at Re = 10 5 . The upper force transducer was mounted on a spherical rod end to rotate freely with the airfoil model. A National Instruments data acquisition card with a 16-bit resolution, NI-DAQ 6212, was connected to each force sensor. The two cards were externally triggered to commence acquisition at the same time. Two 0.9 × 0.9 m acrylic endplates were used to alleviate the free-end effects. The endplates were equipped with a 3D-printed super-elliptic leading edge to reduce the probability of flow separation at the leading edge [27]. The lower plate had a semicircular slot to allow the shaft to pass through it and to mount the Dynamixel motor at the end of the TEF shaft beneath the lower plate. In addition, two vertical 0.9 × 0.9 m acrylic plates were mounted as sidewalls to produce a closed test section, as illustrated in Figure 1c. They were aligned to the sidewalls of the wind tunnel contraction. The force sensors measured the total loads acting on the airfoils, including the inertial, frictional, and aerodynamic forces. Therefore, for every test run, the true aerodynamic responses were deduced by subtracting the force measurements when the wind tunnel was on and off.
The measurements are presented without blockage or other corrections because none are available for dynamic tests. In addition, dynamic stall measurements have low sensitivity to blockage effects and are qualitatively similar after applying any corrections [28,29]. The frequency of the load fluctuation due to DSV shedding strongly depends on V ∞ and the rate of change of the normal distance between the airfoil surface and the DSV. In these experiments, the DSV shedding results in fluctuating C n of the clean airfoil at a frequency of f v = 16 Hz for V ∞ = 10 m/s. Therefore, the natural frequency ( f n ) of the test apparatus in the direction of C n was analyzed by impact loading on the airfoil model and recording the normal-force readings. Figure 2 shows the results of the impact loading on the clean model and the smart one. The peaks in this figure correspond to the natural frequency of the system, with values of f n = 37 and 33.4 Hz for the clean and smart airfoil, respectively. This decrease in f n is attributed to the increase in the shaft length of the model with TEF to accommodate the TEF motor assembly and its corresponding weight, which directly affect the stiffness and weight of the moving parts. Thus, experimentation on the models at V ∞ = 10 m/s achieved the condition of f v < f n for both cases and did not excite the mechanical structure by introducing underdamped oscillations in C n . Therefore, all force readings in this work were filtered using a cut-off frequency of 20 Hz to keep the aerodynamic content and reduce structural interference and environmental noise. A detailed uncertainty analysis was carried out on the steady C l by following the ISO guidelines [30] and the approach described in [31], and the important results are shown in the next section. The measurements were grouped in four sets of data to achieve four main tasks, namely determining the aerodynamic response of the airfoil, tuning the model, validation of the model accuracy, and studying the model effectiveness for closed-loop control. The airfoil was pitched sinusoidally at four mean angles of α m = 5 • , 10 • , 15 • , and 20 • , with amplitude angles of α a = 5 • , 8 • , 8 • , and 8 • , respectively, while the TEF was inactive (β = 0 • ). These load cases were repeated when the TEF was deflected with the same frequency and with a phase shift between the sinusoidal motion of the airfoil and TEF of φ = 0, π/2, π, and 3π/2. Later, a PID controller for the TEF's motor was designed in SIMULINK/MATLAB.

Steady and Unsteady Measurements for the Clean NACA 64 3 -618 Airfoil
In Figure 3, the lift from lifting-line theory with a slope of 2π/rad is compared to the measured C l and the experimental data obtained by Mack et al. [24] at Re = 6.4 × 10 4 for the clean model. The error bars represent uncertainty in the measured C l . To facilitate a comparison with the results obtained by Mack et al., wind tunnel corrections [32] were applied to the present measurements. For instance, the wind tunnel corrections reduced the maximum measured C l by approximately 0.17. C l qualitatively matches the experimental data obtained by Mack et al. [24] in terms of the slope at α < α ss and magnitude at α ss (α ss = 15.6 • for the current study), which validates the selection of the aspect ratio of the airfoil model and 2D flow conditions. However, there is a discrepancy at α = 0 • , owing to the different approaches followed by Mack et al. [24] and the present work in determining the zero angle of attack. Mack et al. [24] considered the α at which C l = 0 to be the zero angle of attack, while the present study assumed α = 0 • when the airfoil was normal to the wind tunnel test section as measured using a special 3D-printed fixture. Note that a small adjustment in α of around 2 • would significantly improve the agreement between the two sets of measurements. Despite the positive camber of this airfoil, C l is negative for 0 • < α < 2.5 • , which can be attributed to the separation of the laminar boundary layers on both sides of the airfoil [33]. As the adverse pressure gradient on the suction surface increases during pitch up, the laminar separation point moves upstream, which results in the lower C l slope at α < 10 • compared to the 2π slope. At elevated values of α, the C l slope increases abruptly due to the sudden reattachment of the separated shear layer and the formation of an LSB to have the plateau in C l shown in Figure 3 for all measurements [24,33,34]. The existence of the LSB on the suction side abruptly increases the effective camber, causing a sudden increase in dC l /dα, which then exceeds the theoretical slope of 2π. This slope comes from the thin airfoil theory, which takes no account of viscous effects. Despite the high sensitivity of the LSB to the testing conditions, which cannot be perfectly duplicated, along with the slightly modified airfoil used by Mack et al. [24], the plateau width is comparable for both studies, Figure 3, before it bursts and C l drops suddenly at α = 17 • . The sudden occurrence of the LSB was found using smoke visualization by Mueller and Batill [33] for the NACA 64 3 -018 airfoil that shows a similar aerodynamic response, and the pressure measurements of Heine et al. [34] for the NACA 64 3 -618 airfoil. At α = 0 • , a laminar separation occurs over the rear portion of both the suction and pressure sides. As α increases, the laminar separation point moves upstream owing to the increase in the adverse pressure gradient over the airfoil suction side. Once the suction pressure downstream from the laminar separation point is approximately equal to the suction pressure of an unseparated turbulent boundary layer without an LSB, the separated turbulent shear layer reattaches [33]. The plateau in C l can be attributed to the approach of the LSB to the leading edge, which increases the airfoil's effective camber at this location. This, in turn, balances the adverse effect of the progression of the trailing edge separation at high values of α. Increasing the effective camber near the leading edge is more effective than increasing it at any other location [35] and modifies the flow substantially to increase the associated suction pressure.
In Figure 4, deep dynamic stall responses (C n and C a ) of the NACA 64 3 -618 airfoil at Re = 5 × 10 4 and 10 5 are compared. The airfoil was sinusoidally pitched by α = 20 • + 10 • sin(kτ − π/2), where k = 0.1 and τ = 2tV ∞ /c. The error bars in Figure 4 represent one standard deviation, σ, from the mean of the measurements over ten cycles, with the maximum deviations of approximately σ = 2.2% and 5.6% of the maximum of C n and C a , respectively. When an airfoil is experiencing deep dynamic stall, multiple overshoots are often observed in C n , as illustrated by points P 4 , P 5 , and P 6 in Figure 4a. The height and number of overshoots are almost comparable for the two Re, confirming that this limited change in Re has a relatively small influence on dynamic stall events [36,37]. α exceeded α ss for a large fraction of the pitching cycle, inducing secondary overshoots in C n at P 5 and P 6 , corresponding to a shedding of secondary and tertiary DSVs, which was observed previously [38][39][40]. Once the primary DSV detaches, C a drops abruptly as the flow separates fully; no fluctuations were observed in this drop [41], as shown by the unsteady C a in Figure 4b. Compared to the steady data for Re = 10 5 , no abrupt increase is seen in steady curves at Re = 5 × 10 4 due to the development of the LSB. This difference is manifested in the dynamic stall loads at Re = 5 × 10 4 and Re = 10 5 by comparing the peaks of the unsteady C a in Figure 4b, and the first quarter of the unsteady C n loop at points P 1 , P 2 , and P 3 in Figure 4a. At Re = 10 5 , the initiation of the LSB is delayed from α = 9 • for the steady case to α = 14 • at point P 1 for the unsteady response in Figure 4a, which is consistent with the results of Lee and Gerontakos [42] for the NACA 0012 airfoil at Re = 1.35 × 10 5 . Carta [43] theoretically demonstrated that the ratio between the unsteady adverse pressure gradient and the steady one is less than unity, which delays the stall. Further, the lower strength of the adverse pressure gradient slows the upstream movement of the laminar separation point. The impact of the separated shear layer reattachment on the unsteady C n is more gradual compared to the steady case, which is likely due to the progression of the trailing edge separation of the turbulent boundary layer downstream. The occurrence and bursting of the LSB have a minimal effect on the development and shedding mechanisms of the DSV when comparing the primary vortex overshoot in C n (the difference between points P 3 and P 4 ) and the time of the peaks P 4 and P 5 of C n loops in Figure 4a. This result is in agreement with previous experimental findings [44,45], which concluded that the bursting of the LSB is not responsible for the dynamic stall phenomenon. However, it probably contributes to the maximum value of C n , which can be seen when comparing the maximum value of C n at Re = 10 5 and 5 × 10 4 .

The Modified ONERA Dynamic Stall Model for Smart Airfoils
The original ONERA semi-empirical model for the unsteady aerodynamic forces in the time domain utilizes differential equations for the aerodynamic forces [15]. The final equations of this model were guided by observing the behavior of the aerodynamic forces of a pitching airfoil at different values of α and k. The version used here is an extension of the modified ONERA dynamic stall model for clean airfoils, developed by Mohamed and Wood [23], to smart airfoils with emphasis on modeling the impact of the DSV development. This modified model is distinguished by its novel approach in determining the onset of the dynamic stall (rapid increase in C n , at the angles shown in Figure 5a) by monitoring the peak of the computed C a in the time-marching simulation (independent of k), which can be utilized efficiently for smart airfoils. As an airfoil pitches up (increasing α), the shedding of the primary DSV is accompanied by a drop in the unsteady C a [23]. Modeling the onset of dynamic stall for a smart airfoil is more complex than for a clean airfoil due to the change in the effective angle of attack, in response to the TEF motion, as shown in Figure 5a. For the same sinusoidal motion of the airfoil, actuating the TEF with different phase angles influences the onset of a dynamic stall, as illustrated by the unsteady C n in Figure 5a. However, the peak C a can follow this change, as shown in Figure 5b, by the angles being lower than the corresponding ones in Figure 5a. The delay time between these angles is around τ D = 2.5. This accords with our earlier observations for the S809 airfoil [23], which showed the peak of C a occurs at time τ D = 2.5 after the initiation of the primary overshoot in C n . Airfoil shape changes can be due to a deflected portion of the airfoil, which can involve a TEF, a variable leading-edge droop, or a variable camber, etc. In the present study, the steady curves for the smart airfoil change with the position of the TEF. To implement the dynamic stall model for a smart airfoil, curves of the steady force coefficients for the range of TEF deflection are required to create a static database that can be used for any arbitrary motion of the airfoil and TEF. This model is mainly driven by the steady load loss (∆C q = C q,lin − C q ), which depends on the airfoil shape. The subscript q can be a or n to represent either C a or C n , respectively. C q,lin is the extrapolated force coefficient of the attached flow regime in the separated flow regime. Full details of these parameters and their determination for a particular airfoil are documented in Reference [23]. Figure 6e,f show examples of interpolating the load losses that can be used for any arbitrary motion of the smart airfoil after calculating the extrapolated linear load regime and the load loss at each β in Figure 6a-d. The steady response of the tested airfoil is discussed later in this section. The motion of the TEF modifies both the force coefficients and the effective angle of attack. C n,lin is translated up and down in Figure 6a,b, according to the direction of the TEF deflection, which implies that C n,lin accounts only for the change in the force coefficients in response to the change in β. ∆C n and ∆C a are translated horizontally in Figure 6c,d, respectively, to account for the change in the effective angle of attack in response to the change in β. However, ∆C n varies slightly between the different curves around α = 12 • in Figure 6c, owing to the nonlinearity introduced by the LSB, which would be negligible at high values of Re. The present approach is versatile and can be used with any morphing airfoil as the steady curves change in response to the actuation of the deployed aerodynamic control device. Compared to the modified model presented in [31], the model equations are presented here as a function of α and β to clarify which variable changes with the change in the TEF position. The final form of the modified ONERA dynamic stall for C n and C a is: a q = a o,q + a 2,q (∆C n (α)) 2 √ r q = r o,q + r 2,q (∆C n (α)) 2 e q = e 2,q (∆C n (α)) 2 where the coefficients s n , k n , σ o , σ 1 can be determined by correlating force measurements for the flat plate or the particular airfoil. The circulation variables, Γ 1 and Γ 2 , represent the circulation in the attached flow regime and the dynamic circulation loss in the separated flow regime, respectively. Considering a purely pitching airfoil, the kinematics of the airfoil motion are described using the induced velocity at the quarter-chord, W o = V ∞ α, and its derivative, W 1 = bα , where b is the half-chord. In Equation (3), the coefficients m n and σ n relate C n to the rate and acceleration of the airfoil pitching, respectively, whereas λ determines the time constant, b/λV ∞ . The parameters r q and a q represent the frequency and damping associated with the dynamic stall, respectively, whereas e q is a phase shift parameter. In Equation (4), it can be observed that the dynamic characteristics (r q , a q , and e q ) of the dynamic stall loads (C n and C a ) are independent of the TEF position and depend only on ∆C n (α). Their coefficients should be correlated to the experimental data of the particular airfoil as they are airfoil dependent. For more details, see References [15,23]. Using C q,lin and ∆C q as functions of α and β in the previous equations would restore the steady curves forα = 0 at the required α and β. Accounting for the delayed separation during pitching-up motion requires forcing ∆C n to have small values after α ss for a time τ d . Consequently, α ss should be identified for all the steady curves of C n as it changes in response to changes in β. The location of the point (α ss , C n,ss ) is illustrated in Figure 7a. For the current smart airfoil, α ss is in the range of 14 ≤ α ≤ 17.5 • when the TEF is deflected in the range of −20 ≤ β ≤ 20 • , as shown in Figure 6a. For all steady curves, α ss was automatically identified in the simulation when the steady C n drops sharply, i.e., dC n /dα < −20 rad −1 andα > 0. The condition dC n /dα < −20 rad −1 is identified based on the slope of the steady C n curve at α ss in Figure 7a, which is almost constant for all the steady curves in Figure 6a as the steady curves are translated with minimal distortion in response to the change in β. It is worth noting that Equation (4) is activated for C n at low α values because of the nonlinearities in C n . The airfoil exhibits static load hysteresis around the stall [46], as shown in Figure 7. The load hysteresis is generated by the dynamics of the laminar boundary layer separation and transition on the suction surface [46,47]. Pitching the airfoil down from high α delays the occurrence of the LSB and extends its existence on the suction side, in contrast to the pitching-up motion, as shown by C n in Figure 7a, for 8 • < α < 17 • . This aerodynamic behavior is in agreement with the experiments conducted by Mack et al. [24] at Re = 1.37 × 10 5 and 6.4 × 10 4 . The delay in the reattachment of the separated boundary layer during pitch down was examined using PIV measurements by Yang et al. [48] for the GA(W)-1 airfoil at Re = 1.6 × 10 5 . It was concluded that the severe reverse flow that occurs at the trailing edge hinders the flow reattachment that begins from the leading edge. Regarding the persistence of the LSB at low α during the pitching-down motion, Mack et al. [24] stated: "this observation can be explained by the fact that once a boundary layer was separated from the surface as for example in the stall region and in the region of laminar separation without reattachment, more energy (entrainment of high momentum fluid) was required to reattach the boundary layer than was required to keep it from separating in the first place", a mechanism that is not fully understood. However, this phenomenon is better explained in terms of the adverse pressure gradient and the type of boundary layer. The adverse pressure gradient at α = 10 • is not sufficient to separate the redeveloped turbulent boundary layer downstream from the LSB during the pitch-down, compared to the laminar boundary layer during the pitch-up. Using the static load hysteresis as an input for the dynamic stall models increases the accuracy of the model, particularly during pitch-down [49,50]. However, the analysis of the present experimental measurements suggested that some conditions should be applied for using the static load hysteresis as an input for the model. The steady curves for both pitch-up and -down should be prepared for every load case, according to the angles of attack α 1 and α 2 defined in Figure 7. In particular, when the airfoil was sinusoidally pitched down from an angle greater than α 1 , the steady curve was used as an input for modeling the corresponding unsteady forces. If the airfoil continued to be pitched to α ≤ α 2 , the steady curve during pitch-up was used as an input for modeling the unsteady forces. This finding is consistent with that of Yang et al. [48], who investigated the static load hysteresis on a NASA low-speed GA(W)-1 airfoil. They suggested that the flow is able to track its history. For smart airfoils, the DSV model depends on the difference C n,onset − C n,ss (β) to relate the influence of the primary DSV on C n [23], starting from the onset of the dynamic stall when τ DS = 0 for a period of non-dimensional time T v : where C n,onset is the computed C n at the onset of the dynamic stall and C n,ss (β) is the static stall value of C n at the given β. If the airfoil is still at an angle greater than α ss , the overshoot of secondary DSV is applied in the form of where the coefficients A n1 , A n2 , and T v can be determined from wind tunnel measurements.
In the simulation, C n,ss is estimated based on the TEF position at the onset of the dynamic stall to correlate the impact of actuating the TEF on the strength of the DSV and the associated load overshoot. A linear fit of C n,ss to β gave C n,ss (β) = 0.017β • + 0.863 The wind tunnel measurements at β = 0 • were used to tune the coefficients of the modified ONERA dynamic stall model, whereas the measurements at β = 0 • were used for validating the approach of utilizing the change in C q,lin and ∆C q to simulate the influence of actuating the TEF on the unsteady loads. The coefficients were identified and tuned in MATLAB, using the Genetic Algorithm Toolbox with an objective function of minimizing the Euclidean norm of the difference between the experimental force coefficients and the model response. The coefficients a 2 , e 2 , and τ d , along with the coefficients of the DSV model (A n1 , A n2 , and T v ), were chosen to tune the model response in the fully separated flow regime following the work in Reference [23]. The tuned values are presented in Table 1, while the remaining coefficients retained their original values: s l = π, k l = π/2, λ = 0.17, m l = 0.53, d 1 = −2.29, σ l = 2π, σ o = −0.45, and σ 1 = −2.29 [15].

Modeling Unsteady Loads on Smart NACA 64 3 -618 Airfoil
In this section, the results of the modified model for the cases used in tuning and validating the model response are presented under several operating conditions. A comparison between the model predictions and the unsteady measurements is shown in Figures 8 and 9 for the cases used in tuning the model when β = 0 • . The phase-averaged cycles of C n , C a , and C d along with one standard deviation of the mean measurement over a total of ten cycles of oscillations are presented. C d is computed from C n and C a by rotating the axes. At a high mean angle of α m = 20 • , the proposed model shows a good correlation with the experimental data, particularly for the overshoots associated with the shedding of the primary DSV and the application time of this effect, as shown in Figure 8 and associated previously with τ D . The peak of the unsteady C a and its location are predicted accurately, which in turn gives a good prediction of the application time of the DSV effect on C n . During pitch-down, the model sometimes poorly correlates the unsteady loops when the LSB starts to form at k = 0.05 and 0.075 for α m = 10 • in Figure 9b,c. By comparing the load cases in Figure 9a,b, increasing the k from 0.025 to 0.05 weakens the LSB excessively while pitching-down, which modifies the apparent camber in a way that the model cannot represent. The worst prediction occurs for C d in Figure 9, affected by the poor correlation of the impact of the LSB on C n and C a .
When the TEF is activated, the proposed model provides accurate predictions, as shown in Figures 10 and 11, and in the supplementary results in Appendix A, Figures A2 and A3. These results confirm the validity of the proposed model for the influence of a TEF on unsteady loads despite the complex aerodynamic response of this airfoil at low Re. The results of modeling the primary overshoot in C n in Figure 10 support the utilization of the peak of the unsteady C a as a robust indicator for the onset of the dynamic stall for smart airfoils. Our model provides excellent correlations to the effect of the LSB on the unsteady loads during the pitch-down at k = 0.025, as provided in Figure 11, compared to the prediction of the load case at the same k in Figure 9a used for tuning the model. Figure 11b demonstrates the importance of using C n,ss as a function of β in Equations (5) and (6) to generate the load overshoot in C n . Using a constant value of C n,ss at β = 0 • would generate no overshoot, as the value of the unsteady C n loop at the onset of dynamic stall is almost C n,ss at β = 0 • due to the motion of the TEF. The worst repeatability of the unsteady measurements occurs for the axial force at α ≤ 10 • , as shown in Figure A3, because of its very small values and the fact that the force transducer could not resolve the tiny unsteady variations of C a . Despite not tuning the model coefficients in the regime α < α ss , the model well predicts the unsteady loops and the change in their direction in response to the TEF deflection for most of the load cases in Figure A3. The accuracy of the modified model suggests that the unsteady measurements of a smart airfoil with an inactive TEF are sufficient to tune the dynamic stall model when the unsteady measurements for the active TEF are unavailable.

Control of Unsteady Loads on Smart NACA 64 3 -618 Airfoil
After modeling the unsteady loads on the smart airfoil, the dynamic stall model was utilized to tune the gains of a PID controller in SIMULINK/MATLAB, which is shown in Figures A4 and A5 in Appendix A. For brevity, we omit details of the SIMULINK computations of the unsteady C n and C a in the model shown in Figure A5. More details can be found in [23]. The steady force coefficients (C n and C a ) were uploaded to the model; ∆C q and C q,lin were interpolated online for the given α and β at every time step. The setpoint for each load case was the mean-time value of C l determined from C n and C a , as shown in Figure A4. The manipulated input (β) was limited to ±20 • , which required applying the anti-windup technique of Reference [51] to minimize the input saturation owing to the integral windup. The tuned gains of the PID controller were K p = 27, K i = 4500, K d = 0.025. The controller derivative was filtered using a first-order filter of T f = 0.007. At this low Re, the control effort required from the TEF is significantly increased because both of the LSB and the DSVs contribute to the load variations. Scaling the airfoil and the setup of the experiment (mounting the airfoil at both ends and directly measuring the aerodynamic forces) was computationally expensive and hindered the execution of real-time control. The aim was to achieve a range of reduced frequencies at which the DSVs develop on the suction surface to investigate the accuracy of the model and the control authority of the TEF to suppress the effect of the DSVs. The experiment required measuring the loads when the wind tunnel was on and off. Consequently, the PID controller was tuned and designed offline, and later the simulated control output (the preset β) was applied to the airfoil for every loading case. The advantage of this approach is investigating the accuracy of the model for the non-sinusoidal motion of the TEF that is required for control. Figures 12-15 present the measurements of implementing the predetermined β compared to the uncontrolled measurements and the output of the control simulation. Additionally, the extent of alleviating the unsteady C l and C d is shown in Figure 16, where δC l and δC d represent the variation of uncontrolled unsteady loops of lift and drag, respectively; and δC contr l and δC contr d represent the variation of controlled unsteady loops of lift and drag, respectively. At α m = 20 • and 15 • , the TEF weakens the DSV structure, which is clear in the overshoot height due to the shedding of the primary DSV for the full range of tested k in the unsteady C l in Figures 12 and 13 and in the unsteady C d in Figure 14. The LSB still significantly contributes to the variation of the aerodynamic loads at k = 0.025, as shown in the results of C l in Figures 12a and 13a, which is not the case at high Re. The measurements of C d in Figure 15 are affected by the limited capabilities of the force transducer to resolve the unsteady variations in C a at low α; however, more than a 40% reduction in the variation of C d was achieved, as shown in Figure 16b. Despite the use of the anti-windup technique, the TEF saturates at β = −20 • -see, for instance, the values of β in Figure 14. This can be attributed to the TEF ineffectiveness in changing the steady aerodynamic loads at β ≤ −10 • , as demonstrated by the curves of C n in Figure 6. Using β < −20 • hardly changes the steady curves of C n . By approaching the fully separated flow regime of the airfoil (load cases at α m = 15 • and 20 • ), the controller performance deteriorates-see the reduction in C l variations at α m = 15 • and 20 • compared to α m = 5 • and 10 • in Figure 16, where the tendency to develop a strong DSV increases. Since the controller acts on unsteady variations in C l only, the reductions in the variation of C d are less than those of C l at the same operating conditions, as shown in Figure 16b. The controller performance here is limited by the prediction of the effect of the LSB, the backlash in the smart motor, and the time delay in the response of the mechanical system and the controller. The results suggest the use of a passive aerodynamic-control device, such as a vortex generator, in conjunction with the TEF to suppress the effects of the DSV at high values of α m . In future investigations of the present model, a model predictive controller could be investigated to control the overshoots of the DSV.

Conclusions
The primary novel contribution of this article is developing and testing experimentally a dynamic stall model for predicting the impact of shedding of multiple DSVs on the unsteady loads of a smart airfoil with a trailing edge flap (TEF). This task was complicated by the low Reynolds number, Re = 10 5 , of the experiment, which caused a substantial laminar separation bubble (LSB) that significantly affects the loads. The model is based on our previous modified ONERA model for the dynamic stall loads on clean airfoils (without modifications). Direct force measurements were obtained in a wind tunnel for a smart 64 3 3-618 airfoil model at Re = 10 5 to tune and validate the model. Steady force coefficients during the pitching-up and -down motion of the airfoil were used as input for the model, and some conditions were applied. The model response was tuned using measurements of the smart airfoil when the TEF was inactive and validated against measurements when the TEF was active. High accuracy was achieved in all regimes except when modeling the LSB impact on the unsteady loads during pitch-down at high reduced frequencies.
In addition, the results for the validating load cases demonstrated that the inclusion of the steady curves change in response to the TEF actuation is sufficient to reproduce the loads when the TEF is active, consequently making it unnecessary to include the unsteady measurements for the active TEF when tuning the model response. Using a PID controller built-in SIMULINK/MATLAB, the TEF was capable of alleviating the variation in the unsteady loads by more than 40% for load cases that did not exceed α ss . At higher α, the TEF becomes less effective as α m increases. The control authority of the TEF was limited in terms of canceling the DSV effect on the unsteady loads, and the LSB was still a major contributor to the variation in the unsteady loads for the tested airfoil at low Re. Funding: The PhD program of the first author was funded by the Egyptian General Administration for Missions, whereas this research was funded by the Natural Science and Engineering Research Council, and the Schulich bequest to the Faculty of Engineering at the University of Calgary.

Data Availability Statement:
The experimental data and the predicted model results can be obtained from the first author.

Conflicts of Interest:
The authors declare no conflict of interest. Figure A3. A comparison between the measured and modeled C n , C a , and C d for the smart NACA 64 3 3-618 airfoil at α = 5 • + 5 • sin(0.1τ − π/2) and β = 20 • sin(0.1τ − φ): (a) φ = 0, (b) φ = π, (c) φ = π/2, and (d) φ = 3π/2. Load cases were not used in tuning the model response.  Figure A5. Figure A5. Subsystem of computing C a and C n using the modified ONERA dynamic stall. Details of the SIMULINK computations of both the unsteady C a and C n are omitted.