Experimental Study of Dynamical Airfoil and Aerodynamic Prediction

: Dynamic stall is a critical limiting factor for airfoil aerodynamics and a challenging problem for active flow control. In this experimental study, dynamic stall was measured by high-fre-quency surface pressure tapes and pressure-sensitive paint (PSP). The influence of the oscillation frequency was examined. Dynamic mode decomposition (DMD) with time-delay embedding was proposed to predict the pressure field on the oscillating airfoil based on scattered pressure measurements. DMD with time-delay embedding was able to reconstruct and predict the dynamic stall based on scattered measurements with much higher accuracy than standard DMD. The reconstruction accuracy of this method increased with the number of delay steps, but this also prolonged the computation time. In summary, using the Koopman operator obtained by DMD with time-delay embedding, the future dynamic pressure on an oscillating airfoil can be accurately predicted. This method provides powerful support for active flow control of dynamic stall.


Introduction
Stall is one of the most critical limiting factors for airfoil aerodynamics. Static stall occurs when the airfoil reaches a certain angle of attack. However, in real-world applications, dynamic stall caused by an unsteady motion is more common. An example of such motion is the rapid pitch-up of the airfoil in rotor blades and wind turbines. Compared with static stall, dynamic stall is a more complicated unsteady phenomenon due to the flow separation caused by the rapid change in the angle of attack and the variable stall frequency. Thus, it leads to an abrupt decay of aerodynamic forces and a fluctuation of the structural load on the airfoil. McCroskey et al. [1] comprehensively studied dynamic stall in an oscillating airfoil and dynamic vortex shedding from the airfoil. Further research [2][3][4][5] has since provided a detailed understanding of the characteristics of dynamic stall with respect to the lift and pitch moment coefficient. Carr et al. [6] reviewed the significant progress in dynamic stall and concluded that the key influencing factors of this process were oscillation amplitude and frequency, Mach number, and the three-dimensional effects of the airfoil.
Dynamic stall is common in practical engineering and detrimental to the performance of air vehicles. Many researchers employed vortex generators to realize effective control of the dynamic stall [7,8]. Vortex generator control is a passive control method, which is easy-implemented and does not require extra energy injection. However, passive control cannot adapt to the change in working conditions. Recently, many researchers focused on active flow control, which is adjustable based on a feedback control loop [9][10][11][12]. Therefore, real-time sensing of the complicated and unsteady flow, e.g., the dynamic stall, is critical in the control loop.
The use of pressure tapes for local measurements is considered a promising method to investigate the flow conditions on an airfoil due to their reliability and ease of implementation. To reduce the complexity of the algorithm or system, some researchers applied some conventional and easy-implemented methods to identify the aerodynamic parameters of the airfoils. Therefore, it is important to predict the flow dynamic based on scattered pressure data. Patel et al. [13] used a pressure sensor array on the upper surface of an airfoil to detect the onset of flow separation. Saini et al. [14] predicted the aerodynamic parameters (i.e., inflow velocity and angle of attack) using discrete surface pressures measured at ports in the vicinity of the leading edge of an airfoil. Provost et al. [15] used a time-shift linear model to predict the roll moment using a sparse set of surface pressure measurements. Juliano et al. [16] found that the pressure field on an oscillating airfoil can have distinct patterns at different locations that are difficult to capture using pressure sensors. Similarly, An et al. [17] develop a method for estimating the instantaneous lift coefficient on a rapidly pitching airfoil that uses a small number of pressure sensors and a measurement of the angle of attack. Gao et al. [18] detected unsteady boundary layer transition on a pitching airfoil using a statistical criterion calculated from thirty surface-pressure transducer measurements. With the development of artificial intelligence, the machine learning method has been successfully applied to sense or predict the aerodynamic parameters [19,20]. However, machine learning methods, such as a blacked box, cannot provide an interpretable physical explanation. Therefore, some researchers employed data-driven techniques to realize the system or parameters identification [21][22][23]. Zhou et al. [24] explained the physical relationship of pressure data and then developed an offlineonline method to sense the real-time aerodynamic parameters. Saini et al. [25] developed the leading-edge flow sensing technique, which uses a few pressures in the airfoil leadingedge to identify and sense various flow events associated with vortex shedding in unsteady airfoils.
Most of the researches mentioned above focused on the flow field around static objects. The associated and available methods typically combine proper orthogonal decomposition (POD) with other algorithms, such as linear stochastic estimation [26] and particle filter [27]. However, it is difficult to predict the evolution of the aerodynamics of dynamic airfoil over time. In this regard, the dynamical system concept should be taken into consideration, and dynamic mode decomposition (DMD) can provide some inspiration, which is an equation-free and data-driven method to model dynamical systems. DMD was initially proposed by Schmid [28] to decompose the flow field into both spatially and temporally coherent structures and has the capability of predicting dynamic systems in the temporal domain. Berger et al. [29] successfully applied DMD to control robotic movement parameters by predicting human-robot interactions based on initial-state input data. However, it is difficult for standard DMD to make predictions based on scattered measurements because the method relies on an infinite-dimensional linear dynamic approximation for a finite-dimensional nonlinear dynamic system. It, therefore, may not provide an accurate result or may even fail when applied to highly nonlinear systems with lowdimensional measurement data [30]. To solve this problem, some data-driven methods have been proposed. The theory of DMD with delay embedding was proposed by Brunton et al. [31] to decompose a three-dimensional chaotic Lorenz system into a linear model in the delay coordinates using a single variable. Mohammad et al. [32] utilized DMD in time delay coordinates to capture the dynamics of a compressible signal containing oscillation perfectly. In addition, Susuki et al. [33] conducted this method to realize the prediction of the voltage dynamics of the rudimentary model directly from the single-bus measurement. The time delay embedding [34] has been investigated in detail on the structure and conditioning of linear models of nonlinear dynamics. Yuan et al. [35] developed DMD with a delay embedding method to predict high temporal and spatial resolved flow fields based on local particle image velocimetry (PIV) measurement. Arbabi et al. [36] applied the delay embedding method to sparse measurements and realized predictive control of nonlinear flow based on a Koopman model framework.
In the current work, the scattered data obtained from pressure tapes on the surface of an airfoil were used to predict the surface pressure of the airfoil during the oscillation process. Pressure-sensitive paint (PSP) surface measurement provided a qualitative visualization to verify the dynamic stall on the oscillating airfoil. Based on the scattered pressure tape measurements, DMD with time-delay embedding was applied to reconstruct the spatiotemporal pressure distribution on the oscillating airfoil and predict the future data in real-time.

Measurement Setup
The experiment was conducted in an open-circuit low-speed wind tunnel at the China Aerodynamics Research and Development Center. The test section of the wind tunnel has a cross-section of 1.8 m × 1.4 m. In the experiment, the freestream velocity was U∞ = 34 m/s. The two-dimensional airfoil model OA309, made of high-strength carbon fiber composite materials, with a chord length of c = 0.4 m and span length of s = 1.8 m, was mounted on a pitching testing platform. The freestream Reynolds number based on the chord length c was Re = 0.9 × 10 6 . The airfoil model was subjected to an unsteady periodic motion. As shown in Figure 1, the pitching motion mechanism was composed of a servo motor, planetary reducer, and encoder. The mechanism allowed control over motion parameters, such as mean angle of attack α0, pitching amplitude α1, and oscillation frequency f. In the experiment, α0 was set as 10°, α1 was set as 10°, and f was varied as 1, 2, 3, and 4 Hz. The reduced frequency was defined as k = πfc/U∞, with values of 0.037, 0.074, 0.111, and 0.148. The formula for the angle of attack α can be expressed as follows: The PSP measurement system and pressure tapes are shown in Figure 1. Part of the upper surface of the airfoil was painted with PSP to evaluate the changes in pressure distribution when the airfoil was oscillating. The PSP was excited by a UV LED (UHP-T-385-EP from Prizmatix Ltd., Givat Shmuel, Israel), and images were captured by a camera (PCO 1600) with a filter. Using a synchronizer (BNC-575) and potentiometer, PSP images were captured during certain phases of the oscillating airfoil. During the PSP measurement process, exposure signals from the camera were acquired simultaneously with measurements from the potentiometer, which contained the angle of attack of the airfoil. Using these two signals, every PSP frame was assigned to a certain angle of attack. Then, the phase-locked pressure distribution of the upper surface of the airfoil could be obtained by PSP measurement. The pressure tape measurements and PSP measurements were taken simultaneously. Forty dynamic differential pressure sensors (8510B, ENDVECO) were arranged chordwise at the midspan of the airfoil model. The measurement range of these pressure sensors is 1 psi, and the measurement inaccuracy is ±0.05% of full scale. Locations of these 40 pressure sensors are detailed in Table 1. These pressure sensors adopted a dynamic sampling frequency to sample 256 data points per oscillation cycle of the dynamic airfoil. The acquisition time was approximately 30 s. The aerodynamic parameters, such as lift and pitching moments, were obtained by the integration of the pressure data. Aerodynamic loads can be calculated based on the integration of the surface pressure coefficient. The pressure coefficients can be calculated from the following expression: where ρ is the air density, Pi is the static pressure obtained from a single pressure tap, P∞ is the static pressure of the incoming flow. The detailed process of the integration can refer to previous material [37].

Dynamic Mode Decomposition with Time-Delay Embedding
DMD is an efficient method of decomposing a dynamical system into specific spatial modes that evolve at a certain frequency and growth-decay rate in time. In DMD, the measurement data obtained by sensors at time j are arranged into a column vector xj = (x1,j, x2,j ,…, xn,j), where xj is a data snapshot and n denotes the number of measurement points (in this case, the number of sensors). All m snapshots were sampled with a constant time interval (sampling time step) Δt. DMD describes a discrete-time system as: where matrix A reflects the dynamical characteristics of the system. The first step of DMD is to arrange all m snapshots into two matrices, X1 = (x1, x2, … , xm−1) and X2 = (x2, x3, … , xm). Then, the linear approximation A may be written in terms of matrices based on Eq. 4 as: Using matrix A, known as the Koopman operator, future data can be predicted from current data. Thus, the core of DMD is to search for the optimal linear approximation A of the nonlinear dynamical system. The crucial step of the DMD algorithm is obtaining the rank-reduced representation of A in terms of a POD-projected matrix , which can be written as follows: where U ∈ n×r , Σ ∈ r×r , and V ∈ m×r are obtained by singular value decomposition (SVD).
Here, r is defined as the truncated number and represents the rank of approximation to X1 using SVD reduction. The next step is to compute the eigen-decomposition of , which is similar to A: where the columns of W are eigenvectors of and Λ is a diagonal matrix containing the eigenvalues λk, which capture the time dynamics of the corresponding DMD modes. These eigenvalues and eigenvectors of can be related back to the similarity-transformed eigenvalues and eigenvectors of A to reconstruct the DMD modes given by: where is the matrix whose columns each represent a DMD mode. Then, a snapshot at any time can be reconstructed from the DMD modes and eigenvalues: where b = Φ −1 × 1 contains the initial values of DMD modes. As shown above, the DMD method relies on the linear approximation of a nonlinear dynamic system. When the subdomains are of limited spatial dimensionality, it is difficult for DMD to provide accurate predictions. In this study, the direct use of discrete sensor measurements as the observables creates a system that is far removed from the infinite-dimensional linear approximation of highly nonlinear systems. To solve this problem, extended DMD [38] and kernel DMD [39] have previously been proposed. These methods use new observables to provide a set of coordinates in which the dynamical systems appear linear. However, a universal and systematic approach needs to be found. The details of time delay embedding methods can refer to the works of Brunton et al. [31]. Time-delay embedding augments the limited spatial observables by embedding future measurements into the current measurements. Therefore, more information can be captured to construct the linear Koopman operator. A requirement of time-delay embedding is to have access to at least nd + 1 consecutive snapshots, where nd is the number of delay steps. The time-delay-embedded matrices X1 and X2 can be transformed into H1 and H2 as follows: where is a new snapshot after delay embedding, and occupies the last n rows . The dimensions of the matrices after delay embedding H1 and H2, are ( × ) × ( − ). When nd = 1, the matrix after delay embedding is the same as the original matrix.
Analogously to Eq. 4, the linear mapping of H1 and H2 can be written as follows: Then, the time-delay-embedded matrices H1 and H2 are used to build the DMD model. The matrix B is the Koopman operator and reflects the dynamical characteristics. Therefore, standard DMD can be applied to find the corresponding eigenvalues and eigenvectors of B. Using the matrix B as follows: The current new snapshot at time + − 1 can be used to predict the future new snapshot , which contains the data at time + .

Flow Visualization
To examine the two-dimensional pressure field on the oscillating airfoil, the images from PSP measurement were phase-averaged over 15 periods. The evolution of the pressure field during the oscillation of the clean airfoil is shown in Figure 2 at two typical AoAs under an oscillation frequency of f = 1 Hz. It should be noted that these PSP measurements only provided a qualitative visualization because the PSP data was influenced by the temperature of the incoming flow. The upper panels represent the downstroke, and the lower panels represent the upstroke. The low-pressure region is indicated in blue, and high pressures are shown in yellow. At α = 9°, the pressure distribution in the downstroke is similar to that in the upstroke. However, when α increases to 13°, an obvious low-pressure region forms at the leading edge of the airfoil in the upstroke, as can be seen in Figure 2d. The color intensity shows that the pressure at the leading edge in the downstroke is higher than that in the upstroke at the same angle of attack, which suggests the onset of flow separation and the decrease in lift. The phenomenon that flow separation occurs at the leading edge in the downstroke but does not occur in the upstroke is consistent with the dynamic stall was also investigated by Heine et al. using PIV measurement [8].

Aerodynamic Loads
The angle of attack of the airfoil was varied to examine the influence of the VGs mounted at the leading edge of the airfoil on the static aerodynamic loads. In our study, the aerodynamic loads of static and dynamical airfoil were measured to verify the phenomenon of dynamic stall. The measurement results of static airfoil are shown in Figure  3. Note that there is no oscillation frequency for this static case. As shown in Figure 3, the curves of the lift and pitch moment show a relatively constant trend until α ≈ 13.5°. With a further increase in the angle of attack, the lift and pitch moment drop abruptly. An increase of 0.5° in the angle of attack leads to a drop of approximately 0.39 in the lift coefficient and 0.08 in the pitch moment coefficient. Therefore, α = 13.5° can be considered as the stall angle with respect to flow separation. These results were in line with our common perception and are consistent with the findings of previous studies [8,40]. The influence of the oscillation frequency f of the dynamical airfoil on the aerodynamic loads was examined in this study. The lift and moment coefficients under different AoA for the oscillating airfoil at different f are shown in Figure 4. It should be reminded that the relevant aerodynamic load data were phase-averaged from 10 measurement periods, and the resolution was improved by interpolation. In the present work, instantaneous pressure data at different periods are similar, so phase-averaged data based on 10 periods is credible. When f = 1 Hz, the lift coefficient of the upstroke is higher than that of the downstroke by up to 0.8 under the same angle of attack. There was a crucial characteristic feature in the curve of loop hysteresis, which was related to a certain α = 6° when f = 1 Hz. When α < 6°, the lift coefficients of the upstroke is lower than that of the downstroke under the same AoA. Once α > 6°, the lift coefficients of the upstroke is higher than that of the downstroke. Therefore, we can say that there are two hysteresis loops during the oscillating process. It can be assumed that the pressure of the lower surface of airfoil can be higher under small AoA due to the downstroke of the airfoil. When increasing the oscillation frequency to 2 Hz (Figure 4b), the certain AoA which decided the two-loop hysteresis became smaller. Continuing to improve the oscillation frequency f to 3 Hz (Figure 4c), only one loop hysteresis is left. However, when α= 0,1,2°, the lift moments under upstroke and downstroke were close. When the oscillation frequency f was increased to 4 Hz (Figure 4d), the lift coefficients under even a small angle of attack during the downstroke were smaller than that during the upstroke. With increasing oscillation frequency f, the loop hysteresis becames more pronounced, and the shapes of the hysteretic loops changed. In detail, the increase of oscillation frequency f induced the slight increase of upstroke and downstroke in the lift coefficient under larger AoAs, but the lift coefficient under smaller AoAs in the downstroke decreased slightly.
In addition to the lift coefficients, the characteristics of pitch moment were also investigated, as shown in Figure 5. For f = 1 Hz, there is no abrupt drop in the moment coefficient such as the static case during the upstroke. The moment coefficient of both upstroke and downstroke was very close, i.e., the difference of moment coefficient between upstroke and downstroke was small. It is worth noting that compared to the static case, the moment coefficient under large AoAs decreased dramatically, which suggested that dynamic stall influenced the moment coefficient. For example, when α = 20°, the moment coefficient in the static case was about −0.055 but it decreased to about −0.085 in the dynamical case. With increasing oscillation frequency f, the curve of moment coefficient formed two-loop hysteresis. The moment coefficient of upstroke and downstroke showed a huge difference, suggesting that high oscillation frequency leads to a large fluctuation in pitch moment. Combined with the change in lift coefficient, it can be inferred that not only did the dynamic stall exert a negative influence on the aerodynamic load, the oscillation frequency f also did. More importantly, high oscillation frequency f makes the aerodynamic performance more intricate (i.e., two-loop hysteresis in a coefficient curve). This problem is detrimental to the performance of real-time active flow control. Therefore, realtime perception, even the prediction of aerodynamic parameters, is challenging and meaningful. In Section Ⅲ.B, we will use a novel data-driven method to reconstruct the time-resolved pressure data and then realize the prediction of future pressure.

Prediction of Future Data
In this section, DMD with time delay embedding was used to reconstruct the history of pressure data of the oscillating airfoil. The high time-resolved pressure data collected by 20 sensors on the upper surface (see Table 1) were used in this method because the leading edge is the region with the most obvious pressure fluctuations. To examine the performance of this method, the error rate e between reconstructed data and true data is defined as: where m is the number of used snapshots, , is the reconstructed pressure field at time t, and , is the true velocity field on the subdomain at time t. POD analysis was applied to the time delay embedded dataset, and the mode energy spectrum is shown in Figure 6. The core of time delay embedding is the improvement of data dimensionality to capture more phase information. Therefore, the first three POD modes contained more than 98% energy. With the increase of the number of delay steps nd, more POD modes were needed to reconstruct the data accurately. It should be noted that when the number of delay steps was increased to 200, the POD energy spectrum tended to be similar. In this experiment, one period contains 256 snapshots. Therefore, when nd > 200, most of the phase information was included in the time delay embedded matrix. The raw and reconstructed pressure data of sensor #1 were used to compare the reconstruction performance of the standard DMD and DMD with delay embedding. Sensor #1 is located at the leading edge, where the most obvious pressure fluctuations happen. The standard DMD (i.e., nd = 0) and DMD with time-delay embedding were then applied to reconstruct the whole-time sequence of the pressure data. As shown in Figure 7a, standard DMD cannot accurately predict the dynamic pressure on sensor #1 even at the initial stage of the oscillation cycle. With time-delay embedding, the prediction performance is clearly improved. With nd = 50 delay steps, DMD can capture the periodic trend, but with an underprediction of the oscillation amplitude and e ≈ 43%. With nd increased to 200, the predicted values agree closely with the ground truth both in periodic trend and amplitude for more than six cycles and e ≈ 7%. Thus, DMD with time-delay embedding demonstrates the ability to reconstruct the periodic single-point pressure data over a long period. With an increasing number of delay steps nd, the accuracy of the reconstruction is improved. However, an excessive number of delay steps will lead to too large a dimensionality of the matrix, which will not only consume computing and storage resources but negatively affect the real-time performance of the prediction algorithm. Therefore, the duration of the prediction cycle should be balanced with the requirement for computation time. In fact, the accurate prediction of just one cycle, that is, the next cycle, is usually helpful for effective decision-making in real-time control systems. Figure 7. Reconstruction of pressure data measured by sensor #1 (see Table 1 Based on the predicted pressures at single measurement points, DMD with time-delay embedding was then used to reconstruct the time-resolved pressure field on the upper surface of the airfoil. The ground truth of the pressure field obtained by 20 pressure sensors is shown in Figure 8a. The values reconstructed using DMD with time-delay embedding are shown in Figure 8b and match well with the measured data and accurately show the shape and peak of the three-dimensional surface. The core of DMD is to find the Koopman operator, which allows the prediction of future data based on past data. The above section has demonstrated the accurate reconstruction of time-resolved pressure data using DMD with time-delay embedding, from which an optimal Koopman operator was obtained. Using the Koopman operator, realtime prediction can be realized. In Figure 9, the black line represents the ground-truth pressure data measured by sensor #1 over one oscillation period. The blue markers represent the past sampled pressure data that were used to build a new snapshot y . The future new snapshot was obtained based on the Koopman operator B provided by DMD with time-delay embedding, which had been verified in the above cases. Then, the future pressure data at the position of sensor #1 could be extracted from . Furthermore, future snapshot can also be calculated with the help of a predicted snapshot in the same way. When this operation is repeated, the data for a long period of time in the future can be predicted. In Figure 9, 256 snapshots (a whole oscillation period) were predicted. The predicted next-step value accurately matches the ground-truth value, indicating the feasibility and potential of DMD with time delay embedding. The computation was implemented in MATLAB, and all the cases are run on an Intel Xeon E5-1620 at 3.50 GHz. Every prediction took 0.09s. Admittedly, it is difficult to meet real-time requirements, but high-performance industrial computers can realize real-time prediction in practical applications. Work is ongoing to further improve the rapidity and adaptability of this method

Conclusions
In this study, dynamic stall of oscillating airfoils was measured by dynamic pressure tapes and PSP visualization. The influence of the oscillation frequency was examined. Dynamic mode decomposition (DMD) with time-delay embedding was proposed to predict the pressure field on the oscillating airfoil based on scattered measurements. The major findings are as follows: (1) Under the dynamic conditions of the oscillating airfoil, with increasing oscillation frequency, although lift coefficient can be improved at large AoAs, the aerodynamic characteristics showed a huge difference between upstroke and downstroke. The specific performance is that the loop hysteresis of both the lift and pitch moment coefficients were enlarged, which is detrimental to aerodynamic performance.
(2) DMD with time-delay embedding was able to reconstruct and predict dynamic stall based on scattered measurements with much higher accuracy than standard DMD. The reconstruction accuracy of this method increased with the number of delay steps, but this also prolonged the computation time. With the help of the Koopman operator obtained by DMD with time-delay embedding, the future dynamic pressure for a long time of an oscillating airfoil can be accurately predicted.