Energy Harvesting Performance of Thick Oscillating Airfoils Using a Discrete Vortex Model

Research on the unsteady aerodynamics of oscillating airfoils has grown in recent years due to their importance in developing high-performance energy harvesters [1–10]. McKinney and De-Laurier first proposed the concept of flow energy harvesting using flapping foils [11]. The oscillating airfoil kinematics include large-amplitude heaving and pitching motions with the optimal range of $0.5c$ to $1.0c$ for heaving amplitude and 70 deg to 85 deg for pitching amplitude [12,13]. A phase shift between heaving and pitching velocities is also required to optimize the flow dynamics such that the largest lift forces occur during maximum heaving velocity to achieve the largest energy output rate [14].

The oscillating foils for energy harvesting operate at extremely large angles of attack, resulting in periodic flow separation with associated dynamic stall characteristics. During dynamic stall, the flow separates near the leading edge of the airfoil surface and subsequently rolls into leading-edge vortices (LEVs). LEVs produce regions of low pressure on the airfoil and generate a large suction force as long as they remain attached to the foil surface [15]. Contrary to traditional wind and hydroturbines, where the flow generally remains attached to the foil surface, to achieve high energy harvesting efficiency, the LEV plays a crucial role in augmenting lift force and enhancing the efficiency of flapping foil energy harvesters. Hence, predicting the onset of separation is critical in understanding the flow physics and developing unsteady low order models.

In addition to the heaving and pitching amplitudes, h₀ and θ₀, a few other parameters are identified to influence energy harvesting, including the reduced frequency, $k=fc/U∞$⁠, the chord-based Reynolds number, $Re=U∞c/ν$⁠, and the airfoil thickness. The optimum range of reduced frequency is $0.12−0.18$ [1,14,16]. The power output tends to be relatively insensitive to the heaving and pitching amplitudes when operating near the optimal reduced frequencies [5,6,9]. The optimum relative phase between heaving and pitching is reported to be near 90 deg [5,9,12,16]. Among the parameters mentioned above, the effect of airfoil thickness on leading-edge separation has been relatively unexplored, which is addressed in this paper.

Several studies, including higher-order computational simulations [1,2,16–18] and experiments [9,19–24], have contributed to understanding the high amplitude oscillating flow dynamics and energy harvesting performance. Unsteady thin airfoil theories assume the flow remains attached and apply to small amplitude kinematics [25,26]. These methods have been extended to solve three-dimensional problems [27,28]. Furthermore, Ramesh et al. [29] proposed the use of a leading-edge suction parameter (LESP), equivalent to the first Fourier coefficient in unsteady thin airfoil theory, to predict the onset of LEV formation. This method relies on CFD wall shear values as input to determine a critical LESP value to predict the onset of leading-edge flow separation.

Studying the effect of airfoil thickness is particularly important for the wind energy industry due to the high thickness-to-chord ratio of the airfoils used in wind turbines and regularly experiencing dynamic stall due to unsteadiness induced by rotation and operation in the atmospheric boundary layer. Sharma and Visbal [30] used high-resolution, wall-resolved large eddy simulations to study the boundary layer flow physics and its effect on the onset of dynamic stall. They showed that airfoil thickness could alter the stall onset mechanism. A similar observation has been reported by Heine et al. [31] when using passive disturbance generators to delay dynamic stall on the OA209 airfoil. Mulleners and Raffel [32] used particle image velocimetry (PIV) and direct pressure measurements to investigate the stall onset in an oscillating OA209 airfoil. They defined the rate of change of angle of attack at the static-stall angle as “effective unsteadiness” and observed that increasing this parameter facilitated dynamic-stall onset. Müller-Vahl et al. [33] experimentally investigated dynamic stall mitigation using active flow control techniques for NACA0018. However, they did not explore the effect of airfoil thickness on boundary layer physics. Larsen et al. [34] investigated the effects of airfoil thickness on dynamic stall using a semi-analytical model developed specifically for wind turbine airfoils with a high thickness-to-chord ratio and observed minor changes to the dynamic stall predictions due to airfoil thickness. Overall, limited work has been done on developing reduced-order models for thick foils. This highlights the need to develop reasonably accurate low order models to better understand the key components that impact the flow separation and the resultant energy harvesting power output.

The growth of the leading edge vortex, or vortex generation in general, has been studied to determine local scaling for the vortex growth dynamics. For example, the optimal vortex formation number has been defined based on the local shear velocity and appropriate length scale [35]. This is based on the idea of a universal vortex formation time [36]. Experimental studies on LEV growth for an oscillating foil show that during early times, the transient growth collapses for reduced frequencies from 0.06 to 0.16 using the time scale $c/U¯shear$⁠, where $U¯shear$ is the cycle-averaged leading-edge shear velocity [37]. This idea suggests defining a timescale based on the shear velocity for the leading-edge separation, which will be discussed later in this paper.

This study aims to better understand the effect of airfoil thickness on the leading-edge separation timing and location along the airfoil chord. A discrete vortex model (DVM) is developed for oscillating airfoil energy harvesters using the leading-edge separation criterion proposed in our previous study, and a trailing-edge correction to account for trailing-edge reversed flow. The model is applied to flow past two-dimensional NACA00XX series airfoils with reduced frequencies ranging from k = 0.06 to k = 0.14. A panel method is used to describe the airfoil geometry. Details of the panel method to formulate the general solution have been provided [38]. Further details of the model, the applied corrections, and the implementation are provided in the Methodology section. The unsteady lift data is compared against CFD results for NACA0015 airfoil. Lastly, the transient and cycle-averaged power output for the NACA00XX airfoil series are calculated, and the energy harvesting performance of different geometries is compared.

The paper is organized as follows. First, an overview of the DVM methodology is presented. Second, the separation criteria results are scaled based on the leading-edge shear velocity and radius of curvature to collapse data over the studied range of reduced frequencies. Third, modifying convection velocity based on the reduced frequency shows some improvement in capturing the LEV dynamics. Fourth, the results of model predictions of energy harvesting performance are presented. The overall goal of this study is to incorporate both timing and location into the leading-edge separation criterion, which can be used for a variety of low order models of thick oscillating foils.

This study evaluates five airfoils with different leading-edge geometries listed in Table 1. All the airfoils have a chord length of $c=0.15(m)$ and have a pitching axis about the midchord of the airfoil. The NACA airfoils have no camber, the foil is symmetric across the x-axis, and the thickest part of the airfoil occurs at 30% of the chord length from the tip of the leading-edge. As part of the nomenclature, the NACA00XX has a maximum thickness of XX% of the chord length. The radius of curvature is measured at the leading-edge.

The airfoils are modeled using inviscid theory, based on the surface distribution of single vortex points on panels with vortices being shed near the leading-edge and at the trailing-edge. Compared to a more analytical method such as thin airfoil theory, this method is more flexible to model complexities such as different airfoil cambers and relative motions at the leading and trailing-edges. Figure 1 shows the panel discretization of the airfoil surface and the motion kinematics. The airfoils are discretized into flat panels around the surface. The number of panels shown in Fig. 1 is not the actual number and only represents how the airfoil is discretized. The number of panels was tested to show that further spatial resolution with increased number of panels does not change the results. The number of panels was increased until a converged solution was found for all values of the reduced frequency, as is shown in Fig. 2. Results show that increasing the panel number from 100 to 200 significantly changes the lift coefficient during most parts of the cycle. However, changing the panel number from 200 to 400 resulted in very minimal change in results during the entire cycle. Consequently, 200 panels were used for all simulations for all reduced frequencies.

where n refers to the normal direction on each panel i, $h˙n,i$ is the normal heaving velocity, $(θṗ× ri)n$ is the normal pitching velocity, and $(∂ΦLEV∂n)$ and $(∂ΦTEV∂n)$ are the normal velocities induced by the previously shed vortices near the leading-edge and at the trailing-edge, respectively, and $Φ$ is the perturbation potential from the shed leading-edge and trailing-edge vortices. The term $∑AijΓj$ is the bound vortices velocity contribution on the foil where A_ij is the influence coefficient for vortex j on collocation point of panel i, and Γ_j is the bound vortex strength which is being calculated for panel i. More details on the calculation of the influence coefficients are presented in Ref. [38].

In this study, each cycle included 1000 time steps. At each time-step, a vortex is shed from the trailing-edge; however, vortex shedding from the leading-edge depends on the leading-edge separation criterion, presented later. The shed vortices are positioned at 1/3 of the distance of the previously shed vortex to approximate the shape of the shear layer [39]. The vortex strength is determined iteratively using the 1-D Newton's method coupled with Kelvin's theorem given in Eq. (4) [38].

In a previous study [40], the authors proposed a criterion for the initiation of leading-edge separation for a thin airfoil based on the wall shear stress distribution near the leading-edge. This criterion requires the local shear stress to pass through zero (a positive-to-negative transition) during the airfoil oscillation. As mentioned earlier, LEV shedding only happens if flow separation triggers near the leading-edge. At this occurrence, a measure of the value of the leading-edge shear velocity is set as the critical parameter. LEV shedding is only active at the time steps when the leading-edge shear velocity exceeds the critical value. This study uses a similar approach to find the leading-edge separation instant. In addition, the exact leading-edge separation location (zero wall-stress location near the leading-edge) is found from the CFD simulations, and point vortices are shed from this location rather than the airfoil leading-edge.

In general, energy harvesting occurs due to both heaving and pitching contributions; however, the pitching contribution to the power output has been shown to be very small for the range of reduced frequencies in this study and for pitching about the midchord [1,12,21].

where $h˙$ is the linear heaving velocity, F_Y is the lift force in the heaving direction, and b is the airfoil span.

where the subscript LE refers to the pitching contribution evaluated at the leading-edge. This is not strictly an effective angle of attack since the pitching contribution is only valid at the leading-edge. However, it does provide a measure that relates to the leading-edge shear velocity by accounting for both heaving and pitching contributions at the leading-edge.

where $CN,non$ is the noncirculatory component of the lift force normal to the chord. A comparison of using $αeff,LE$ versus α_eff to find the fictitious separation point (f^sep) along the chord is presented in Fig. 3 for NACA0021. The biggest difference is seen at t/T between 0.1 and 0.25, where $fLEsep$ is smaller and better accounts for the over-prediction of the transient lift force, as is shown later in the force results. Furthermore, as the reduced frequency, k, increases, this difference increases, as is shown in the Results section, to better predict the lift force transient predictions.

This study uses fluent 2019 R2 as a high-resolution two-dimensional simulation solver and is used to determine the separation time and location conditions that are then applied to the DVM model. The CFD model is based on RANS modeling with a dynamic mesh of the turbulent flow field to account for the oscillating motion of the foil. The turbulence model used for this study is the Spalart–Allmaras (SA) model, which has been designed specifically for external aerodynamic flows. This model has also been tested and used for previous oscillating foil simulations and found to provide excellent agreement with experimental results [12,40]. This is not to say that alternative turbulence models may be helpful in resolving some of the flow characteristics.

Uniform time steps of 5000th of a cycle were used (the time interval changes as the oscillating frequency changes). This time was found to converge results based on evaluating results using time steps ranging from 1000 to 5000 per cycle. A multilayer dynamic grid allowing a mesh to remain in a fixed reference frame relative to the airfoil is illustrated in Fig. 4. The grid contained 330,000 cells. Multiple cycles were simulated with the initial cycles not included in the results to mitigate initial startup variations. The total computational run time was approximately 70 h per case. Grid independence studies were conducted based on the calculated coefficient of lift. First, the required time-step size for convergence was found, as stated above, which was then followed by increasing the cell count within the mesh from approximately 100,000 to 330,000 cells until the results showed grid convergence [49]. Validation was obtained in a previously published study by the authors [40] comparing transient cycle force results for the identical airfoil geometry and conditions [12].

Wall shear stress along the chord was evaluated to determine the time and location of separation, which is defined as the critical condition of the sign reversal of the local wall shear stress. This is discussed in greater detail in Refs. [40] and [49]. NACA00XX geometries are studies for Reynolds number range of Re = 15,000–40,000 and reduced frequency range of $k=0.06−0.14$⁠. The dynamic mesh consists of two mesh regions connected through a nonconforming sliding interface. The first mesh has a fixed frame of reference, while the second (smaller) mesh region is enclosed and heaves vertically and pitches rotationally at the same rate as the airfoil. Additional details on the CFD model, including discretization and validation, are presented in Ref. [49].

This velocity component, which is normal to the chord, represents the tendency for the flow to accelerate around the leading edge; the greater the acceleration, the greater the local suction pressure occurring at the leading edge. Considering the given kinematic motion of heaving and pitching, the linear heaving velocity and the angular pitching velocity both contribute to the shear velocity. As the reduced frequency increases, by either the freestream velocity decreasing or the oscillating frequency increasing, the local shear velocity will decrease.

It is shown in the authors' thin flat plate study [40], the critical separation time for a thin flat airfoil can be scaled using the time scale $τν=c/U¯shear$⁠. However, attempts to collapse the critical separation time data for the NACA thicker airfoils indicate that the temporal scaling used for the thin flat airfoil has a strong dependence on the foil thickness when using the chord length, c, as the length scale. Using an alternate timescale based on the radius of curvature (r_c) as $rc/U¯shear$ better collapses the data, as shown in Fig. 5(a). In this figure is shown the time dependence of $Ushear/U∞$ for a range of reduced frequencies, k, and foil radius of curvature associated with the different foil geometries during the cycle. Note that the negative velocities indicate the shear flow is downward when the heaving and pitching result in a downward leading-edge motion, a consequence of the pitching angle effect on the freestream as indicated in Eq. (15). When the shear velocity is zero, there is a balance of contributions from the freestream and the heaving and pitching motions. The critical separation times, corresponding to the time of zero wall-shear stress during the cycle, are shown with black markers on each curve. The critical times are shown to occur at earlier scaled times for larger radius of curvature leading edge foils.

The CFD results of the separation location based on the wall shear stress criteria are shown in Fig. 6. As the radius of curvature increases, the flow remains attached for a longer distance along the leading-edge; hence the leading-edge separation location shifts further from the leading-edge. While the reduced frequency impacts the time of separation (as shown in Fig. 5(a), the relative speed of the airfoil motion does not impact the location of the separation point near the leading-edge, and the separation location appears to be independent of the reduced frequency while linearly increasing with the radius of curvature.

Dynamic stall conditions have been shown to be influenced by airfoil thickness [30]. One possible reason for this is the timing associated with the motion of the leading-edge vortex from initial formation through its final separation from the airfoil surface. In addition, the LEV convection has been shown to strongly influence the transient lift force on the airfoil [37,50]. The CFD results illustrating the LEV growth slightly after initial formation for three different values of k are shown in Fig. 7. The model results are also shown in the center series of images. It can be seen that the vortex growth is not modeled well, even though the CFD indicated separation time has been applied. This causes the force predictions of the model to be in error. To remedy this in the model, adjustments were made to the convection velocity applied to the discrete vertices.

The growth and convection of the LEV are strongly affected by the reduced frequency [37]. An empirical approach was used to better match the strength and size of the LEV to the CFD results during the cycle. A series of attempts were made to adjust the local convection velocity to better match the observed LEV growth. Once the LEV growth and convection were matched, the transient force data were compared between the DVM and CFD results. Table 2 shows the adjustment made to the convection velocity as a function of the reduced frequency, indicating the need for an increasing convection velocity away from the airfoil (chord-normal direction) with increasing k. The improvement to the LEV prediction is shown in Fig. 7 by comparing the DVM adjusted results in the right side column to the CFD results in the left hand column.

The size, shape, and location of the LEV are much better represented in the DVM results compared with the CFD generated vorticity fields when using the convection velocity adjustments. Also, the convection velocity modified force predictions shown in Fig. 8 show a significant improvement. As shown in Fig. 7, the vorticity distribution near the trailing-edge does not compare well with that obtained in the CFD. This discrepancy is evidence that the trailing-edge correction is in fact needed to obtain improved force estimates. It should be noted that the correction used in this study adjusts the force prediction “globally” to account for the separation and consequently does not specifically modify the vorticity distribution in the model. It is recommended that further studies address a more local correction directly to the vorticity field to better capture the trailing-edge effects on the vorticity and force distributions. This may improve the ability to better model the effects of separation at higher reduced frequencies.

Figure 8 shows the transient lift coefficient profiles for NACA0015 generated from the low order model and compared with CFD results for the same condition. The leading-edge separation is based on the use of the critical value of $Ushear/U∞$ at the associated critical time of $tcrit/τcrit≈−8$ (shown earlier in Fig. 9). The leading-edge separation location is also set based on the zero wall-shear location along the chord identified from CFD, as shown in Fig. 6 for the range of k values in this study. Overall the DVM results compare reasonably well with the CFD results with transient lift force coefficient root-mean-square deviations under 10% for the lower k values and approximately 15–20% for the larger k values. The low order model captures the same trends, despite some local quantitative differences that depend on the reduced frequency. The general description of the occurrence of the primary and secondary peaks is dependent on the reduced frequency [37] and is captured well in these results; however, for k = 0.12 and k = 0.14, the secondary peak is not as pronounced, and the model is not able to capture the peaks at precisely the same times as indicated in the CFD results. Nonetheless, the model results are still very close. The trailing-edge separation correction is also applied to the force calculation (Impulse-Corrected). The correction uses $αeff,LE$ instead of α_eff and is shown to reduce any difference between the model and CFD results, which is more noticeable in higher reduced frequencies. The most significant improvement with the correction is seen for k = 0.12 and k = 0.14. Contrary to using α_eff, which would lead to an under-estimation of the primary and secondary peaks for reduced frequencies greater than k = 0.12 [40], using $αeff,LE$ takes into account the dependency of the empirical separation correction as a function of reduced frequency. Overall, the low order model, with the application of the leading-edge separation criteria for timing and location, captures the general cycle variations and performance characteristics well. Results show that agreement between the CFD and DVM decreases somewhat for the largest reduced frequency cases in this study. In particular, the transient force results for k = 0.12 and k = 0.14 shown in Fig. 8 indicate the development of the secondary peak happens around the half-cycle (near $t/T=0.48$⁠), which is a little later than CFD results. If the viscous effects become dominant on the vortex dynamics near the surface, it is expected that the DVM would not accurately predict the resultant force and yield larger loads, as is shown in Fig. 8. Despite improvements with adjusting the convection velocity for higher reduced frequencies, this speculation requires further investigation and implementing other improvements at the very high reduced frequencies.

The results of the calculated transient power coefficients using the low order model are shown in Fig. 10 for both lower (Fig. 10(a), k = 0.08) and higher (Fig. 10(b), k = 0.14) reduced frequencies. The transient power is based on the heaving contribution identified in Eq. (9). Compared to the flat plate, all NACA geometries have higher power output. Furthermore, as the airfoil thickness increases, the transient power peak moves later in the cycle, and its amplitude also increases, regardless of reduced frequency. Compared to k = 0.08, the power output is significantly larger for k = 0.14 for all of the geometries. This is consistent with the results reported in the literature [1,21,51] that suggests the optimal reduced frequency to be around $k=0.12−0.15$⁠.

As shown in Fig. 11, the cycle-averaged power coefficient increases with increasing reduced frequency, which agrees well with the results reported in the literature [1,12]. Improved energy extraction performance for thicker airfoils has been previously reported for high Renumbers (similar to the range used in this study) [17]. These DVM results also show that the trend persists for all of the thicker foil results presented here, with the thickest foils performing best. For the thicker airfoils, the average power output is higher, with greater improvements seen for the higher reduced frequencies, where the power output is highest. For instance, NACA0021 gives $10−20%$ higher heaving power coefficients compared to the other NACA geometries and the flat plate. In general, increasing the radius of curvature greatly affects leading-edge separation by reducing the acceleration around the leading edge and suffering from a lower suction pressure near the leading edge. Since the airfoil undergoes high amplitude heaving and pitching motions, separation and dynamic stall are still anticipated for the relatively thick airfoils. Furthermore, LEV strength and convection after detachment are key in power output results. LEV detached convection is determined by the external flow experienced by the airfoil and is influenced by the TEV dynamics associated with the vortex roll-up along the suction side of the foil [52]. Airfoil thickness also affects the viscous effects near the foil surface; hence inviscid models need to account for the relevant effects to simulate the power output more accurately.

This study investigates the effect of airfoil thickness on the onset of leading-edge separation. Leading-edge separation timing and location are identified using CFD simulations for five airfoil geometries, including a flat plate and four NACA00XX airfoils. All the airfoils undergo high amplitude heaving and pitching motions with reduced frequencies ranging from k = 0.06 to k = 0.14. This separation condition is utilized to develop low order models and predicts the aerodynamic force and power output. This critical separation time is shown to collapse to a single nondimensional value using a time scale based on the leading-edge shear velocity and the radius of curvature. This time scale and critical time provide a means to predict separation based on the foil geometry and motion kinematics in a given freestream flow and is advantageous for providing separation criteria for a wide range of low order models. A discrete vortex model is developed to predict transient lift forces for the range of reduced frequencies and to evaluate the effect of airfoil thickness on the energy harvesting performance. The point vortex convection velocities are also modified based on the reduced frequency and enhance the shape of LEV. A trailing-edge separation correction is also applied to the model results, which uses the leading edge effective angle of attack, and decreases the slight difference between the model prediction of the lift force and CFD, especially around the peaks. This study further shows the effect of airfoil thickness on the leading-edge separation during dynamic stalls. The use of a low order DVM with corrections for leading-edge separation time and location, as well as a modified trailing-edge separation criteria using the pitching affected angle of attack, and introducing a reduced frequency modified convection velocity is shown to provide very good predictions of energy harvesting of an oscillating foil with nonzero thickness.

• National Science Foundation CBET, Directorate for Engineering (Award No. 1804964; Funder ID: 10.13039/100000084).

A_ij =

influence coefficient for vortex j on collocation point i

c =

chord length, m

$CN, circ$ =

unadjusted circulatory normal force coefficient

$CN,non$ =

non-circulatory normal force coefficient

$CN,sep$ =

adjusted circulatory normal force coefficient

C_s =

unadjusted suction force coefficient

$Cs,sep$ =

adjusted suction force coefficient

C_Y =

lift coefficient

$CP,h$ =

heaving power output coefficient

DVM =

discrete vortex model

f =

heaving/pitching frequency, Hz

$f0sep$ =

steady-state separation point

f_sep =

fictitious separation point

$FC,Y$ =

circulatory lift force, N

$FNC,Y$ =

circulatory lift force, N

$Fs,Y$ =

suction force, N

F_Y =

lift force, N

h =

heaving motion, m

h₀ =

heaving amplitude, m

k =

reduced frequency

LEV =

leading-edge vortex

P_h =

instantaneous heaving power, W

r_c =

leading-edge radius of curvature, m

Re =

Reynolds number

t =

time, s

TEV =

trailing-edge vortex

U_shear =

leading-edge shear velocity, m/s

$U∞$ =

freestream velocity, m/s

α₁ =

critical angle of static separation point, deg

α_eff =

effective angle of attack, deg

$αeff,LE$ =

leading-edge effective angle of attack, deg

θ₀ =

pitching amplitude, deg

θ_p =

angular pitching velocity, rad

ν =

kinematic viscosity, m²/s

ρ =

density, kg/m³

$τ1,τ2$ =

empirical relaxation time coefficients

Γ =

circulation, m²/s

$Φ$ =

velocity potential