Unsteady aerodynamic effects in smallamplitude pitch oscillations of an airfoil
Abstract
Highfidelity wallresolved largeeddy simulations (LES) are utilized to investigate the flowphysics of smallamplitude pitch oscillations of an airfoil at . The investigation of the unsteady phenomenon is done in the context of natural laminar flow airfoils, which can display sensitive dependence of the aerodynamic forces on the angle of attack in certain “offdesign” conditions. The dynamic range of the pitch oscillations is chosen to be in this sensitive region. Large variations of the transition point on the suctionside of the airfoil are observed throughout the pitch cycle resulting in a dynamically rich flow response. Changes in the stability characteristics of a leadingedge laminar separation bubble has a dominating influence on the boundary layer dynamics and causes an abrupt change in the transition location over the airfoil. The LES procedure is based on a relaxationterm which models the dissipation of the smallest unresolved scales. The validation of the procedure is provided for channel flows and for a stationary wing at .
Keywords  unsteady aerodynamics, dynamicresponse, transition, wallresolved les, laminar separation bubble, local stability
1 Introduction
A large focus of the research on unsteady wings tends towards stall dynamics such as the earlier works of McCroskey (1981); McCroskey et al. (1982); McCroskey (1973); McCroskey et al. (1976); Carr et al. (1977); Ericsson and Reding (1988a, b), etc. More recent works by Dunne and McKeon (2015); Rival and Tropea (2010); Choudhry et al. (2014); Visbal (2011, 2014); Visbal and Garmann (2017); Alferez et al. (2013); Rosti et al. (2016) etc. continue the investigation which appears far from complete. The review by McCroskey (1982) and a more recent one by Coorke and Thomas (2015) provide an overview of the dynamic stall process. Studies on the aerodynamic behavior of smallamplitude pitch oscillations are typically done from the perspective of aeroelasticity where investigations tend to focus on the time varying aerodynamic forces on the airfoil with much less attention paid to the boundarylayer characteristics. Some studies focusing on the time dependent boundary layer in small pitch amplitudes include the work done by Pascazio et al. (1996) which shows a time delay in laminarturbulent transition during pitching. Nati et al. (2015) analyzed the effect of small amplitude pitching on a laminar separation bubble at low Reynolds numbers. Mai and Hebler (2011) and Hebler et al. (2013) investigate the boundarylayer changes in an oscillating natural laminar flow airfoil in transonic conditions.
Such cases qualitatively represent small changes in operating conditions, such as the changes due to structural deformations or small trailingedge flap deflections. The understanding of flow response to such changes can be crucial in cases where small perturbations induce large variations in aerodynamic forces. Such sensitive dependence of aerodynamic forces may be found in the static characteristics of natural laminar flow (NLF) airfoils for a certain range of angle of attack. Their performance critically depends on maintaining laminar flow over the suction side of the airfoil and a loss of laminar flow over the airfoil causes large variations of the aerodynamic forces. In addition to such sensitive dependence of the aerodynamic characteristics, recent transonic unsteady experiments using NLF airfoils have brought to light a peculiar property of these airfoils. The unsteady aerodynamic coefficients for laminar airfoils exhibit a nonlinear dynamic response to simple harmonic pitch motions (Mai and Hebler, 2011; Hebler et al., 2013). Such a nonlinear response is inconsistent with the predictions of classical unsteady aerodynamic models (Theodorsen, 1935) which are typically based on inviscid assumptions. Similar experiments within the subsonic range have been performed by Lokatt (2017) who also found strongly nonlinear behavior of the normal force coefficient. These nonlinearities occur only for oscillations within a certain range of angle of attack () and have been strongly linked to the free movement of transition over the suction side of the airfoil. They seem to be nearly absent when suction side transition is fixed at the leadingedge (Mai and Hebler, 2011; Lokatt, 2017). =While the above experiments were performed at a Reynolds number range of , simulations and experiments at lower transitional Reynolds numbers on a NACA 0012 airfoil have also been shown to exhibit nonlinear unsteady aerodynamic characteristics (Poirel et al., 2008; Poirel and Yuan, 2010; Barnes and Visbal, 2016). At these transitional Reynolds numbers, the laminar separation bubble has been shown to play a key role in the unsteady dynamics (Yuan et al., 2013; Barnes and Visbal, 2018). In all the above cases, the changes in boundarylayer characteristics lead to dynamic phenomena which can not be modeled using linear theories. With these effects in mind, the current work seeks to investigate the flow dynamics of an unsteady airfoil within the aerodynamic regime where large changes in the boundarylayer characteristics are observed. The analysis of the unsteady laminar separation bubble which develops near the leadingedge is done using local linear stability analysis in order to explain some of the observed unsteady dynamics.
The present work investigates the flow physics of smallamplitude pitch oscillations of a laminar airfoil (figure 1). The airfoil was designed at the Aeronautical and Vehicle Engineering department of KTH, where it has been used in previous experimental and numerical works (Lokatt and Eller, 2017). The same airfoil was used in the unsteady experiments of Lokatt (2017). The simulations are performed at a chordbased Reynolds number of . The angle of attack range for the oscillation was chosen from the static characteristics of the airfoil. The static characteristics were calculated using an integral boundary layer code XFOIL (Drela, 1989), which predicted sharp changes in the coefficient of moment () and suctionside transition location (figure 2) above an angle of attack . The steep slope of the coefficient of moment curve indicates a region where aerodynamic forces are sensitive to small changes in . It is important to note that the results obtained from XFOIL are only used to approximately identify the above mentioned region of sensitive dependence for the airfoil at a low computational cost. Simulations with a static airfoil are performed to ensure that such sensitive dependence on angle of attack is indeed observed with the highfidelity LES. The exact dynamic range of the pitching motion is prescribed based on the results of the static simulations (described in section 3).
In recent works, wallresolved largeeddy simulations have proven to be an effective tool for studying flow physics at high Reynolds numbers with a computational cost which is much lower than that of direct numerical simulations (DNS). Some of the works to utilize this method include spatially evolving boundary layers (EitelAmor et al., 2014), turbulent channel flows (Schlatter et al., 2006a, b), pipe flows (Chin et al., 2015) and flow over wings (Uzun and Hussaini, 2010; Lombard et al., 2016). Successful application of the approach has motivated its use in the present work, which aims to gain insight into the flowphysics of unsteady airfoils undergoing small amplitude pitch oscillations.
The remainder of the paper is divided into 3 sections. Section 2 describes the numerical setup and presents the results of the validation of the LES. Results of both the stationary and pitching simulations are discussed in Section 3. The conclusions of the study are presented in Section 4.
2 Computational setup
2.1 Numerical method
The computational code used for the simulations is Nek5000, which is an opensource incompressible Navier–Stokes solver developed by Fischer et al. (2008) at Argonne National Laboratory. It is a based on a spectralelement method which allows the mapping of elements to complex geometries along with a highorder spatial discretization within the elements. The method uses Lagrange interpolants as basis functions and utilizes Gauss–Lobatto–Legendre (GLL) quadrature for the distribution of points within the elements. The spatial discretization is done by means of the Galerkin approximation, following the  formulation (Maday and Patera, 1989). An order polynomial approximation is used for the velocity with a order approximation for pressure. The nonlinear terms are treated explicitly by thirdorder extrapolation (EXT3), whereas the viscous terms are treated implicitly by a thirdorder backward differentiation scheme (BDF3). Aliasing errors are removed with the use of overintegration. The ArbitraryLagrangianEulerian (ALE) formulation is used to account for the mesh deformation due to the motion of the pitching airfoil (Ho and Patera, 1990, 1991). All equations are solved in nondimensional units with the velocities normalized by the reference freestream velocity and the length scales in all directions are normalized by the chord length . The resultant nondimensional time unit is given by .
2.2 Relaxationterm largeeddy simulation (RTLES)
The LES method is based on the RT3D variant of the ADMRT approach first used by Schlatter et al. (2004). The method supplements the governing equations with a dissipative term (). The equations for the resolved velocity and pressure thus read as
(1a)  
(1b) 
where is a highpass spectral filter and is a model parameter. Together the two parameters determine the strength of the dissipative term. The method has been used in earlier studies of spatially developing boundary layers (EitelAmor et al., 2014) and channel flows (Schlatter et al., 2006b). The method has been shown to be reliable in predicting transition location and also preserving the characteristic structures which are seen in the DNS of transitional flows (Schlatter et al., 2006b).
A number of tests were carried out in a channel flow at a friction Reynolds number of , and the results are compared with the DNS data of Moser et al. (1999). The final mesh was set up such that the streamwise resolution was and the spanwise resolution was . The first point in the wallnormal direction was set at and the wallnormal resolution near the boundary layer edge was . The superscript indicates normalization in inner units. A comparison of the results for the turbulent channel flow is shown for the mean velocity in figure 2(a), and for the turbulent kinetic energy (TKE) budget in figure 2(b). The dissipation profile shown in the figure is the sum of resolved dissipation and the added dissipation by the relaxation term. A good agreement with the DNS is found for the mean velocity and all the kinetic energy budget terms (including the total dissipation). The resolution used this study is much finer than the typical LES and closer to the coarse DNS resolutions used in the studies of turbulent flows. A very similar resolution is used in the simulation of spatially developing boundarylayer over a flatplate by EitelAmor et al. (2014) where the ADMRT model is used. Similarly, LES cases of pipe flows at by Chin et al. (2015) uses slightly coarser resolutions than the one used in the present study.
2.3 Mesh generation
The optimum mesh resolution (in inner units) obtained in the channel flow results is then used to design the mesh around the airfoil. Wallshear stress data is obtained using XFOIL to estimate the grid spacing on the airfoil. A trip is introduced in XFOIL at to obtain turbulent wallshear values on both the suction and pressure sides of the airfoil. Here denotes the chord length. The grid design around the airfoil uses the following criteria:

For , , and , using the local wallshear () values on the airfoil. Since the flow is expected to be laminar on the pressure side, the streamwise resolution is slightly relaxed to while keeping the same wallnormal resolution.

For , the peak value over the suction side of the airfoil is used to estimate the grid spacing.

For , the suction side experiences a large adverse pressure gradient which significantly reduces values. Therefore, the values from the pressure side are used for both the suction and pressure sides.

A structured mesh is used, which is extruded in the spanwise direction. Hence the spanwise resolution is constant throughout the domain. The resolution is set to , where the the peak value from the suction side is used.
A different criterion is needed for defining the resolution in the wake where the wallbased criteria are not valid. Accordingly, Reynoldsaveraged Navier–Stokes (RANS) simulations were performed using the transition  SST model (Langtry and Menter, 2009) with ANSYS^{®} FLUENT, to estimate the Kolmogorov length scale () in the wake region. The RANS is setup with domain boundaries at a distance of 100 chords from the airfoil. The grid in the wake region for the LES is designed such that the average grid spacing between the GLL points follows the criteria: . The computational domain is set up such that the far field boundaries of the computational domain are two chords away from the airfoil leading edge in either direction. The outflow boundary is four chords downstream from the airfoil leading edge and the inlet is designed as a curved inflow boundary with a constant radial distance of two chords from the leading edge of the airfoil. The spanwise width of the domain is chords. The domain can be visualized in figure 4(a) and a closeup view of the spectralelements is shown in (figure 4(b)). Each of the spectralelements are further discretized by grid points in 3D, corresponding to an order spectral discretization. Periodic boundary conditions are imposed on the spanwise boundaries, while the energystabilized outflow condition suggested by Dong et al. (2014) is imposed on the outflow boundary. Velocity field data for locations corresponding to the boundaries of the LES computational domain is extracted from an unsteady RANS simulation. The extracted data is imposed as a Dirichlet boundary condition on these boundaries for the LES. The method is very similar to the one used by Hosseini et al. (2016) in their DNS of flow around a wing section. In order to simulate low turbulence flight conditions, freestream turbulence of intensity is superimposed on the Dirichlet boundary conditions. The freestream turbulence is generated using Fourier modes with a von Kármán spectrum. The procedure is similar to the one described in Schlatter (2001); Brandt et al. (2004) and Schlatter et al. (2008) and has been used for the study of transition in flat plate boundary layers under the influence of freestream turbulence. The same method has also been used for generating gridturbulence in simulations of windturbines (Kleusberg, 2017).
A validation of the above methodology for complex geometries such as a wing section was performed at a chord based Reynolds number of for NACA4412 airfoil. The LES grid resolution was setup with the same grid criteria as described above. The domain boundaries and boundary conditions are identical to the setup in Hosseini et al. (2016). A comparison of the wallnormal profiles of the normalized kinetic energy budget is shown in figure 4. The profiles are extracted at a streamwise location of on the suction side of the airfoil. The LES profiles (lines) match well with the DNS data (circles) of Hosseini et al. (2016), signifying the high accuracy of the LES with the current resolution.
3 Results and discussion
3.1 Steady results
Simulations with a stationary airfoil were performed to investigate the location of transition without pitching motion. The simulations were performed for at two different angles of attack ( and ). As observed in figure 6, the isocontours of coherent structures, identified by negative method (Jeong and Hussain, 1995), show a substantial change in transition location for a small . For the strong pressure gradient effects near the trailing edge cause transition at . While for , a leadingedge laminar separation bubble forms, causing flow transition much closer to the leading edge at . Such a leadingedge laminar separation bubble is not observed for the case. The results are consistent with the trends obtained from XFOIL calculations, showing a large variation in the transition point within a small change (figure 2).
3.2 Unsteady boundary layer characteristics
Once the static characteristics of the airfoil are obtained, the dynamic effects on the boundary layer are investigated by pitching the airfoil about a mean angle with an amplitude of . The reduced frequency of oscillation is and the pitch axis is located at . The reduced frequency is defined as , where is the angular frequency of oscillation, is the semichord length and is the freestream velocity. The motion of the airfoil is prescribed by equation 2. The pitching motion corresponds to an oscillation time period of .
(2) 
The time variation of the coefficient of lift () is shown in figure 7. The initial phase of pitching motion is carried out using a lower resolution (polynomial order ) to simulate the initial transient period of the flow at a lower computational cost. The polynomial order is then smoothly raised to before the fourth pitch cycle. Due to the fairly large separation at the trailing edge, effects of transition movement and turbulence, successive pitch cycles are not expected to have identical behavior, however some of qualitatively repeating trends can be observed. The behavior of the lift coefficient shows a chaotic but qualitatively repeating pattern where shows a smooth increase during the pitchup motion, with strong secondary effects occurring near the maxima of the pitch cycles. Similarly in the pitchdown phase the lift decreases smoothly with secondary effects again becoming important at the minima of the pitchcycles.
Note that the lift appears to be marginally ahead in phase as compared to the instantaneous angle of attack. This is in fact a feature of pitch oscillations at slightly high reduced frequencies. The linear flutter theory by Theodorsen (1935) divides the unsteady response of a pitching airfoil into addedmass and circulatory components. While the circulatory component lags the instantaneous angle of attack, the addedmass component exhibits a phase gain (for pure pitch oscillations). At higher reduced frequencies the addedmass components dominate the unsteady response and thus lead to a net phase gain in the lift response. Variation of the theoretical unsteady lift phase with reduced frequencies as well as comparisons with experimental data from Halfman (1952) and Rainey (1957) can be found in Leishman (2000).
In order to understand the time variation of the spatially developing boundary layer on the airfoil, we analyze the spacetime evolution of the instantaneous spanwise averaged wallshear stress. The spacetime surface plot is shown in figure 7(a), which spans the fourth and fifth pitch cycles. The color specifies the value of wallshear stress on the suction side of the airfoil. Regions with color intensity strongly towards red are indicative of high shear and thus turbulent flow. The exception to the rule being the region close to the leading edge where the flow is laminar but a high shear region exists due to the extremely thin boundary layer close to the stagnation point. The same spacetime surface is shown again as a binary colored surface plot in figure 7(b), where black colored regions indicate negative wallshear stress and hence separated flow, while the white region corresponds to locations with attached flow ().
It is obvious from the two plots in figure 8, that the developing boundary layer on the airfoil exhibits a dynamically rich response to smallamplitude pitch oscillations, with different key boundary layer characteristics controlling the dynamics of the flow in different phases of the pitch cycle. We identify some of the key boundarylayer characteristics to paint an overall picture of the dynamics. A persistent trailingedge separation can be identified in figure 7(b) beyond . The trailingedge separation does not exhibit reverse flow of the time, as can be seen from the white patches dispersed between largely black colored regions. An isolated separated region (distinct from the trailing edge separation) is observed at at times and . This is identified as a trailingedge LSB. This LSB is short lived in time, existing for slightly less than a quarter of the pitchcycle. A large separated region near the leading edge is a leadingedge LSB, which persists much longer in time, spanning nearly half a pitch cycle. Evident from figure 7(a) is that the transition point changes substantially throughout the pitch cycle. Interestingly, the flow over the airfoil differs significantly during the pitchup and the pitchdown phases for the same angle of attack. For example, when the instantaneous angle of attack is at phase , which represents mean angle of attack but in the pitchup phase, the flow over the airfoil is mostly laminar up to . On the other hand, for a phase of , representing the airfoil at the mean angle of attack but in the pitchdown phase, the flow is almost entirely turbulent with the start of the turbulent region approximately at .
3.3 Transition location
In the present work we focus on the variation of transition location throughout the pitch cycles. Since the flow case is unsteady, and the transition location does not remain fixed, a criteria based on the instantaneous state of the flow is needed to determine the transition location. The details of the evaluation of the transition location can be found in A. The temporal variation and the phase portrait of the calculated transition location are shown in figure 9. The magenta line in figure 8 shows the calculated transition locations superposed on the spacetime plot of the wallshear stress. The empirical transition locations are consistent with the picture of wallshear stress, with transition marginally preceding regions of turbulent flow.
The superposed plots in figure 8 indicate that the two LSBs play a dominating role in the flow dynamics and that transition location is governed by the characteristics of these LSBs. Figure 10 (left) shows the isocontours of instantaneous vortical structures, identified by the method (Jeong and Hussain, 1995), at four different times during the pitchup cycle when the transition is moving upstream. This phase is marked by the appearance of a leadingedge LSB which grows in size. The top figure shows the flow state near the mean angle of attack () during the pitchup phase. The flow is mostly laminar across the airfoil with no structures observed prior to flow transition (close to trailing edge) and there is no leadingedge LSB. The high adverse pressure gradient near the trailing edge causes the laminar flow to easily separate, forming a LSB and flow transitions over this separated shear layer. Figure 10 (left) from top to bottom shows the part of the oscillation cycle when transition is moving upstream at time instants of and , which corresponds to instantaneous angles of and respectively. Note that pitchup phase completes at when the airfoil is at the highest angle of attack. Thus upstream movement of transition starts nearly at the end of the pitchup phase and continues to move upstream even during the pitchdown cycle. The laminar separation bubble close to the trailing edge ceases to exist as transition moves upstream. At the flow transition is seen to occur on the separated shear layer of the leadingedge LSB (bottom left in figure 10).
As the airfoil progresses through the pitchdown cycle the leadingedge LSB shrinks in size and the transition point then starts moving downstream, initiating the process of flow relaminarization (figure 10 right, top to bottom). The leadingedge LSB eventually disappears, although the transition point moves downstream of the LSB before it completely disappears. The flow over the airfoil is not completely relaminarized even when the airfoil is at the lowest angle of attack and the relaminarization process continues into the pitchup phase. There is a marked asymmetry between the upstream and downstream movement of the transition point. An approximate velocity for both the upstream and downstream motion of the transition point is calculated across the points marked by circles in figure 8(a) which correspond to transition movement with near constant velocity. The velocity of upstream transition movement is calculated across the black circles and is equal to while the velocity of the downstream motion of transition is calculated across the green circles and is equal to . Thus the upstream spread of turbulent regions is much faster than flow relaminarization.
Specifying a velocity of transition movement implies that the transition location changes smoothly with time. This is true for the downstream movement of transition, however the final stages of the upstream movement appear to be more complex. Figure 11 shows the instantaneous vortical structures at two time instants during the upstream movement phase. In the left figure () a single connected region of turbulence can be observed which is preceded by a laminar region identified by the near absence of small vortical structures. This region starts at about and the entire boundary layer downstream is turbulent. On the right figure () there is a similar large connected region of turbulence, spreading from approximately of the chord. However there is also a nascent region of turbulence starting at . These two regions are separated by a patch of laminar flow. The figure on the right then has two distinct locations where transition to turbulence takes place. After a short while, turbulence spreads downstream from this newly formed transition location and eventually the entire boundary layer downstream is turbulent. This flow state, where two distinct turbulent regions can be identified exists only for a short duration and by no intermediate regions of laminar flow can be identified. However the emergence of two distinct transition locations indicates that at least two competing mechanisms exist for the growth of disturbances in the boundary layer and that their relative strength changes during the pitch cycle. Given that this new transition occurs at the separated shear layer of the leadingedge LSB, the stability properties of the LSB are likely to play a role in the emergence of this new transition point.
3.4 Stability characteristics of the laminar separation bubble
Stability characteristics of steady laminar separation bubbles have been studied with different perspectives by several previous authors like Hammond and Redekopp (1998); Alam and Sandham (2000); Marxen et al. (2003); Lang et al. (2004); Marxen and Rist (2010); Jones et al. (2008, 2010); Marxen et al. (2012, 2013) etc. In unsteady cases, the stability characteristics of the leadingedge LSB can help shed some light on the changing transition locations throughout the pitch cycle and in particular on the existence of competing mechanisms for transition. The competition between convective and absolute instability characteristics may provide a possible explanation for the transient existence of two distinct points of transition. The change of flow state from laminar to turbulent is often governed by the linear amplification of small disturbances within the boundary layer. For flows with streamwise and spanwise homogeneity the evolution of small perturbations within the boundary layer is governed by the OrrSommerfeld equation
(3) 
where is the wallnormal component of the velocity fluctuations and is the second derivative in the wallnormal direction of the streamwise velocity . Due to Squire’s theorem, analyzing the twodimensional perturbations is sufficient to study the modal stability characteristics. Following Schmid and Henningson (2001) and assuming an asantz function for the wallnormal perturbations as
(4) 
results in a dispersion relation between the complex frequency and the streamwise wavenumber given by
(5) 
Here represents the derivative in the wallnormal () direction. While strictly valid for homogeneous flows, the OrrSommerfeld equation has often been used for flows that exhibit weak inhomogeneity such as the Blasius boundary layer and also for the case of laminar separation bubbles (Alam and Sandham, 2000; Hammond and Redekopp, 1998; Häggmark et al., 2001); Boutilier and Yarusevych (2012); Marxen et al. (2012).
A temporal stability analysis using a real spatial wavenumber results in an eigenvalue problem for the frequency . Resulting complex frequencies with a positive imaginary part signify that the boundary layer is unstable and that small perturbations within the boundary layer would grow in time and cause transition.
To explore the time varying stability properties of the LSB, temporal stability analysis of the OrrSommerfeld equations is performed using the instantaneous wallnormal profiles of tangential velocity, calculated as per equation 6. Several velocity profiles can be considered for local analysis. Alam and Sandham (2000) and Hammond and Redekopp (1998) in their local analysis associate the stability characteristics of the LSB with the maximum reverse flow intensities. In accordance with the previous studies, we focus on the wallnormal profiles of tangential velocity which exhibit the maximum reverse flow intensity relative to the boundarylayer edge velocity. The edge velocity of the local profiles was determined using the criterion of vanishing spanwise vorticity i.e. . Since there is a very small but finite amount of vorticity in the farfield due to the freestream turbulence, the criterion for boundary layer height is set as the wallnormal distance where the vorticity decays to of its maximum value in the boundary layer. Figure 12 shows the wallnormal profiles of tangential velocity and spanwise vorticity along with the evaluated height of the boundary layer.
Figure 13 shows the unstable complex frequencies obtained from the temporal stability analysis (with varying ) for instantaneous profiles at several different time instants in the fourth pitch cycle. The reverseflow intensity continues to increase with time until flow transition occurs at the LSB. The flow is unstable for all the analyzed velocity profiles as shown by the existence of complex frequencies with a positive imaginary part. However the highest amplification rate (frequency with maximum imaginary part) does not monotonically increase with reverseflow intensity. At first the maximum amplification rate increases in time, but later it is seen to decrease. This changing characteristics can be associated with structural changes in the LSB, where at first the region of maximum reverseflow is found near the center of the LSB, but as the LSB grows, this region of strong reverse flow moves closer to the downstream end of the LSB. Such qualitative changes in the shape of the LSB can be observed in figure 14.
While the local stability analysis shows that the boundary layer is unstable, an additional distinction needs to be made with regards to the nature of the instability which may be classified as either convective or absolute (Briggs, 1964). The instability characteristics are usually elucidated using the simple concept of group velocity of perturbations, . Growing perturbations that travel with a positive group velocity are deemed convectively unstable since they move away from the source of disturbance. On the other hand, perturbations with zero group velocity are referred to as absolutely unstable, since they do not convect away from the origin of instability and create a selfsustaining process of perturbation growth. In the context of LSBs, high reverseflow intensity has been associated with the presence of an absolute instability. Alam and Sandham (2000) with their local stability analysis of a twoparameter family of reverse flow profiles indicated that reverse flow intensities above may cause the flow to be locally absolutely unstable. With a similar analysis on a three parameter family of profiles Hammond and Redekopp (1998) obtained onset of absolute instabilities at reverse flow velocities. In the same study the authors also performed global stability analysis on a synthetically created boundary layer with a symmetric separation bubble and found the flow to be globally unstable for reverse flow velocities. Figure 15 shows the absolute value of the maximum reverse flow intensity in the leadingedge LSB found in the present study. The reverse flow velocities are found to be higher than in both the fourth and fifth pitch cycles, which is higher than the thresholds indicated by earlier investigators for absolute instability.
In such a case it is likely that the leadingedge LSB changes character to become absolutely unstable during the pitch cycle. The change of characteristics may explain the simultaneous existence of two different transition points in the boundary layer. Initially when the reverseflow intensity in the LSB is small, the flow is convectively unstable and perturbations grow while traveling downstream. Thus transition occurs downstream of the LSB. As the reverseflow intensity becomes larger the region of absolute instability may exist within the LSB which would cause perturbations to grow in time without being convected away. When these perturbations become large enough they would cause transition over the LSB. Thus momentarily, there would be two distinct transition points, one due to the growth of absolute instabilities in the LSB and one due to convectively amplifying disturbances which cause transition downstream of the LSB.
To explore such a possibility, from the local stability analysis, one needs to identify unstable modes with zero group velocity. Briggs (1964) proposed a general method for the identification of absolute instabilities in a system which is commonly referred to as the Briggs method. The method involves solving the spatial stability problem for a range of complex , thus mapping contours on the complex frequency plane on to the complex wavenumber plane through the dispersion relation. Saddle points obtained in the complex wavenumber plane are the points where the dispersion relation has a double root. These points are known as the “pinchpoints” in the complex wavenumber plane where the branches corresponding to the forward and backward propagating solutions of the dispersion relation meet via a double root. These pinch points correspond to perturbation modes that have zero group velocity. The mapping of the saddlepoint on the frequency plane gives the absolute frequency . If this absolute frequency lies in the unstable half of the plane, then the system exhibits an absolute instability. Kupfer et al. (1987) proposed the inverse method called the cuspmap method, which maps contours in the complex wavenumber plane on to the complex frequency plane via the dispersion relation and identified the absolute instability by locating the cusp in the complex frequency plane. This allowed for solving the simpler linear eigenvalue problem for rather than the nonlinear eigenvalue problem for in the Briggs method (Briggs, 1964). Several contours need to be mapped from the wavenumber to the frequency plane for the location of the cusp. Kupfer et al. (1987) proposed mapping contours along lines with constant real part of . Figure 16 and 17 show the cusps obtained in the complex frequency plane for the profiles with the highest reverseflow intensity of tangential velocity. In the fourth cycle at an unstable cusp is obtained with a positive absolute amplification rate of . On the other hand the cusp found in the fifth cycle is just marginally stable with an absolute amplification rate of at . Since the velocity profiles considered here are averaged values instead of stationary baseflow profiles, it is possible that small temporal variation and/or nonlinear modification of the flow may hide the absolute instability properties when an averaged flowfield is considered. Nonetheless the results provide a strong indication that the LSB changes character to become absolutely unstable during the pitch cycle.
According to Chomaz et al. (1991) and Huerre and Monkewitz (1990), the existence of a finite region of absolute instability is a necessary criterion for flows to exhibit a global hydrodynamic instability. Local stability analysis is therefore performed at multiple chordwise locations at . Figure 18 shows the spatial variation of the absolute growth rate along with the contours of negative streamwise velocity in the LSB. Clearly a small finite region of absolute instability exists close to the rear end of the LSB, where the highest intensities of reverse flow exist.
Caution needs to be taken in the interpretation of the stability analysis. The local stability analysis using the OrrSommerfeld equations assumes stationary homogeneous base state. Neither of which are strictly fulfilled when using the instantaneous, spanwise averaged velocity profiles. One can compare the properties of the LSB and the instability timescales to judge the validity of the such an analysis. The ratio of the timeperiod of the absolute instability ( cycle) to the time period of oscillation is , suggesting that the boundary layer would appear nearly stationary to the amplifying disturbances. For the spatial inhomogeneity, the ratio of the maximum boundary layer height to the length of the separation bubble can be used as an indicator. This ratio is equal to which suggests a weak spatial inhomogeneity. Here is the maximum value of the boundary layer height over the LSB and is the spatial extent of the LSB. Both the above quantities indicate the quasisteady, homogeneous flow assumptions may be used to obtain qualitative features of the flow case. The analysis then suggests that the convective instability properties of the boundary layer initially become stronger as the LSB grows. The LSB also changes in shape, with regions of high reverseflow moving to the end of the bubble. At high reverseflow intensities the LSB changes character and exhibits a region of absolute instability. This can potentially explain the emergence of two distinct transition locations. The upstream transition would be caused by the temporally growing instabilities which amplify within the region of absolute instability without being convected downstream. On the other hand spatially growing waves associated with convective instabilities would trigger transition further downstream of the LSB. The emergence of the second transition point associated with absolute instabilities would cause abrupt changes in the boundarylayer characteristics.
As mentioned earlier, due to conditions of weak spatial inhomogeneity and slow timevariation, the analysis of local stability analysis remains indicative rather than conclusive. Global instabilities may be triggered by mechanisms different from those arising due to an absolute instability (Huerre and Monkewitz, 1990). Studies have claimed alternate mechanisms of global instability in laminar separation bubbles. Cherubini et al. (2010); Theofilis et al. (2000); Rodríguez and Theofilis (2010) propose topological changes in the separated flow region as the cause of global instability. Jones et al. (2008); Marxen et al. (2013) provide some evidence that secondary instabilities of the vortices shed from the LSB can lead to a selfsustaining transition. Jones et al. (2008) claim that a combination of hyperbolic instability in the braid region of vortices and elliptic instability of the vortex core causes the flow to behave like an oscillator. Marxen et al. (2013) also find evidence of elliptic and hyperbolic instabilities in their numerical study. In the current work we do not make a rigorous attempt to rule out all other mechanism of global instability. However some comments can be made on the possible existence of these other mechanisms in the current case. Topological changes, as proposed by Theofilis et al. (2000); Rodríguez and Theofilis (2010) or the one by Cherubini et al. (2010) cause the LSB to get divided into two regions. Velocity contours inside the LSB for our case (figure 14) show no such division of the LSB. Therefore such topological changes may not be the cause of transition in the present case. The flow cases studied by Jones et al. (2008); Marxen et al. (2013) rely on the secondary instability for sustained turbulence. However in both the cases the reverse flow intensities are small () which ruled out the possibility of absolute instabilities in those cases. In contrast in the present case the LSB exhibits reverse flow intensities of which makes it substantially more likely for an absolute instability mechanism to be active.
4 Conclusion
A relaxationterm filtering procedure is used for wallresolved LES of flow over a pitching airfoil. Validation of the LES procedure is done in a channel flow at and for a wing section at and the results show a good agreement with available DNS data sets.
Flow over an airfoil is simulated using the LES procedure at a chord based Reynolds number of within an angle of attack range where the aerodynamic forces on the airfoil exhibit sensitive dependence on the angle of attack. This sensitive dependence is captured in the steady simulations at different angles of attack with large changes in transition location within a small variation of .
Pitch oscillations of the airfoil within this range of sensitive dependence displays a rich variety of unsteady flow phenomena. The flow goes through alternating periods of fully turbulent and laminar flow over the suctionside of the airfoil with different governing mechanisms for transition through the oscillating phases. When the flow is mostly laminar over the airfoil surface it separates easily under the effect of adverse pressure gradient, forming an LSB near the trailingedge. Flow transition occurs over this separated shear layer. As the angle of attack increases, a leadingedge LSB is formed which first excites spatially growing waves (convective instability) causing transition downstream of the LSB. Initially the amplification rate of these spatially growing waves increases as the size of the LSB grows causing transition to move upstream. Eventually a region of absolute instability develops within the LSB and flow transition occurs abruptly on the separated shear layer. When transition is first triggered by this absolute instability mechanism the flow exhibits two distinct transition locations and abrupt changes in the boundary layer follow.
In the pitchdown cycle, the reverse phenomenon occurs where the leadingedge LSB shrinks in size and the region of absolute instability ceases to exit. The transition is then again governed by spatially amplifying waves. The spatial amplification rate now reduces as the LSB shrinks and transition moves downstream. The flow thus goes through states of convective and absolute instability, leading to constantly changing transition location. The upstream and downstream velocities of the transition point movement however are vastly different, with an average upstream velocity being around and a much slower downstream velocity of . This asymmetry is yet to be investigated, but may be an important parameter in unsteady turbulence modeling.
Acknowledgement
The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High Performance Computing at the Royal Institute of Technology (KTH). Simulation have also been performed at the Barcelona Supercomputing Center, Barcelona, with computer time provided by the PRACE Project Access Call (number 2015133182). The work was partially funded by European Research Council under grant agreement 694452TRANSEPERC2015AdG. The work was also partially funded by Vinnova through the NFFP project UMTAPS, with grant number 201400933. We would like to thank Dr. David Eller and Dr. Mikaela Lokatt for providing us with the NLF design and the numerous discussions on different aerodynamic aspects of the project.
Appendix A Empirical transition location
In order to quantify the transition location, we calculate statistical quantities which are averaged in the homogeneous spanwise direction, and also averaged for a short duration () in time. This averaged value is considered as the “instantaneous” value for that quantity. Thus this instantaneous value of any statistical quantity can be evaluated as in equation 6:
(6) 
Here is the spanwise extent of the computational domain. In order for such a quantity to be representative of the instantaneous state of the flow, the time duration of the averaging must be small. For the current case we use , which amounts to of the oscillation time period during which the flow can be assumed to remain in approximately the same state. Using this procedure we evaluate the fluctuating Reynolds stress, , in order to determine the instantaneous transition location. The first streamwise location (on the suctionside) where the fluctuating Reynolds stress becomes larger than a certain threshold is denoted as the transition point. In order to prescribe a suitable threshold, the maximum value of across the entire boundary layer is evaluated for all times. This maximum value does not exhibit large variations, remaining within the same order of magnitude for all times with its mean value being . The threshold for determining transition is set to of this value. The transition point is thus the first point where . Since we use an adhoc criterion for transition location, this is crosschecked by evaluating the variance of the spanwise velocity fluctuations , and following an identical procedure as described above. In this case the transition criterion is prescribed as , since the peak spanwise fluctuation intensities are found to be nearly twice the peak Reynolds stress . Physically, growing spanwise velocity fluctuations indicate the onset of threedimensionality in the boundary layer, which are associated with secondary instabilities. While the thresholds specified may be considered adhoc, the qualitative picture of transition movement is not very sensitive to small changes in the threshold. Changing the thresholds by a factor of 2 still produces the same qualitative trends as seen in figure 19. Moreover, the time variation of transition point determined by two different physical quantities agree very well with each other. Thus we consider the quantities as a good representations of instantaneous flow characteristics.
References
 Alam and Sandham (2000) Alam, M. and Sandham, N. D. (2000). Direct numerical simulation of ’short’ laminar separation bubbles with turbulent reattachment. Journal of Fluid Mechanics, 410:1–28.
 Alferez et al. (2013) Alferez, N., Mary, I., and Lamballais, E. (2013). Study of stall development around an airfoil by means of high fidelity large eddy simulation. Flow, Turbulence and Combustion, 91(3):623–641.
 Barnes and Visbal (2016) Barnes, C. J. and Visbal, M. (2016). Highfidelity les simulations of selfsustained pitching oscillations on a naca0012 airfoil at transitional reynolds numbers. In 54th AIAA Aerospace Sciences Meeting, AIAA SciTech Forum.
 Barnes and Visbal (2018) Barnes, C. J. and Visbal, M. R. (2018). On the role of flow transition in laminar separation flutter. Journal of Fluids and Structures, 77:213 – 230.
 Boutilier and Yarusevych (2012) Boutilier, M. S. H. and Yarusevych, S. (2012). Separated shear layer transition over an airfoil at a low reynolds number. Physics of Fluids, 24(8):084105.
 Brandt et al. (2004) Brandt, L., Schlatter, P., and Henningson, D. S. (2004). Transition in boundary layers subject to freestream turbulence. Journal of Fluid Mechanics, 517.
 Briggs (1964) Briggs, R. J. (1964). ElectronStream Interaction with Plasmas. The MIT Press, Cambridge, Massachusetts.
 Carr et al. (1977) Carr, L. W., McAlister, K. W., and McCroskey, W. J. (1977). Analysis of the development of dynamic stall based on oscillating airfoil experiments. Technical report, NASA Ames Research Center; Moffett Field, CA, United States.
 Cherubini et al. (2010) Cherubini, S., Robinet, J.C., and Palma, P. D. (2010). The effects of nonnormality and nonlinearity of the navier–stokes operator on the dynamics of a large laminar separation bubble. Physics of Fluids, 22(1):014102.
 Chin et al. (2015) Chin, C., Ng, H., Blackburn, H., Monty, J., and Ooi, A. (2015). Turbulent pipe flow at : A comparison of wallresolved largeeddy simulation, direct numerical simulation and hotwire experiment. Computers and Fluids, 122:26 – 33.
 Chomaz et al. (1991) Chomaz, J.M., Huerre, P., and Redekopp, L. G. (1991). A frequency selection criterion in spatially developing flows. Studies in Applied Mathematics, 84(2):119–144.
 Choudhry et al. (2014) Choudhry, A., Leknys, R., Arjomandi, M., and Kelso, R. (2014). An insight into the dynamic stall lift characteristics. Experimental Thermal and Fluid Science, 58:188 – 208.
 Coorke and Thomas (2015) Coorke, T. C. and Thomas, F. O. (2015). Dynamic stall in pitching airfoils: Aerodynamic damping and compressibility effects. Annual Review of Fluid Mechanics, 47(1):479–505.
 Dong et al. (2014) Dong, S., Karniadakis, G. E., and Chryssostomidis, C. (2014). A robust and accurate outflow boundary condition for incompressible flow simulations on severelytruncated unbounded domains. Journal of Computational Physics, 261:83–105.
 Drela (1989) Drela, M. (1989). XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils, pages 1–12. Springer Berlin Heidelberg, Berlin, Heidelberg.
 Dunne and McKeon (2015) Dunne, R. and McKeon, B. J. (2015). Dynamic stall on a pitching and surging airfoil. Experiments in Fluids, 56(8):157.
 EitelAmor et al. (2014) EitelAmor, G., Örlü R., and Schlatter, P. (2014). Simulation and validation of a spatially evolving turbulent boundary layer up to . International Journal of Heat and Fluid Flow, 47:57–69.
 Ericsson and Reding (1988a) Ericsson, L. and Reding, J. (1988a). Fluid mechanics of dynamic stall part i. unsteady flow concepts. Journal of Fluids and Structures, 2(1):1 – 33.
 Ericsson and Reding (1988b) Ericsson, L. and Reding, J. (1988b). Fluid mechanics of dynamic stall part ii. prediction of full scale characteristics. Journal of Fluids and Structures, 2(2):113 – 143.
 Fischer et al. (2008) Fischer, P. F., Lottes, J. W., and Kerkemeier, S. G. (2008). Nek5000 web page. http://nek5000.mcs.anl.gov.
 Häggmark et al. (2001) Häggmark, C. P., Hildings, C., and Henningson, D. S. (2001). A numerical and experimental study of a transitional separation bubble. Aerospace Science and Technology, 5(5):317–328.
 Halfman (1952) Halfman, R. L. (1952). Experimental aerodynamic derivatives of a sinusoidally oscillating airfoil in twodimensional flow. Technical report, National Advisory Committee for Aeronautics; Washington, DC, United States.
 Hammond and Redekopp (1998) Hammond, D. and Redekopp, L. (1998). Local and global instability properties of separation bubbles. European Journal of Mechanics  B/Fluids, 17(2):145 – 164.
 Hebler et al. (2013) Hebler, A., Schojda, L., and Mai, H. (2013). Experimental investigation of the aeroelastic behavior of a laminar airfoil in transonic flow. In Proceedings IFASD.
 Ho and Patera (1990) Ho, L.W. and Patera, A. T. (1990). A Legendre spectral element method for simulation of unsteady incompressible viscous freesurface flows. Computer Methods in Applied Mechanics and Engineering, 80(1):355 – 366.
 Ho and Patera (1991) Ho, L.W. and Patera, A. T. (1991). Variational formulation of threedimensional viscous freesurface flows: Natural imposition of surface tension boundary conditions. International Journal for Numerical Methods in Fluids, 13(6):691–698.
 Hosseini et al. (2016) Hosseini, S. M., Vinuesa, R., Schlatter, P., Hanifi, A., and Henningson, D. S. (2016). Direct numerical simulation of the flow around a wing section at moderate Reynolds number. International Journal of Heat and Fluid Flow, 61:117 – 128.
 Huerre and Monkewitz (1990) Huerre, P. and Monkewitz, P. A. (1990). Local and global instabilities in spatially developing flows. Annual Review of Fluid Mechanics, 22(1):473–537.
 Jeong and Hussain (1995) Jeong, J. and Hussain, F. (1995). On the identification of a vortex. Journal of Fluid Mechanics, 285.
 Jones et al. (2008) Jones, L. E., Sandberg, R. D., and Sandham, N. D. (2008). Direct numerical simulations of forced and unforced separation bubbles on an airfoil at incidence. Journal of Fluid Mechanics, 602:175–207.
 Jones et al. (2010) Jones, L. E., Sandberg, R. D., and Sandham, N. D. (2010). Stability and receptivity characteristics of a laminar separation bubble on an aerofoil. Journal of Fluid Mechanics, 648:257–296.
 Kleusberg (2017) Kleusberg, E. (2017). Wind turbine simulations using spectral elements. Licentiate thesis, Royal Institute of Technology (KTH), Stockholm, Sweden.
 Kupfer et al. (1987) Kupfer, K., Bers, A., and Ram, A. K. (1987). The cusp map in the complexfrequency plane for absolute instabilities. The Physics of Fluids, 30(10):3075–3082.
 Lang et al. (2004) Lang, M., Rist, U., and Wagner, S. (2004). Investigations on controlled transition development in a laminar separation bubble by means of lda and piv. Experiments in Fluids, 36(1):43–52.
 Langtry and Menter (2009) Langtry, R. B. and Menter, F. R. (2009). CorrelationBased Transition Modeling for Unstructured Parallelized Computational Fluid Dynamics Codes. AIAA Journal, 47:2894–2906.
 Leishman (2000) Leishman, J. G. (2000). Principles of Helicopter Aerodynamics. Cambridge University Press.
 Lokatt (2017) Lokatt, M. (2017). On Aerodynamic and Aeroelastic Modeling for Aircraft Design. Doctoral thesis, KTH Royal Institute of Technology.
 Lokatt and Eller (2017) Lokatt, M. and Eller, D. (2017). Robust viscousinviscid interaction scheme for application on unstructured meshes. Computers Fluids, 145:37 – 51.
 Lombard et al. (2016) Lombard, J.E. W., Moxey, D., Sherwin, S. J., Hoessler, J. F. A., Dhandapani, S., and Taylor, M. J. (2016). Implicit LargeEddy Simulation of a Wingtip Vortex. AIAA Journal, 54:506–518.
 Maday and Patera (1989) Maday, Y. and Patera, A. T. (1989). Spectral element methods for the incompressible NavierStokes equations. In Stateoftheart surveys on computational mechanics (A9047176 2164). New York, American Society of Mechanical Engineers, 1989, p. 71143. Research supported by DARPA., pages 71–143.
 Mai and Hebler (2011) Mai, H. and Hebler, A. (2011). Aeroelasticity of a laminar wing. In Proceedings IFASD, Paris.
 Marxen et al. (2012) Marxen, O., Lang, M., and Rist, U. (2012). Discrete linear local eigenmodes in a separating laminar boundary layer. Journal of Fluid Mechanics, 711:1–26.
 Marxen et al. (2013) Marxen, O., Lang, M., and Rist, U. (2013). Vortex formation and vortex breakup in a laminar separation bubble. Journal of Fluid Mechanics, 728:58–90.
 Marxen et al. (2003) Marxen, O., Lang, M., Rist, U., and Wagner, S. (2003). A combined experimental/numerical study of unsteady phenomena in a laminar separation bubble. Flow, Turbulence and Combustion, 71(1):133–146.
 Marxen and Rist (2010) Marxen, O. and Rist, U. (2010). Mean flow deformation in a laminar separation bubble: separation and stability characteristics. Journal of Fluid Mechanics, 660:37–54.
 McCroskey (1973) McCroskey, W. J. (1973). Inviscid flowfield of an unsteady airfoil. AIAA Journal, 11(8):1130 – 1137.
 McCroskey (1981) McCroskey, W. J. (1981). Phenomenon of dynamic stall. Technical report, NASA Ames Research Center; Moffett Field, CA, United States.
 McCroskey (1982) McCroskey, W. J. (1982). Unsteady airfoils. Annual Review of Fluid Mechanics, 14(1):285–311.
 McCroskey et al. (1976) McCroskey, W. J., Carr, L. W., and McAlister, K. W. (1976). Dynamic stall experiments on oscillating airfoils. AIAA Journal, 14(1):57 – 63.
 McCroskey et al. (1982) McCroskey, W. J., McAlister, K. W., Carr, L. W., and Pucci, S. L. (1982). An experimental study of dynamic stall on advanced airfoil sections. volume 1: Summary of the experiment. Technical report, NASA Ames Research Center, Moffett Field, CA, United States.
 Moser et al. (1999) Moser, R. D., Kim, J., and Mansour, N. N. (1999). Direct numerical simulation of turbulent channel flow up to . Physics of Fluids, 11(4):943–945.
 Nati et al. (2015) Nati, A., De Kat, R., Scarano, F., and Van Oudheusden, B. W. (2015). Dynamic pitching effect on a laminar separation bubble. Experiments in Fluids, 56(9):172.
 Pascazio et al. (1996) Pascazio, M., Autric, J., Favier, D., and Maresca, C. (1996). Unsteady boundarylayer measurement on oscillating airfoilstransition and separation phenomena in pitching motion. In 34th Aerospace Sciences Meeting and Exhibit, page 35.
 Poirel et al. (2008) Poirel, D., Harris, Y., and Benaissa, A. (2008). Selfsustained aeroelastic oscillations of a naca0012 airfoil at lowtomoderate reynolds numbers. Journal of Fluids and Structures, 24(5):700 – 719.
 Poirel and Yuan (2010) Poirel, D. and Yuan, W. (2010). Aerodynamics of laminar separation flutter at a transitional reynolds number. Journal of Fluids and Structures, 26(7):1174 – 1194.
 Rainey (1957) Rainey, A. G. (1957). Measurement of aerodynamic forces for various mean angles of attack on an airfoil oscillating in pitch and on two finitespan wings oscillating in bending with emphasis on damping in the stall. Technical report, National Advisory Committee for Aeronautics. Langley Aeronautical Lab.; Langley Field, VA, United States.
 Rival and Tropea (2010) Rival, D. and Tropea, C. (2010). Characteristics of pitching and plunging airfoils under dynamicstall conditions. Journal of Aircraft, 47(1):80–86.
 Rodríguez and Theofilis (2010) Rodríguez, D. and Theofilis, V. (2010). Structural changes of laminar separation bubbles induced by global linear instability. Journal of Fluid Mechanics, 655:280–305.
 Rosti et al. (2016) Rosti, M. E., Omidyeganeh, M., and Pinelli, A. (2016). Direct numerical simulation of the flow around an aerofoil in rampup motion. Physics of Fluids, 28(2):025106.
 Schlatter (2001) Schlatter, P. (2001). Direct numerical simulation of laminarturbulent transition in boundary layer subject to freestream turbulence. Diploma thesis, Royal Institute of Technology (KTH), Stockholm, Sweden.
 Schlatter et al. (2008) Schlatter, P., Brandt, L., de Lange, H. C., and Henningson, D. S. (2008). On streak breakdown in bypass transition. Physics of Fluids, 20(10):101505.
 Schlatter et al. (2004) Schlatter, P., Stolz, S., and Kleiser, L. (2004). LES of transitional flows using the approximate deconvolution model. International Journal of Heat and Fluid Flow, 25(3):549 – 558. Turbulence and Shear Flow Phenomena (TSFP3).
 Schlatter et al. (2006a) Schlatter, P., Stolz, S., and Kleiser, L. (2006a). Analysis of the SGS energy budget for deconvolution and relaxationbased models in channel flow, pages 135–142. Springer Netherlands, Dordrecht.
 Schlatter et al. (2006b) Schlatter, P., Stolz, S., and Kleiser, L. (2006b). Largeeddy simulation of spatial transition in plane channel flow. Journal of Turbulence, 7:N33.
 Schmid and Henningson (2001) Schmid, P. J. and Henningson, D. S. (2001). Stability and Transition in Shear Flows. Springer.
 Theodorsen (1935) Theodorsen, T. (1935). General theory of aerodynamic instability and the mechanism of flutter. Technical report, National Advisory Committee for Aeronautics; Langley Aeronautical Lab.; Langley Field, VA, United States.
 Theofilis et al. (2000) Theofilis, V., Hein, S., and Dallmann, U. (2000). On the origins of unsteadiness and threedimensionality in a laminar separation bubble. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 358(1777):3229–3246.
 Uzun and Hussaini (2010) Uzun, A. and Hussaini, M. Y. (2010). Simulations of vortex formation around a blunt wing tip. AIAA Journal, 48:1221–1234.
 Visbal (2011) Visbal, M. R. (2011). Numerical investigation of deep dynamic stall of a plunging airfoil. AIAA Journal, 49(10):2152 – 2170.
 Visbal (2014) Visbal, M. R. (2014). Analysis of the onset of dynamic stall using highfidelity largeeddy simulations. In 52nd Aerospace Sciences Meeting, AIAA SciTech Forum. AIAA.
 Visbal and Garmann (2017) Visbal, M. R. and Garmann, D. J. (2017). Numerical investigation of spanwise end effects on dynamic stall of a pitching naca 0012 wing. In 55th AIAA Aerospace Sciences Meeting, AIAA SciTech Forum. AIAA.
 Yuan et al. (2013) Yuan, W., Poirel, D., and Wang, B. (2013). Simulations of pitchheave limitcycle oscillations at a transitional reynolds number. AIAA Journal, 51:1716 – 1732.