# Vortex Shedding Lock-In due to Pitching Oscillation of a Wind Turbine Blade Section at High Angles of Attack

Authors: Craig Meskell; Alberto Pellegrino
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

The unsteady flow around a pitching two-dimensional airfoil section (NREL S809) has been simulated using unsteady RANS with the transition SST turbulence model. This geometry is chosen to represent a wind turbine blade in a standstill configuration. The Reynolds number isRe=106 based on a chord length of 1 m. A prescribed sinusoidal pitching motion has been applied at a fixed amplitude of 7° for a range of high angles of attack 30°<α<150°. At these incidences, the airfoil will behave more like a bluff body and may experience periodic vortex shedding. It is well known that, in bluff body flows, oscillations can lead to a lock-in (lock-in) of the vortex shedding frequency, fv, with the body’s motion frequency, fp. In order to investigate the susceptibility of airfoil to lock-in, the frequency ratio r (r=fp/fv0) has been varied around r=1. The lock-in region boundaries have been proposed, and an analysis of the effect of the oscillation amplitude has been conducted. The lock-in map obtained suggests that, for the vibration amplitude considered, the risk of vortex-induced vibration is more significant in the regions of α≈40° and α≈140°, i.e., for shallower characteristic lengths. Finally, a lumped parameter wake oscillator model has been proposed for pitching airfoils. This simple model is in qualitative agreement with the CFD results.


---

## Body

## 1. Introduction

### 1.1. Context

In 2012, there were a total of 284 GW of wind generating capacity installed globally [1]. The majority of the machines installed are horizontal axis design (HAWT). As with any technology that harvests ambient energy, wind turbines must be permanently exposed to environmental extremes. As a result, the minimum design requirement for wind turbines as defined by the IEC Standard [2] includes specifications for extreme wind loading conditions (see Section 6.3 of this IEC Standard). For example, the structure must be designed to withstand steady, turbulent, and deterministic gusts which occur on average only every 50 years (i.e., a 50-year return period). As part of the transient wind (i.e., gust) environment, a sudden change of wind direction of up to 180° must be considered. In the extreme change of direction scenario, the blade will be suddenly inclined at high angle of attack to the local flow. Hence, the section will behave more like a bluff body than a streamlined airfoil, and so periodic vortex shedding might be expected. To avoid the blades acting as bluff bodies, the rotor is always allowed to rotate in high winds, with wind shed from the blades so as to reduce the power production. The motion of the blades ensures a relatively low angle of incidence regardless of wind direction. However, in a standstill condition, due perhaps to mechanical failure of the transmission, the possibility high local angles of attack increases. Furthermore, as more turbines are installed offshore, the calendar window for maintenance is reduced and so the possibility of a wind turbine in standstill for extended periods increases. This paper will focus on this scenario: a large wind turbine in standstill at high angles of incidence. As there is no rotation, a 2D model is sufficient. The flexibility of the blade is modeled as a prescribed 2 degrees of motion.


### 1.2. Background to Vortex Shedding and Wake Oscillator Model

Vortex shedding induced vibration (VIV) can be a problem for a wide range of structures exposed to a cross flow. A major feature of VIV is “lock-in,” a resonance phenomenon that occurs when the vortex shedding frequency of an elastic bluff body is close to its structural frequency. At this condition, amplification of motion can lead to fatigue or even catastrophic failure. In particular, VIV can be expected in flexible slender bodies of any cross section, and therefore, this field is still an active area of research, in a wide range of engineering fields of interest such as civil constructions (bridge decks, bridge suspension cables, and tall buildings), industrial applications (heat-exchanger tube arrays in nuclear plants, ocean riser pipes, tall chimneys, and marine piles), and turbomachinery applications. The literature is rich in comprehensive reviews on VIV for bluff bodies (e.g., [3–6]). The most common lock-in problem studied is an elastically mounted circular cylinder constrained to move transverse to the freestream (see for example experiments by Bishop and Hassan [7] or Williamson and Roshko [8]).The pitching motion of airfoils has been studied for many years, focusing on the dynamic stall phenomenon, e.g., [9, 10] for rotorcraft or turbomachinery blades. Bekhti et al. [11] has examined the aeroelastic behavior of the NREL S809 airfoil. In recent years, the dynamic stall problem study has been extended to wind turbine rotors as well, e.g., [12, 13]. Also, several experimental studies, e.g., [14, 15] and numerical studies, e.g., [16, 17] investigate the vortex shedding characteristics for a pitching airfoil. These studies have been confined to low Reynolds flow regimes, perhaps due to the experimental facility limits. In addition, most studies refer to symmetrical airfoils (e.g., NACA0012 or elliptic shaped airfoils). High Reynolds flow is considered by Raveh and Dowell [18] who conducted a study on the frequency lock-in phenomenon for oscillating airfoils in transonic buffeting flows. In the turbomachinery field, Zhu and Kamemoto [19] determined experimentally the lock-in boundaries, for certain angles of attack, of an axial pump blade subjected to pitching motion. An interest in frequency lock-in of airfoils at relatively high Reynolds flow regimes is demonstrated by some work done recently. Tang and Dowell [20] conducted experiments on a NACA0012 at high angles of attack and Re=1×105→5×105. The lock-in behavior of the airfoil oscillating with a prescribed pitching motion about the quarter-chord was determined (“V”-shaped lock-in region). Vortex-induced vibrations on the same airfoil in pitching motion at the same Reynolds number regime have been studied by Besem et al. [21], both experimentally and numerically. Wind tunnel experiments were conducted for incidences α=25°→90°. A computational analysis was also carried out in order to determine the lock-in region at α=40°, which is observed to be not perfectly symmetric about the natural shedding frequency. At this angle of attack, they calculate a Strouhal number of about 0.16 (static airfoil configuration), for a Reynolds number of 2×105. In addition, they report that the airfoil does not shed numerically for α=30°, while in the experiments, the vortex shedding starts to appear at α=25°. This may be due to substantial differences between experimental and numerical studies, e.g., surface roughness (present in the experiments) and turbulence model (used in the CFD model). Lind and Jones [22] have reported experimentally determined Strouhal numbers and aerodynamic loads for a static NACA0012 at Reynolds numbers of up to 6.6×105 and they note that the behavior of thicker sections is qualitatively similar. Pellegrino and Meskell [23] have examined the effect of streamwise translational oscillation on vortex shedding from the S809 airfoil section. The goal of the first part of the current work is to determine the extent of lock-in, for a pitching cambered airfoil (NREL S809) at a given oscillation amplitude for different positive high angles of attack, using uRANS. It should be noted that it in general it is artificial to consider flow-induced vibration mechanisms in isolation. One might expect that stall flutter may interact with vortex shedding. However, as the motion is prescribed, the effect of self-excited vibration mechanisms will be irrelevant in the flow field as the changes to magnitude and phase of the fluid force are exactly counteracted by whatever is required to yield a pure sinusoidal motion at that amplitude. In effect, this study is concerned with the stability of the wake for a given airfoil motion. Some experimental work has been done by Ehrmann et al. [24] on transverse VIV (i.e., perpendicular to the flow direction) of airfoils at α=90°, comparing different shapes to circular cylinder results. Furthermore, in the confined flow of a ducted axial compressor disk, the so-called nonsynchronous vibration (NSV) may occur. The physical mechanics of this self-excited vibration is still not fully understood. It has been attributed to blade tip vortex shedding (e.g., [25]). Clark et al. [26–28] have shown that it is possible to exhibit NSV-like behavior in a 2D simulation, suggesting that vortex shedding may be responsible and propose using CFD data to identify a lumped parameter model based on van der Pol-type representation of NSV.The aim of the second part of this paper is to build a lumped parameter VIV model for a pitching airfoil in an unconfined flow field. An extensive review on the mathematical models used to investigate VIV of circular cylinders has been presented by Gabbai and Benaroya [29]. The idea of using a van der Pol-type oscillator to represent the oscillating force on a cylinder due to vortex shedding is credited to Bishop and Hassan [30]. A representation of fluctuating lift behavior in terms of a proposed wake oscillator is presented by Hartlen and Currie [31], demonstrating the potential of the wake oscillator concept for future VIV studies. This wake oscillator is self-excited and it is subjected to a coupling term proportional to the transverse cylinder velocity. They succeeded in qualitatively capturing experimental behavior, for example, a large cylinder oscillation amplitude resonance region when the vortex shedding frequency is near the natural frequency of the cylinder. Facchinetti et al. [32] examined three different types of coupling (displacement, velocity, and acceleration), i.e., the action of the structure on the fluid wake oscillator. It was found that the acceleration coupling best succeeded in qualitatively modeling the VIV features considered. The VIV wake oscillator model proposed in the present work for a pitching airfoil is based on this new model by Facchinetti et al. [32]. The goal is to qualitatively reproduce the VIV behavior observed with uRANS simulations.To summarize, while there is a considerable body of research in the literature on vortex induced vibration, there is very limited information on the behavior of camber airfoils at high angles, especially at high Reynolds numbers. Furthermore, while there is several mathematical models which explain the phenomenon of lock-in well, it is not clear if these can be applied to airfoils to obtain quantitative results. This paper seeks to contribute to filling this gap.


## 1.1. Context

In 2012, there were a total of 284 GW of wind generating capacity installed globally [1]. The majority of the machines installed are horizontal axis design (HAWT). As with any technology that harvests ambient energy, wind turbines must be permanently exposed to environmental extremes. As a result, the minimum design requirement for wind turbines as defined by the IEC Standard [2] includes specifications for extreme wind loading conditions (see Section 6.3 of this IEC Standard). For example, the structure must be designed to withstand steady, turbulent, and deterministic gusts which occur on average only every 50 years (i.e., a 50-year return period). As part of the transient wind (i.e., gust) environment, a sudden change of wind direction of up to 180° must be considered. In the extreme change of direction scenario, the blade will be suddenly inclined at high angle of attack to the local flow. Hence, the section will behave more like a bluff body than a streamlined airfoil, and so periodic vortex shedding might be expected. To avoid the blades acting as bluff bodies, the rotor is always allowed to rotate in high winds, with wind shed from the blades so as to reduce the power production. The motion of the blades ensures a relatively low angle of incidence regardless of wind direction. However, in a standstill condition, due perhaps to mechanical failure of the transmission, the possibility high local angles of attack increases. Furthermore, as more turbines are installed offshore, the calendar window for maintenance is reduced and so the possibility of a wind turbine in standstill for extended periods increases. This paper will focus on this scenario: a large wind turbine in standstill at high angles of incidence. As there is no rotation, a 2D model is sufficient. The flexibility of the blade is modeled as a prescribed 2 degrees of motion.


## 1.2. Background to Vortex Shedding and Wake Oscillator Model

Vortex shedding induced vibration (VIV) can be a problem for a wide range of structures exposed to a cross flow. A major feature of VIV is “lock-in,” a resonance phenomenon that occurs when the vortex shedding frequency of an elastic bluff body is close to its structural frequency. At this condition, amplification of motion can lead to fatigue or even catastrophic failure. In particular, VIV can be expected in flexible slender bodies of any cross section, and therefore, this field is still an active area of research, in a wide range of engineering fields of interest such as civil constructions (bridge decks, bridge suspension cables, and tall buildings), industrial applications (heat-exchanger tube arrays in nuclear plants, ocean riser pipes, tall chimneys, and marine piles), and turbomachinery applications. The literature is rich in comprehensive reviews on VIV for bluff bodies (e.g., [3–6]). The most common lock-in problem studied is an elastically mounted circular cylinder constrained to move transverse to the freestream (see for example experiments by Bishop and Hassan [7] or Williamson and Roshko [8]).The pitching motion of airfoils has been studied for many years, focusing on the dynamic stall phenomenon, e.g., [9, 10] for rotorcraft or turbomachinery blades. Bekhti et al. [11] has examined the aeroelastic behavior of the NREL S809 airfoil. In recent years, the dynamic stall problem study has been extended to wind turbine rotors as well, e.g., [12, 13]. Also, several experimental studies, e.g., [14, 15] and numerical studies, e.g., [16, 17] investigate the vortex shedding characteristics for a pitching airfoil. These studies have been confined to low Reynolds flow regimes, perhaps due to the experimental facility limits. In addition, most studies refer to symmetrical airfoils (e.g., NACA0012 or elliptic shaped airfoils). High Reynolds flow is considered by Raveh and Dowell [18] who conducted a study on the frequency lock-in phenomenon for oscillating airfoils in transonic buffeting flows. In the turbomachinery field, Zhu and Kamemoto [19] determined experimentally the lock-in boundaries, for certain angles of attack, of an axial pump blade subjected to pitching motion. An interest in frequency lock-in of airfoils at relatively high Reynolds flow regimes is demonstrated by some work done recently. Tang and Dowell [20] conducted experiments on a NACA0012 at high angles of attack and Re=1×105→5×105. The lock-in behavior of the airfoil oscillating with a prescribed pitching motion about the quarter-chord was determined (“V”-shaped lock-in region). Vortex-induced vibrations on the same airfoil in pitching motion at the same Reynolds number regime have been studied by Besem et al. [21], both experimentally and numerically. Wind tunnel experiments were conducted for incidences α=25°→90°. A computational analysis was also carried out in order to determine the lock-in region at α=40°, which is observed to be not perfectly symmetric about the natural shedding frequency. At this angle of attack, they calculate a Strouhal number of about 0.16 (static airfoil configuration), for a Reynolds number of 2×105. In addition, they report that the airfoil does not shed numerically for α=30°, while in the experiments, the vortex shedding starts to appear at α=25°. This may be due to substantial differences between experimental and numerical studies, e.g., surface roughness (present in the experiments) and turbulence model (used in the CFD model). Lind and Jones [22] have reported experimentally determined Strouhal numbers and aerodynamic loads for a static NACA0012 at Reynolds numbers of up to 6.6×105 and they note that the behavior of thicker sections is qualitatively similar. Pellegrino and Meskell [23] have examined the effect of streamwise translational oscillation on vortex shedding from the S809 airfoil section. The goal of the first part of the current work is to determine the extent of lock-in, for a pitching cambered airfoil (NREL S809) at a given oscillation amplitude for different positive high angles of attack, using uRANS. It should be noted that it in general it is artificial to consider flow-induced vibration mechanisms in isolation. One might expect that stall flutter may interact with vortex shedding. However, as the motion is prescribed, the effect of self-excited vibration mechanisms will be irrelevant in the flow field as the changes to magnitude and phase of the fluid force are exactly counteracted by whatever is required to yield a pure sinusoidal motion at that amplitude. In effect, this study is concerned with the stability of the wake for a given airfoil motion. Some experimental work has been done by Ehrmann et al. [24] on transverse VIV (i.e., perpendicular to the flow direction) of airfoils at α=90°, comparing different shapes to circular cylinder results. Furthermore, in the confined flow of a ducted axial compressor disk, the so-called nonsynchronous vibration (NSV) may occur. The physical mechanics of this self-excited vibration is still not fully understood. It has been attributed to blade tip vortex shedding (e.g., [25]). Clark et al. [26–28] have shown that it is possible to exhibit NSV-like behavior in a 2D simulation, suggesting that vortex shedding may be responsible and propose using CFD data to identify a lumped parameter model based on van der Pol-type representation of NSV.The aim of the second part of this paper is to build a lumped parameter VIV model for a pitching airfoil in an unconfined flow field. An extensive review on the mathematical models used to investigate VIV of circular cylinders has been presented by Gabbai and Benaroya [29]. The idea of using a van der Pol-type oscillator to represent the oscillating force on a cylinder due to vortex shedding is credited to Bishop and Hassan [30]. A representation of fluctuating lift behavior in terms of a proposed wake oscillator is presented by Hartlen and Currie [31], demonstrating the potential of the wake oscillator concept for future VIV studies. This wake oscillator is self-excited and it is subjected to a coupling term proportional to the transverse cylinder velocity. They succeeded in qualitatively capturing experimental behavior, for example, a large cylinder oscillation amplitude resonance region when the vortex shedding frequency is near the natural frequency of the cylinder. Facchinetti et al. [32] examined three different types of coupling (displacement, velocity, and acceleration), i.e., the action of the structure on the fluid wake oscillator. It was found that the acceleration coupling best succeeded in qualitatively modeling the VIV features considered. The VIV wake oscillator model proposed in the present work for a pitching airfoil is based on this new model by Facchinetti et al. [32]. The goal is to qualitatively reproduce the VIV behavior observed with uRANS simulations.To summarize, while there is a considerable body of research in the literature on vortex induced vibration, there is very limited information on the behavior of camber airfoils at high angles, especially at high Reynolds numbers. Furthermore, while there is several mathematical models which explain the phenomenon of lock-in well, it is not clear if these can be applied to airfoils to obtain quantitative results. This paper seeks to contribute to filling this gap.


## 2. CFD Model

### 2.1. Methodology

A 2D blade section needs three independent degrees of freedom to describe any general motion. The conventional motion decomposition is shown in Figure1: two translational (transverse and streamwise, i.e., perpendicular and in-line with the flow direction, respectively) and one rotational (pitching).Figure 1
Schematic of a 2D airfoil degrees of freedom at a positive incidence with respect to the flow direction. The airfoil is the NREL S809. The highlighted motion is the only considered in this study.In this work, only the pure pitching motion has been considered. A series of two-dimensional unsteady simulations of the flow around the pitching NREL S809 airfoil at different angles of attack were performed. Figure1 also shows the positive angular direction considered for the airfoil inclination with respect to the freestream direction. All calculations have been made at a Reynolds number Re=106, referred to the airfoil’s chord length. The software used is the commercially available uRANS solver ANSYS Fluent. In order to determine the lock-in region for the pitching airfoil considered at a given mean angle of attack, the vortex shedding characteristics for the airfoil in a static configuration at that angle of attack must be known; see Pellegrino and Meskell [33]. The vortex shedding frequency depends on the mean angle of attack. With the static vortex shedding frequency fv0 (i.e., shedding frequency with no airfoil motion) at each mean angle of attack considered, αmean, and fp being the frequency of pitching oscillation, the frequency ratio, r, can be defined as follows:
(1)r=fpfv0.In this study, the frequency ratio,r, is varied for a sinusoidal pitching amplitude of 7°. As the static vortex shedding frequency is fixed, the pitching frequency, fp, is varied to obtain each value of r. The range of mean incidences investigated is (35°≤α≤145°), as no periodic vortex shedding has been observed in the ranges (0°≤α≤30°) and (150°≤α≤180°) at the Reynolds number considered; see [33]. In the ranges (35°≤α≤50°) and (130°≤α≤145°), a αmean increment of 5° was used, whereas in the range (50°≤α≤130°), the αmean increment was 10°, as in this region, no rapid change in the lock-in characteristics has been observed. Furthermore, only the positive range of angles of attack (0°≤α≤180°) has been considered; as for the negative range (0°≥α≥−180°), it is assumed that the behavior will be qualitatively the same. This is based on the vortex shedding behavior for the static airfoil [33]. The assessment of the lock-in condition is obtained with an analysis of the autospectrum of the lift force coefficient, since the frequencies present in the wake are expected to be equal to the frequencies observed in the fluctuating lift force coefficient [34]. Therefore, if the vortex shedding frequency locks in to the structural frequency, a unique dominant frequency and its harmonics will be observed in the fluctuating lift force coefficient [7]. The lift force coefficient results are obtained for a simulated time of 30 s.


### 2.2. Model Set-Up

The fluid domain has been created with the ANSYS default meshing tool. The domain, shown in Figure2, extends to 10 c (airfoil chord) length from the body in the upstream direction, 40 c in the downstream direction, and 20 c in the cross-flow direction. The reference chord length is 1 m. The model is similar to that which was implemented by Pellegrino and Meskell [23] for the case of vortex shedding from a translating S809 airfoil.Figure 2
Flow field domain: dimensions and boundary conditions.The fluid domain has been divided into different zones in which the mesh density is different, keeping it finer in the vicinity of the airfoil and fading in the far wake downstream. Note that these zones are used for mesh generation and are not boundaries in the flow field. The inner zone has been modeled as a circle in order to facilitate the motion of the airfoil with respect to the incoming flow. That is, for simplicity, a sliding mesh technique has been used rather than either a remeshing or deforming mesh scheme in order to reduce the computation time. The pivot point has been set at the middle of chordc/2 (i.e., the point on the chord line is equidistant from leading and trailing edge). One could argue that it is common to use the pitching moment around the quarter of chord (25% of the chord from the leading edge) as this is theoretically the aerodynamic center in a subsonic regime. The reason for choosing the midchord point is that this airfoil is at such a high mean angle of attack that the theoretical aerodynamic center is not relevant and the behavior as a bluff body is of most interest here. Furthermore, for a vibrating blade in torsional motion (i.e., pitching), the elastic axis is likely to be aft of the quarter-chord point. For a thick airfoil section such at the one used here, the elastic axis will be closed to the midsection.The pitching motion has been applied with a compiled user defined function (UDF). The motion is described as follows:(2)αt=αmean+αamplsin2πfpt,where αmean is the mean angle of attack, αampl is the motion angular amplitude, and fp is the pitching motion frequency.The mesh is unstructured, with triangular elements throughout except the quadratic elements in the near-wall region. The near-wall mesh has 512 cells in the chordwise direction, and the overall mesh is composed in total of approximately 100,000 nodes. The no-slip condition was used at the airfoil surface. Inlet flow was specified at the upstream boundary of the domain, with the inflow conditions shown in Table1. Flow velocity has been set based on a Reynolds number value of Re=106 with air as the fluid. The outlet condition was specified in the downstream boundary of the domain with a reference pressure pgauge=0.Table 1
Inlet flow properties.
Inlet flow properties (Re=106)
Velocity
U∞
14.607
m/s
Density
ρ
1.225
kg/m3
Viscosity
μ
1.7894×10−5
kg/msTurbulence closure is modeled by the shear stress transport (SST) eddy viscosity model. This transition turbulence model requires a fine mesh in the vicinity of the airfoil. For that purpose, the height of the first cell has been tuned in order to obtain ay+ value below 2 everywhere on the airfoil surface. Computations have been made using the PISO scheme, with second-order upwind discretization for the momentum and a time step of 2.5 ms. It is also useful to express time in a dimensionless fashion as follows:
(3)τ=tU∞c.Therefore, the dimensionless time step used isΔτ=0.0365.A sensitivity analysis of the model (airfoil static case) has been conducted considering the Strouhal number of the fixed airfoil atα=40° (St=0.148). The aim is to assess if a sufficient grid resolution is applied, as well as a reasonable time step value. This implicitly assumes that the Strouhal number is well simulated, and then the lock-in range will also be properly captured. The size of the grid elements has been doubled and halved; the finer mesh does not give a different result of the Strouhal number, while it does the coarser one (see Table 2). Hence, the used grid can be considered appropriate. Then, four time steps (nΔt, with n=4,2,1,1/2) have been used with the same grid (see Table 3). Higher time steps give a different Strouhal number value, whereas a smaller time step gives the same result. Hence, the choice of Δt=0.0025s.Table 2
Grid sensitivity analysis.
Grid
St
Used
0.148
Coarser
0.145
Finer
0.148Table 3
Time step sensitivity analysis.
Δt (s)
St
0.01
0.144
0.005
0.146
0.0025
0.148
0.00125
0.148


### 2.3. Results

The unsteady flow has been simulated at first for one oscillation amplitude only (αampl=7°) over the range at which vortex shedding has been observed, (35°≤α≤145°), for a range of frequency ratios around r=1. The principal objective was to find the extent of the lock-in range (i.e., where lock-in will occur). Then, an analysis on the oscillation amplitude effect has been conducted for angular positions at which the lock-in appears to be more significant.In order to illustrate the effect lock-in has on the vortex shedding excitation that would be experienced by the structure, the behavior of the lift coefficient atα=140° is considered. Figure 3 shows a portion of the time record of the lift coefficient for two frequency ratios. For the sake of clarity, only 2 seconds of simulated time is shown, after the wake has well developed; dimensionless time is used. At a frequency ratio r=0.8, shown in Figure 3(a), the vortex shedding is not locked-in to the pitching frequency, while Figure 3(b) shows the equivalent record for the locked-in case (r=1.1). It is clear that when the vortex shedding has locked-in, the lift coefficient is periodic, with a single dominant frequency, albeit with harmonics. In contrast, in Figure 3(a), the trace is clearly not periodic, characteristic of the locked-out condition. This behavior is consistent with the results of an oscillating circular cylinder studied among others by Bishop and Hassan [7]. Although this appears to be narrowband random, in fact the behavior is a pseudo-periodic response typical of a signal with two distinct, incommensurate but similar frequencies.Figure 3
Detail of the time record of the lift force coefficient atα=140°, for different frequency ratios showing the effect of lock-out (a) and lock-in (b).
(a)
r=0:8 (locked-out)
(b)
r=1:1 (locked-in)Figure4 shows the autospectra associated with the two cases seen in Figure 3. Frequencies have been normalized with the static vortex shedding frequency, giving a nondimensional frequency: f∗=f/fv0. In the locked-out case, for r=0.8 (Figure 4(a)), the pitching and the vortex shedding frequencies can be seen at f∗=0.8 and in the vicinity of f∗=1, respectively. A vortex shedding frequency at f∗=1 (highlighted by the dash-dot line) would occur if the airfoil was static. There is also a smaller peak at f∗=1.6, second harmonic of the excitation frequency, fp. Note that representing this pseudo-periodic response series in the frequency domain based on a relatively short time series gives rise to a broadening of the vortex shedding peak. In contrast, the spectrum obtained for the lock-in cases shown in Figure 4(b) presents clearly the unique dominating frequency (and its second harmonic). It can be clearly seen that as the motion frequency moves closer to the vortex shedding frequency, eventually these two frequencies coalesce (i.e., the vortex shedding frequency jumps onto the motion frequency).Figure 4
Spectra of the lift force coefficient atα=140°, for different frequency ratios showing the effect of lock-out (a) and lock-in (b). PSD arbitrarily scaled. Dash-dot line indicates f∗=1, corresponding to the static case.
(a)
r=0:8 (locked-out)
(b)
r=1:1 (locked-in)The lock-in behavior of vortex shedding in the presence of pitching motion can be summarized over the entire angle of attack and frequency ranges, shown in Figure5. The lock-in condition is present for the entire angle of attack range, but the range of frequencies depends strongly on angle. As can be clearly seen, the lock-in range is in general more significant when the airfoil presents a smaller characteristic length (vertical projected chord), i.e., in the regions around α=40° and α=140°.Figure 5
Lock-in map for positive angles of attack, pitching motion atαampl=7°: ○, lock-in; ×, lock-out and lock-in region (in green).The effect of the oscillation amplitude on the lock-in aroundr=1 has been assessed for angles α=40° and α=140°, for which the lock-in regions are shown, respectively, in Figures 6(a) and 6(b). As expected, the extent of the frequency lock-in range decreases if the oscillation amplitude is reduced. Moreover, it can be seen that the lock-in region is broadly smaller for α=140°: this suggests that a sharp trailing edge in the upstream position seems to slightly reduce the risk of lock-in. It must be noted that for α=40° and r>1, the lock-in boundary is represented by a vertical line (Figure 6(a)). This may not be indicative of the real behavior, being due to the resolution of the frequency ratio (the used step is r=0.05), meaning that this boundary might not be vertical using a higher resolution in the parametric study.Figure 6
Lock-in region for pitching amplitude analysis at two incidences:○, lock-in; ×, lock-out and lock-in region (in-green).
(a)
α=40°
(b)
α=140°Generally, it is not clear why the wake should be more susceptible to synchronizing at angles of attack with shallower characteristic length (Figure5), but it may be related to the overall structure of the wake before lock-in occurs. Figure 7 shows the vortex street at different angles of attack for static airfoil condition. It can be seen that as α is increased towards a vertical position of the airfoil to the flow (Figures 7(a), 7(c), and 7(d)), the vortex shape becomes bigger and more elongated. Also, the vortex spacing increases both in the x and y directions, consistent with results found in the literature about the Kármán structure [35]. Figures 7(a) and 7(b) show that the wake structure for α=40° and α=140° is very similar, seen in a specular way as the trailing edge will be either downwind or upwind. Figure 8 shows the vortex street (i.e., instantaneous vorticity field) behind the pitching airfoil for three different cases (two of them already considered above in Figures 3 and 4). Consider the instantaneous vorticity field of two locked-out cases. Figure 8(c) shows the vorticity for an angle of attack, α=70°, and a frequency ratio, r=0.97. The wake has arranged itself into an antisymmetric configuration similar to the Kármán street often associated with periodic vortex shedding. This configuration is relatively stable and so would be more resilient against external forcing and hence exhibits a narrower lock-in range. In this situation, the separation points (the leading and trailing edge of the airfoil) are almost at the same streamwise location, as the blade is almost normal to the flow. In Figure 8(a), the wake behind the airfoil at α=140° and r=0.8 shows a different structure. The vortices have been formed and shed periodically, but they have partially paired, in a pattern reminiscent of the 2P mode. This may be due to fact that upstream separation point (i.e., the sharp “trailing edge”) is now 0.76 chord lengths upstream of the other separation point (i.e., the blunt “leading edge”). This configuration is less stable than the 2S mode. Once lock-in is established (e.g., r=1.1), the wake structure becomes antisymmetric as shown in Figure 8(b).Figure 7
Instantaneous vorticity magnitude at different incidences for static airfoil,t=10s, showing different vortex street structures, counter-clockwise (positive) vorticity (in pink), and clockwise (negative) vorticity (in blue).
(a)
α=40°
(b)
α=140°
(c)
α=50°
(d)
α=70°Figure 8
Instantaneous vorticity field for pitching motion at different conditions:αmean=140° (a, b) and αmean=70° (c), counter-clockwise (positive) vorticity (in pink), and clockwise (negative) vorticity (in blue).
(a)
α=140°; frequency ratio, r=0:8. Vortex shedding locked-out
(b)
α=140°; frequency ratio, r=1:1. Vortex shedding locked-in
(c)
α=50°; frequency ratio, r=0:97. Vortex shedding locked-outThe point made above, about the significantly different behavior whether the airfoil position presents a smaller or larger characteristic length, can be explained with the help of Figure9. Two incidences, αmean=40° with a smaller characteristic length and αmean=70° with a larger one, are used to contrast how the flow differently reattaches on the leading edge during the downstroke pitching motion. Figure 9(a) shows the instantaneous horizontal velocity for αmean=40° at its maximum pitching position, while Figure 9(b) at its minimum one (for αampl=7°). It can be seen that the downstroke motion causes the reattachment point of the flow to move substantially on the upper surface. This leads to a reduction of the effective characteristic length of the airfoil. Conversely, for αmean=70°, Figures 9(c) and 9(d) show that the reattachment point does not move significantly comparing the upper (Figure 9(c)) and the lower (Figure 9(d)) pitching positions. This is consistent with what was observed in Figures 5 and 8, i.e., the more resilient wake for αmean=70° as well as its two separation points being almost at the same streamwise location.Figure 9
Instantaneous horizontal velocity showing the flow reattachment point on the leading edge forαmean=40° (a, b) and αmean=70° (c, d) on their max and min pitching position (αampl=7°); flow velocity left-to-right (in pink) and flow velocity right-to-left (in blue).
(a)
αmean=40°; upper position
(b)
αmean=40°; lower position
(c)
αmean=70°; upper position
(d)
αmean=70°; lower position


## 2.1. Methodology

A 2D blade section needs three independent degrees of freedom to describe any general motion. The conventional motion decomposition is shown in Figure1: two translational (transverse and streamwise, i.e., perpendicular and in-line with the flow direction, respectively) and one rotational (pitching).Figure 1
Schematic of a 2D airfoil degrees of freedom at a positive incidence with respect to the flow direction. The airfoil is the NREL S809. The highlighted motion is the only considered in this study.In this work, only the pure pitching motion has been considered. A series of two-dimensional unsteady simulations of the flow around the pitching NREL S809 airfoil at different angles of attack were performed. Figure1 also shows the positive angular direction considered for the airfoil inclination with respect to the freestream direction. All calculations have been made at a Reynolds number Re=106, referred to the airfoil’s chord length. The software used is the commercially available uRANS solver ANSYS Fluent. In order to determine the lock-in region for the pitching airfoil considered at a given mean angle of attack, the vortex shedding characteristics for the airfoil in a static configuration at that angle of attack must be known; see Pellegrino and Meskell [33]. The vortex shedding frequency depends on the mean angle of attack. With the static vortex shedding frequency fv0 (i.e., shedding frequency with no airfoil motion) at each mean angle of attack considered, αmean, and fp being the frequency of pitching oscillation, the frequency ratio, r, can be defined as follows:
(1)r=fpfv0.In this study, the frequency ratio,r, is varied for a sinusoidal pitching amplitude of 7°. As the static vortex shedding frequency is fixed, the pitching frequency, fp, is varied to obtain each value of r. The range of mean incidences investigated is (35°≤α≤145°), as no periodic vortex shedding has been observed in the ranges (0°≤α≤30°) and (150°≤α≤180°) at the Reynolds number considered; see [33]. In the ranges (35°≤α≤50°) and (130°≤α≤145°), a αmean increment of 5° was used, whereas in the range (50°≤α≤130°), the αmean increment was 10°, as in this region, no rapid change in the lock-in characteristics has been observed. Furthermore, only the positive range of angles of attack (0°≤α≤180°) has been considered; as for the negative range (0°≥α≥−180°), it is assumed that the behavior will be qualitatively the same. This is based on the vortex shedding behavior for the static airfoil [33]. The assessment of the lock-in condition is obtained with an analysis of the autospectrum of the lift force coefficient, since the frequencies present in the wake are expected to be equal to the frequencies observed in the fluctuating lift force coefficient [34]. Therefore, if the vortex shedding frequency locks in to the structural frequency, a unique dominant frequency and its harmonics will be observed in the fluctuating lift force coefficient [7]. The lift force coefficient results are obtained for a simulated time of 30 s.


## 2.2. Model Set-Up

The fluid domain has been created with the ANSYS default meshing tool. The domain, shown in Figure2, extends to 10 c (airfoil chord) length from the body in the upstream direction, 40 c in the downstream direction, and 20 c in the cross-flow direction. The reference chord length is 1 m. The model is similar to that which was implemented by Pellegrino and Meskell [23] for the case of vortex shedding from a translating S809 airfoil.Figure 2
Flow field domain: dimensions and boundary conditions.The fluid domain has been divided into different zones in which the mesh density is different, keeping it finer in the vicinity of the airfoil and fading in the far wake downstream. Note that these zones are used for mesh generation and are not boundaries in the flow field. The inner zone has been modeled as a circle in order to facilitate the motion of the airfoil with respect to the incoming flow. That is, for simplicity, a sliding mesh technique has been used rather than either a remeshing or deforming mesh scheme in order to reduce the computation time. The pivot point has been set at the middle of chordc/2 (i.e., the point on the chord line is equidistant from leading and trailing edge). One could argue that it is common to use the pitching moment around the quarter of chord (25% of the chord from the leading edge) as this is theoretically the aerodynamic center in a subsonic regime. The reason for choosing the midchord point is that this airfoil is at such a high mean angle of attack that the theoretical aerodynamic center is not relevant and the behavior as a bluff body is of most interest here. Furthermore, for a vibrating blade in torsional motion (i.e., pitching), the elastic axis is likely to be aft of the quarter-chord point. For a thick airfoil section such at the one used here, the elastic axis will be closed to the midsection.The pitching motion has been applied with a compiled user defined function (UDF). The motion is described as follows:(2)αt=αmean+αamplsin2πfpt,where αmean is the mean angle of attack, αampl is the motion angular amplitude, and fp is the pitching motion frequency.The mesh is unstructured, with triangular elements throughout except the quadratic elements in the near-wall region. The near-wall mesh has 512 cells in the chordwise direction, and the overall mesh is composed in total of approximately 100,000 nodes. The no-slip condition was used at the airfoil surface. Inlet flow was specified at the upstream boundary of the domain, with the inflow conditions shown in Table1. Flow velocity has been set based on a Reynolds number value of Re=106 with air as the fluid. The outlet condition was specified in the downstream boundary of the domain with a reference pressure pgauge=0.Table 1
Inlet flow properties.
Inlet flow properties (Re=106)
Velocity
U∞
14.607
m/s
Density
ρ
1.225
kg/m3
Viscosity
μ
1.7894×10−5
kg/msTurbulence closure is modeled by the shear stress transport (SST) eddy viscosity model. This transition turbulence model requires a fine mesh in the vicinity of the airfoil. For that purpose, the height of the first cell has been tuned in order to obtain ay+ value below 2 everywhere on the airfoil surface. Computations have been made using the PISO scheme, with second-order upwind discretization for the momentum and a time step of 2.5 ms. It is also useful to express time in a dimensionless fashion as follows:
(3)τ=tU∞c.Therefore, the dimensionless time step used isΔτ=0.0365.A sensitivity analysis of the model (airfoil static case) has been conducted considering the Strouhal number of the fixed airfoil atα=40° (St=0.148). The aim is to assess if a sufficient grid resolution is applied, as well as a reasonable time step value. This implicitly assumes that the Strouhal number is well simulated, and then the lock-in range will also be properly captured. The size of the grid elements has been doubled and halved; the finer mesh does not give a different result of the Strouhal number, while it does the coarser one (see Table 2). Hence, the used grid can be considered appropriate. Then, four time steps (nΔt, with n=4,2,1,1/2) have been used with the same grid (see Table 3). Higher time steps give a different Strouhal number value, whereas a smaller time step gives the same result. Hence, the choice of Δt=0.0025s.Table 2
Grid sensitivity analysis.
Grid
St
Used
0.148
Coarser
0.145
Finer
0.148Table 3
Time step sensitivity analysis.
Δt (s)
St
0.01
0.144
0.005
0.146
0.0025
0.148
0.00125
0.148


## 2.3. Results

The unsteady flow has been simulated at first for one oscillation amplitude only (αampl=7°) over the range at which vortex shedding has been observed, (35°≤α≤145°), for a range of frequency ratios around r=1. The principal objective was to find the extent of the lock-in range (i.e., where lock-in will occur). Then, an analysis on the oscillation amplitude effect has been conducted for angular positions at which the lock-in appears to be more significant.In order to illustrate the effect lock-in has on the vortex shedding excitation that would be experienced by the structure, the behavior of the lift coefficient atα=140° is considered. Figure 3 shows a portion of the time record of the lift coefficient for two frequency ratios. For the sake of clarity, only 2 seconds of simulated time is shown, after the wake has well developed; dimensionless time is used. At a frequency ratio r=0.8, shown in Figure 3(a), the vortex shedding is not locked-in to the pitching frequency, while Figure 3(b) shows the equivalent record for the locked-in case (r=1.1). It is clear that when the vortex shedding has locked-in, the lift coefficient is periodic, with a single dominant frequency, albeit with harmonics. In contrast, in Figure 3(a), the trace is clearly not periodic, characteristic of the locked-out condition. This behavior is consistent with the results of an oscillating circular cylinder studied among others by Bishop and Hassan [7]. Although this appears to be narrowband random, in fact the behavior is a pseudo-periodic response typical of a signal with two distinct, incommensurate but similar frequencies.Figure 3
Detail of the time record of the lift force coefficient atα=140°, for different frequency ratios showing the effect of lock-out (a) and lock-in (b).
(a)
r=0:8 (locked-out)
(b)
r=1:1 (locked-in)Figure4 shows the autospectra associated with the two cases seen in Figure 3. Frequencies have been normalized with the static vortex shedding frequency, giving a nondimensional frequency: f∗=f/fv0. In the locked-out case, for r=0.8 (Figure 4(a)), the pitching and the vortex shedding frequencies can be seen at f∗=0.8 and in the vicinity of f∗=1, respectively. A vortex shedding frequency at f∗=1 (highlighted by the dash-dot line) would occur if the airfoil was static. There is also a smaller peak at f∗=1.6, second harmonic of the excitation frequency, fp. Note that representing this pseudo-periodic response series in the frequency domain based on a relatively short time series gives rise to a broadening of the vortex shedding peak. In contrast, the spectrum obtained for the lock-in cases shown in Figure 4(b) presents clearly the unique dominating frequency (and its second harmonic). It can be clearly seen that as the motion frequency moves closer to the vortex shedding frequency, eventually these two frequencies coalesce (i.e., the vortex shedding frequency jumps onto the motion frequency).Figure 4
Spectra of the lift force coefficient atα=140°, for different frequency ratios showing the effect of lock-out (a) and lock-in (b). PSD arbitrarily scaled. Dash-dot line indicates f∗=1, corresponding to the static case.
(a)
r=0:8 (locked-out)
(b)
r=1:1 (locked-in)The lock-in behavior of vortex shedding in the presence of pitching motion can be summarized over the entire angle of attack and frequency ranges, shown in Figure5. The lock-in condition is present for the entire angle of attack range, but the range of frequencies depends strongly on angle. As can be clearly seen, the lock-in range is in general more significant when the airfoil presents a smaller characteristic length (vertical projected chord), i.e., in the regions around α=40° and α=140°.Figure 5
Lock-in map for positive angles of attack, pitching motion atαampl=7°: ○, lock-in; ×, lock-out and lock-in region (in green).The effect of the oscillation amplitude on the lock-in aroundr=1 has been assessed for angles α=40° and α=140°, for which the lock-in regions are shown, respectively, in Figures 6(a) and 6(b). As expected, the extent of the frequency lock-in range decreases if the oscillation amplitude is reduced. Moreover, it can be seen that the lock-in region is broadly smaller for α=140°: this suggests that a sharp trailing edge in the upstream position seems to slightly reduce the risk of lock-in. It must be noted that for α=40° and r>1, the lock-in boundary is represented by a vertical line (Figure 6(a)). This may not be indicative of the real behavior, being due to the resolution of the frequency ratio (the used step is r=0.05), meaning that this boundary might not be vertical using a higher resolution in the parametric study.Figure 6
Lock-in region for pitching amplitude analysis at two incidences:○, lock-in; ×, lock-out and lock-in region (in-green).
(a)
α=40°
(b)
α=140°Generally, it is not clear why the wake should be more susceptible to synchronizing at angles of attack with shallower characteristic length (Figure5), but it may be related to the overall structure of the wake before lock-in occurs. Figure 7 shows the vortex street at different angles of attack for static airfoil condition. It can be seen that as α is increased towards a vertical position of the airfoil to the flow (Figures 7(a), 7(c), and 7(d)), the vortex shape becomes bigger and more elongated. Also, the vortex spacing increases both in the x and y directions, consistent with results found in the literature about the Kármán structure [35]. Figures 7(a) and 7(b) show that the wake structure for α=40° and α=140° is very similar, seen in a specular way as the trailing edge will be either downwind or upwind. Figure 8 shows the vortex street (i.e., instantaneous vorticity field) behind the pitching airfoil for three different cases (two of them already considered above in Figures 3 and 4). Consider the instantaneous vorticity field of two locked-out cases. Figure 8(c) shows the vorticity for an angle of attack, α=70°, and a frequency ratio, r=0.97. The wake has arranged itself into an antisymmetric configuration similar to the Kármán street often associated with periodic vortex shedding. This configuration is relatively stable and so would be more resilient against external forcing and hence exhibits a narrower lock-in range. In this situation, the separation points (the leading and trailing edge of the airfoil) are almost at the same streamwise location, as the blade is almost normal to the flow. In Figure 8(a), the wake behind the airfoil at α=140° and r=0.8 shows a different structure. The vortices have been formed and shed periodically, but they have partially paired, in a pattern reminiscent of the 2P mode. This may be due to fact that upstream separation point (i.e., the sharp “trailing edge”) is now 0.76 chord lengths upstream of the other separation point (i.e., the blunt “leading edge”). This configuration is less stable than the 2S mode. Once lock-in is established (e.g., r=1.1), the wake structure becomes antisymmetric as shown in Figure 8(b).Figure 7
Instantaneous vorticity magnitude at different incidences for static airfoil,t=10s, showing different vortex street structures, counter-clockwise (positive) vorticity (in pink), and clockwise (negative) vorticity (in blue).
(a)
α=40°
(b)
α=140°
(c)
α=50°
(d)
α=70°Figure 8
Instantaneous vorticity field for pitching motion at different conditions:αmean=140° (a, b) and αmean=70° (c), counter-clockwise (positive) vorticity (in pink), and clockwise (negative) vorticity (in blue).
(a)
α=140°; frequency ratio, r=0:8. Vortex shedding locked-out
(b)
α=140°; frequency ratio, r=1:1. Vortex shedding locked-in
(c)
α=50°; frequency ratio, r=0:97. Vortex shedding locked-outThe point made above, about the significantly different behavior whether the airfoil position presents a smaller or larger characteristic length, can be explained with the help of Figure9. Two incidences, αmean=40° with a smaller characteristic length and αmean=70° with a larger one, are used to contrast how the flow differently reattaches on the leading edge during the downstroke pitching motion. Figure 9(a) shows the instantaneous horizontal velocity for αmean=40° at its maximum pitching position, while Figure 9(b) at its minimum one (for αampl=7°). It can be seen that the downstroke motion causes the reattachment point of the flow to move substantially on the upper surface. This leads to a reduction of the effective characteristic length of the airfoil. Conversely, for αmean=70°, Figures 9(c) and 9(d) show that the reattachment point does not move significantly comparing the upper (Figure 9(c)) and the lower (Figure 9(d)) pitching positions. This is consistent with what was observed in Figures 5 and 8, i.e., the more resilient wake for αmean=70° as well as its two separation points being almost at the same streamwise location.Figure 9
Instantaneous horizontal velocity showing the flow reattachment point on the leading edge forαmean=40° (a, b) and αmean=70° (c, d) on their max and min pitching position (αampl=7°); flow velocity left-to-right (in pink) and flow velocity right-to-left (in blue).
(a)
αmean=40°; upper position
(b)
αmean=40°; lower position
(c)
αmean=70°; upper position
(d)
αmean=70°; lower position


## 3. VIV Lumped Parameter Model

### 3.1. Methodology

In this section, a lumped parameter model for VIV of a pitching airfoil is presented. It is based on the wake oscillator proposed by Facchinetti et al. [32], which in turn is based on the classical van der Pol oscillator. A nonlinear model is chosen as the amplitude of oscillation in the simulations is by necessity nonzero. So even though limit cycle oscillations or force magnitudes are not of concern, the interaction between oscillation amplitude and resulting fluid force with lock-in must be accounted for. The wake is modeled as a nonlinear oscillator as follows:
(4)q¨+εΩfq2−1q̇+Ωf2q=F,where q is the dimensionless wake variable (overdots mean derivative with respect to the dimensionless time), ε is effectively a damping parameter, Ωf is the vortex shedding angular frequency, and F models the effect of the structure motion on the near wake by providing coupling with structural dynamics. The wake variable, q, is here associated with the fluctuating lift coefficient of the structure. The coupling term is based on the acceleration of the structure as recommended by Facchinetti et al. [32] as follows:
(5)F=A/Ds¨,where s¨ is a generic acceleration of the structure, A is the coupling parameter, and D is the characteristic length (e.g., cylinder diameter). Note that Facchinetti et al. [32] presented a coupled 2 degrees of freedom system, i.e., the structure was free to vibrate. In this paper, only the wake oscillator is considered, since the acceleration is imposed (the structure motion is prescribed as per equation 2).Remembering that the focus of this paper is the airfoil pitching motion, the wake oscillator equation has been reformulated as follows:(6)q¨+ε2πfv0q2−1q̇+2πfv02q=A/Lα¨,where q,ε,A are defined in the same way as in equation (4); fv0=StU∞/L is the static vortex shedding frequency, L is the characteristic length of the airfoil, and α¨ is the angular acceleration. Note that if the mean angle of attack varies, so does the characteristic length as it is defined to be the vertical projection of the chord with respect to the flow direction [33].The pitching motion acceleration can be derived by equation (2), yielding as follows:
(7)α¨t=−2πfp2αamplsin2πfpt.


### 3.2. Verification and Parameter Setting

Figure10 shows the lock-in region for a circular cylinder compared with the parameters specified by Facchinetti et al. [32]. The lock-in behavior obtained from the current model implementation compares well with the lock-in boundaries presented by Facchinetti et al. [32], offering confidence in the current implementation.Figure 10
Comparison of predicted lock-in boundaries for a circular cylinder:□, upper boundary (current implementation); ∇, lower boundary (current implementation); lock-in region (current implementation) (in green) and [32] (in black line).As the wake oscillator model is a semiempirical model, some of its parameters require empirical data to be quantified. In the current study, due to the apparent lack of experimental data on VIV of pitching airfoils at high angle of attack and high Reynolds number, the “empirical” data used are the lock-in boundaries presented in Figures5 and 6 above from CFD simulations. The parameters to be defined are ε and A. A harmonic balance method was used initially for two sample cases (α=40° and α=140°), using the harmonic coefficients obtained from the study of vortex shedding from the static NREL S809 [33]. However, the prediction of lock-in boundaries did not compare well with the dynamic CFD results. This suggests that maybe it is not possible to obtain a universal set of parameters for a VIV model of airfoil independent of angle of attack. An attempt has been made to renormalize the coupling factor based on the projected chord, in a similar manner to that used by the authors [33] to reduce the dependency of the Strouhal number on angle of attack. However, this was not successful.Another approach has been used: estimatingA for each angle of attack and keeping ε constant. The rationale for this approach is the sensitivity of the model predictions to these parameters. A sensitivity analysis of ε and A has been conducted using the model proposed by Facchinetti et al. [32] for a circular cylinder (m★=4). It can be clearly seen in Figure 11(a) that the reduced velocity range for lock-in is not affected by changing the ε value (at least for the considered variation of ±100%). In contrast, Figure 11(b) shows that a variation of A significantly affects the system behavior. Note that the apparent stepping behavior is simply due to the resolution of reduced velocity chosen—only integer values are used. A smoother variation could be achieved, but it would not change the core point: the choice of coupling parameter significantly changes the lock-in range.Figure 11
Sensitivity of the VIV model of a circular cylinder by [32] to parameter values. Lock-in region: □, upper boundary; ∇, lower boundary; lock-in region (in green).
(a)
Sensitivity analysis of the damping parameterε with reduced velocity, Ur, for m∗=4
(b)
Sensitivity analysis of the coupling parameterA with reduced velocity, Ur, for m∗=4For every angle of attack, two estimates ofA are obtained: one for the postcoincidence boundary and the other for the precoincidence boundary. The proposed values of A is the mean of these estimates. The ε value used is 0.45 based on the mean of ε obtained with the initial harmonic balance approach for the sample cases of α=40° and α=140°. For comparison, the ε value used by Facchinetti et al. [32] for a circular cylinder is 0.3.


### 3.3. Results

The dynamic behavior of the wake system has been investigated for the whole range of angles of attack considered for the CFD analysis in Section2. The aim is to match the lock-in results obtained from the CFD model, in order to propose a set of parameters (i.e., ε,A) to be used for a preliminary VIV risk assessment at each angle of attack.Figure12 shows the variation of the coupling parameter with angle of attack over the whole range considered, for a constant value of pitching amplitude, αampl=7°. For comparison, the value proposed by Facchinetti et al. [32] for a circular cylinder (A=12) is indicated with a dash-dot line. It can be concluded that a generic value of A≈1 might be used in the range of (45°≤α≤140°), from a qualitative point of view. At incidences α=35°,40°,145°, higher values are present, correlated to the larger frequency lock-in range observed in the CFD study.Figure 12
Variation of the coupling factor,A, with angle of attack, for αampl=7°. Dash-dot line indicates the value used by [32] for a circular cylinder.Figure13 shows the proposed lock-in boundaries obtained from the wake oscillator model, overlapped with the lock-in region obtained from the CFD model (Figure 5), referring to αampl=7°. It can be seen that, qualitatively speaking, the VIV model reproduces the CFD results. Moreover, in the range of (45°≤α≤140°), a good match is achieved quantitatively. For incidences α=35°,40°,145°, the wake oscillator model overestimates the lock-in frequency range, due to a higher value of A.Figure 13
Proposed lock-in boundaries atαampl=7°: VIV model (in black line) and CFD model (in green).An analysis of the effect of the pitching amplitude has also been conducted, at the same incidences considered for the CFD model in Section2, shown in Figure 14 using the value of A from Figure 12. The proposed lock-in boundaries are overlapped with the lock-in region obtained from the CFD model (shown in detail in Figure 6). Results indicate that the wake oscillator model follows the expected behavior (i.e., lock-in range decreasing with amplitude reduced). However, while there is reasonable qualitative agreement, a dedicated tuning of the A value would be needed for every amplitude to achieve a close quantitative match with the CFD results.Figure 14
Proposed lock-in boundaries for pitching amplitude analysis at two incidences: VIV model (in black line) and CFD model (in green).
(a)
α=40°
(b)
α=140°


## 3.1. Methodology

In this section, a lumped parameter model for VIV of a pitching airfoil is presented. It is based on the wake oscillator proposed by Facchinetti et al. [32], which in turn is based on the classical van der Pol oscillator. A nonlinear model is chosen as the amplitude of oscillation in the simulations is by necessity nonzero. So even though limit cycle oscillations or force magnitudes are not of concern, the interaction between oscillation amplitude and resulting fluid force with lock-in must be accounted for. The wake is modeled as a nonlinear oscillator as follows:
(4)q¨+εΩfq2−1q̇+Ωf2q=F,where q is the dimensionless wake variable (overdots mean derivative with respect to the dimensionless time), ε is effectively a damping parameter, Ωf is the vortex shedding angular frequency, and F models the effect of the structure motion on the near wake by providing coupling with structural dynamics. The wake variable, q, is here associated with the fluctuating lift coefficient of the structure. The coupling term is based on the acceleration of the structure as recommended by Facchinetti et al. [32] as follows:
(5)F=A/Ds¨,where s¨ is a generic acceleration of the structure, A is the coupling parameter, and D is the characteristic length (e.g., cylinder diameter). Note that Facchinetti et al. [32] presented a coupled 2 degrees of freedom system, i.e., the structure was free to vibrate. In this paper, only the wake oscillator is considered, since the acceleration is imposed (the structure motion is prescribed as per equation 2).Remembering that the focus of this paper is the airfoil pitching motion, the wake oscillator equation has been reformulated as follows:(6)q¨+ε2πfv0q2−1q̇+2πfv02q=A/Lα¨,where q,ε,A are defined in the same way as in equation (4); fv0=StU∞/L is the static vortex shedding frequency, L is the characteristic length of the airfoil, and α¨ is the angular acceleration. Note that if the mean angle of attack varies, so does the characteristic length as it is defined to be the vertical projection of the chord with respect to the flow direction [33].The pitching motion acceleration can be derived by equation (2), yielding as follows:
(7)α¨t=−2πfp2αamplsin2πfpt.


## 3.2. Verification and Parameter Setting

Figure10 shows the lock-in region for a circular cylinder compared with the parameters specified by Facchinetti et al. [32]. The lock-in behavior obtained from the current model implementation compares well with the lock-in boundaries presented by Facchinetti et al. [32], offering confidence in the current implementation.Figure 10
Comparison of predicted lock-in boundaries for a circular cylinder:□, upper boundary (current implementation); ∇, lower boundary (current implementation); lock-in region (current implementation) (in green) and [32] (in black line).As the wake oscillator model is a semiempirical model, some of its parameters require empirical data to be quantified. In the current study, due to the apparent lack of experimental data on VIV of pitching airfoils at high angle of attack and high Reynolds number, the “empirical” data used are the lock-in boundaries presented in Figures5 and 6 above from CFD simulations. The parameters to be defined are ε and A. A harmonic balance method was used initially for two sample cases (α=40° and α=140°), using the harmonic coefficients obtained from the study of vortex shedding from the static NREL S809 [33]. However, the prediction of lock-in boundaries did not compare well with the dynamic CFD results. This suggests that maybe it is not possible to obtain a universal set of parameters for a VIV model of airfoil independent of angle of attack. An attempt has been made to renormalize the coupling factor based on the projected chord, in a similar manner to that used by the authors [33] to reduce the dependency of the Strouhal number on angle of attack. However, this was not successful.Another approach has been used: estimatingA for each angle of attack and keeping ε constant. The rationale for this approach is the sensitivity of the model predictions to these parameters. A sensitivity analysis of ε and A has been conducted using the model proposed by Facchinetti et al. [32] for a circular cylinder (m★=4). It can be clearly seen in Figure 11(a) that the reduced velocity range for lock-in is not affected by changing the ε value (at least for the considered variation of ±100%). In contrast, Figure 11(b) shows that a variation of A significantly affects the system behavior. Note that the apparent stepping behavior is simply due to the resolution of reduced velocity chosen—only integer values are used. A smoother variation could be achieved, but it would not change the core point: the choice of coupling parameter significantly changes the lock-in range.Figure 11
Sensitivity of the VIV model of a circular cylinder by [32] to parameter values. Lock-in region: □, upper boundary; ∇, lower boundary; lock-in region (in green).
(a)
Sensitivity analysis of the damping parameterε with reduced velocity, Ur, for m∗=4
(b)
Sensitivity analysis of the coupling parameterA with reduced velocity, Ur, for m∗=4For every angle of attack, two estimates ofA are obtained: one for the postcoincidence boundary and the other for the precoincidence boundary. The proposed values of A is the mean of these estimates. The ε value used is 0.45 based on the mean of ε obtained with the initial harmonic balance approach for the sample cases of α=40° and α=140°. For comparison, the ε value used by Facchinetti et al. [32] for a circular cylinder is 0.3.


## 3.3. Results

The dynamic behavior of the wake system has been investigated for the whole range of angles of attack considered for the CFD analysis in Section2. The aim is to match the lock-in results obtained from the CFD model, in order to propose a set of parameters (i.e., ε,A) to be used for a preliminary VIV risk assessment at each angle of attack.Figure12 shows the variation of the coupling parameter with angle of attack over the whole range considered, for a constant value of pitching amplitude, αampl=7°. For comparison, the value proposed by Facchinetti et al. [32] for a circular cylinder (A=12) is indicated with a dash-dot line. It can be concluded that a generic value of A≈1 might be used in the range of (45°≤α≤140°), from a qualitative point of view. At incidences α=35°,40°,145°, higher values are present, correlated to the larger frequency lock-in range observed in the CFD study.Figure 12
Variation of the coupling factor,A, with angle of attack, for αampl=7°. Dash-dot line indicates the value used by [32] for a circular cylinder.Figure13 shows the proposed lock-in boundaries obtained from the wake oscillator model, overlapped with the lock-in region obtained from the CFD model (Figure 5), referring to αampl=7°. It can be seen that, qualitatively speaking, the VIV model reproduces the CFD results. Moreover, in the range of (45°≤α≤140°), a good match is achieved quantitatively. For incidences α=35°,40°,145°, the wake oscillator model overestimates the lock-in frequency range, due to a higher value of A.Figure 13
Proposed lock-in boundaries atαampl=7°: VIV model (in black line) and CFD model (in green).An analysis of the effect of the pitching amplitude has also been conducted, at the same incidences considered for the CFD model in Section2, shown in Figure 14 using the value of A from Figure 12. The proposed lock-in boundaries are overlapped with the lock-in region obtained from the CFD model (shown in detail in Figure 6). Results indicate that the wake oscillator model follows the expected behavior (i.e., lock-in range decreasing with amplitude reduced). However, while there is reasonable qualitative agreement, a dedicated tuning of the A value would be needed for every amplitude to achieve a close quantitative match with the CFD results.Figure 14
Proposed lock-in boundaries for pitching amplitude analysis at two incidences: VIV model (in black line) and CFD model (in green).
(a)
α=40°
(b)
α=140°


## 4. Conclusions

It has been found that a NREL S809 wind turbine blade section subject to periodic pitching at high mean angles of attack may exhibit vortex shedding at a frequency other than the nominal vortex shedding frequency for a fixed geometry. This behavior is similar to that reported for other prismatic sections, most notably a circular cylinder. The extent of the frequency ratio range required for lock-in for a pitching NREL S809 airfoil has been estimated using uRANS over the incidence range at which periodic vortex shedding has been observed (35°≤α≤145°). It has been found that the wake is more susceptible to lock-in when the blade profile is at incidences presenting shallower characteristic length (α≈40°,140°). This is true whether the sharp trailing edge is upstream or downstream of the blunt leading edge.The effect of the oscillation amplitude has also been investigated. The behavior observed was expected, i.e., the extent of the lock-in frequency range decreases as the motion amplitude is reduced. Moreover, it has been observed that a sharp trailing edge in upstream position seems to slightly reduce the risk of lock-in.Finally, a lumped parameter VIV model (wake oscillator model) has been proposed for pitching airfoils. The lock-in boundaries qualitatively match the overall results obtained with the CFD model. It is necessary to tune the coupling factor for each angle of attack and perhaps each amplitude.It can be concluded that the wake oscillator model proposed by Facchinetti et al. [32] for a circular cylinder can be manipulated and adapted for an airfoil. However, it seems that a universal coupling parameter can not be used, as it is strictly dependent on the geometric configuration. This clearly limits the usefulness of the approach. Further work is need to investigate if the coupling parameter can be made independent of angle of attack, by choosing an appropriate characteristic length.Further work is needed to investigate the effect of other modes of oscillation, such as pure translational (i.e., transverse and streamwise) as well as coupled motions, and the effect of hysteresis on the wake behind the NREL S809.



---

*Source: 6919505-2019-03-11.xml*