Abstract

The energy harvesting performance of thick oscillating airfoils is predicted using an inviscid discrete vortex model (DVM). National Advisory Committee for Aeronautics (NACA) airfoils with different leading-edge geometries are modeled that undergo sinusoidal heaving and pitching with reduced frequencies, k=fc/U, in the range 0.060.14, where f is the heaving frequency of the foil, c is the chord length, and U is the freestream velocity. The airfoil pitches about the midchord with heaving and pitching amplitudes of h0=0.5c and θ0=70deg, respectively, known to be in the range of peak energy harvesting efficiencies. A vortex shedding initiation criteria is proposed based on the transient local wall stress distribution determined from computational fluid dynamics (CFD) simulations and incorporates both timing and location of leading-edge separation. The scaled shedding times are shown to be predicted over the range of reduced frequencies using a timescale based on the leading-edge shear velocity and radius of curvature. The convection velocity of the shed vortices is also modeled based on the reduced frequency to better capture the dynamics of the leading-edge vortex. An empirical trailing-edge separation correction is applied to the transient force results using the effective angle of attack modified to include the pitching component. Impulse theory is applied to the DVM to calculate the transient lift force and compares well with the CFD simulations. Results show that the power output increases with increasing airfoil thickness and is most notable at higher reduced frequencies where the power output efficiency is highest.

Introduction

Research on the unsteady aerodynamics of oscillating airfoils has grown in recent years due to their importance in developing high-performance energy harvesters [110]. 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, h0 and θ0, a few other parameters are identified to influence energy harvesting, including the reduced frequency, k=fc/U, the chord-based Reynolds number, Re=Uc/ν, and the airfoil thickness. The optimum range of reduced frequency is 0.120.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,1618] and experiments [9,1924], 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.

Methodology

Airfoil Geometry Definition.

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.

Table 1

Airfoil geometries, all airfoils are symmetric with a chord length of 150 mm

Airfoil nameAirfoil shapeRadiusofCurvaturechord
Flat plate0.0036
NACA00120.0120
NACA00150.0188
NACA00180.0270
NACA00210.0368
Airfoil nameAirfoil shapeRadiusofCurvaturechord
Flat plate0.0036
NACA00120.0120
NACA00150.0188
NACA00180.0270
NACA00210.0368

Low Order Model Formulation.

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.

Fig. 1
Inertial frame of reference (X, Y), chord frame of reference (x, y), and airfoil with discretization into N panels during downstroke. θp˙ and h˙ are angular pitching velocity and linear heaving velocity, respectively. Panel “i  ” is magnified as an example to show the placement of collocation and vortex points on each panel.
Fig. 1
Inertial frame of reference (X, Y), chord frame of reference (x, y), and airfoil with discretization into N panels during downstroke. θp˙ and h˙ are angular pitching velocity and linear heaving velocity, respectively. Panel “i  ” is magnified as an example to show the placement of collocation and vortex points on each panel.
Close modal
Fig. 2
Effect of the number of panels used in the low order model. Lift coefficient results are for NACA0015, k = 0.06 during downstroke.
Fig. 2
Effect of the number of panels used in the low order model. Lift coefficient results are for NACA0015, k = 0.06 during downstroke.
Close modal
The heaving and pitching motions with a 90 deg phase shift are defined according to the following:
h(t)=h0cos(2πft)θ(t)=θ0sin(2πft)
(1)
A phase shift of 90 deg between heaving and pitching optimizes the flow dynamics to take advantage of the dynamic stall. This way, the largest lift forces occur near the maximum heaving velocity to achieve the largest energy output rate [14]. The airflow is approaching from the left side with a uniform velocity of U. The airfoil heaves normal to the incoming flow with a linear heaving velocity of h˙, and simultaneously pitches about the midchord with angular pitching velocity of θ˙p, both given in the following equation:
h˙=2πfh0sin(2πft)θp˙=2πfθ0cos(2πft)
(2)
The point vortices and collocation points are placed on each panel at the 1/4 and 3/4 panel length, respectively. Boundary conditions include (i) surface impermeability and (ii) Kutta condition. Surface impermeability is enforced at each collocation point of each panel and leads to N linear equations shown in Eq. (3). The Kutta condition states zero vorticity at the trailing-edge, which is implemented by adding a wake panel along the chord at the trailing-edge, and setting the pressure difference across this panel to zero. The resulting equation becomes
AijΓj=U,n,i+h˙n,i+(θṗ×ri)n(ΦLEVn)(ΦTEVn)
(3)

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 (ΦLEVn) and (ΦTEVn) 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 Aij 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 order to satisfy the boundary conditions in the unsteady flow, the bound circulation changes in a way to be balanced by the shed vortices from the leading and trailing-edges. This concept is expressed in Kelvin's theorem and is presented in the following equation:
ΣiΓi+ΣΓTEV+ΣΓLEV=0
(4)

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].

Leading-Edge Separation Timing and Location.

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.

Depending on the established leading-edge vortex shedding criterion, if a vortex was shed near the leading-edge at the previous time-step, the current shed vortex is positioned similarly to the shed vortices from the trailing-edge presented above. Otherwise, the vortex is placed based on the induced velocity normal to the leading-edge from the freestream, heave, and pitch contributions. Depending on the reduced frequency, this induced velocity is modified to better represent the viscous effects near the airfoil surface, and is discussed in the Results. To solve for the strengths of the vortices shed near the leading-edge and at the trailing-edge, a 2-D Newton's method is used with the updated expression for Kelvin's theorem shown in the following equation:
ΣiΓi+Σj=1m1Γj,TEV+Γm,k+Σj=1q1Γj,LEV+Γq,k=0
(5)
where subscript k refers to the iteration number, m is the vortex shed at the trailing-edge, and q is the vortex shed near the leading-edge. As mentioned earlier, all shed vortices are convected with the inviscid velocity times an adjustment factor (shown in the Results section) which is a function of the reduced frequency. Therefore, the induced velocity from point vortices on each panel as well as vortices shed from near the leading-edge and trailing-edge is calculated using influence coefficients and added to the freestream velocity. Then, this adjusted inviscid velocity is multiplied by the time-step to find the position increment. To avoid artificially large velocities on the point vortices induced by the others in their proximity, vortex blobs with finite core radii are defined to give more realistic induced velocities as well as vorticity distributions [41]. In the present work, the core radius is taken to be 1.3 times the average spacing between the shed vortices [42] as presented in the following equation:
rcore=1.3UΔt
(6)

Transient Force and Power Calculation.

The unsteady lift force on the airfoil is calculated using the momentum impulse formulation [43] shown in Eq. (7)
FY=ρΓb(Udxbdt)ρxbdΓbdt+ρΓLEV(UdxLEVdt)ρxLEVdΓLEVdt+ρΓTEV(UdxTEVdt)ρxTEVdΓTEVdt
(7)
where b refers to the bound circulation and xi is the X-position of each vortex in the inertial frame. The circulatory plus suction force and the noncirculatory force in the Y direction are given in the following equation:
FC,Y+Fs,Y=ρΓb(Udxbdt)+ρΓLEV(UdxLEVdt)+ρΓTEV(UdxTEVdt)FNC,Y=ρxbdΓbdtρxLEVdΓLEVdtρxTEVdΓTEVdt
(8)

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].

The heaving power output (Ph), as well as the power output coefficient (Cp,h), are defined in the following equation:
Ph=h˙FYCP,h=Ph12ρU3cb
(9)

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

Trailing-Edge Separation Correction.

At higher angles of attack, flow separation near the airfoil trailing-edge most likely occurs, resulting in reduced lift forces. Since the low order model does not forecast this, a correction can be applied based on a modified Kirchhoff flow approximation to overcome this over-prediction of lift force [4446]. This correction calculates a fictitious separation point along the chord given in Eq. (10) and adjusts the circulatory part of the force to account for the bound circulation losses
f0sep={10.3exp(|αeff,LE|α1S1)|αeff,LE|<α10.04+0.66exp(α1|αeff,LE|S2)|αeff,LE|α1
(10)
where S1=3.0,S2=2.3,α1=15.25deg define the static stall characteristics and are determined from experimental data [47]. Originally, this equation is based on the effective angle of attack, αeff, for a heaving airfoil, which does not include pitching motion. To remedy this, a new effective angle is determined, which includes the pitching effect on the relative flow at the leading edge and is designated as
αeff,LE=θparctan(h˙+c2(θp˙cosθp)Uc2(θp˙cosθp))
(11)

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.

In order to consider the separation point movement along the airfoil and capture the transient dynamic effects for unsteady flow conditions, a first order differential equation is applied as follows [45,48]:
τ1dfsepdt+fsep=f0sep(αeff,LEτ2α˙eff,LE)
(12)
where τ1=0.52CU and τ2=4.5CU are empirical relaxation time constants. This empirical correction is applied to the circulatory and suction components of the lift force in the direction normal to the chord as follows:
CN,sep=CN,circ(1+fsep122)2Cs,sep=Csfsep12
(13)
Finally, the corrected lift force is then calculated in the following equation:
CY=(CN,sep+CN,non)cosθp+Cs,sepsinθp
(14)

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 (fsep) 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.

Fig. 3
Comparison of fictitious separation point used in the trailing-edge separation correction during the cycle using αeff (fsep) versus using αeff,LE (fsep−LE) for NACA0021 at reduced frequencies of (a) k = 0.08, (b) k = 0.12, and (c) k = 0.14
Fig. 3
Comparison of fictitious separation point used in the trailing-edge separation correction during the cycle using αeff (fsep) versus using αeff,LE (fsep−LE) for NACA0021 at reduced frequencies of (a) k = 0.08, (b) k = 0.12, and (c) k = 0.14
Close modal

Computational Fluid Dynamics Formulation.

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].

Fig. 4
Example of dynamic double mesh used to build the CFD model, showing the sliding interface that enables a mesh to remain in a fixed reference frame relative to the airfoil
Fig. 4
Example of dynamic double mesh used to build the CFD model, showing the sliding interface that enables a mesh to remain in a fixed reference frame relative to the airfoil
Close modal

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.060.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].

Results

Effect of Airfoil Thickness and Reduced Frequency on Separation Onset.

To help interpret the results, a scaling velocity is introduced that better defines the flow in the region near the leading edge. The shear velocity, Ushear, is defined as the relative velocity component normal to the airfoil chord line at the leading edge [13]. This velocity represents the induced velocity at the leading-edge due to the airfoil motion and the freestream velocity and is given as
UshearU=sin(θp)h˙Ucos(θp)θ˙c2U
(15)

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 (rc) 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.

Fig. 5
(a) Relative shear velocity versus the nondimensional time using the radius of curvature and mean shear velocity; black markers indicate the critical separation times. (b) Critical shear velocity of all tested NACA airfoils versus the critical scaled time with the indicated linear fit.
Fig. 5
(a) Relative shear velocity versus the nondimensional time using the radius of curvature and mean shear velocity; black markers indicate the critical separation times. (b) Critical shear velocity of all tested NACA airfoils versus the critical scaled time with the indicated linear fit.
Close modal
The results for the scaled critical shear velocity versus the scaled critical times of separation for the range of NACA leading-edge radii of curvature are shown in Fig. 5(b). These data follow a nearly linear trend of decreasing shear velocity amplitude, tending toward zero, as the radius of curvature decreases. The linear fit to the critical points shown in Fig. 5(b) is presented in the following equation:
Ushear,critU=0.00075(tcritU¯shearrc)0.50
(16)
The effect of varying the reduced frequency, k, on the critical times when scaled with a convection velocity of rc/U¯shear as shown in Eq. (16), results in a collapse of data associated with each reduced frequency. Since the value of U¯shear/U depends only on the motion kinematics, it can be determined for a given heaving and pitching amplitude and phase of pitching. In this study, the heaving and pitching amplitudes and the phase difference are all held constant. So for this study, U¯shear/U is constant for all k values (approximately 0.6). Using these scaling arguments, it is possible to rescale the critical separation times using the critical shear velocity, which can be expressed as
tcritτcrit=(tcritU¯shearrc)(Ushear,critU)=(tcritUshear,critrc)(U¯shearU)
(17)

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.

Fig. 6
Location along the leading edge at which the wall shear stress changes direction versus the leading edge radius of curvature. The separation location appears to be independent of the reduced frequency, k.
Fig. 6
Location along the leading edge at which the wall shear stress changes direction versus the leading edge radius of curvature. The separation location appears to be independent of the reduced frequency, k.
Close modal

Discrete Vortex Model Convection Velocity.

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.

Fig. 7
Effect of modifying convective velocity for NACA0015 airfoil at reduced frequencies of (a) k = 0.06, (b) k = 0.10, and (c) k = 0.14; the vorticity fields are from the instant slightly after leading-edge separation
Fig. 7
Effect of modifying convective velocity for NACA0015 airfoil at reduced frequencies of (a) k = 0.06, (b) k = 0.10, and (c) k = 0.14; the vorticity fields are from the instant slightly after leading-edge separation
Close modal

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.

Table 2

Adjustment factor multiplied by convection velocity of point vortices versus reduced frequency

Ky-Component convection velocity multiplier
0.061.02
0.081.05
0.101.10
0.121.15
0.141.20
Ky-Component convection velocity multiplier
0.061.02
0.081.05
0.101.10
0.121.15
0.141.20

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.

Fig. 8
Comparison of the transient lift coefficient (CL) for NACA0015 airfoil predicted by the low order model using the impulse method compared with CFD results for (a) k = 0.06, (b) k = 0.08, (c) k = 0.10, (d) k = 0.12, and (e)k = 0.14; “Corrected” refers to the use of the trailing-edge separation correction
Fig. 8
Comparison of the transient lift coefficient (CL) for NACA0015 airfoil predicted by the low order model using the impulse method compared with CFD results for (a) k = 0.06, (b) k = 0.08, (c) k = 0.10, (d) k = 0.12, and (e)k = 0.14; “Corrected” refers to the use of the trailing-edge separation correction
Close modal

Aerodynamic Lift Force Results.

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/τcrit8 (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.

Fig. 9
Critical separation times scaled using Eq. (17) versus reduced frequency, k
Fig. 9
Critical separation times scaled using Eq. (17) versus reduced frequency, k
Close modal

Energy Harvesting Performance of NACA Airfoils.

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.120.15.

Fig. 10
Transient power coefficient for NACAXX and a flat plate at (a) k = 0.08 and (b) k = 0.14
Fig. 10
Transient power coefficient for NACAXX and a flat plate at (a) k = 0.08 and (b) k = 0.14
Close modal

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 1020% 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.

Fig. 11
Cycle-averaged power coefficient versus k for NACAXX geometries and the flat plate
Fig. 11
Cycle-averaged power coefficient versus k for NACAXX geometries and the flat plate
Close modal

Conclusion

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.

Funding Data

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

Nomenclature

Aij =

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

Cs =

unadjusted suction force coefficient

Cs,sep =

adjusted suction force coefficient

CY =

lift coefficient

CP,h =

heaving power output coefficient

DVM =

discrete vortex model

f =

heaving/pitching frequency, Hz

f0sep =

steady-state separation point

fsep =

fictitious separation point

FC,Y =

circulatory lift force, N

FNC,Y =

circulatory lift force, N

Fs,Y =

suction force, N

FY =

lift force, N

h =

heaving motion, m

h0 =

heaving amplitude, m

k =

reduced frequency

LEV =

leading-edge vortex

Ph =

instantaneous heaving power, W

rc =

leading-edge radius of curvature, m

Re =

Reynolds number

t =

time, s

TEV =

trailing-edge vortex

Ushear =

leading-edge shear velocity, m/s

U =

freestream velocity, m/s

α1 =

critical angle of static separation point, deg

αeff =

effective angle of attack, deg

αeff,LE =

leading-edge effective angle of attack, deg

θ0 =

pitching amplitude, deg

θp =

angular pitching velocity, rad

ν =

kinematic viscosity, m2/s

ρ =

density, kg/m3

τ1,τ2 =

empirical relaxation time coefficients

Γ =

circulation, m2/s

Φ =

velocity potential

References

1.
Zhu
,
Q.
,
2011
, “
Optimal Frequency for Flow Energy Harvesting of a Flapping Foil
,”
J. Fluid Mech.
,
675
, pp.
495
517
.10.1017/S0022112011000334
2.
Kinsey
,
T.
, and
Dumas
,
G.
,
2012
, “
Computational Fluid Dynamics Analysis of a Hydrokinetic Turbine Based on Oscillating Hydrofoils
,”
ASME J. Fluids Eng.
,
134
(
2
), p.
021104
.10.1115/1.4005841
3.
Xiao
,
Q.
,
Liao
,
W.
,
Yang
,
S.
, and
Peng
,
Y.
,
2012
, “
How Motion Trajectory Affects Energy Extraction Performance of a Biomimic Energy Generator With an Oscillating Foil?
,”
Renewable Energy
,
37
(
1
), pp.
61
75
.10.1016/j.renene.2011.05.029
4.
Liu
,
W.
,
Xiao
,
Q.
, and
Cheng
,
F.
,
2013
, “
A Bio-Inspired Study on Tidal Energy Extraction With Flexible Flapping Wings
,”
Bioinspiration Biomimetics
,
8
(
3
), p.
036011
.10.1088/1748-3182/8/3/036011
5.
Young
,
J.
,
Lai
,
J. C.
, and
Platzer
,
M. F.
,
2014
, “
A Review of Progress and Challenges in Flapping Foil Power Generation
,”
Prog. Aerosp. Sci.
,
67
, pp.
2
28
.10.1016/j.paerosci.2013.11.001
6.
Xiao
,
Q.
, and
Zhu
,
Q.
,
2014
, “
A Review on Flow Energy Harvesters Based on Flapping Foils
,”
J. Fluids Struct.
,
46
, pp.
174
191
.10.1016/j.jfluidstructs.2014.01.002
7.
Ramesh
,
K.
,
Granlund
,
K.
,
Ol
,
M. V.
,
Gopalarathnam
,
A.
, and
Edwards
,
J. R.
,
2018
, “
Leading-Edge Flow Criticality as a Governing Factor in Leading-Edge Vortex Initiation in Unsteady Airfoil Flows
,”
Theor. Comput. Fluid Dyn.
,
32
(
2
), pp.
109
136
.10.1007/s00162-017-0442-0
8.
Wu
,
J.
,
Chen
,
Y.
, and
Zhao
,
N.
,
2015
, “
Role of Induced Vortex Interaction in a Semi-Active Flapping Foil Based Energy Harvester
,”
Phys. Fluids
,
27
(
9
), p.
093601
.10.1063/1.4930028
9.
Kim
,
D.
,
Strom
,
B.
,
Mandre
,
S.
, and
Breuer
,
K.
,
2017
, “
Energy Harvesting Performance and Flow Structure of an Oscillating Hydrofoil With Finite Span
,”
J. Fluids Struct.
,
70
, pp.
314
326
.10.1016/j.jfluidstructs.2017.02.004
10.
Su
,
Y.
, and
Breuer
,
K.
,
2019
, “
Resonant Response and Optimal Energy Harvesting of an Elastically Mounted Pitching and Heaving Hydrofoil
,”
Phys. Rev. Fluids
,
4
(
6
), p.
064701
.10.1103/PhysRevFluids.4.064701
11.
McKinney
,
W.
, and
DeLaurier
,
J.
,
1981
, “
Wingmill: An Oscillating-Wing Windmill
,”
J. Energy
,
5
(
2
), pp.
109
115
.10.2514/3.62510
12.
Kinsey
,
T.
, and
Dumas
,
G.
,
2008
, “
Parametric Study of an Oscillating Airfoil in a Power-Extraction Regime
,”
AIAA J.
,
46
(
6
), pp.
1318
1330
.10.2514/1.26253
13.
Onoue
,
K.
, and
Breuer
,
K. S.
,
2016
, “
Vortex Formation and Shedding From a Cyber-Physical Pitching Plate
,”
J. Fluid Mech.
,
793
, pp.
229
247
.10.1017/jfm.2016.134
14.
Xie
,
Y.
,
Lu
,
K.
, and
Zhang
,
D.
,
2014
, “
Investigation on Energy Extraction Performance of an Oscillating Foil With Modified Flapping Motion
,”
Renewable Energy
,
63
, pp.
550
557
.10.1016/j.renene.2013.10.029
15.
Polhamus
,
E. C.
,
1971
, “
Predictions of Vortex-Lift Characteristics by a Leading-Edge Suctionanalogy
,”
J. Aircr.
,
8
(
4
), pp.
193
199
.10.2514/3.44254
16.
Veilleux
,
J.-C.
, and
Dumas
,
G.
,
2017
, “
Numerical Optimization of a Fully-Passive Flapping-Airfoil Turbine
,”
J. Fluids Struct.
,
70
, pp.
102
130
.10.1016/j.jfluidstructs.2017.01.019
17.
Ashraf
,
M.
,
Young
,
J.
,
Lai
,
J. S.
, and
Platzer
,
M.
,
2011
, “
Numerical Analysis of an Oscillating-Wing Wind and Hydropower Generator
,”
AIAA J.
,
49
(
7
), pp.
1374
1386
.10.2514/1.J050577
18.
Wang
,
Z.
,
Du
,
L.
,
Zhao
,
J.
,
Thompson
,
M.
, and
Sun
,
X.
,
2020
, “
Flow-Induced Vibrations of a Pitching and Plunging Airfoil
,”
J. Fluid Mech.
,
885
, p. A36.10.1017/jfm.2019.996
19.
Kinsey
,
T.
,
Dumas
,
G.
,
Lalande
,
G.
,
Ruel
,
J.
,
Mehut
,
A.
,
Viarouge
,
P.
,
Lemay
,
J.
, and
Jean
,
Y.
,
2011
, “
Prototype Testing of a Hydrokinetic Turbine Based on Oscillating Hydrofoils
,”
Renewable Energy
,
36
(
6
), pp.
1710
1718
.10.1016/j.renene.2010.11.037
20.
Siala
,
F.
, and
Liburdy
,
J. A.
,
2015
, “
Energy Harvesting of a Heaving and Forward Pitching Wing With a Passively Actuated Trailing Edge
,”
J. Fluids Struct.
,
57
, pp.
1
14
.10.1016/j.jfluidstructs.2015.05.007
21.
Siala
,
F. F.
, and
Liburdy
,
J. A.
,
2020
, “
Power Estimation of Flapping Foil Energy Harvesters Using Vortex Impulse Theory
,”
Renewable Energy
,
154
, pp.
894
902
.10.1016/j.renene.2020.03.067
22.
Siala
,
F. F.
,
2019
, “
Experimental Investigations of the Unsteady Aerodynamics of Oscillating Airfoils Operating in the Energy Harvesting Regime
,” Ph.D. thesis,
Oregon State University
, Corvallis, OR.
23.
Totpal
,
A. D.
,
Siala
,
F. F.
, and
Liburdy
,
J. A.
,
2018
, “
Energy Harvesting of an Oscillating Foil at Low Reduced Frequencies With Rigid and Passively Deforming Leading Edge
,”
J. Fluids Struct.
,
82
, pp.
329
342
.10.1016/j.jfluidstructs.2018.04.022
24.
Siala
,
F. F.
,
Kamrani Fard
,
K.
, and
Liburdy
,
J. A.
,
2020
, “
Experimental Study of Inertia-Based Passive Flexibility of a Heaving and Pitching Airfoil Operating in the Energy Harvesting Regime
,”
Phys. Fluids
,
32
(
1
), p.
017101
.10.1063/1.5119700
25.
Wagner
,
H.
,
1924
, “
Über Die Entstehung Des Dynamischen Auftriebes Von Tragflügeln
,” Dynamischer Auftrieb von Tragfftigeln, Band 5, Heft, pp.
17
34
.
26.
Theodorsen
,
T.
, and
Mutchler
,
W.
,
1935
, “
General Theory of Aerodynamic Instability and the Mechanism of Flutter
,” NACA Report No. 496.
27.
Boutet
,
J.
, and
Dimitriadis
,
G.
,
2018
, “
Unsteady Lifting Line Theory Using the Wagner Function for the Aerodynamic and Aeroelastic Modeling of 3D Wings
,”
Aerospace
,
5
(
3
), p.
92
.10.3390/aerospace5030092
28.
Bird
,
H. J.
,
Otomo
,
S.
,
Ramesh
,
K. K.
, and
Viola
,
I. M.
,
2019
, “
A Geometrically Non-Linear Time-Domain Unsteady Lifting-Line Theory
,”
AIAA Scitech 2019 Forum
, San Diego, CA, Jan. 7–11, p.
1377
.
29.
Ramesh
,
K.
,
Gopalarathnam
,
A.
,
Granlund
,
K.
,
Ol
,
M. V.
, and
Edwards
,
J. R.
,
2014
, “
Discrete-Vortex Method With Novel Shedding Criterion for Unsteady Aerofoil Flows With Intermittent Leading-Edge Vortex Shedding
,”
J. Fluid Mech.
,
751
, pp.
500
538
.10.1017/jfm.2014.297
30.
Sharma
,
A.
, and
Visbal
,
M.
,
2019
, “
Numerical Investigation of the Effect of Airfoil Thickness on Onset of Dynamic Stall
,”
J. Fluid Mech.
,
870
, pp.
870
900
.10.1017/jfm.2019.235
31.
Heine
,
B.
,
Mulleners
,
K.
,
Joubert
,
G.
, and
Raffel
,
M.
,
2013
, “
Dynamic Stall Control by Passive Disturbance Generators
,”
AIAA J.
,
51
(
9
), pp.
2086
2097
.10.2514/1.J051525
32.
Mulleners
,
K.
, and
Raffel
,
M.
,
2012
, “
The Onset of Dynamic Stall Revisited
,”
Exp. Fluids
,
52
(
3
), pp.
779
793
.10.1007/s00348-011-1118-y
33.
Müller-Vahl
,
H. F.
,
Nayeri
,
C. N.
,
Paschereit
,
C. O.
, and
Greenblatt
,
D.
,
2016
, “
Dynamic Stall Control Via Adaptive Blowing
,”
Renewable Energy
,
97
, pp.
47
64
.10.1016/j.renene.2016.05.053
34.
Larsen
,
J. W.
,
Nielsen
,
S. R.
, and
Krenk
,
S.
,
2007
, “
Dynamic Stall Model for Wind Turbine Airfoils
,”
J. Fluids Struct.
,
23
(
7
), pp.
959
982
.10.1016/j.jfluidstructs.2007.02.005
35.
Dabiri
,
J. O.
,
2009
, “
Optimal Vortex Formation as a Unifying Principle in Biological Propulsion
,”
Annu. Rev. Fluid Mech.
,
41
(
1
), pp.
17
33
.10.1146/annurev.fluid.010908.165232
36.
Gharib
,
M.
,
Rambod
,
E.
, and
Shariff
,
K.
,
1998
, “
A Universal Time Scale for Vortex Ring Formation
,”
J. Fluid Mech.
,
360
, pp.
121
140
.10.1017/S0022112097008410
37.
Siala
,
F. F.
, and
Liburdy
,
J. A.
,
2019
, “
Leading-Edge Vortex Dynamics and Impulse-Based Lift Force Analysis of Oscillating Airfoils
,”
Exp. Fluids
,
60
(
10
), p.
157
.10.1007/s00348-019-2803-5
38.
Katz
,
J.
, and
Plotkin
,
A.
,
2001
,
Low-Speed Aerodynamics
, Vol.
13
,
Cambridge University Press
, New York.
39.
Ansari
,
S.
,
Żbikowski
,
R.
, and
Knowles
,
K.
,
2006
, “
Non-Linear Unsteady Aerodynamic Model for Insect-Like Flapping Wings in the Hover Part 1: Methodology and Analysis
,”
Proc. Inst. Mech. Eng., Part G
,
220
(
2
), pp.
61
83
.10.1243/09544100JAERO49
40.
Kamrani Fard
,
K.
,
Ngo
,
V.
, and
Liburdy
,
J. A.
,
2021
, “
A Leading-Edge Vortex Initiation Criteria for Large Amplitude Foil Oscillations Using a Discrete Vortex Model
,”
Phys. Fluids
,
33
(
11
), p.
115123
.10.1063/5.0065097
41.
Chorin
,
A. J.
,
1973
, “
Numerical Study of Slightly Viscous Flow
,”
J. Fluid Mech.
,
57
(
4
), pp.
785
796
.10.1017/S0022112073002016
42.
Leonard
,
A.
,
1980
, “
Vortex Methods for Flow Simulation
,”
J. Comput. Phys.
,
37
(
3
), pp.
289
335
.10.1016/0021-9991(80)90040-6
43.
Li
,
J.
, and
Wu
,
Z.-N.
,
2016
, “
A Vortex Force Study for a Flat Plate at High Angle of Attack
,”
J. Fluid Mech.
,
801
, pp.
222
249
.10.1017/jfm.2016.349
44.
Leishman
,
J. G.
, and
Beddoes
,
T.
,
1989
, “
A Semi-Empirical Model for Dynamic Stall
,”
J. Am. Helicopter Soc.
,
34
(
3
), pp.
3
17
.10.4050/JAHS.34.3.3
45.
Fan
,
Y.
,
1997
, “
Identification of an Unsteady Aerodynamic Model Up to High Angle of Attack Regime
,” Ph.D. thesis,
Virginia Polytechnic Institute and State University
, Blacksburg, VA.
46.
Liu
,
Z.
,
Lai
,
J. C.
,
Young
,
J.
, and
Tian
,
F.-B.
,
2017
, “
Discrete Vortex Method With Flow Separation Corrections for Flapping-Foil Power Generators
,”
AIAA J.
,
55
(
2
), pp.
410
418
.10.2514/1.J055267
47.
Bauchau
,
O. A.
,
2007
,
Dymore User's Manual
,
Georgia Institute of Technology
,
Atlanta, GA
.
48.
Goman
,
M.
, and
Khrabrov
,
A.
,
1994
, “
State-Space Representation of Aerodynamic Characteristics of an Aircraft at High Angles of Attack
,”
J. Aircr.
,
31
(
5
), pp.
1109
1115
.10.2514/3.46618
49.
Ngo
,
V.
,
2021
, “
Leading Edge Separation Analysis of Oscillating Airfoil Energy Harvesters
,” M.S. thesis,
Oregon State University
,
Corvallis, OR
.
50.
Stevens
,
P.
, and
Babinsky
,
H.
,
2017
, “
Experiments to Investigate Lift Production Mechanisms on Pitching Flat Plates
,”
Exp. Fluids
,
58
(
1
), pp.
1
17
.10.1007/s00348-016-2290-x
51.
Kamrani Fard
,
K.
, and
Liburdy
,
J. A.
,
2020
, “
Discrete Vortex Modeling of a Flapping Foil Energy Harvester With Lev Shedding Criterion
,”
ASME
Paper No. IMECE2020-24216. 10.1115/IMECE2020-24216
52.
Babinsky
,
H.
,
Stevens
,
R. J.
,
Jones
,
A. R.
,
Bernal
,
L. P.
, and
Ol
,
M. V.
,
2016
, “
Low Order Modelling of Lift Forces for Unsteady Pitching and Surging Wings
,”
AIAA
Paper No. 2016-0290.10.2514/6.2016-0290