A leading-edge vortex initiation criteria for large amplitude foil oscillations using a discrete vortex model

There is growing interest in the ability to provide a reasonably accurate low order model for oscillating airfoils where oscillations include both large-amplitude heaving and pitching motions,^1–3 with a review focusing on leading-edge vortex (LEV) modeling.⁴ These flows are of interest for both energy harvesting and propulsion applications. The distinction between energy harvesting and propulsion can be made based on the operating reduced frequency, $k=fc/U∞$ where f is the heaving frequency of the foil, c is the chord length, and $U∞$ is the freestream velocity. The feathering parameter, χ, defined as the ratio of maximum pitching amplitude to the angle of attack based on the maximum heaving velocity, has been used for an oscillating foil to distinguish between propulsion (χ less than one) and power generation (χ greater than one).⁵ 

Application of a heaving and pitching foil for energy harvesting was proposed by McKinney and DeLaurier;⁶ one of their main aerodynamic differences when compared with rotary turbines is the dynamic stall condition imposed during large heaving and pitching amplitudes. The occurrence of dynamic stall takes advantage of the leading-edge separation and the resulting enhanced lift forces during a cycle. These systems use a pitching and heaving foil and have energy extraction efficiencies ranging from 16.8% to 40%,^5–8 making them competitive with traditional aero-hydrokinetic turbines. The oscillating foil systems have a number of advantages in that they are versatile, making available additional resource sites, and are easy to scale and they operate with lower tip speeds, which reduces overall noise, blade stress, and wildlife mortality.^9–15 To illustrate the basic flow configuration, Fig. 1 shows a thin flat airfoil undergoing combined pitching and heaving motions defined in this study as

where θ₀ and h₀ are the heaving and pitching amplitudes, and the freestream velocity is $U∞$⁠.

A number of studies have examined several kinematic motions in an attempt to optimize flow conditions affecting forces on a foil surface. Sinusoidal motions for both heaving and pitching have received the most attention in the studies reported in the literature. However, several non-sinusoidal motions such as trapezoidal pitch profiles have shown improvements to energy harvesting, particularly at low pitch amplitudes.^14,16–18 The heaving and pitching wave forms have a strong influence on energy harvesting¹⁷ with the optimal range for heaving and pitching amplitudes found to be 0.5c to 1.0c and $70°$ to $85°$⁠, respectively.^5,13,19 These large heaving and pitching amplitudes lead to dynamic stall conditions during the cycle, which needs to be considered in modeling the transient force amplitudes on the foil.

In general, energy harvesting occurs due to both heaving and pitching contributions. A phase shift between heaving and pitching is required to optimize the flow dynamics and to take advantage of the dynamic stall, such that the largest lift forces occur during maximum heaving velocity to achieve the largest energy output rate.^15,20 The total power output (P) and the total power output coefficient (C_P) are defined in Eq. (2),

where $ḣ$ is the linear heaving velocity, $θṗ$ is the angular pitching velocity, F_L is the lift force in the heaving direction (normal to $U∞$⁠), M is the airfoil pitching moment, and b is the airfoil span.

The power output is linked to the cycle dynamics associated with the leading-edge separation in that enhanced forces occur when compared to static separation conditions. Beyond this, however, the leading-edge vortex dynamics, including its growth, advection, and possible interactions with the trailing-edge flow conditions, play a large role in overall energy harvesting performance.^21,22 Although there are still unknown aspects of the flow dynamics with both leading and trailing-edge vortex flows contributing to the force distributions and power potential of oscillating foils, the ability to predict separation is critical in generating low order models of oscillating energy harvesters.

The important parameters that have been identified influencing energy harvesting are the reduced frequency, $k=fc/U∞$⁠, the heaving amplitude per chord length, $h0/c$⁠, pitching amplitude, θ₀, the airfoil geometry, chord-based Reynolds number, $Re=ρUc/μ$⁠, and oscillating wave form. It has been shown that the optimum range of reduced frequency is 0.12–0.18, which has been confirmed by a number of studies.^8,10,20,23 The power output tends to be relatively insensitive to the heaving and pitching amplitudes when operating near the optimal reduced frequencies.^13,14,24 The highest power output occurs at a relative phase between heaving and pitching to be near $90°$⁠.^5,13,23,24

A number of studies, including experiments^7,24–28 as well as higher order computational simulations,^{10,23,29–31} have contributed to the understanding of the high amplitude oscillating flow dynamics associated with energy harvesting. The application of classical unsteady airfoil theories is typically limited to small amplitudes and mostly attached leading-edge flows.^32,33 The more recent extension to larger amplitude oscillations has been somewhat successful in predicting force transients and overall power output.^1,2,34,35 Although a theoretical approach to predicting flow separation during high amplitude oscillations has not been developed, further details of unsteady boundary layer separation during dynamic stall have been evaluated based on the operating reduced frequency³⁶ along with leading-edge flow separation prediction using empirical data for an unsteady airfoil.³⁷ In addition, more recent studies have identified vortex shedding using angle of attack as a critical parameter during pitching.³⁸ 

Some of the early efforts to address leading-edge separation have been based on local pressure and velocity conditions, such as the occurrence of strong adverse pressure gradients after a local suction peak³⁹ and the magnitude of the leading-edge acceleration.⁴⁰ For oscillating foils, the pitch rate was shown as an important parameter correlating angle of attack with separation.⁴¹ Unsteady thin airfoil theory provides a means to help predict the leading-edge separation.^1,42 This approach uses the Fourier coefficients of the vortex sheet developed along the camber line with the first coefficient, A₀, representative of the distribution of the downwash velocity relative to the freestream velocity that is associated with the airfoil dynamics. This includes the Theodorsen function to account for wake effects on the foil vortex distribution. Further modifications better account for large angles of attack.¹ Using this approach results in a leading-edge suction parameter, or LESP, based on A₀, which is determined from the low order model using thin airfoil theory. It is then hypothesized that a critical value of A₀ is associated with the occurrence of leading-edge separation. The coefficient A₀ accounts for the instantaneous angle of attack and its critical value is said to be independent of motion kinematics.⁴² Previously, however, it has been shown that varying kinematic motions play a role in leading-edge vortex formation.^16,43 Despite this, the use of a critical LESP has been shown to provide useful predictions of the transient forces on thin foils. However, the critical value of LESP must be determined by some alternative means. For example, use of detailed computational fluid dynamics (CFD) analysis to identify a separation criterion, such as an unsteady value of the skin friction coefficient near the leading-edge,³ or the use of a least root mean square (RMS) error matching of force predictions of the low order model with experimental or CFD results, to identify an optimal critical LESP that matches experimental or computational results.^1,3,28 More recently, the role of a critical leading-edge suction parameter has been shown to influence leading edge separation in unsteady flows; however, the existence of trailing edge reverse flow is suggested to be a cause for non-universal critical values.⁴⁴ Consequently, there is a need to predict leading-edge separation conditions that can be applied directly to a low order model and ideally is determined by the operating characteristics of the energy harvester.

Since the flow dynamics are cyclic, the time scales associated with the flow can be expected to be related to the cycle period or operating frequency. The timescale associated with leading-edge vortex formation has not been established, in general. However, there are some studies indicating local scaling concepts for shear generated vortex growth. For instance, it has been shown that the vortex strength can be used to define a vortex formation number based on the local shear velocity and appropriate length scale.⁴⁵ This has been developed from the more general approach to defining a universal vortex formation time.⁴⁶ Experiments on the leading-edge vortex growth in both size and strength for an oscillating foil show that, for a range of reduced frequencies from 0.06–0.16, the transient growth collapses using the timescale (⁠$cU¯SL$⁠), where $U¯SL$ is the cycle-averaged leading-edge shear velocity.⁴⁷ However, during later times in the vortex dynamics the growth rates are seen to be dependent on the reduced frequency, with lower reduced frequencies obtaining larger relative strength overall. This early time universal collapse of growth suggests that a local shear layer velocity dependent timescale may be appropriate to define the oscillating foil separation leading to vortex formation, as discussed later in this paper.

A low order model based on the discrete vortex model has been applied to the dynamic stall conditions of an oscillating foil and provides a means to determine the transient forces and overall power generation. For general details, see Refs. 37 and 48–52. More recent applications to oscillating foils include.^3,53–56 Details of the formulation of panel methods to formulate the general solution have been provided.⁵⁷ Corrections are needed in this modeling approach for leading-edge and trailing-edge separation effects. Further details of the model used in this study, the applied corrections, and implementation are provided in Sec. II. In this paper, we develop a leading-edge vortex initiation criterion applicable for low order models of an oscillating foil energy harvester derived from CFD simulations. A discrete vortex model (DVM) for a thin airfoil is used to evaluate the unsteady aerodynamics within the energy harvesting regime using this criterion. The influence of reduced frequency and Reynolds number on energy harvesting output is evaluated and compared with CFD results. The DVM in this study illustrates the use of the proposed leading-edge separation criteria. The model is based on a discrete vortex method using a panel method to describe the airfoil surface. Low order models can provide a computationally inexpensive means of analyzing the basic flow features of rather complex flows. This comes at a cost of spatial resolution, in this case the panel size used to define the geometry. Also, the model is inviscid (although a vortex blob is used to account for vortex diffusion) and lacks the ability to provide physics-based predictions of leading and trailing-edge separation. This study is focused on the implementation of leading-edge separation criteria, with a trailing-edge flow correction provided to account for trailing-edge reversed flow. Details of the model are given in Sec. II.

The paper is organized as follows. First, an overview of the low order model methodology is presented. Second, the criteria for leading-edge separation are presented and discussed. Third, the separation criteria results are scaled based on a leading-edge shear velocity to collapse data over the wide range of reduced frequencies studied. Fourth, the results of model predictions of energy harvesting performance are presented. The overall outcome is a leading-edge separation criterion that is determined based on the kinematic motion and can be used for a variety of low order models of oscillating foils.

In this study, a thin flat two dimensional airfoil (actual dimensions are 150 mm chord length, 6.5 mm width and elliptical leading and trailing-edges with a major to minor axis ratio of 6:1) is modeled using inviscid theory using a surface distribution of single vortex points on panels with vortices being shed at the leading and trailing-edges. Compared to a more analytical method such as thin airfoil theory, this method is more flexible to model complexities such as the airfoil chamber and relative motion at the leading-edge. Figure 2 shows the panel discretization and the airfoil motion kinematics. The air flow, approaching from the left, has uniform velocity of $U∞$⁠. The airfoil heaves normal to the incoming flow with linear heaving velocity $ḣ$ and pitches about the mid-chord with angular pitching velocity $θṗ$ based on the sinusoidal functions given in Eq. (1). The point vortices and collocation points are placed on each panel at 1/4 and 3/4 panel length, respectively. Boundary conditions include (i) surface impermeability which is enforced at each collocation point on every panel using N linear equations given in Eq. (3) below and (ii) the Kutta condition (zero vorticity at the trailing-edge) which is implemented by adding an additional wake panel beyond the last specified airfoil panel at the trailing-edge, and setting the pressure difference across this panel to zero,

In this equation, n refers to the normal direction on each panel i, $hn,i̇$ is the heaving motion contribution, $(θp→.×ri→)n$ is the pitching contribution, and $(∂ΦLEV∂n)$ and $(∂ΦTEV∂n)$ are the contributions from previously shed vortices at the leading-edge and the trailing-edge, respectively, and $Φ$ is the perturbation potential from the shed vortices from the leading and trailing-edges. The term $ΣAijΓj$ is the velocity contribution from the bound vortices 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 solved for. Further details of the calculation of the influence coefficients can be found in Ref. 57.

During leading and trailing-edge vortex separation, the instantaneous strength of a given leading and trailing-edge vortex is found from Kelvin's theorem, presented in Eq. (4), which states that the total circulation must be zero for a converged solution:⁵⁷ 

Vortices are discretely shed at the trailing-edge, and only at the trailing-edge when the leading-edge separation criterion is not met. The shed trailing-edge vortex (TEV) is positioned in a way to approximate the shape of the shear layer and placed at 1/3 of the distance of the previously shed vortex.⁴⁹ The strength of this shed vortex is determined iteratively using the 1D Newton's method coupled with the inviscid conservation of vorticity requirement given in Eq. (4).⁵⁷ 

Unlike trailing-edge vortices that are shed at each time step, leading-edge vortices are shed only if flow separation initiates at the leading-edge. Hence, the onset criterion of leading-edge separation needs to be accurately quantified in order to be implemented into the model. In this study, we propose the initiation of leading-edge separation based on the wall shear stress distribution near the leading-edge requiring the local shear stress to pass through zero (a positive to negative transition) during the airfoil oscillation. At this occurrence, a measure of the value of the leading-edge velocity is set as the critical parameter. If the leading-edge velocity exceeds the critical value, vortices are shed from the leading-edge. In addition, LEV shedding terminates when this measure falls below this critical value. FLUENT 2019 R2 was used to perform highly resolved CFD flow simulations to determine the time and location of separation. This includes evaluation of the wall shear stress to find the critical conditions for separation. The CFD model is based on RANS modeling with a dynamic mesh of the turbulent flow field to account for the oscillating foil, with further details provided in Sec. III B. The low order model then uses this time in the cycle to define the leading-edge measure based on the leading-edge flow velocity. Consequently, this criterion represents the local leading-edge shear velocity which occurs when the wall shear stress in the leading-edge region of the flow passes through zero. The local leading-edge shear velocity is determined by the kinematic motion within a given freestream velocity. Further details of determining the onset of leading-edge flow separation from the CFD results, as well as defining the leading-edge shear velocity as the critical parameter in the low order model, are discussed below in Sec. III along with the use of proper scaling of the separation time to collapse results over a wide range of reduced frequencies.

Applying the leading-edge vortex shedding criterion, if a vortex was shed from the leading-edge at the previous time step, the current shed vortex is positioned using the same method applied for the shed vortices from the trailing-edge, which was presented above. Otherwise, the vortex is placed based off of the induced velocity normal to the leading-edge from the freestream, heave, and pitch. To solve for the strengths of the vortices shed at the leading-edge and trailing-edge, a 2D Newton's method is used with the updated expression for Kelvin's theorem presented in Eq. (5),

where the subscript k refers to the iteration number, q is the vortex shed at the leading-edge, and m is the vortex shed at the trailing-edge. All shed vortices are convected with the low order model inviscid velocity. Therefore, the induced velocity from point vortices on each panel as well as vortices shed from the leading and trailing-edges is calculated using influence coefficients and added to the freestream velocity. Then, this inviscid velocity is multiplied by the time step to find the position increment.

Typically, discrete vortex models use point vortices with zero radii to model the flow. However, this method leads to artificially large velocities on the point vortices induced by the others in their close proximity. To overcome this issue, vortex blobs with finite core radii are used to give more realistic induced velocities as well as vorticity distributions.⁵⁸ In the present work, a vortex-core model is used for all shed vortices.⁵⁹ The core radius is taken to be 1.3 times the average spacing between the shed vortices⁴⁸ as presented in Eq. (6),

The lift force on the airfoil is calculated in the low order model using two different methods: unsteady Bernoulli's equation and the momentum impulse formulation as described below.

The pressure difference between the lower and upper surfaces of the thin airfoil can be found from the transient form of Bernoulli's equation⁵⁷ presented in Eq. (7),

where s denotes the tangential direction in the airfoil frame, $Φ$ is the perturbation potential that includes the influence of vortices on the foil as well as shed vortices, and Γ is the bound circulation per unit length. Further details of the application of Eq. (7) to oscillating foils are available.⁶⁰ For thin airfoils, the suction force acting at the leading-edge in the chord-wise direction contributes to the total lift force and can be expressed as^1,57

The total lift force becomes

The lift force can be broken up into a circulatory component F_C (the instantaneous circulation effect) and non-circulatory component F_NC (the time dependency effect) as

The force can also be calculated from the vortex impulse³⁴ using Eq. (11),

where x_i is the X-position of each vortex in inertial frame and b refers to the bound circulation. The circulatory plus suction force and the non-circulatory force in the Y direction are given by Eq. (12),

At appropriate angles of attack, an airfoil trailing-edge separation my occur at the top surface, rather than the trailing-edge, resulting in reduced lift forces. To overcome this over-prediction of lift force, a correction can be applied based on a modified Kirchhoff flow approximation.^40,61,62 This correction calculates a fictitious separation point along the chord and adjusts the circulatory part of the force to account for the bound circulation losses. This separation point for a static airfoil is given by Eq. (13),⁶³ 

where constants $S1=3.0, S2=2.3, α1=15.25°$ are determined from static experimental data.⁶⁴ The function above gives a stationary separation point; for the unsteady flow conditions, a first order differential equation is applied to describe the separation point movement and capture the transient dynamic effects,^61,65

where $τ1=0.52CU$ and $τ2=4.5CU$ are empirical relaxation time constants. This empirical correction, only applied to the circulatory and suction components of the lift force in the chord-normal direction, is then

The corrected lift force is then calculated in Eq. (16),

where $CN,non$ is the non-circulatory component of the lift force in the chord- normal direction.

The computational simulations for this study use FLUENT 2019 R2 as a high-resolution, two-dimensional solver for a reduced frequency range of $k=0.08–0.16$ and Reynolds number range of 1000–45 000. A PISO scheme was used for velocity-pressure coupling. A second-order accurate scheme was used for pressure terms, while a second order upwind and first order upwind schemes were used for discretizing momentum and modified turbulent viscosity, respectively. A dynamic mesh⁵ was applied to the computational model using two mesh regions connected through a non-conforming sliding interface. The first mesh was a fixed reference frame, while the smaller second mesh region enclosed and moved vertically heaving and rotationally pitching at the same rate of the airfoil. Shown in Fig. 3 are compared with Kinsey and Dumas⁵ for a full cycle with t/T being the cycle time divided by the cycle period, T. Slight differences occur slightly after the top and bottom heave positions which may be a result of the application of different turbulence models as well as dissimilar parameters for chord length, freestream velocity, and heaving frequency. Despite these slight differences, the key flow physics between the two results are very close throughout the cycle.

The leading-edge separation criterion uses the condition of wall shear stress sign reversal along the foil surface (top surface during the upstroke and bottom surface during downstroke). In addition, as can be shown in the momentum equation evaluated very close to the surface, a sign change of the stress gradient away from the surface is approximately balanced by the pressure gradient along the surface.⁶⁶ Consequently, a change of sign of the stress at the wall would be consistent with a change of sign of the pressure gradient along the surface. The low-pressure spike (large suction) formed near the leading-edge results in a pressure gradient change of sign as well. Therefore, it is proposed that using the wall stress sign reversal will be a strong indicator of separation consistent with a large suction event occurring during the cycle. As such, the wall shear stress along the foil surface is monitored and the first instance of wall shear stress of zero coupled with a pressure gradient minimum near the leading-edge is used as the separation criteria. Figure 4 illustrates the wall shear stress and pressure gradient distributions along the surface near the leading-edge as well as contour plots of vorticity and pressure for the case of k = 0.08 for the sinusoidal heaving and pitching motion. The wall shear stress is seen to pass through a minimum of zero at the critical time with the vorticity plots showing a small concentration of clockwise vorticity on the bottom surface as the foil cycles downward. Similarly, the pressure gradient forms a peak at a slightly earlier time and slightly closer to the leading-edge, with the pressure contours showing the beginning of a large suction pressure close to the leading-edge. It is not expected that the critical times of stress and pressure gradient coincide, but rather the pressure gradient forms initially that then causes the reduction of stress leading to separation.⁶⁷ It is chosen to use the wall stress criteria time, rather than the pressure gradient, due to the time lag required to reduce and reverse the flow momentum needed for separation.

The same approach has been applied for the entire range of reduced frequencies, $k=0.06–0.16$⁠, and the critical times corresponding to leading-edge separation are shown in Fig. 5. As can be seen, the critical time for separation is delayed in the cycle with increasing reduced frequency. Since increased reduced frequency can occur with either decreased freestream velocity (and lower suction pressure) and/or increased oscillating frequency (increased velocity component downward toward the surface), this delay of separation is expected and shown to deviate from a linear trend at the higher reduced frequencies.

Previous studies have shown that the onset of leading-edge flow separation can be characterized by a critical flow parameter at the leading-edge such as the leading-edge pressure gradient³⁹ and flow velocity.⁴⁰ More recently, studies show that the leading-edge suction is directly related to the first Fourier coefficient, A₀, from thin airfoil theory.⁴² A₀ can be found in terms of downwash velocity, W(x, t), as

where the downwash velocity is given by

A₀ accounts for the instantaneous angle of attack resulting from motion kinematics. It has been hypothesized that a critical value of A₀ would correspond to the onset of LEV formation for a given airfoil and Reynolds number and found A₀ to be a good measure of the flow behavior at the leading-edge since it is the only term in thin airfoil theory that results in a singularity in the vorticity distribution at the leading-edge.⁴² Furthermore, a few researchers have found that for oscillating airfoils, the pitching rate and angle of attack correspond to initiation of dynamic stall and laminar flow separation at the leading-edge.^41,68 Inspired by these findings, we modify the effective angle of attack to include the pitching component at the leading-edge. This is not strictly an effective angle of attack since the pitching contribution is only valid at the leading-edge and not further along the chord. However, it does provide a measure that relates to the leading-edge shear velocity by accounting for both heaving and pitching at the leading-edge.

The modified effective angle of attack at the leading-edge of the airfoil is defined here as

where the subscript LE refers to evaluating the pitching contribution at the leading-edge. It is anticipated that this parameter will provide a means to account for the leading-edge separation through a critical value for $αeff,LE$ during the cycle when separation occurs as determined using the wall shear stress indicator presented previously. This parameter is only dependent of motion kinematics and freestream velocity, a benefit to modeling since leading-edge separation may be predicted based on a critical value within the cycle without having to run the flow model. The results for $αeff,LE$ as well as A₀ and their associated critical values predicted based on the wall shear stress criteria are presented in Sec. III.

The separation criteria are shown for a range of reduced frequencies $k=0.06–0.16$⁠, shown to produce higher energy harvesting effciencies.¹⁰ A combined sinusoidal heaving/pitching motion is used, with a heaving amplitude of $h0/c=0.5$⁠, a pitching amplitude of $θ0=70°, 90°$ phase shift, and pitching axes located at the mid-chord. The Reynolds number ranges from approximately 15 400 to 41 000, respectively for the highest to lowest reduced frequencies.

The separation times, scaled using the leading-edge shear layer velocity (U_SL), are used to find corresponding values of the pitching modified effective angle of attack parameter, $αeff,LE$⁠, given in Eq. (19). The shear layer velocity is approximated as the vector sum of the leading-edge local velocity and the component of the freestream velocity (⁠$U∞$⁠) in the direction of the airfoil motion¹⁹ as follows:

The leading-edge shear velocity is a combination of the freestream, heaving, and pitching contributions to the flow at the leading-edge. In addition, results of the low order model calculations of the lift force as well as the transient vorticity distributions are presented and compared with the CFD results.

The timescale associated with the cycle frequency is the cycle period, T. In addition, an associated timescale for the vortex formation can be formed as $τv∼LvU¯SL$ where L_v is a length scale related to the leading-edge and $U¯SL$ is the cycle averaged shear layer velocity generated at the leading-edge presented in Eq. (20). The shear layer velocity can be thought of as being proportional to the relative strength of the local vorticity, if divided by an appropriate length scale, such as possibly the thickness of the shear layer. Assuming that the shear layer thickness is proportional to the leading-edge radius of curvature, r_LE, then one could define L_v based on this radius. However, for a thin airfoil this becomes problematic, so assuming that the radius is proportional to the chord length, c, we use L_v = c as a first approximation and define the timescale as

This timescale is expected to be rather large, since the length scale used is c, resulting in small values of $t/τv$ at the critical condition. When applied to thick foils, this timescale could be modified by multiplying by the ratio: $rLE/c$⁠. For a circular leading-edge, the radius of curvature of the airfoil tip would be the circle radius, and for an elliptically shaped leading-edge, the radius of curvature would be $b2/a$ (where b is the semi-minor axis and a is the semi-major axis). However, since separation may occur further along the foil surface than at the tip, the appropriate radius of curvature is not fully obvious. Results are shown below using an assumed circular radius for the leading-edge with the radius then set to one half of the thickness of the flat foil.

The results of identifying the separation cycle time, t_critical, scaled with the cycle period, T, using the analysis of the wall shear stress zero crossing sign change, as illustrated in Fig. 4, for all of the reduced frequencies are shown in Fig. 6. This plot shows the values of A₀ obtained from the low order model as a function of time with the critical times indicated in the plot. These results show that as the reduced frequency increases, the separation time increases as does the (negative) value of the critical A₀. Therefore, for increasing k resulting from either reduced $U∞$ or increased oscillating frequency, separation occurs later within the cycle period. Also shown for each reduced frequency are the results of the scaled critical separation time criteria based on the definition of τ_v in Eq. (21). As can be seen, these critical times collapse to a single value over the entire range of reduced frequencies. This indicates that the non-dimensional time for separation is independent of reduced frequency. These results are contrary to the assertion that a single value of critical A₀ is valid for all reduced frequencies. It implies that the separation, as denoted by the wall stress criteria given, has a relatively constant value of approximately $tcritical/τv=0.4$ for the large amplitude oscillations used in this study. This value is, as mentioned, a consequence of using the large length scale, c in the definition of τ_v. Estimating the radius of curvature for the current airfoil to be equal to the distance along the chord from the tip to the point of maximum thickness (which assumes a circular geometry radius of curvature) then the non-dimensional critical time becomes $tcritical/τv=3.1$⁠. This value approaches the “optimal vortex formation number” proposed for shear layer driven vortex formation.^45,46 In both cases, using the definition of τ_v as the ratio of a foil length scale and $U¯SL$ results in a constant non-dimensional critical cycle time for leading-edge separation for the range of reduced frequencies evaluated. Importantly, given this constant, the actual cycle time can be found based on the given kinematic motion and airfoil geometry for a given freestream velocity. Further details need to be worked out depending on airfoil geometry and separation location.

Consistent with these results, it has also been shown that the leading-edge vortex strength and size growth vs time, when the time is scaled using Eq. (21), collapse into a single curve for the same range of reduced frequencies as used in this study.⁴⁷ However, this collapse onto a single curve for all values of k is only valid at early growth times, up to $t/τv∼1.4$⁠. For later times, low and high values of k (k < 0.1 or k > 0.1) result in different scaled growth rates, which indicates that other flow dynamics, not necessarily attributed to the leading-edge shear velocity, become important such as interactions of vortices with trailing-edge flow effects.

The relative leading-edge shear velocity, $USL/U∞$⁠, as well as the proposed separation criteria parameter, $αeff,LE$⁠, is shown in Fig. 7 vs the cycle time t/T as well as $t/τv$ for the full range of reduced frequencies. As shown, as k increases the amplitude of the relative shear velocity decreases, with values in this study near the freestream value. There is a slight shift to later cycle times for the peak shear velocity (from t/T = 0.30 to 0.36 as k increases from 0.06 to 0.16). This is reflective of the increase in both heaving and pitching velocities relative to the freestream velocity. Also, shown in Fig. 7 is the time variation of $αeff,LE$ vs both t/T and $t/τv$⁠. The trends are similar to the relative shear velocity based on the given definition of $αeff,LE$⁠. The critical values of $αeff,LE$ appear to collapse to a small value near zero. This indicates, according to the definition of U_SL, Eq. (20), that the critical condition occurs when the heaving and pitching velocity contribution of the foil nearly balance the chord-normal freestream velocity component. This condition can be defined based on the kinematic motion of the foil.

The effect of changing Reynolds number while maintaining the same reduced frequency was investigated in terms of its effect on $αeff,LE$⁠. This was done by maintaining the same oscillating frequency while changing the chord length and freestream velocity with the constraint of a specified Reynolds number. For instance, if k = 0.08, with a frequency of 1.6, and Re = 10 000 in air, then the chord length is 0.086 m and $U∞$ is 1.719 (m/s). The CFD simulations were run for k = 0.08 and 0.14 for Re values of 1000, 10 000, 30 455, and 45 000. For each case, the critical time based on the wall shear stress criteria was determined, and at this time, the value of the relative shear velocity $USL/U∞$⁠, as well as $αeff,LE$ was calculated. The results are shown in Fig. 8, indicating the change of both the relative shear velocity and $αeff,LE$ at the critical time. These results indicate that these two critical values trend toward an asymptotic value toward zero with increasing Re, or a critical value largely independent of Re at sufficiently large Re. However, for low Re operation there is a much stronger dependence on Re and this would need to be explored when operating in this regime. Also, the higher reduced frequency case, k = 0.14, shows a trend toward an asymptote for lower Re values compared with the lower reduced frequency case. All of the results presented in this paper have a range of Re between 1000 and 45 000 which seems to be within the asymptotic range. It is speculated that in the limit of very large Re the results would approach zero indicating the critical condition being $αeff,LE=0$⁠. Further study is required of these trends which would also include variations in oscillating wave form and airfoil geometry.

The non-dimensional span-wise vorticity evolution found from the low order model is presented in Fig. 9 and is compared with CFD results for (a) k = 0.08 and (b) k = 0.14. The flow approaches from the left and the non-dimensional times $t/T=0$ and $t/T=0.5$ correspond to the top and bottom heaving positions, respectively. Results are the same for the upstroke and downstroke so only the downstroke half cycle is shown. Clockwise rotation is shown in blue and counterclockwise in red. Low order model results of the vorticity fields for both reduced frequencies show good agreement with CFD results, with the model capturing the growth and advection of both leading-edge vortices (LEV) and trailing edge vortices (TEV). Since the convective timescale of the flow is relatively larger for the higher reduced frequencies, the formation and advection of both LEV and TEV occur at a slower rate for k = 0.14 compared with k = 0.08. Moreover, the size of LEV is significantly smaller for the larger reduced frequency. This is the result of a decrease in the shear layer velocity (U_SL) for higher reduced frequencies due to either decreased freestream velocity or increased oscillating frequency as can be seen in Eq. (20). An additional major difference between the low and high reduced frequency cases is that the trailing-edge shear layer is not strong enough to generate roll-up into a TEV. This is shown by examining the two reduced frequency cases at later times, where for k = 0.08 there is a strong clockwise rotation at the trailing edge, whereas for k = 0.14 this does not occur. This is a result of the LEV being advected later in the cycle so that it is present near the trailing edge, coinciding with a smaller angle of attack, which suppresses TEV growth. A consequence of this is that the LEV growth rate and its advection influence the transient force distribution on the foil during the cycle which seems to be well modeled using the discrete vortex method over this range of reduced frequencies.

The results of the low order model predictions of energy harvesting performance in terms of the transient forces (⁠$CL=FL/12ρU∞2cb$⁠) generated during a cycle are presented in Fig. 10 and compared with CFD results for the same conditions. The leading-edge separation is based on the use of the critical value of $αeff,LE$ at the associated critical time of $t/τv=0.4$⁠. The model force results are shown for both the application of the transient Bernoulli equation and the impulse method. The results compare well with the CFD results which include viscous effects. The low order model captures the same trends, albeit with some local quantitative differences that depend on the reduced frequency. The two force methods match almost identically, with only small differences near the primary and secondary peaks for higher reduced frequencies. The general description of the occurrence of the primary and secondary peaks has been shown to be dependent on the reduced frequency,⁴⁷ as is also shown in these results. The trailing-edge separation correction was also applied to both force calculation methods (Bernoulli-Corrected and Impulse-Corrected). The correction slightly decreases the difference between the model and CFD results. The greatest improvement with the correction is seen for k = 0.10 and k = 0.12. For reduced frequencies greater than 0.12, the correction leads to an under-estimation of the primary and secondary peaks. Since this correction is empirically based using the effective angle of attack, further work is required to determine the nature of the corrections needed for oscillating energy harvesting foils, especially as a function of reduced frequency. Overall, however, the low order model with the application of the leading-edge separation criteria captures the general cycle variations and performance characteristics.

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 shown in Fig. 10 show the development of a secondary peak late in the half cycle (near $t/T=0.45$⁠) which is more pronounced than that given for with CFD model results. Comparing low and high reduced frequency cases in Fig. 9 shows that the leading-edge vortex occurs later in the high reduced frequency case and prevents a large trailing edge vortex from forming near the end of the half cycle. However, there is a significant backflow along the foil shown as blue vorticity indicative of reverse up toward the leading-edge. Although this vorticity distribution is identified in the low order model its viscous interactions with the foil are most likely poorly evaluated. This difference in vortex dynamics for the low and high reduced frequency conditions are thought to play a significant role in the observed changes in the transient loads on the foil. If the viscous effects become important on the vortex dynamics near the surface it is expected that the DVM model would not accurately predict the resultant force and yield larger loads, as is shown by the secondary peaks in Fig. 10. This speculation requires further investigation and possible low order model improvements at the very high reduced frequencies.

The results of the calculated total transient power output using the low order model, expressed as the power coefficient, are shown in Fig. 11(a). The transient power is based on both the heaving and pitching contributions identified in Eq. (2). Experimental details of the vortex dynamic effects on the transient power production for a range of reduced frequencies show how vortex strength and convection influence the power production with a peak near k = 0.12.²⁶ The peak power determined by the low order model follows this trend as is more evident based on the cycle averaged power coefficient given in Fig. 11(b) vs the reduced frequency.

This study examines the separation criteria associated with a large amplitude heaving and pitching airfoil operating in the energy harvesting regime. The separation condition is required to further develop low order models to predict aerodynamic loads and power output of oscillating energy harvester. Using CFD simulations, a separation criterion based on the cycle time of the zero crossing of the wall shear stress near the leading-edge is identified for a wide range of reduced frequencies for a thin airfoil. The corresponding leading-edge pressure gradient minimum is also evaluated to be consistent with the generation of leading-edge suction. This criterion results in a critical time for separation that is shown to collapse to a single non-dimensional value using a timescale based on the leading-edge shear velocity and chord length. This timescale and critical time show consistency with the vortex formation time and provide a means of predicting separation for energy harvesters that can be used in lower order models. It is also shown that the separation criterion can be expressed in terms of a modified effective leading-edge angle of attack that accounts for pitching and heaving contributions. This provides a means to predict separation solely based on the foil kinematic motion in a given freestream flow. This has the added benefit of providing separation criteria for a wide range of low order models. This modified effective angle of attack is shown to approach zero at the critical separation condition. A Reynolds number dependence is shown at low operating Reynolds numbers but an asymptotic value of zero for the effective leading-edge angle of attack seems appropriate at high Reynolds numbers. A low order model based on the discrete vortex method is used to assess the ability to use the separation criterion to predict transient lift forces for a range of reduced frequencies. A trailing-edge separation correction is also applied to the model results, which is shown to decrease the slight discrepancy between the model and CFD for reduced frequencies smaller than k = 0.12, especially around the primary and secondary peaks. Although the focus of this study is on oscillating foil energy harvesting, it may also contribute to the general understanding of oscillating flow separation criteria.

The authors would like to acknowledge the financial support provided by the National Science Foundation CBET Award No. 1804964.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.