Stiffness effects on laminar separation flutter

Journal of Fluids and Structures 91 (2019) 102767

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Stiffness effects on laminar separation flutter ∗

Caleb J. Barnes , Miguel R. Visbal Air Force Research Laboratory, Wright-Patterson AFB, OH 45433, United States

article

info

Article history: Received 4 June 2019 Received in revised form 18 September 2019 Accepted 9 October 2019 Available online xxxx Keywords: Fluid–structure interaction Laminar separation flutter Transitional flow Limit-cycle oscillation Flutter Low-to-moderate Reynolds number

a b s t r a c t This work explores self-sustained pitching oscillations of a NACA0012 airfoil operating at low-to-moderate Reynolds numbers in which the aerodynamic flow is in a transitional regime. One-degree of freedom (1-DOF) pitching oscillations are explored for chordbased Reynolds numbers of Rec = 77, 000 and 110,000, which fall on the lower end of the flutter regime, over a large range of torsional spring constants. Laminar separation flutter (LSF) is shown to persist over several orders of magnitude of structural rigidity. At the Reynolds numbers tested, the effect of changing structural rigidity at a fixed Reynolds number exhibits remarkable similarity to the effect of varying Reynolds number at a fixed stiffness. Although the two effects are not perfectly analogous, both parameters influence the timing of transition events relative to the limit-cycle. Raising stiffness values at either Reynolds number causes the structural response to outpace boundary layer transition which in turn yields an increasingly nonlinear aeroelastic response. Oscillations are no longer sustained when the structural frequency merges with the frequency of the aeroelastic response. At this point, the slow relaxation time of boundary layer transition fails to keep pace with the structural frequency and the processes that sustain LSF are inhibited. Oscillations terminate at a lower stiffness in the lower Reynolds number case—an apparent effect of slower boundary layer transition. Interestingly, hysteresis patterns and flow topologies coalesce into nearly identical topologies for both Reynolds numbers when stiffness is very low. Published by Elsevier Ltd.

1. Introduction Extended regions of laminar flow are feasible in a number of realistic scenarios such as high-altitude long-endurance (HALE) platforms, small unmanned aerial vehicles, or in the use of laminar flow control devices to reduce skin friction drag. All of these scenarios exploit energy efficiency initiatives to meet endurance and sustainability goals. Creative design of light-weight multi-functional airframes can further contribute to satisfying the aforementioned design criteria. However, light-weight laminar wings are not without risk as aeroelastic interactions in the transitional flow regime are not well explored to date. Multidisciplinary interactions may arise due to disruptions in transition location, unsteady separation, or complicating factors posed by flight-path disturbances or atmospheric turbulence. Recent wind tunnel experiments by Poirel et al. (2008) and Poirel and Mendes (2014) report evidence of pitchdominated self-sustained oscillations in the transitional flow regime. This phenomenon was attributed to nonlinearity in the aerodynamic loads provided by laminar boundary layer separation (Poirel et al., 2008; Poirel and Yuan, 2010; Yuan et al., 2013) leading to the so-called laminar separation flutter (LSF). Several computational studies have confirmed nonlinear aerodynamic loading and negative aerodynamic damping in the transitional flow regime (Poirel and Yuan, 2010; ∗ Corresponding author. E-mail address: [email protected] (C.J. Barnes). https://doi.org/10.1016/j.jfluidstructs.2019.102767 0889-9746/Published by Elsevier Ltd.

2

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

Yuan et al., 2013). The unsteady processes that appear to drive the sustained oscillations are sensitive to a number of external factors. For instance, tripping the laminar boundary layer or introducing free stream turbulence can suppress pitching oscillations (Poirel et al., 2008, 2011; Yuan et al., 2015). This phenomenon is distinct from other well-known limit-cycle-oscillations (LCO). Occurring at low speed and low angle of attack, it lacks the large shock motions of transonic flutter or massive separation inherent to stall flutter. Gibbs et al. (2014) have reported low Reynolds number LCO on a stiff flag configuration in which no separation was observed. However, LSF appears to require both laminar separation and the subsequent flow transition processes (Poirel et al., 2011; Barnes and Visbal, 2018b). Barnes and Visbal (2018c) explored the roles of flow transition in LSF using high-fidelity implicit large-eddy simulation (ILES) coupled with structural dynamics over a range of Reynolds numbers. Several flow transition processes were found to be characteristic at all Reynolds numbers studied. Laminar separation forms on the rear portion of the wing. Eventually, the detached shear layer becomes unstable, transitions to turbulence, and reattaches at the trailing edge. This process imparts a brief spike in moment coefficient which tends to damp the LCO. Transitional flow is briefly alleviated following this event. The long laminar separation bubble (LSB) initially formed by turbulent reattachment is then maintained by laminar reattachment due to momentum transfer of spanwise coherent boundary layer vortices. Later, a second excursion of transitional flow emerges and eventually meets the LSB which is further modified and shortened by turbulent reattachment. This flow pattern alternates between the upper and lower surfaces in subsequent strokes. Reynolds number tends to influence the timing of the aforementioned transitional flow events (Barnes and Visbal, 2018c). Turbulent reattachment of open separation occurs earlier as Reynolds number is increased. This effect causes positive aerodynamic damping to initiate at lower angles of incidence. Enhanced second-stage transition shortens the separation bubble and delays its approach to the trailing-edge. Barnes and Visbal (2018c) found that higher Reynolds number cases, Rec = 150,000 and 200,000, required a disturbance to initiate LCO which could be realistically achieved by an external flow disruption (Barnes and Visbal, 2018a). While some progress has been made in understanding self-sustained oscillations in the transitional flow regime, knowledge of this new phenomenon is largely constrained to a limited set of parameters. Poirel et al. (2008) considered the effects of stiffness and elastic axis location on amplitude and frequency over a range of Reynolds numbers but the specific changes to unsteady flow structure were not explored in detail. It remains to be demonstrated whether LSF is specific to the limited configurations tested thus far or prevails over a broad range of parameters—as this article demonstrates. The precursory works by the authors (Barnes and Visbal, 2018c, 2016a) report the effects of Reynolds number and preliminary efforts toward understanding effects of airfoil compliance on LSF. The present study seeks to address this gap in the literature by considering a wide range of torsional rigidity, 0.003 N m/rad ≤ Kθ ≤ 30.0 N m/rad at Reynolds numbers on the lower end of the flutter regime. Values of Rec = 77,000 and Rec = 110,000 are considered. The lower number produced the most unique and nonlinear response at the baseline stiffness of Kθ = 0.3 N m/rad from Barnes and Visbal (2018c) while the higher value yielded the largest amplitude. High-fidelity ILES coupled with linear structural equations of motion is used to explore the unsteady evolution of fine-scale flow structure for the present aeroelastic cases. 2. Computational framework 2.1. Aerodynamics The high-order ILES solver FDL3DI (Visbal and Gaitonde, 1999; Gaitonde and Visbal, 1998) is used for all simulations in the present study. This computational framework solves the full, unfiltered, compressible Navier–Stokes equations cast in strong conservation form on a general time-dependent curvilinear coordinate system. The system of equations are integrated in time using the implicit, approximate factorization of Beam and Warming (1978) and simplified through the diagonalization of Pulliam and Chaussee (1981). The time-integration scheme is augmented through a Newton-like sub-iteration procedure to maintain temporal accuracy (Rai and Chakravarthy, 1986; Rai, 1989). Fourth-order, nonlinear dissipation terms (Jameson et al., 1981; Pulliam, 1986) are appended to the implicit operator to improve stability. The explicit operator of the implicit time-integration scheme represents the numerical approximation and dictates the formal order of accuracy. Spatial derivatives in the explicit operator are discretized along a coordinate line in the computational domain using the implicit, 6th-order, formulation of compact-differencing (Lele, 1992). High-order onesided formulas, designed to retain the tri-diagonal form of the system of equations, are applied at the computational boundaries (Visbal and Gaitonde, 1999; Gaitonde and Visbal, 1998). The solution approach for the Navier–Stokes equations described above is used to solve laminar, transitional, and turbulent flow regions without change using an ILES procedure. The ILES approach does not require sub-grid-scale (SGS) models or additional heat flux terms required by standard large-eddy-simulations (LES). Alternatively, a high-order, lowpass Padé-type filter, based on the templates proposed by Lele (1992) and Alpert (1981), is applied to eliminate spurious components. The filter is applied to the conserved variables along each transformed coordinate direction once after each time step or sub-iteration. An 8th-order filter is used for the interior points in the present work which selectively damps only the poorly resolved high-wavenumber content. The one-sided filtering strategies described by Visbal and Gaitonde (1999) and Gaitonde and Visbal (1999) are applied at near-boundary points. The high-order low-pass filtering used in conjunction with high-order spatial discretization provides an effective alternative to the use of SGS models as shown by Visbal and Rizzetta (2002), Visbal et al. (2003), and more recently by Garmann et al. (2013). A reinterpretation of this ILES approach in the context of an Approximate Deconvolution Model (Stolz and Adams, 1999) has been presented by Mathew et al. (2003). As the grid resolution increases or Reynolds number decreases, the ILES approach becomes effectively direct numerical simulation (DNS).

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

3

2.2. Structural dynamics and coupling The structural component solves the generic linear system of equations of motion,

¨ + Cd˙ + Kd = r Md

(1)

where d is the vector of degrees of freedom, M is the mass matrix, C is the damping matrix, K is the stiffness matrix, and r is the external load vector. The specific values contained in M, C and K in Eq. (1) are problem dependent and introduced for each case as needed. Eq. (1) is expanded to a 2N system in state-space form as a series of 1st-order ordinary differential equations by introducing a new set of variables a = {d1 , d˙ 1 , . . . , dN , d˙ N }T where N is the degrees of freedom of the original system. The resulting system of equations is integrated in time using the 2nd-order implicit scheme of Beam and Warming (1978) augmented with Newton-like sub-iterations. Greater detail on the structural solution procedure can be found in Morton et al. (1998) and Barnes and Visbal (2016b). A similar approach including cubic non-linearity terms has also been previously implemented and verified against experimental data by Golubev et al. (2014) and Nguyen et al. (2015). Loose temporal coupling between the aerodynamic and structural models can adversely affect solution integrity (Gordnier and Visbal, 2002). The time lag in the present system of equations is eliminated by implicitly coupling the two physics models through the global sub-iteration procedure. Because both physics models are cast in iterative form, a fully implicit coupling between the aerodynamic and structural models can be obtained. Within each sub-iteration, aerodynamic forces are integrated on the airfoil surface and then passed to the structural model. The resulting displacements are then returned to the aerodynamics solver and used to move the fluid dynamics mesh. This interchange is repeated within each time step thereby synchronizing the two procedures and preserving second-order temporal accuracy for both the fluid and structural dynamics models. 3. Details of the computations 3.1. Configuration A NACA0012 airfoil section elastically mounted in pitch is considered in this work as a representative model for torsional deformations. The equation of motion in terms of pitch-angle, θ , is given in dimensional form as Iea

d2 θ dt 2

+ Dθ

dθ dt

+ Kθ (θ − α) = Mea

(2)

where θ is the pitching displacement from the reference angle of attack, α . A positive value for θ corresponds to a pitch-up displacement. Iea is the moment of inertia about the elastic axis, Kθ is the torsional stiffness, Dθ is the damping, and Mea is the pitching moment about the elastic axis where a pitch-up moment is positive. Eq. (2) is cast in non-dimensional form using the total pitching moment coefficient given as 2 CM = 2Mea /ρ∞ U∞ sc 2

where s is the spanwise width of the airfoil section and c is the airfoil chord length, as well as the quantities

ωθ =

√

Kθ /Iea ,

u=

2U∞ c ωθ

,

Dθ

ζθ = √

2 Iθ Kθ

,

µθ =

2Iea

ρ∞ sc 4

,

defined as the undamped natural frequency, reduced velocity, damping ratio, and inertia parameter, respectively. These substitutions yield the final form of the structural dynamics given as d2 θ dτ 2

+ 2ζθ

2

dθ

u

dτ

( )

( )2 +

2 u

(θ − α) =

CM

µθ

(3)

where the non-dimensional time is defined as τ = tc /U∞ . The pitching moment coefficient is taken about the elastic axis in this article which is located at xea /c = 0.186 downstream from the leading edge. The resulting Eq. (3) is used to construct the matrices in Eq. (1). The structural properties, elastic axis, and reference angle of attack listed below in Table 1. These values are the same as used by Barnes and Visbal (2018c) with the exception of torsional rigidity, Kθ which is varied in this work. 3.2. Computational mesh and boundary conditions An O-mesh topology is chosen to discretize the computational domain around the airfoil section as shown in Barnes and Visbal (2016b). Grid sensitivity was evaluated for a static NACA0012 airfoil at a Reynolds number of Rec = 200,000 and angle of attack of α = 4◦ in Barnes and Visbal (2016b). The static case provides a conservative estimate for this problem as transition to turbulence in the aeroelastic cases fails to develop to the scales experienced in the static case. The mid-level mesh proved a reasonable choice for this more restrictive scenario and therefore can be extended to the current Reynolds numbers of Rec = 110,000 and 77,000 with confidence.

4

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767 Table 1 Physical properties for elastically mounted airfoil. Quantity

Value

Units

U∞ c s Iea Kθ Dθ xea

7.5, 10.7 0.156 1.0 0.00135 0.003, 0.30, 1.0, 3.0, 30.0 0.002 0.186c 0◦

m/s m m kg m2 N m/rad N m s

α

Table 2 Computation parameters, (Rec = 110,000). Case Kθ Kθ Kθ Kθ Kθ

= 0.003 = 0.300 = 1.000 = 3.000 = 30.00

(N (N (N (N (N

m/rad) m/rad) m/rad) m/rad) m/rad)

Run time

Time step

Span extent

Phase-avg. cycles

Disturbance

τ τ τ τ τ

∆τ ∆τ ∆τ ∆τ ∆τ

s/c s/c s/c s/c s/c

10 10 15 15 N/A

No No No No Yes

= 344 = 300 = 250 = 250 = 150

= 1 × 10−4 = 1 × 10−4 = 1 × 10−4 = 1 × 10−4 = 1 × 10−4

= 0.2 = 0.2 = 0.2 = 0.2 = 0.2

The no-slip condition is enforced at the airfoil surface where grid velocities are directly translated from the structural equations of motion. Third-order adiabatic (∂ T /∂ n = 0), and zero normal pressure gradient (∂ P /∂ n = 0), conditions are also employed at the surface. Spatial periodicity is imposed at both the O-grid cut line and spanwise boundaries through a 5-point coincident overlap. The spanwise extent of the computational domain is 0.2c and is sufficiently wide to allow the spatial evolution of transitional flow features. Freestream conditions are specified at the far field boundary located 100 chord lengths away from the wing surface. Mesh stretching in conjunction with high-order low-pass spatial filtering provides a buffer zone to eliminate startup transients and flow acoustics as they propagate away from the airfoil. This buffer effect is reinforced by the application of an absorbing sponge zone (ASZ) near the far field boundary and described in more detail in Barnes and Visbal (2016b). 4. Results This section explores the effect of stiffness on self-sustained pitching oscillations at Reynolds numbers of Rec = 110,000 and 77,000. Both flow speeds were previously considered in more detail using a spring constant of Kθ = 0.3 (N m/rad) in Barnes and Visbal (2018c) which will be considered as the baseline for comparison in the present discussion. The higher Reynolds number was examined because it provided the largest oscillation amplitude in the prior study while the lower value exhibited uniquely nonlinear behavior compared to all other cases. The torsional stiffnesses considered are recorded above in Table 1. The softest configuration is two orders of magnitude below the baseline value largely removing the spring’s restoring forces while the stiffest case is two orders of magnitude larger than the baseline rigidity. The damping, Dθ , and moment of inertia, Iθ , are unchanged from those employed in Barnes and Visbal (2018c). Each of the simulations were initiated from the previously completed static solution at α = 0◦ with the wing section undeflected and at rest; θ = 0◦ , θ˙ = 0. Cases which do not oscillate on their own were subjected to a disturbance in the form of a smooth square-pulse torque described by 1 CM ,max (tanh(2a(τ − t1 )) − tanh(2a(τ − t2 ))) (4) 2 where a = 10, t1 and t2 are the start and end times of the applied disturbance and CM ,max is the amplitude of the applied load. Values for t1 , t2 , and CM ,max are case dependent. CM ,applied (τ ) =

4.1. Rec = 110,000 Effects of stiffness on self-sustained oscillation at a Reynolds number of Rec = 110,000 are presented first. Details regarding the run time, number of cycles for phase-averaging, and computational parameters used are recorded in Table 2. Fig. 1 shows the dynamic response over the first τ = 200 characteristic times at all five spring constants. Harmonic, self-sustained oscillations emerged at the four lowest stiffness values of Kθ = 0.003, 0.3, 1.0, and 3.0 (N m/rad). The Kθ = 30.0 (N m/rad) case was subjected to a disturbance in the form of Eq. (4) where t1 = 1.0, t2 = 2.0 and CM ,max = 1.0 which produced a θo > 5◦ initial deflection. Following the disturbance, oscillations decayed to an effectively static state, Fig. 1(e), demonstrating flutter is not sustained in the stiffest configuration. Figs. 2(a) and (b) compare the pitching amplitude and dimensional first-harmonic frequency, respectively, for each spring constant. The results of the Rec = 110,000 simulations are shown alongside those of the Rec = 77,000 cases

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

5

Fig. 1. Time history of pitching oscillations for Rec = 110,000; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), (c) Kθ = 1.0 (N m/rad), (d) Kθ = 3.0 (N m/rad) and (e) Kθ = 30.0 (N m/rad).

Fig. 2. Effect of stiffness on time-mean (a) pitching amplitude, (b) dimensional first-harmonic frequency, and (c) non-dimensional first-harmonic frequency. The horizontal axis uses a logarithmic scale.

which are addressed in the next section. In general, amplitude decays and frequency rises as stiffness is increased. The dimensional frequency of the oscillating system is compared against the undamped natural frequency of the spring at each value of stiffness, portrayed by a gray line in Fig. 2(b). At low values of stiffness, Kθ = 0.003 and 0.3 (N m/rad), the response of the coupled system is far removed from the structure’s natural frequency and appears to be converging on a lower limit as stiffness goes to 0. This finding suggests the LCO is dominated by purely fluidic effects on softer structures as suggested by Poirel et al. (2008). This is an interesting observation in that small, aerodynamically driven oscillations on the control surface of a low Reynolds number or laminar flow wing may exacerbate the occurrence of free-play flutter (Dowell and Tang, 2002). At all spring constants, the coupled system is stiffer than that of the structure alone as indicated in Fig. 2(b) which shows that the dimensional frequency of the aeroelastic configuration always falls above the natural frequency of the structure. The frequency of the coupled system slowly merges toward the elastic mount’s natural frequency as the spring constant is increased. At Kθ = 3.0 (N m/rad) the coupled frequency is only slightly above the spring’s natural frequency. Oscillations are no longer sustained at stiffnesses beyond where these two curves merge. Fig. 2(c) shows the non-dimensional frequency (f + = fc /U∞ ) plotted against reduced velocity as a comparison of the various cases in a non-dimensional context. The values for the first-harmonic frequency range from f + = 0.04 in the softest case to 0.11 in the stiffest case. In comparison, the dominant frequency in the moment frequency spectra of the rigid, Rec = 110,000, α = 0◦ case reported in Barnes and Visbal (2016b) was f + = 5.6. This value, attributed to vortex shedding from the static wing, is nearly two orders of magnitude larger than the present LCO frequencies. The short time-scales associated with transitional flow structures are unlikely to provide a meaningful coupling effect with the slow time-scale of the structural response. Rather, as discussed in the previous paper (Barnes and Visbal, 2018c), the longer time-scales associated with boundary layer response and transition excursions on the airfoil surface are more closely related to the present phenomenon. The moment coefficient, taken about the elastic axis, is compared next. These measurements were processed in the same manner described in Barnes and Visbal (2016b,a) using both phase-averaging and an explicit 4th-order spatial low-pass filter. High-frequency content on the order of transitional/turbulent scales is not phase-locked and therefore intentionally removed by the phase-average and filter process. The remaining signal unmasks the broader trends of primary consequence to self-sustained oscillations from the high-frequency noise apparent in a single cycle. Fig. 3

6

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

Fig. 3. Effect of stiffness on moment coefficient vs. t/T at Rec = 110,000. The gray line shows a qualitative representation of the pitch deflection.

compares the phase-averaged and filtered moment coefficient against cycle-normalized time, t /T , for the four cases that exhibit flutter. For reference, t /T = 0 corresponds with the equilibrium point, θ = 0◦ , during pitch-up. Previously, Barnes and Visbal (2018c) described two nonlinear events common to all Reynolds numbers considered in that study for a spring constant of Kθ = 0.3 (N m/rad). The dominant feature was a brief but strong spike in moment coefficient that occurred when open separation transitioned and reattached to the surface near the end of each pitch stroke. This event tended to impose a short-lived structural damping effect, hence it is referred to as transition damping. Turbulent reattachment during this process formed a long LSB responsible for negative aerodynamic damping in the subsequent stroke. The second nonlinear event was a less severe moment peak that crested at intermediate angles of incidence. This feature tended to emerge as the LSB neared the trailing edge. Suction beneath the LSB in conjunction with a lengthening moment arm appeared to produce this brief excursion in the moment coefficient. At high Reynolds numbers, second-stage transition diminished the LSB which tempered its influence on the moment, whereas at low Reynolds numbers, reduced second stage transition permitted a stronger influence of the LSB. These two events were referred to as the ‘primary’ and ‘secondary’ peaks, respectively, in Barnes and Visbal (2018c). Both of the aforementioned features are clearly visible in the pitching moment of all four cases in Fig. 3. The dominant peak is annotated as (1) and the secondary peak is labeled as (2). The effect of stiffness on these salient events is similar to but inverse the effect of Reynolds number presented in Barnes and Visbal (2018c). As Reynolds number was decreased in Barnes and Visbal (2018c), the primary peak emerged later in the cycle and diminished in amplitude while the secondary peak became more pronounced and crested earlier in the cycle. Likewise, in Fig. 3 the primary peak is increasingly delayed and attenuates as stiffness increases to Kθ = 3.0 (N m/rad). Also like the effect of decreasing Reynolds number, the secondary peak becomes increasingly pronounced as stiffness is increased. Unlike the effect of diminishing Reynolds number, the secondary peak crests later in the cycle at high stiffness. Although the effects of structural rigidity and Reynolds number portray remarkable similarity, the latter observation indicates the two trends are not perfectly analogous. Regardless, the relation between timing of flow transition events appears to be paramount in both parameters. The frequency spectra of CM are compared between levels of rigidity in Fig. 4. The spectra for the soft and baseline airfoils take on the same general form. Namely, they exhibit a dominant peak which coincides with the frequency of oscillation and additional clearly defined peaks at odd-superharmonics of the pitching frequency. The stiffer cases, Kθ = 1.0 and 3.0 (N m/rad), take on a notably different form. As before, the dominant peak coincides with the oscillation frequency, but its amplitude is nearly matched by that of the 3rd superharmonic. At Kθ = 3.0 (N m/rad) only the 1st and 3rd harmonic frequencies are clearly defined; consistent with emergence of the secondary peak in Fig. 3 as a dominant feature in the stiffer cases. This trend is again similar to that observed in Barnes and Visbal (2018c) which showed the third superharmonic becomes increasingly relevant at low Reynolds numbers. A physical understanding for the aforementioned shifts in aeroelastic behavior is elucidated by observing the transient flow topology for a representative cycle of the LCO. Fig. 5 provides a quantitative comparison of the flow structure between cases using a space–time distribution of two quantities on the airfoil upper surface. The horizontal axis is chordwise position beginning at the leading edge and the vertical axis is cycle normalized time beginning at θ = 0◦ during pitch-up. The blue contour line represents zero skin friction coefficient, Cf = 0, filtered in the same manner as in Barnes and Visbal (2018c) and helps to identify laminar separation and approximate turbulent reattachment. The colored contours represent the magnitude of streamwise vorticity averaged across the span and then integrated along coordinate lines normal to the ∫ surface, ⟨|ωx |⟩ dn. This quantity helps identify the onset of spanwise instability and transition to turbulence. Together, these values help to compare the flow topology between cases and correlate timing of flow transition events. The threedimensional flow structure is shown using an iso-surface of Q-criterion at selected instances of time in Fig. 6 for further clarification. The lower surface flow topology portrays the same trends although out of phase with the upper surface and therefore is not presented for the sake of brevity.

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

7

Fig. 4. Effect of stiffness on pitching moment spectra sampled at a rate of ∆τ = 1 × 10−4 over an interval of τ = 150 for Rec = 110,000; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), (c) Kθ = 1.0 (N m/rad) and (d) Kθ = 3.0 (N m/rad).

∫

⟨|ωx |⟩ dn overlaid by a contour line of filtered Cf = 0 on the top surface of Rec = 110,000; (a) Kθ = Fig. 5. Space–time distributions of 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), (c) Kθ = 1.0 (N m/rad), and (d) Kθ = 3.0 (N m/rad). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The baseline case, Fig. 5(b), presents a typical flow topology for LSF. First, open laminar separation forms on the upper surface during early phases of pitch-up motion, t /T ∗16 > 12. This detached shear-layer becomes unstable, abruptly transitions to turbulence, and reattaches between t /T ∗16 = 2 and 3 and is visible in Fig. 6(b) at t /T = 3/16. This first excursion of turbulent flow, labeled as ‘T1’ in Figs. 5 and 6, was described as first-stage transition in Barnes and Visbal (2018c). This event is responsible the transition damping effect previously described. Flow transition recedes following its initial excursion and the newly formed long LSB is sustained by laminar reattachment, 2 < t /T ∗16 < 6, due to momentum transfer of spanwise-coherent structures, again visible at the corresponding phases of motion in Fig. 6(b).

8

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

Fig. 6. Snapshots of instantaneous Q-criterion for the Rec = 110,000 cases; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), (c) Kθ = 1.0 (N m/rad), and (d) Kθ = 3.0 (N m/rad).

Eventually, t /T ∗16 > 6, a second transition front intersects and interacts with the LSB shortening the bubble in the later phases of the pitch-down stroke. This larger excursion of transitional flow, labeled as ‘T2’ in Figs. 5 and 6, was referred to as second-stage transition in Barnes and Visbal (2018c). Juxtaposition of the LSB with the trailing edge appears to be responsible for the secondary peak described previously. Toward the bottom of the pitch-down stroke, t /T ∗16 > 9, laminar separation and flow transition wash off the airfoil and the upper surface briefly returns to attached laminar flow. The softer case, Kθ = 0.003 (N m/rad) in Figs. 5(a) and 6(a), follows a very similar pattern with a few modifications. First, transition and reattachment emerges a bit earlier in the softer case, consistent with the behavior of the dominant peak in Fig. 3. Second, the separation front penetrates further upstream as a consequence of the larger deflection achieved by this case. The transition front is more aggressive than in the baseline case due to both a larger peak deflection and longer timescale of oscillation. Consequently, second-stage transition encounters and truncates the LSB at an earlier phase of the cycle than for the baseline torsional stiffness. Flow patterns for the stiffer cases, Kθ = 1.0 and 3.0 (N m/rad) in Figs. 5(c, d) and 6(c, d), differ more significantly from those of the softer cases. Turbulent reattachment of open laminar separation emerges at increasingly later phases of the oscillation cycle as stiffness is increased. Emergence of the second-stage flow transition is largely inhibited in the stiffer cases and laminar reattachment prevails during the pitch-down stroke. The lower amplitude and shorter timescale of oscillations on these two cases keep the LSB closer to the trailing edge. This enhanced moment arm enforces a stronger influence on the moment coefficient. Actual effects of the aforementioned flow topologies are better conveyed by relating the moment coefficient to the pitch angle. Fig. 7 compares the phase-averaged and filtered pitching moment plotted against θ . The area between the pitch-up and pitch-down strokes represents net-cycle damping. Counterclockwise oriented regions depict net-positive damping and clockwise regions portray net-negative damping. The hysteresis curve is similar in nature between the softest and baseline cases, Kθ = 0.003 and 0.3 (N m/rad), again indicating the physical processes are similar between these two levels of rigidity. Both are dominated by a central clockwise-oriented region indicative of negative aerodynamic damping and capped at large angles of incidence by counterclockwise regions indicative of the transition damping phenomenon. Below a certain threshold for structural rigidity, hysteresis patterns appear to be fairly regularized. The responses of the stiffer cases, Kθ = 1.0 and 3.0 (N m/rad), are quite different. As before, the main features include primary peaks near the maximum and minimum angles of incidence and clear secondary peaks at intermediate angles of incidence. Contrary to the stiffer configurations, the secondary peaks are substantially pronounced. This secondary feature becomes so prominent that it results in a locally negative effective aerodynamic stiffness at small angles of incidence; that is, the general slope of the hysteresis loop (negative) locally reverses (becomes positive) in the central portion. Two factors appear to influence this modification. First, as described for the flow topology in Fig. 5, proximity of the LSB to the trailing edge enhances its influence on the moment. Second, both cases lack an appreciable second excursion of transition

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

9

Fig. 7. Pitching moment vs. pitch angle for Rec = 110,000; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), (c) Kθ = 1.0 (N m/rad), and (d) Kθ = 3.0 (N m/rad). Table 3 Computation parameters, (Rec = 77,000) Case Kθ Kθ Kθ Kθ

= 0.003 = 0.300 = 1.000 = 3.000

(N (N (N (N

m/rad) m/rad) m/rad) m/rad)

Run time

Time step

Span extent

Phase-avg. cycles

Disturbance

τ τ τ τ

∆τ ∆τ ∆τ ∆τ

s/c s/c s/c s/c

5 10 15 N/A

No No Yes Yes

= 300 = 300 = 300 = 100

= 1 × 10−4 = 1 × 10−4 = 1 × 10−4 = 1 × 10−4

= 0.2 = 0.2 = 0.2 = 0.2

in Fig. 5. Therefore, the LSB is largely preserved by laminar reattachment after its initial formation. This behavior indicates that a prevalence of laminar reattachment may promote an outsized influence of the secondary peak. The stiffest configuration, Kθ = 3.0 (N m/rad), presents further modifications to the aforementioned trends. Orientation of the hysteresis curve is entirely clockwise but very narrow at small angles of incidence. In this case, the primary peak actually contributes to a positive net energy transfer to the structure. This effect is a consequence of the actual timing of this event relative to the cycle as evinced by Fig. 3. The dominant peaks are delayed into the subsequent strokes in which their contribution reinforces the direction of travel rather than provide a damping effect. Both of these cases are more similar to the Rec = 77,000 case shown in Barnes and Visbal (2018c)—also depicted in the following section—than they are to the baseline case at Rec = 110,000. 4.2. Rec = 77,000 This section reinforces the observations and conclusions of the previous section by analyzing the effects of stiffness for the lower Reynolds number of Rec = 77,000. The four lowest torsional spring constants were considered in this section as presented in Table 3. Time histories of pitch deflections for all four cases are shown in Fig. 8. Aeroelastic oscillations were sustained for the smallest three values of Kθ = 0.003, 0.3, and Kθ = 1.0 (N m/rad). LCO for the Kθ = 1.0 (N m/rad) case only emerged after exposure to a disturbance in the form of Eq. (4). Oscillations were not sustained for the stiffest case, Kθ = 3.0 (N m/rad), regardless of the disturbance supplied. Pitch amplitude and frequency of oscillation are presented alongside those of the Rec = 110,000 cases for comparison in Fig. 2. Generally, the same trends are apparent for both Reynolds numbers with a few notable distinctions. First, it is interesting to note that amplitude and non-dimensional frequency are very similar for the most compliant airfoil, Kθ = 0.003 (N m/rad). Amplitude decays faster with increasing stiffness at the lower Reynolds number. Consistent with the previous section, the dimensional frequency of the coupled system at Rec = 77,000 is stiffer than that of the elastic mount but lies below the coupled frequencies of the Rec = 110,000 cases. The coupled frequency of the lower Reynolds number case merges with the structural natural frequency at a lower stiffness. Consequently, LCO is eliminated on a less rigid structure in lower Reynolds number flows. This trend would also suggest that oscillations did not occur at higher stiffnesses as Reynolds number is increased. However, there appears to be a limit to this effect. Preliminary simulations at Rec = 200,000 found that oscillations were not sustained indefinitely at Kθ = 3.0 (N m/rad) which was the largest stiffness to sustain oscillations at Rec = 110,000. Fig. 2(c) shows the non-dimensional first-harmonic frequency against reduced velocity. It is interesting to note that the relation between these two non-dimensional parameters is remarkably similar for both Reynolds numbers shown.

10

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

Fig. 8. Time history of pitching oscillations for the first τ = 200 at Rec = 77,000; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), (b) Kθ = 1.0 (N m/rad), and (d) Kθ = 3.0 (N m/rad).

Fig. 9. Amplitude vs. non-dimensional frequency f + = fc /U∞ for all cases.

Fig. 9 shows the mean amplitude plotted against the non-dimensional frequency, f + , for the Rec = 77,000 and 110,000 cases as well as the values for other Reynolds numbers at Kθ = 0.3 (N m/rad) from Barnes and Visbal (2018c). It is easy to notice that oscillation amplitude and non-dimensional frequency roughly share an inverse power relation in the form of |θ| = 0.3/f + as indicated by the wide gray band in Fig. 9 which was obtained from a least squares fit of the available data. If a minimum survivable amplitude could be identified, this relation might serve as a guideline for structural design. For example, if the minimum survivable amplitude for LSF were later identified to be θmin = 2.5◦ —this choice is purely arbitrary—a structural stiffness of Kθ = 12.0 (N m/rad) at a Reynolds number of Rec = 200,000 to match the predicted aeroelastic frequency and eliminate LSF. LCO at each Reynolds number appears to terminate at a different minimum amplitude which suggests the this criteria would be sufficient to inhibit LCO but overly conservative in most cases. The phase-averaged and filtered moment coefficient is plotted against cycle normalized time in Fig. 10(a) for the Rec = 77,000 cases. Moment coefficient exhibits a similar evolution with increasing stiffness as observed previously for the higher Reynolds number. The dominant spikes in moment are increasingly delayed as stiffness rises, secondary peaks become more pronounced, and superharmonic content, Fig. 11, becomes increasingly relevant as stiffness is increased. Perhaps most interesting is the comparison between the Rec = 77,000 and 110,000 cases in Fig. 10(b) for the softest elastic mount, Kθ = 0.003 (N m/rad). As noted previously, the amplitude and frequency for both Reynolds numbers are very similar. Likewise, the moment coefficient takes on an almost identical form except that the secondary peak is more prominent and emerges earlier at the lower Reynolds number. This similarity further indicates some degree of saturation of LSF when stiffness is sufficiently small. Frequency spectra for the moment coefficient are presented in Fig. 11 for each stiffness. As before, the softest case exhibits a dominant peak at the oscillation frequency and multiple smaller clearly defined peaks at odd-superharmonics. As stiffness is increased, the relative significance of these odd-superharmonics increases such that they nearly match the amplitude of the fundamental frequency. This effect is the same as observed in the previous section. Fig. 12 shows space–time distributions of the upper surface flow topology in the same manner as in Fig. 5 as well as the three-dimensional flow structure in Fig. 13. As a further indication of saturation in flutter behavior, the separation and transition patterns on the most compliant airfoil, Kθ = 0.003 (N m/rad), are exceptionally similar to those of the Rec = 110,000 case at the same stiffness. At Kθ = 0.3 and 1.0 (N m/rad), the flow topologies strongly resemble those of the stiffer cases at Rec = 110,000. Namely, the Kθ = 0.3 and 1.0 (N m/rad) cases experience a lack of second-stage

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

11

Fig. 10. Pitching moment vs. time; (a) effect of stiffness at Rec = 77,000 and (b) comparison between Rec = 77,000 and Rec = 110,000 at Kθ = 0.003 (N m/rad). The gray line shows a qualitative representation of the pitch deflection.

Fig. 11. Pitching moment vs. pitch angle; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad) at Rec = 77,000.

Fig. 12. Space–time distributions of ⟨|ωx |⟩ dn overlaid by a contour line of filtered Cf = 0 for Rec = 77,000; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), and (c) Kθ = 1.0 (N m/rad). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

∫

12

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

Fig. 13. Snapshots of instantaneous Q-criterion at Rec = 77,000; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad), and (c) Kθ = 1.0 (N m/rad).

Fig. 14. Pitching moment vs. pitch angle at Rec = 77,000; (a) Kθ = 0.003 (N m/rad), (b) Kθ = 0.3 (N m/rad) and (c) Kθ = 1.0 (N m/rad).

transition and prevalence of laminar reattachment during the later stages of pitch-down. Turbulent reattachment during first-stage transition occurs progressively later in the cycle. These same effects occurred for Rec = 110,000 at higher spring constants of Kθ = 1.0 and 3.0 (N m/rad). Similarity in the effect of stiffness is further indicated in the hysteresis curves presented in Fig. 14. The lowest torsional spring constant yields a typical LSF hysteresis consisting of a central clockwise region capped by counterclockwise portions just before the end of each stroke. The clockwise portion differs from that of the higher Reynolds number in Fig. 7(a) in that it survives a smaller range of deflections and the moment coefficients are of higher magnitude. As stiffness is increased, Fig. 14(b), the central clockwise region skews toward a positive slope due to proximity of the LSB to the trailing edge. Further, the damping phases of motion (primary peaks) subsequently diminish due to later onset of transition damping. This trend continues in Fig. 14(c) in that the effect of transition and reattachment contributes to negative aerodynamic damping as was also the case in Fig. 7(d). This self-similarity is likely limited to the lower half of the flutter regime. Preliminary tests at Re = 200,000 exhibit hysteresis patterns that are saturated across a range of stiffness values but do not converge onto that shown for the lower Reynolds number cases.

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

13

5. Conclusions In this work, high-fidelity ILES was implicitly coupled with structural dynamics to explore limit-cycle-oscillations on an airfoil operating at transitional Reynolds numbers. This phenomenon is a 1-DOF limit cycle oscillation not well understood to date and distinct from other well known flutter phenomena. A NACA0012 wing section elastically mounted in pitch at x/c = 0.186 was considered over a wide range of torsional stiffnesses, 0.003 N m/rad ≤ Kθ ≤ 30.0 N m/rad at two Reynolds numbers of Rec = 77,000 and 110,000 which fall on the lower end of the flutter regime. LSF was observed to occur over a wide range of torsional rigidity for both Reynolds numbers. From a non-dimensional standpoint, LSF appeared to saturate into a common pattern when stiffness was sufficiently small. The coupled aeroelastic system was stiffer than the elastic mount as indicated by the dimensional frequency. As stiffness of the elastic mount was increased, the aeroelastic and structural frequencies began to merge. Once the aeroelastic and structural frequencies merged for either Reynolds number, oscillations were no longer sustained. The effects of stiffness on flow structure and loading for the cases presented were remarkably similar to but inverse the effects of Reynolds number reported in the previous work. At high stiffness, turbulent reattachment occurred at later phases of the cycle. The subsequent LSB became increasingly dominated by laminar – rather than turbulent – reattachment. This shift in flow structure resulted in an increasingly non-linear loading cycle. Likewise, the same effects were observed to occur when Reynolds number was decreased at a fixed stiffness in the preceding work. At either Reynolds number, LSF exhibited a fairly regularized pattern of laminar separation, turbulent reattachment, and transient LSB at very low stiffness The relative timing of flow transition events appears to be the common factor between stiffness and Reynolds number. When the natural frequency of the structure was higher than that driven by LSF, the flow transition processes began to lag behind the structural response. Because flow transition was slower at the lower Reynolds number, as indicated by the lower dimensional frequency of the aeroelastic response, it became sensitive to the effects of stiffness at a lower structural frequency accounting for an earlier termination of self-sustained oscillations. The higher Reynolds number boundary layer adjusted more rapidly allowing aerodynamic effects to dominate and sustain oscillations at a higher torsional stiffness. This observation would suggest that oscillations at even higher Reynolds numbers should survive on stiffer structures. However, preliminary tests at Rec = 200,000 indicate there is a limit to this effect as this Reynolds number did not sustain oscillations at the maximum stiffness achieved at Rec = 110,000. The effect of stiffness at Reynolds numbers on the higher end of the flutter regime (those beyond the maximum amplitude) remains a topic in need of further investigation but outside the scope of the present work. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This material is based upon work supported by the Air Force Office of Scientific Research, United States under a Lab Task monitored by Dr. G. Abate and by a grant of HPC time from the DoD HPC Major Shared Resource Centers at AFRL, United States and ERDC, United States. References Alpert, P., 1981. Implicit filtering in conjunction with explicit filtering. J. Comput. Phys. 44 (1), 212–219. http://dx.doi.org/10.1016/0021-9991(81) 90047-4. Barnes, C.J., Visbal, M.R., 2016a. Aeroelastic response of an airfoil at transitional Reynolds numbers, AIAA paper 2016-3634. AIAA http://dx.doi.org/ 10.2514/6.2016-3634. Barnes, C.J., Visbal, M.R., 2016b. High fidelity LES simulations of self-sustained pitching oscillations on a NACA0012 airfoil at transitional Reynolds numbers, AIAA paper 2016-1353. AIAA http://dx.doi.org/10.2514/6.2016-1353. Barnes, C.J., Visbal, M.R., 2018a. Gust response of rigid and elastically mounted airfoils at a transitional Reynolds number. Aerosp. Sci. Technol. 74, 112–119. http://dx.doi.org/10.1016/j.ast.2017.12.025. Barnes, C.J., Visbal, M.R., 2018b. Mitigation of laminar spearation flutter using plasma-based actuators, AIAA paper 2018-3523. AIAA http://dx.doi. org/10.2514/6.2018-3523. Barnes, C.J., Visbal, M.R., 2018c. On the role of flow transition in laminar separation flutter. J. Fluids Struct. 77, 213–230. http://dx.doi.org/10.1016/j. jfluidstructs.2017.12.009. Beam, R., Warming, R., 1978. An implicit factored scheme for the compressible Navier-Stokes equations. AIAA J. 16 (4), 393–402. http://dx.doi.org/ 10.2514/3.60901. Dowell, E.H., Tang, D., 2002. Nonlinear aeroelasticity and unsteady aerodynamics. AIAA J. 40 (9), 1697–1707. http://dx.doi.org/10.2514/2.1853. Gaitonde, D.V., Visbal, M.R., 1998. High-Order Schemes for Navier-Stokes Equations: Algorithm and Implementation Into FDL3DI. Technical report AFRL-VI-WP-TR-1998-3060, Air Force Research Laboratory, Wright-Patterson AFB. Gaitonde, D., Visbal, M., 1999. Further development of a navier-stokes solution procedure based on higher-order formulas, AIAA paper 1999-0557. AIAA http://dx.doi.org/10.2514/6.1999-557. Garmann, D.J., Visbal, M.R., Orkwis, P.D., 2013. Comparative study of implicit and subgrid-scale model large-eddy simulation techniques for low-Reynolds number airfoil applications. Internat. J. Numer. Methods Fluids 71 (12), 1546–1565. http://dx.doi.org/10.1002/fld.3725.

14

C.J. Barnes and M.R. Visbal / Journal of Fluids and Structures 91 (2019) 102767

Gibbs, S.C., Fichera, S., Zanotti, A., Ricci, S., Dowell, E.H., 2014. Flow field around the flapping flag. J. Fluids Struct. 48, 507–513. http://dx.doi.org/10. 1016/j.jfluidstructs.2014.02.011. Golubev, V.V., Nguyen, L., Counts, M., MacKunis, W., Pedroza, N.R., Pasiliao, C.L., 2014. Robust nonlinear control of airfoil gust-induced limit cycle oscillations using synthetic jet actuators, AIAA paper 2014-2936. AIAA http://dx.doi.org/10.2514/6.2014-2936. Gordnier, R., Visbal, M., 2002. Development of a three-dimensional viscous aeroelastic solver for nonlinear panel flutter. J. Fluids Struct. 16 (4), 497–527. http://dx.doi.org/10.1006/jfls.2000.0434. Jameson, A., Schmidt, W., Turkel, E., 1981. Numerical solutions of the Euler equations by finite volume methods using Ruge-Kutta time stepping schemes, AIAA paper 1981-1259. AIAA http://dx.doi.org/10.2514/6.1981-1259. Lele, S., 1992. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 103 (1), 16–42. http://dx.doi.org/10.1016/00219991(92)90324-R. Mathew, J., Lechner, R., Foysi, H., Sesterhenn, J., Friedrich, R., 2003. An explicit filtering method for large eddy simulation of compressible flows. Phys. Fluids 15 (8), 2279–2289. http://dx.doi.org/10.1063/1.1586271. Morton, S.A., Melville, R.B., Visbal, M.R., 1998. Accuracy and coupling issues of aeroelastic Navier-Stokes solutions on deforming meshes. J. Aircr. 35 (5), 798–805. http://dx.doi.org/10.2514/2.2372. Nguyen, L., Golubev, V.V., MacKunis, W., Pedroza, N.R., Pasiliao, C.L., 2015. High-accuracy simulations of robust LCO control using synthetic jet actuatorsx, AIAA paper 2015-0809. AIAA http://dx.doi.org/10.2514/6.2015-0809. Poirel, D., Harris, Y., Benaissa, A., 2008. Self-sustained aeroelastic oscillations of a NACA0012 airfoil at low-to-moderate Reynolds numbers. J. Fluids Struct. 24 (5), 700–719. http://dx.doi.org/10.1016/j.jfluidstructs.2007.11.005. Poirel, D., Mendes, F., 2014. Experimental small-amplitude self-sustained pitch-heave oscillations at transitional Reynolds numbers. AIAA J. 52 (8), 1581–1590. http://dx.doi.org/10.2514/1.J052541. Poirel, D., Métivier, V., Dumas, G., 2011. Computational aeroelastic simulations of self-sustained pitch oscillations of a NACA0012 at transitional Reynolds numbers. J. Fluids Struct. 27, 1262–1277. http://dx.doi.org/10.1016/j.jfluidstructs.2011.05.009. Poirel, D., Yuan, W., 2010. Aerodynamics of laminar separation flutter at a transitional Reynolds number. J. Fluids Struct. 26, 1174–1194. http: //dx.doi.org/10.1016/j.jfluidstructs.2010.06.005. Pulliam, T., 1986. Artificial dissipation models for the Euler equations. AIAA J. 24 (12), 1931–1940. http://dx.doi.org/10.2514/3.9550. Pulliam, T., Chaussee, D., 1981. A diagonal form of an implicit approximate-factorization algorithm. J. Comput. Phys. 39 (2), 347–363. http: //dx.doi.org/10.1016/0021-9991(81)90156-X. Rai, M.M., 1989. Three-dimensional Navier-Stokes simulations of turbine rotor-stator interaction; part i - methodology. J. Propul. Power 5 (3), 305–311. http://dx.doi.org/10.2514/3.23154. Rai, M.M., Chakravarthy, S.R., 1986. An implicit form for the Osher upwind scheme. AIAA J. 24 (5), 735–743. http://dx.doi.org/10.2514/3.9340. Stolz, S., Adams, N., 1999. An approximate deconvolution procedure for large-eddy simulation. Phys. Fluids 11 (7), 1699–1701. http://dx.doi.org/10. 1063/1.869867. Visbal, M.R., Gaitonde, D.V., 1999. High-order-accurate methods for complex unsteady subsonic flows. AIAA J. 37 (10), 1231–1239. http://dx.doi.org/ 10.2514/2.591. Visbal, M.R., Morgan, P.E., Rizzetta, D.P., 2003. An implicit LES approach based on high-order compact differencing and filtering schemes, AIAA paper 2003-4098. AIAA http://dx.doi.org/10.2514/6.2003-4098. Visbal, M.R., Rizzetta, D.P., 2002. Large-eddy simulation on curvilinear grids using compact differencing and filtering schemes. J. Fluids Eng. 124 (4), 836–847. http://dx.doi.org/10.1115/1.1517564. Yuan, W., Poirel, D., Wang, B., 2013. Simulations of pitch-heave limit-cycle oscillations at a transitional Reynolds number. AIAA J. 51 (7), 1716–1732. http://dx.doi.org/10.2514/1.J052225. Yuan, W., Poirel, D., Wang, B., Benaissa, A., 2015. Effect of freestream turbulence on airfoil limit-cycle-oscillations at transitional Reynolds numbers. J. Aircr. 52 (4), 1214–1225. http://dx.doi.org/10.2514/1.C032807.