Unsteady aerodynamic effects in small-amplitude pitch oscillations of an airfoil

TL;DR

This study uses high-fidelity simulations to explore how small-amplitude pitch oscillations affect the unsteady aerodynamics of an airfoil, revealing complex flow behaviors influenced by laminar separation bubbles and transition points.

Contribution

It provides new insights into the flow physics of oscillating airfoils in sensitive regions, highlighting the impact of laminar separation bubbles on unsteady aerodynamics.

Findings

Large variations in transition points during pitch cycle

Dynamic changes in laminar separation bubble stability

Abrupt transition location shifts affecting flow response

Abstract

High-fidelity wall-resolved large-eddy simulations (LES) are utilized to investigate the flow-physics of small-amplitude pitch oscillations of an airfoil at Re = 100,000. 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 "off-design" conditions. The dynamic range of the pitch oscillations is chosen to be in this sensitive region. Large variations of the transition point on the suction-side of the airfoil are observed throughout the pitch cycle resulting in a dynamically rich flow response. Changes in the stability characteristics of a leading-edge 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 relaxation-term 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 Re = 400,000.

Equations

$$\frac{\partial u}{\partial t}+u\cdot \nabla u=-\frac{1}{\rho }\nabla p+\frac{1}{Re}{\nabla }^{2}u-\chi \mathcal{H}(u),$$

$$\nabla \cdot u=0,$$

$$\alpha ={\alpha }_0+\Delta \alpha \sin (\Omega t).$$

$$\displaystyle\bigg{[}\big{(}\frac{\partial}{\partial t}+U\frac{\partial}{\partial x}\big{)}\nabla^{2}-U^{\prime\prime}\frac{\partial}{\partial x}-\frac{1}{Re}\nabla^{4}\bigg{]}v=0,$$

$$v=\tilde{v}(y)e^{i(k_{x}x-\lambda t)},$$

$$\displaystyle\bigg{[}(-i\lambda+ik_{x}U)(\mathcal{D}^{2}-k_{x}^{2})-ik_{x}U^{\prime\prime}-\frac{1}{Re}(\mathcal{D}^{2}-k_{x}^{2})^{2}\bigg{]}\tilde{v}=0.$$

$$\overline{q}(x,y,t)=\left(\frac{1}{z_{max}-z_{min}}\right)\left(\frac{1}{\Delta t}\right){\int }_{t^{\prime }=t}^{t^{\prime }=t+\Delta t}{\int }_{z=z_{min}}^{z=z_{max}}q(x,y,z,t^{\prime })dzd{t}^{\prime }.$$

Full text

Keywords - unsteady aerodynamics, dynamic-response, transition, wall-resolved 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 small-amplitude 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 boundary-layer 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 laminar-turbulent 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 boundary-layer 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 trailing-edge 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 non-linear dynamic response to simple harmonic pitch motions (Mai and Hebler, 2011; Hebler et al., 2013). Such a non-linear 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 non-linear behavior of the normal force coefficient. These non-linearities occur only for oscillations within a certain range of angle of attack ($\alpha$) 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 leading-edge (Mai and Hebler, 2011; Lokatt, 2017). =While the above experiments were performed at a Reynolds number range of $O(10^{6})$, simulations and experiments at lower transitional Reynolds numbers on a NACA 0012 airfoil have also been shown to exhibit non-linear 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 boundary-layer 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 boundary-layer characteristics are observed. The analysis of the unsteady laminar separation bubble which develops near the leading-edge 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 small-amplitude 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 chord-based Reynolds number of $Re_{c}=100,000$. 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 ($C_{m}$) and suction-side transition location (figure 2) above an angle of attack $\alpha>6^{\circ}$. The steep slope of the coefficient of moment curve indicates a region where aerodynamic forces are sensitive to small changes in $\alpha$. 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 high-fidelity 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, wall-resolved large-eddy 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 (Eitel-Amor 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 flow-physics 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 open-source incompressible Navier–Stokes solver developed by Fischer et al. (2008) at Argonne National Laboratory. It is a based on a spectral-element method which allows the mapping of elements to complex geometries along with a high-order 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 $P_{N}$-$P_{N-2}$ formulation (Maday and Patera, 1989). An $11^{th}$ order polynomial approximation is used for the velocity with a $9^{th}$ order approximation for pressure. The nonlinear terms are treated explicitly by third-order extrapolation (EXT3), whereas the viscous terms are treated implicitly by a third-order backward differentiation scheme (BDF3). Aliasing errors are removed with the use of over-integration. The Arbitrary-Lagrangian-Eulerian (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 non-dimensional units with the velocities normalized by the reference free-stream velocity $U_{0}$ and the length scales in all directions are normalized by the chord length $c$. The resultant non-dimensional time unit is given by $c/U_{0}$.

2.2 Relaxation-term large-eddy simulation (RT-LES)

The LES method is based on the RT3D variant of the ADM-RT approach first used by Schlatter et al. (2004). The method supplements the governing equations with a dissipative term ($\chi\mathcal{H}(u)$). The equations for the resolved velocity and pressure thus read as

[TABLE]

where $\mathcal{H}$ is a high-pass spectral filter and $\chi$ 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 (Eitel-Amor 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 $Re_{\tau}=395$, 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 $\Delta x^{+}=18$ and the spanwise resolution was $\Delta z^{+}=9$. The first point in the wall-normal direction was set at $\Delta y_{w}^{+}=0.64$ and the wall-normal resolution near the boundary layer edge was $y_{max}^{+}=11$. 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 3(a), and for the turbulent kinetic energy (TKE) budget in figure 3(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 boundary-layer over a flat-plate by Eitel-Amor et al. (2014) where the ADM-RT model is used. Similarly, LES cases of pipe flows at $Re_{\tau}\approx 1000$ 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. Wall-shear stress data is obtained using XFOIL to estimate the grid spacing on the airfoil. A trip is introduced in XFOIL at $x/c\approx 0.1$ to obtain turbulent wall-shear values on both the suction and pressure sides of the airfoil. Here $c$ denotes the chord length. The grid design around the airfoil uses the following criteria:

• $\bullet$

For $0.1<x/c<0.6$, $\Delta x^{+}=18$, $\Delta y_{wall}^{+}=0.64$ and $\Delta y_{max}^{+}=11$, using the local wall-shear ($\tau_{w}$) values on the airfoil. Since the flow is expected to be laminar on the pressure side, the stream-wise resolution is slightly relaxed to $\Delta x^{+}=25$ while keeping the same wall-normal resolution.

• $\bullet$

For $x/c<0.1$, the peak $\tau_{w}$ value over the suction side of the airfoil is used to estimate the grid spacing.

• $\bullet$

For $x/c>0.6$, the suction side experiences a large adverse pressure gradient which significantly reduces $\tau_{w}$ values. Therefore, the $\tau_{w}$ values from the pressure side are used for both the suction and pressure sides.

• $\bullet$

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 $\Delta z^{+}=9$, where the the peak $\tau_{w}$ value from the suction side is used.

A different criterion is needed for defining the resolution in the wake where the wall-based criteria are not valid. Accordingly, Reynolds-averaged Navier–Stokes (RANS) simulations were performed using the transition $k$-$\Omega$ SST model (Langtry and Menter, 2009) with ANSYS® FLUENT, to estimate the Kolmogorov length scale ($\eta$) 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: $\Delta x/\eta<9$. 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 $l_{z}=0.25$ chords. The domain can be visualized in figure 5(a) and a close-up view of the spectral-elements is shown in (figure 5(b)). Each of the spectral-elements are further discretized by $12\times 12\times 12$ grid points in 3D, corresponding to an $11^{th}$ order spectral discretization. Periodic boundary conditions are imposed on the spanwise boundaries, while the energy-stabilized 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, free-stream turbulence of intensity $Ti=0.1\%$ is superimposed on the Dirichlet boundary conditions. The free-stream 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 free-stream turbulence. The same method has also been used for generating grid-turbulence in simulations of wind-turbines (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 $Re_{c}=400,000$ 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 wall-normal profiles of the normalized kinetic energy budget is shown in figure 4. The profiles are extracted at a streamwise location of $x/c=0.7$ 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 $Re_{c}=100,000$ at two different angles of attack ($\alpha=6.7^{\circ}$ and $\alpha=8.0^{\circ}$). As observed in figure 6, the iso-contours of coherent structures, identified by negative $\lambda_{2}$ method (Jeong and Hussain, 1995), show a substantial change in transition location for a small $\Delta\alpha=1.3^{\circ}$. For $\alpha=6.7^{\circ}$ the strong pressure gradient effects near the trailing edge cause transition at $x/c\approx 0.7$. While for $\alpha=8.0^{\circ}$, a leading-edge laminar separation bubble forms, causing flow transition much closer to the leading edge at $x/c\approx 0.2$. Such a leading-edge laminar separation bubble is not observed for the $\alpha=6.7^{\circ}$ case. The results are consistent with the trends obtained from XFOIL calculations, showing a large variation in the transition point within a small $\alpha$ 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 $\alpha_{0}=6.7^{\circ}$ with an amplitude of $\Delta\alpha=1.3^{\circ}$. The reduced frequency of oscillation is $k=0.5$ and the pitch axis is located at $(x_{0},y_{0})=(0.35,0.034)$. The reduced frequency is defined as $k=\Omega b/U_{0}$, where $\Omega$ is the angular frequency of oscillation, $b$ is the semi-chord length and $U_{0}$ is the free-stream velocity. The motion of the airfoil is prescribed by equation 2. The pitching motion corresponds to an oscillation time period of $T_{osc}=2\pi$.

[TABLE]

The time variation of the coefficient of lift ($C_{L}$) is shown in figure 7. The initial phase of pitching motion is carried out using a lower resolution (polynomial order $N=5$) to simulate the initial transient period of the flow at a lower computational cost. The polynomial order is then smoothly raised to $N=11$ 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 $C_{L}$ shows a smooth increase during the pitch-up motion, with strong secondary effects occurring near the maxima of the pitch cycles. Similarly in the pitch-down phase the lift decreases smoothly with secondary effects again becoming important at the minima of the pitch-cycles.

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 added-mass and circulatory components. While the circulatory component lags the instantaneous angle of attack, the added-mass component exhibits a phase gain (for pure pitch oscillations). At higher reduced frequencies the added-mass 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 space-time evolution of the instantaneous spanwise averaged wall-shear stress. The space-time surface plot is shown in figure 8(a), which spans the fourth and fifth pitch cycles. The color specifies the value of wall-shear 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 space-time surface is shown again as a binary colored surface plot in figure 8(b), where black colored regions indicate negative wall-shear stress and hence separated flow, while the white region corresponds to locations with attached flow ($\tau_{w}>0$).

It is obvious from the two plots in figure 8, that the developing boundary layer on the airfoil exhibits a dynamically rich response to small-amplitude 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 boundary-layer characteristics to paint an over-all picture of the dynamics. A persistent trailing-edge separation can be identified in figure 8(b) beyond $x/c>0.8$. The trailing-edge separation does not exhibit reverse flow $100\%$ 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 $x/c\approx 0.6$ at times $t/T_{osc}\approx 3$ and $t/T_{osc}\approx 4$. This is identified as a trailing-edge LSB. This LSB is short lived in time, existing for slightly less than a quarter of the pitch-cycle. A large separated region near the leading edge is a leading-edge LSB, which persists much longer in time, spanning nearly half a pitch cycle. Evident from figure 8(a) is that the transition point changes substantially throughout the pitch cycle. Interestingly, the flow over the airfoil differs significantly during the pitch-up and the pitch-down phases for the same angle of attack. For example, when the instantaneous angle of attack is at phase $\phi=0\ (t/T_{osc}=3,\ 4)$, which represents mean angle of attack but in the pitch-up phase, the flow over the airfoil is mostly laminar up to $x/c\approx 0.7$. On the other hand, for a phase of $\phi=\pi\ (t/T_{osc}=3.5,\ 4.5)$, representing the airfoil at the mean angle of attack but in the pitch-down phase, the flow is almost entirely turbulent with the start of the turbulent region approximately at $x/c\approx 0.22$.

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 space-time plot of the wall-shear stress. The empirical transition locations are consistent with the picture of wall-shear 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 iso-contours of instantaneous vortical structures, identified by the $\lambda_{2}$ method (Jeong and Hussain, 1995), at four different times during the pitch-up cycle when the transition is moving upstream. This phase is marked by the appearance of a leading-edge LSB which grows in size. The top figure shows the flow state near the mean angle of attack ($t/T_{osc}=3.09$) during the pitch-up 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 leading-edge 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 $t/T_{osc}=3.09,\ 3.2,\ 3.3$ and $3.47$, which corresponds to instantaneous angles of $\alpha=7.4^{o},7.94^{o},7.94^{o}$ and $6.94^{o}$ respectively. Note that pitch-up phase completes at $t/T_{osc}=3.25$ when the airfoil is at the highest angle of attack. Thus upstream movement of transition starts nearly at the end of the pitch-up phase and continues to move upstream even during the pitch-down cycle. The laminar separation bubble close to the trailing edge ceases to exist as transition moves upstream. At $t/T_{osc}=3.47$ the flow transition is seen to occur on the separated shear layer of the leading-edge LSB (bottom left in figure 10).

As the airfoil progresses through the pitch-down cycle the leading-edge LSB shrinks in size and the transition point then starts moving downstream, initiating the process of flow re-laminarization (figure 10 right, top to bottom). The leading-edge LSB eventually disappears, although the transition point moves downstream of the LSB before it completely disappears. The flow over the airfoil is not completely re-laminarized even when the airfoil is at the lowest angle of attack and the re-laminarization process continues into the pitch-up 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 9(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 $V^{tr}_{u}=-0.60$ while the velocity of the downstream motion of transition is calculated across the green circles and is equal to $V^{tr}_{d}=0.17$. Thus the upstream spread of turbulent regions is much faster than flow re-laminarization.

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 ($t/T_{osc}=3.33$) 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 $40\%$ and the entire boundary layer downstream is turbulent. On the right figure ($t/T_{osc}=3.35$) there is a similar large connected region of turbulence, spreading from approximately $40\%$ of the chord. However there is also a nascent region of turbulence starting at $x/c\approx 0.2$. 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 $t/T_{osc}=3.4$ 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 leading-edge 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 leading-edge 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 Orr-Sommerfeld equation

[TABLE]

where $v$ is the wall-normal component of the velocity fluctuations and $U^{\prime\prime}$ is the second derivative in the wall-normal direction of the streamwise velocity $U$. Due to Squire’s theorem, analyzing the two-dimensional perturbations is sufficient to study the modal stability characteristics. Following Schmid and Henningson (2001) and assuming an asantz function for the wall-normal perturbations as

[TABLE]

results in a dispersion relation between the complex frequency $\lambda$ and the streamwise wavenumber $k_{x}$ given by

[TABLE]

Here $\mathcal{D}$ represents the derivative in the wall-normal ($y$) direction. While strictly valid for homogeneous flows, the Orr-Sommerfeld 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 $k_{x}$ results in an eigenvalue problem for the frequency $\lambda$. 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 Orr-Sommerfeld equations is performed using the instantaneous wall-normal 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 wall-normal profiles of tangential velocity which exhibit the maximum reverse flow intensity relative to the boundary-layer edge velocity. The edge velocity of the local profiles was determined using the criterion of vanishing spanwise vorticity i.e. $\overline{\omega}_{z}\approx 0$. Since there is a very small but finite amount of vorticity in the far-field due to the free-stream turbulence, the criterion for boundary layer height is set as the wall-normal distance where the vorticity decays to $1\%$ of its maximum value in the boundary layer. Figure 12 shows the wall-normal 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 $k_{x}$) for instantaneous profiles at several different time instants in the fourth pitch cycle. The reverse-flow 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 reverse-flow 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 reverse-flow 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, $C_{g}$. 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 self-sustaining process of perturbation growth. In the context of LSBs, high reverse-flow intensity has been associated with the presence of an absolute instability. Alam and Sandham (2000) with their local stability analysis of a two-parameter family of reverse flow profiles indicated that reverse flow intensities above $15\%$ 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 $20\%$ 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 $30\%$ reverse flow velocities. Figure 15 shows the absolute value of the maximum reverse flow intensity in the leading-edge LSB found in the present study. The reverse flow velocities are found to be higher than $30\%$ 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 leading-edge 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 reverse-flow 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 reverse-flow 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 $\lambda$, 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 “pinch-points” 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 saddle-point on the frequency plane gives the absolute frequency $\lambda^{0}$. 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 cusp-map 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 $\lambda$ rather than the non-linear eigenvalue problem for $k_{x}$ 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 $k_{x}$. Figure 16 and 17 show the cusps obtained in the complex frequency plane for the profiles with the highest reverse-flow intensity of tangential velocity. In the fourth cycle at $t/T_{osc}=3.32$ an unstable cusp is obtained with a positive absolute amplification rate of $\lambda^{0}_{i}=2.46$. On the other hand the cusp found in the fifth cycle is just marginally stable with an absolute amplification rate of $\lambda^{0}_{i}=-0.5$ at $t/T_{osc}=4.4$. Since the velocity profiles considered here are averaged values instead of stationary base-flow profiles, it is possible that small temporal variation and/or non-linear modification of the flow may hide the absolute instability properties when an averaged flow-field 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 chord-wise locations at $t/T_{osc}=3.32$. 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 Orr-Sommerfeld 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 time-scales to judge the validity of the such an analysis. The ratio of the time-period of the absolute instability ($4^{th}$ cycle) to the time period of oscillation is $0.05$, 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 $\delta^{max}/L_{x}=0.03$ which suggests a weak spatial inhomogeneity. Here $\delta^{max}$ is the maximum value of the boundary layer height over the LSB and $L_{x}$ is the spatial extent of the LSB. Both the above quantities indicate the quasi-steady, 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 reverse-flow moving to the end of the bubble. At high reverse-flow 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 boundary-layer characteristics.

As mentioned earlier, due to conditions of weak spatial inhomogeneity and slow time-variation, 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 self-sustaining 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 ($\sim 12\%$) which ruled out the possibility of absolute instabilities in those cases. In contrast in the present case the LSB exhibits reverse flow intensities of $30-40\%$ which makes it substantially more likely for an absolute instability mechanism to be active.

4 Conclusion

A relaxation-term filtering procedure is used for wall-resolved LES of flow over a pitching airfoil. Validation of the LES procedure is done in a channel flow at $Re_{\tau}=395$ and for a wing section at $Re_{c}=400,000$ 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 $Re_{c}=100,000$ 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 $\alpha$.

Pitch oscillations of the airfoil within this $\alpha$ 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 suction-side 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 trailing-edge. Flow transition occurs over this separated shear layer. As the angle of attack increases, a leading-edge 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 pitch-down cycle, the reverse phenomenon occurs where the leading-edge 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 $V^{tr}_{u}\approx-0.60$ and a much slower downstream velocity of $V^{tr}_{d}\approx 0.17$. 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 $12^{th}$ PRACE Project Access Call (number 2015133182). The work was partially funded by European Research Council under grant agreement 694452-TRANSEP-ERC-2015-AdG. The work was also partially funded by Vinnova through the NFFP project UMTAPS, with grant number 2014-00933. 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 ($\Delta t$) in time. This averaged value is considered as the “instantaneous” value for that quantity. Thus this instantaneous value of any statistical quantity $\overline{q}(x,y,t)$ can be evaluated as in equation 6:

[TABLE]

Here $(z_{max}-z_{min})$ 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 $\Delta t=4\times 10^{-3}$, which amounts to $0.06\%$ 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, $\overline{u^{\prime}v^{\prime}}(x,y,t)$, in order to determine the instantaneous transition location. The first streamwise location (on the suction-side) 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 $|\overline{u^{\prime}v^{\prime}}|$ 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 $|\overline{u^{\prime}v^{\prime}}|_{max}\approx 0.05$. The threshold for determining transition is set to $5\%$ of this value. The transition point is thus the first point where $|\overline{u^{\prime}v^{\prime}}|>Tol_{uv}=0.0025$. Since we use an ad-hoc criterion for transition location, this is cross-checked by evaluating the variance of the spanwise velocity fluctuations $\overline{w^{\prime}w^{\prime}}$, and following an identical procedure as described above. In this case the transition criterion is prescribed as $|\overline{w^{\prime}w^{\prime}}|>Tol_{ww}=0.005$, since the peak spanwise fluctuation intensities are found to be nearly twice the peak Reynolds stress $|\overline{u^{\prime}v^{\prime}}|$. Physically, growing spanwise velocity fluctuations indicate the onset of three-dimensionality in the boundary layer, which are associated with secondary instabilities. While the thresholds specified may be considered ad-hoc, 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.