Assessment of a numerical design tool for pitching airfoils

Abstract

The aim of this study was to assess the applicability of a numerical design tool for pitching airfoils; by evaluating their aerodynamic performance comparing experimental and numerical data at different levels of complexity. Experimental findings of non-harmonically pitching airfoil configurations in a water tunnel at medium Reynolds numbers are compared to solutions of a modified double-wake vortex-panel method with boundary-layer formulation and unsteady Reynolds-averaged Navier–Stokes simulations. For steady-state airfoil configurations at high angle of attack, large eddy simulation data are also considered. Rigorous analysis of the sources of error and uncertainty in numerical methods and measurements proves that from an engineering point of view, a low-cost double-wake model is well suited for the design process of, for example, wind-turbine blades. That is, the overall first-order error and uncertainty of the double-wake panel method with a strong viscous-inviscid interaction model are in a reasonable range against the measurements and the much more expensive numerical methods. Thus, the low computational costs plus the accuracy make the unsteady panel-boundary-layer method a useful and reliable tool to analyze the aerodynamics of pitching airfoils in the design process.

Keywords

Pitching airfoil high-speed pressure-probe measurement double-wake panel method unsteady Reynolds-averaged Navier–Stokes simulation uncertainty analysis

Introduction

The precise prediction and understanding of dynamic stall and unsteady flow separation is of major importance to improve the performance of numerous technical applications such as wind turbines, helicopters, or flapping wings. For helicopters operating at very high tip-speed ratios, local dynamically separating and reattaching flow regions constitute a major source of vibration loads (Leishman, 2006). For flapping wings with a combined pitch-plunge motion, the description of the temporal lift development is determined by the dynamically separating flow near the leading edge of the wing (Hedenström et al., 2007). This strong stall is characterized by the formation of large vortical structures shed by the separated boundary layer on the wing surface. Horizontal axis wind turbines (HAWTs) and vertical axis wind turbines (VAWTs) may operate in flow regimes where dynamic stall and spectral loads occur. Moreover, in VAWT applications, the dynamic stall occurring at low-to-moderate tip-speed ratios significantly affects the near-wake development whose understanding is crucial for the overall performance and reliable operation of the turbine. Complex three-dimensional vortex systems are formed affecting the circulation distribution and pressure gradient on the blade mainly in the downwind passage. However, the effects of these vortex systems are still felt in the upwind passage determining the maximum tangential force (Ferreira et al., 2014). Due to the rapid change of wind conditions in conjunction with the comparatively slow response time of pitch-regulated HAWT and VAWT machines, dynamic stall is inevitable (Danao et al., 2014; Lee and Wu, 2013).

In numerous studies, the processes that characterize weak and strong dynamic stall regimes, as well as the influence of parameters such as airfoil geometry, reduced frequency, and Reynolds and Mach numbers, are described (Corke and Thomas, 2015; Lee and Gerontakos, 2004; Leishman, 2006; McCroskey, 1982). However, a unified understanding of the physics of the pitching airfoil flow field does not exist. In fact, the understanding of many aspects of dynamic stall inception even in incompressible flow remains incomplete (Corke and Thomas, 2015).

The computationally least expensive method to describe dynamic stall on a wing is the lifting-line model (Leishman and Beddoes, 1986, 1989). Although low in computational costs, the semi-empirically parameterized dynamic stall model of Leishman and Beddoes provides useful first results of the dynamic integral loads if standard conditions and profiles, for example, National Advisory Committee for Aeronautics (NACA) profiles, are chosen or if the coefficient matrix is tuned a priori to match reference findings. The lack of generality, for example, the uncertainty of the model coefficients for a wide range of applications, requires a careful application of the lifting-line method. Unsteady Reynolds-averaged Navier–Stokes (URANS) simulations are widely accepted as a standard computational fluid dynamic (CFD) approach in industrial applications. For the investigation of dynamic stall, however, the dynamic loads are highly susceptible to the choice of the turbulence model (Ekaterinaris and Platzer, 1998) since it is a challenge to properly describe incipient trailing-edge separation even at steady flow conditions (Celić and Hirschel, 2006). Even two or higher equation models do not necessarily increase the confidence in the quantitative prediction of the resulting dynamic loads (Ekaterinaris and Platzer, 1998; Richter et al., 2011). That is, a major uncertainty of the outcome is due to the applied turbulence closure in the RANS and URANS simulations, which has primarily been developed for attached boundary layers or generic mixing processes and due to the error caused by the applied numerical scheme and grid quality. On top of that, for a first-stage design procedure, where a vast number of flow cases must be investigated, the URANS method is still computationally too demanding. The detached eddy simulation (DES) concept (Deck and Laraufie, 2013) is neither attractive as a design tool since the underlying RANS formulation suffers from the same disadvantages as those of pure URANS simulations at even higher computational costs. Large eddy simulations (LES) representing a higher order turbulence modeling approach show a very low-level uncertainty at much higher computational costs compared to the aforementioned methods. Therefore, the LES approach is not appropriate in the (pre-)design stage to tackle large-scale unsteady flow problems.

In this study, the focus is on panel-boundary-layer models which fall in the class of vortex methods. Their applicability to the simulation of the aforementioned unsteady flow phenomena at a sufficient accuracy and low computational costs in the early design stage is evaluated. These models have successfully been applied to steady aerodynamic problems (Drela and Giles, 1987) and have been extended to unsteady problems (Riziotis and Voutsinas, 2008; Voutsinas, 2006; Voutsinas and Riziotis, 1999). Relying on a potential flow description, the meshless vortex or particle methods are able to describe the flow conditions around a pitching airfoil using source/sink and vorticity distributions. If wall-bounded and separating shear layers are to be described by the particle method, the computational costs are comparable to that of URANS methods due to the required number of particles near the wall to satisfy the no-slip condition (Koumoutsakos and Leonard, 1995). Using a semi-empirical turbulence model to locate the flow separation, the vortex model, which is then a panel-boundary-layer method, retains its computational efficiency compared to URANS models (Riziotis and Voutsinas, 2008; Zanon et al., 2013). In terms of load prediction, the order of the turbulence closure, that is, using zero-, one-, or two-equation models to describe viscous flow effects, is not of major importance if vortex methods are used (Riziotis and Voutsinas, 2008). Note that such methods have been successfully applied to describe the flow around pitching airfoils (Riziotis and Voutsinas, 2008; Zanon et al., 2013).

The target of this work is to substantiate the confidence in the application of the methods being computationally less expensive compared to URANS, DES, or LES, such as the panel-boundary-layer formulation, by analyzing the errors and uncertainties of this method and by discussing the findings against URANS results. The focus is on the design process of airfoils under pitching conditions. This process should be reliably executed with a panel-boundary-layer method without reducing the confidence in resulting aerodynamical loads compared to computationally more expensive methods.

In the following, a modified panel-boundary-layer code, a standard RANS model, and an LES are used to first investigate the steady flow separation at high angles of attack. The results are compared to experimental data and discussed in detail. Then, the panel-boundary-layer formulation is applied to a set of pitching airfoil problems and the solutions are analyzed using the corresponding experimental data and URANS. Since the necessary details of the experiments and the numerical configurations are known, first-order errors and uncertainties of each configuration can be assessed and thoroughly compared. Furthermore, a juxtaposition with URANS and panel-boundary-layer simulations in terms of computational costs is given.

Methods

Experimental setup

Numerous measurements were performed in a water tunnel to obtain reliable experimental data to validate the numerical results. The test section has an area of 1.0 m × 1.5 m. The water tunnel operates at a maximum flow velocity of 1.7 m/s which leads to a maximum Reynolds number per meter of about 1.7 × 10⁶. The turbulence intensity measured by particle image velocimetry (PIV) is about 0.4%. In the experiments, the bulk velocity was set as 1.1 m/s such that a Reynolds number per meter of 1.4 × 10⁶ was obtained. The measurements were performed on the upper and lower sides of a 0.4-m chord NACA0018 airfoil which led to a two-dimensional (2D) geometric blockage of about 5.0%. The spanwise extension of the airfoil was 1 m. No boundary-layer tripping was used. A total of 64 Qlite sensors which are clustered near the leading and the trailing edge on the upper and lower sides of the airfoil guaranteed a sufficiently high spatial resolution of the pressure distribution. The measurement frequency was 1000 Hz which is about two orders of magnitude higher than the frequency of the time scales of the relevant force coefficients. No shear-stress measurements were conducted since the focus of this work was on the experimental investigation of the aerodynamics of pressure-dominated separated flows. The data acquisition and processing were performed using the commercial software LabView™ and MATLAB™. The experimental setup is presented in Figure 1. The boundary conditions of the measurement campaign in the water tunnel are precisely known such that the error sources and uncertainties can be determined.

Figure 1.

Schematic diagram of the experimental setup.

Numerical setups

Next, the numerical methods to investigate the steady and unsteady flow fields around airfoils are concisely discussed. First, the double-wake panel method, denoted as 2D-AISE (aeroelastic integrated simulation environment), is described. Then, the RANS/URANS and LES methods are presented.

Panel method

The superposition of elementary solutions is an effective approach for the simulation of the incompressible irrotational flow past an airfoil (Katz and Plotkin, 2001). In the present model, several singularities are used: panels of length l with constant source strengths on the airfoil surface and constant vorticity strength on the airfoil, the near wake, and the separation panel. In the far wake, the vortex panel is changed into a vortex blob which is characterized by vortex strength $Γ_{F V M}$ . The details of the model will be given for the unsteady separated flow condition since the attached flow condition can be regarded as a particular case of the separated flow. Figure 2 shows the distribution of singularities in the attached and separated flow regions. The surface of the airfoil is approximated by N panels with source distributions of constant strength per panel $σ_{i}$ and two circulations per unit length of strength $γ_{1}$ and $γ_{2}$ . The wakes originating at the trailing edge and the separation point are modeled as a collection of vortex blobs, except for the parts close to the airfoil that are discretized by vortex segments described by trailing edge and separation panels (near-wake panels) with length $l_{t e}$ and $l_{s p} .$ The near-wake panel at the blades’ trailing edges is initially aligned to coincide with the local flow determined by the current inviscid solution (outer loop). Its initial orientation is defined by the local velocity at the airfoil trailing edge ${\vec{v}}_{t e, i n v i s c i d},$ while its length $l_{t e}$ equals $| Δ t \times {\vec{v}}_{t e, i n v i s c i d} . |$ After the first iteration step of the boundary-layer method (center loop), a second panel, referred to as separation panel, is introduced starting at the calculated separation point. Since the velocity at the separation point ${\vec{v}}_{s p}$ is constrained to be zero, an auxiliary velocity ${\vec{v}}_{s p^{'}}$ is constructed by neglecting the influence of the singularities on the airfoil surface. The quantity ${\vec{v}}_{s p^{'}}$ is used to determine the initial orientation and length of the separation panel, whose circulation equals that of the upstream panel next to the separation point location. Using an iterative scheme (inner loop), the orientation and length of both near-wake panels and the strength of the singularities on the airfoil are adjusted to match the viscous solution of the boundary-layer method. For a detailed description of this interaction, the reader can refer to the study of Zanon (2011). Once the boundary-layer method is converged, the circulations around the wake panel and the separation panel are assigned to the corresponding vortex blobs. Subsequently, the displacement of each vortex blob is determined, which equals the product of the local velocity, that is, the sum of the freestream and the induced velocities at its position, and the time step. Further details on the panel method and the double-wake formulation can be found in the studies of Basu and Hancock (1978), Riziotis and Voutsinas (2008), and Zanon (2011).

Figure 2.

Panel source/vortex distribution and coordinate systems for double-wake configuration.

Boundary-layer method and viscous inviscid coupling

To determine the boundary layer, the integral formulation for steady flows developed by Drela (1989) was chosen. As stated by Zanon (2011), due to the comparatively simple implementation of this steady model, this has rather been chosen than the fully unsteady boundary-layer formulation of Riziotis and Voutsinas (2008) since the latter increases the complexity of the calculation considerably. It has been shown that the present model based on the quasi-steady boundary-layer approach constitutes a reasonable approximation to describe a comparatively slow airfoil motion. Similar to the results of Zanon (2011), it will be shown in section “Results” that the present boundary-layer approach leads to a good agreement with experimental results in the investigated range of reduced frequencies.

This approach for steady incompressible flow is based on von Kármán’s integral equation and the kinetic energy shape parameter equation (Kármán, 1921; Whitfield, 1978); see equations (1) and (2), which are solved together with different empirical relations depending on the flow condition, that is, laminar, turbulent, or free-wake boundary layer (Drela and Giles, 1987). Falkner–Skan one-parameter velocity profiles are used for laminar closure. The turbulent closure equations include an additional parameter, that is, the shear-stress coefficient, which is computed by integrating the so-called shear-lag equation. For the laminar-turbulent transition, the present model uses the simplified version of the $e^{n}$ method (Drela and Giles, 1987). For the free wake, the skin-friction coefficient is set to zero to close the equations. To efficiently couple the solution of the panel method and that of the boundary layer, the viscous-inviscid coupling method developed by Zanon (2011) is applied. In this method, the transpiration velocity which is imposed as boundary condition in the panel scheme is determined by a sub-iterating scheme which compares the parallel flow velocities at each panel of the inviscid panel method and the boundary-layer solution

\frac{d Θ}{d ξ} + (2 + H + H^{*} - M_{e}^{2}) \frac{Θ}{u_{e}} \frac{d u_{e}}{d ξ} = \frac{C_{f}}{2}

(1)

Θ \frac{d H^{*}}{d ξ} + (2 + H^{* *} + H^{*} (1 - H)) \frac{Θ}{u_{e}} \frac{d u_{e}}{d ξ} = 2 C_{D} - H^{*} \frac{C_{f}}{2}

(2)

In the above equations, the quantities $Θ$ , $H$ , $H^{*}$ , and $H^{* *}$ denote the momentum thickness, shape parameter, kinetic energy shape parameter, and the density flux shape parameter, repectively; $ξ$ , $M_{e}$ , and $M_{e}$ define the streamwise coordinate, edge Mach number, and edge velocity, respectively; and $c_{f}$ and $c_{D}$ are the skin-friction and dissipation coefficients.

Dynamic stall

At dynamic separation, for example, induced by a pitching motion of the airfoil, the integral boundary-layer method has to be extended to meet the increased complexity of the recirculating unsteady flow field. The computation of the skin-friction coefficient $c_{f}$ is extended by taking into account convective velocities of shed vortices to avoid artificial oscillations of the point of separation induced by single vortex elements that are trapped in the recirculation zone. Furthermore, the non-physical high-frequency oscillations of the separation point are damped using a Butterworth low-pass filter on the time evolution of the skin-friction distribution. Considering the uncertainty in determining the separation locations by RANS turbulence models and experimental investigations (Celić and Hirschel, 2006; Roidl et al., 2013), the separation criterion $c_{f, s e p} = 0$ is slightly relaxed by applying $c_{f, s e p} = C_{s e p} \cdot c_{f, e q},$ where $C_{s e p} < < 1$ is a constant. The quantity $c_{f, e q}$ is the skin-friction coefficient of the virtual equilibrium boundary layer computed by the Reynolds number $R e_{Θ}$ based on the momentum thickness $Θ$ . In this study, $C_{s e p}$ was set to be 0.05. This modification is reflected in the sources for the quantification of the uncertainty of the method.

URANS and LES methods

The 2D-URANS simulations are computed using Star-CD^® plus the shear-stress transport (SST) k − $\underset{.}{ω}$ model for turbulent closure (Menter, 1994). A second-order accurate Monotonic Upwind Scheme for Conservation Laws approach and three-stage Runge–Kutta method for the spatial and temporal discretization are used. For steady airfoil configurations, the numerical reference data are computed by LES. A detailed description of the fundamental LES solver is given by Meinke et al. (2002). The code was successfully applied to analyze in detail separating subcritical and supercritical boundary layers at positive and negative pressure gradients (Roidl et al., 2013, 2014).

Problem definition

Next, the essential parameters of the various numerical formulations are defined and summarized in Tables 1 and 2. For the panel method, the number of panels for the steady flow case and the pitching motion was 180. They were clustered near the leading and trailing edge and the time step was 0.05 $c / U$ _∞ A merging curtain is implemented at $x / c = 10$ to control the overall number of vortex cores in the computational domain.

Table 1.

Computational domain, grid resolution, time step, and total number of mesh points for LES and RANS computations at $R e_{c} = 5.5 \times 10^{5} .$

	Domain size		Resolution				No. of grid points
	$L_{f a r f i e l d}$	$L_{s p a n w i s e}$	$Δ x^{+}$	$Δ y_{w a l l}^{+}$	$Δ z^{+}$	∆t (c / U_∞)	Total
LES	25c	0.091c	25	2	15	10⁻⁶	35 × 10⁻⁶
RANS, 2D, steady	25c	0.01	200	1.0	100	10⁻³	70 × 10⁻³
URANS, 2D, pitching	25c	0.01	200	1.0	100	10⁻⁴	70 × 10⁻⁴

LES: large eddy simulation; RANS: Reynolds-averaged Navier–Stokes simulation; 2D: two-dimensional; URANS: unsteady Reynolds-averaged Navier–Stokes similation.

Table 2.

Spatio-temporal configuration for panel method.

	Merging curtain	Temp. resolution	Spatial resolution		No. of panels
	$D_{m e r g . c u r t a i n}$	$t_{p a n e l}$	$Δ s_{m a x}^{+}$	$Δ s_{m i n}^{+}$	Total
2D-AISE, steady	10c	0.05 $c / U_{\infty}$	10	180	150
2D-AISE, pitching	10c	0.025 $c / U_{\infty}$ c / U_∞	10	180	150

AISE: aeroelastic integrated simulation environment.

A c-shaped grid was applied in the URANS formulation. The extension in the streamwise, the wall-normal, and the spanwise direction was $25 c$ , $25 c$ , and $0.01 c$ , respectively, where $c$ defines the chord length. The grid resolution on the surface was $Δ x^{+} = 100$ and $Δ y_{w a l l}^{+} = 1$ which led to an overall number of cells for the NACA0018 airfoil at the chord-based Reynolds number $R e_{c} = 5.5 \times 10^{5}$ of approximately $5 \times 10^{5}$ and for the NACA0015 airfoil at $R e_{c} = 1.5 \times 10^{6}$ of approximately $7 \times 10^{5} .$ The non-dimensional time step is $δ t {= 10}^{- 4}$ which is approximately two orders of magnitude smaller than the time scale of the shed flow recirculations and unsteady trailing vortex systems.

Similar to the RANS configuration, the LES method consists of a c-shaped grid with extensions in the streamwise, the wall-normal, and the spanwise direction of $30 c,$ $25 c,$ and $0.1 c,$ respectively. The grid resolution on the surface was $Δ x^{+} = 25$ , $Δ y_{w a l l}^{+} = 2$ , and $Δ z^{+} = 15$ which led to an overall number of cells for the NACA0018 airfoil at the chord-based Reynolds number $R e_{c} = 5.5 \times 10^{5}$ of approximately $35 \times 10^{6} .$ The non-dimensional time step was $δ t {= 10}^{- 6}$ to ensure numerical stability of the fully explicit scheme.

The NACA0018 airfoil has been chosen for the investigation since its susceptibility to full stall at high angles of attack is lower compared to that of the NACA0012. Hence, trailing-edge flow separation phenomena can be studied over a wider range of angles of attack.

Results

This section is structured as follows. First, the results of the 2D-AISE are compared to the outcome of the measurements, RANS, and LES findings to discuss more fundamentally the flow separation on stationary airfoils using different solutions of various computational approaches. Then, pitching airfoils are investigated numerically via 2D-AISE and URANS and compared to measurements covering a wide range of flow separation mechanisms. That is, increasingly, intricate flow features are analyzed using the measurements, 2D-AISE, and URANS solutions. Based on the analysis, the numerical technique which depicts best the underlying physics of separated flows is determined. Subsequently, the errors and uncertainties of the measurements and the panel-boundary-layer method are analyzed in detail in section “Error estimation and uncertainty quantification.”

Stationary airfoil

The panel, that is, the AISE method, RANS, and LES-based findings of the NACA0018 flow at $α {= [0, 20]}^{\circ}$ and $R e_{c} = 5.5 \times 10^{5}$ are analyzed and compared to experimental data. To ensure a flow transition that can be juxtaposed to experimental flow conditions, the amplification factor in the $e^{n}$ method is $N = 5$ for the panel method and the RANS simulations. The laminar-turbulent transition in the LES simulations has been triggered at the location given by the RANS results.

The distributions of the lift coefficient are shown in Figure 3 for the 2D-AISE, RANS, and LES solution and the experimental data. A good agreement of the 2D-AISE result with the experimental data is obtained over the entire range of angle of attack. The LES results for $α {= 16}^{\circ}$ and $α {= 20}^{\circ}$ also show a comparable quality. Only the RANS solutions using the SST $k - ω$ turbulence model predict separation already at $α \geq 14^{\circ} .$

Figure 3.

(a) Lift coefficient $c_{l}$ and (b) pressure drag coefficient $c_{d}$ of the NACA0018 at $R e = 5.5 \times 10^{5}$ computed by 2D-AISE, RANS, and LES and compared to experimental data.

In Figure 4, the flow field is shown by the contours of $α {= 20}^{\circ}$ for the 2D-AISE, RANS, and LES solutions. The RANS solution predicts a massively separated flow region, whereas the 2D-AISE and LES distributions show hardly any detached flow.

Figure 4.

Contours and streamlines of the NACA0018 flow at $α {= 20}^{\circ}$ and $R e = 5.5 \times 10^{5}$ computed by (a) 2D-AISE, (b) RANS, and (c) LES to illustrate the mean flow field.

Hence, the numerical solutions and the RANS data show good agreement with the experimental data. Although the flow does not fully separate, the RANS solution predicts massive stall at $α {= 20}^{\circ} .$ The high susceptibility to separation is known for all kinds of standard RANS turbulence models (Celić and Hirschel, 2006). The solution is very sensitive to the grid resolution and the spatial and temporal discretization. The deviation between the mean 2D-AISE results and the experimental findings is approximately 10% to 15% based on the lift coefficients at stalled flow and less than 5% at attached flow.

The pressure coefficient $c_{p}$ and skin-friction coefficient $c_{f}$ distributions shown in Figure 5 emphasize that while the pressure shows hardly any sensitivity to the flow field, the wall shear stress differs considerably. This sensitivity is captured by the 2D-AISE and LES solutions, which are quite similar to each other, whereas the RANS solution differs massively for $α {= 16}^{\circ}$ and $α {= 20}^{\circ} .$

Figure 5.

Pressure coefficient $c_{p}$ and skin-friction coefficient $c_{f}$ distributions of the NACA0018 flow at (a, c) $α {= 16}^{\circ}$ and (b, d) $α {= 20}^{\circ}$ at $R e = 5.5 \times 10^{5},$ 2D-AISE, RANS, LES, and experimental data.

Pitching airfoil

Since the computational costs for an LES of a pitching airfoil at a Reynolds number of $R e_{c} = 1.5 \times 10^{6}$ or $R e_{c} = 5.5 \times 10^{5}$ is about two orders of magnitude higher compared to the stationary airfoil results in section “Stationary airfoil,” only URANS and 2D-AISE solutions are compared to experimental data for this problem. In the following, three flow fields differing in angle of attack, reduced frequency, and deviation from harmonic pitching motion are analyzed in detail. First, a well-established pitching airfoil flow case for a NACA0015 at $R e_{c} = 1.5 \times 10^{6},$ $k = 0.15,$ and $3^{\circ} \leq α \leq 19^{\circ}$ (Galbraith et al., 1992) denoted as configuration A is analyzed.

Using the experimental setup described in section “Experimental setup,” a wide range of reduced frequencies $k \in [0, 1.0],$ amplitudes $| α | \in {[5}^{\circ} {, 25}^{\circ}],$ and Reynolds numbers $R e_{c} \in [2.5 \times 10^{5}, 5.5 \times 10^{5}]$ was measured. Two-flow configurations of these experiments using the NACA0018 airfoil are numerically investigated, that is, $- 24^{\circ} \leq α \leq 24^{\circ},$ $k = 0.1$ and $- 25^{\circ} \leq α \leq 25^{\circ},$ $k = 0.07,$ both at $R e_{c} = 5.5 \times 10^{5}$ . These configurations are denoted as B and C. Note that in the latter case, non-harmonic pitching motion is investigated. The details of the numerical method match those given in section “Stationary airfoil.” The results are compared to experimental data. The flow cases B and C undergo symmetric pitching motions about $α_{0} = 0$ such that attached and detached flow conditions occur. The boundary layer is not tripped. Since the freestream turbulence is about 0.4% which corresponds to the conditions in the water tunnel, the critical amplification factor for free transition is set to $N = 5$ in the 2D-AISE computation. The same transition model was used for the URANS computations for cases B and C.

Harmonic pitching: case A

To validate the 2D-AISE code, the flow field over the NACA0015 airfoil at $R e_{c} = 1.5 \times 10^{6},$ $k = 0.15,$ and $3^{\circ} \leq α \leq 19^{\circ}$ is analyzed. The lift, drag, and moment coefficient hystereses are shown in Figure 6 and are compared to experimental results of Galbraith et al. (1992) and results of Riziotis and Voutsinas (2008), who used a coupled viscous-inviscid panel code. A satisfactory agreement is obtained. At high angles of attack, where the separation location and the detached flow are governed by the non-linearity of the Navier–Stokes equations, the results deviate slightly.

Figure 6.

(a) Lift $c_{l},$ (b) drag $c_{d},$ and (c) moment coefficient $c_{m}$ as a function of the angle of attack $α$ of the NACA0015 flow at $R e_{c} = 1.5 \times 10^{6},$ $k = 0.15,$ and $3^{\circ} \geq α \geq 19^{\circ}$ computed by 2D-AISE and compared to experimental data (Galbraith et al., 1992) and numerical data (Riziotis and Voutsinas, 2008).

The transition from incipient trailing-edge separation to full stall appears rather abrupt, that is, the point of separation quickly propagates upstream. Due to the moderately reduced frequency, the subsequent reattachment process of the flow occurs over the entire upstroke cycle. The details of both flow features are very difficult to reproduce by any URANS simulation (Ekaterinaris and Platzer, 1998; Richter et al., 2011). Nevertheless, it can be stated that the current combination of the double-wake panel method and the integral boundary-layer formulation captures important flow features that are responsible for the intricate dynamical process of a pitching motion at considerable reduced frequencies at deep-stall conditions.

The moment coefficient hystereses show some peaks for the 2D-AISE solution when the flow reattaches. This can be explained by the time step at which the vortices are shed from the airfoil which influences the monotonic translation of the point of separation. Using a smaller time step, the peaks can be reduced significantly.

The comparison with the experimental findings shows that the results of the non-symmetric flow case proved the 2D-AISE code to be capable to reproduce the physically demanding reattaching flow from deep stall.

Harmonic pitching: case B

Next, the flow field over the NACA0018 at $R e_{c} = 5.5 \times 10^{5},$ $- 24^{\circ} \leq α \leq 24^{\circ},$ and a reduced frequency $k = 0.1$ is determined using the URANS and 2D-AISE methods. The numerical data are compared with the measurements, and the effect of a regular pitching motion on the aerodynamic performance is discussed.

The hystereses of the lift, drag, and moment coefficients are juxtaposed in Figure 7. For the 2D-AISE solution, a satisfactory agreement with the experimental data is obtained, whereas the URANS results are not convincing. The mismatch of the drag force coefficient distribution between the numerical and the experimental results during the upstroke motion substantiates the complex flow physics of reattaching flow to be associated with a significant uncertainty. The 2D-AISE computation relies on a simple vortex flow model that is coupled to a steady boundary-layer formulation. From Figure 8, it can be concluded that this coupling seems to be sufficient to predict the lift and drag force coefficients at a quality comparable to the experimental results. Note that the configuration of the turbulence model used in the URANS computation, that is, the SST $k - ω$ model, is not explicitly tuned for time-dependent reattaching turbulent boundary layers. The non-linear nature of the coupling of the motion with the boundary-layer formulation is difficult to predict, that is, the current results of the 2D-AISE formulation are to be interpreted to show the limits of such a design method. Nevertheless, the agreement of the time-averaged results with the experimental results is good from a design perspective. This also shows that the major physical effects, that is, the dynamic flow separation and reattachment are reproduced.

Figure 7.

(a) Lift $c_{l},$ (b) drag $c_{d},$ and (c) moment coefficient $c_{m}$ as a function of the angle of attack $α$ of the NACA0018 flow at $R e_{c} = 5.5 \times 10^{5},$ $k = 0.1,$ and $- 24^{\circ} \geq α \geq 24^{\circ}$ computed by 2D-AISE and URANS and compared to experimental data.

Figure 8.

Pressure coefficient distribution $c_{p}$ of the NACA0018 flow at $R e_{c} = 5.5 \times 10^{5},$ $k = 0.1,$ and $- 24^{\circ} \geq α \geq 24^{\circ}$ computed by 2D-AISE and compared to experimental data and URANS solutions: (a) $α = - 22^{\circ}$ (downstroke), (b) $α = - 17^{\circ}$ (upstroke), (c) $α {= 17}^{\circ}$ (upstroke), and (d) $α {= 22}^{\circ}$ (downstroke).

Similar to the steady airfoil flow, the deviation between the mean 2D-AISE results and the experimental findings is approximately 10% to 15% based on the maximum force coefficients at stalled flow and less than 5% at attached flow.

Non-harmonic pitching: case C

A non-harmonic pitching motion at $- 25^{\circ} \leq α \leq 25^{\circ}$ and a reduced frequency $k = 0.07$ is analyzed which is defined by

α (t) = - atan (\frac{\sin (ω t)}{\cos (ω t) + λ})

(3)

where $λ = 2.4$ is a scaling parameter for the time-dependent angle of attack (Paraschivoiu, 2002). The 2D-AISE and URANS solutions are compared to the experimental results to analyze and evaluate the effect of a non-harmonic pitching motion on the distributions of the aerodynamic coefficients. Note that such a non-harmonic motion characterizes the flow for vertical axis wind-turbine blades.

A satisfactory agreement is obtained for the lift, drag, and moment coefficient hystereses shown in Figure 9. The rate of change of the angle of attack is high where the lift coefficient is low and vice versa. That is, the flow conditions alter rapidly from deep stall to reattached flow conditions. Due to the double-wake formulation which determines the dynamic vortex formation near the point of separation, the major contributor to the complex flow system, that is, the dynamically shed vortex system, is described sufficiently accurate.

Figure 9.

(a) Lift $c_{l},$ (b) drag $c_{d},$ and (c) moment coefficient $c_{m}$ as a function of the angle of attack $α$ of the NACA0018 flow at $R e_{c} = 5.5 \times 10^{5},$ $k = 0.07,$ and $- 24^{\circ} \geq α \geq 24^{\circ}$ computed by 2D-AISE and compared to experimental data and URANS solutions.

In Figure 10, the pressure coefficient distributions at four angles of attack are shown. The good agreement between the results of the 2D-AISE solutions and the experimental data indicates that the important flow features such as the dynamic movement of the separation point and the induced velocities by the complex vortex system due to the double-wake formulation are well reproduced. The URANS computation results in a reattachment of the boundary layer in the upstroke movement that is significantly retarded compared to the experimental and 2D-AISE findings.

Figure 10.

Pressure coefficient $c_{p}$ distribution of the NACA0018 flow at $R e_{c} = 5.5 \times 10^{5},$ $k = 0.07,$ and $- 24^{\circ} \geq α \geq 24^{\circ}$ computed by 2D-AISE and compared to experimental data and URANS solutions: (a) $α = - 22^{\circ}$ (downstroke), (b) $α = - 17^{\circ}$ (upstroke), (c) $α {= 17}^{\circ}$ (upstroke), and (d) $α {= 22}^{\circ}$ (downstroke).

The non-harmonic pitching airfoil at high angles of attack is characterized by a global non-constant temporal gradient of the separation point and the resulting wake structure which interacts non-linearly with the boundary-layer formulation. This flow physics is also well reproduced by the 2D-AISE code.

Although this flow structure is even more complicated than that discussed for case B, the maximum deviation of the time-averaged results between the measurements and the computed values is approximately 10% to 15% based on the maximum values of the force coefficients. Note that the agreement of the results for pitching airfoil configurations obtained by the 2D-AISE code with the experimental results is similar to that of the URANS results using standard turbulence models (Celić and Hirschel, 2006; Ekaterinaris and Platzer, 1998).

Before a conclusive summary can be drawn, it has to be emphasized that the flow cases are characterized by high sensitivity of the flow field to various perturbations. On one hand, considering the measurements, the uncertainty of the experimental results at high angles of attack must be taken into account. On the other hand, the numerical simulation also includes errors and uncertainties due to methodological and modeling issues. This issue is further addressed in the following section.

Error estimation and uncertainty quantification

Following the definition given by Iaccarino (2009), the errors, aleatoric (stochastic), and epistemic (deterministic) uncertainties are discussed for the 2D-AISE code and experiments.

To quantify the uncertainty that is connected to the expectation value and the standard deviation of the force coefficients obtained from numerical simulations or measurements, the potential sources are discussed in the following sections. First, the numerical and experimental errors are described in detail followed by a description of the sources of uncertainties in the numerical procedure and the measurement campaign in the water tunnel. Note that the evaluation of errors and uncertainties itself is prone to inaccuracies due to their simplified determination. A significant range of reduced frequencies was investigated during the error estimation process since the impact on the flow physics and, thus, the error is rather pronounced. A complete Monte-Carlo simulation of the parameter space would have required 10⁴ 2D-AISE computations which is an extreme effort. Instead, the impact of the isolated errors and uncertainties is estimated using a first-order sensitivity analysis. Checking the results by a random parameter recombination showed that the estimated upper and lower bounds of the computed overall error and uncertainty were not exceeded.

The normally distributed input parameter space for the numerical simulations to determine the error and uncertainties included approximately 100 parameters. The input parameters were uniformly distributed within the defined limits and the simulations were executed for a range of reduced frequencies $k \in [0.1, 1.0]$ and amplitudes $| α | \in {[5}^{\circ} {, 25}^{\circ}]$ . In the following, the sources of errors and uncertainties for the numerical analysis and the experiments are briefly discussed and their influence on the force coefficients is described.

Numerical errors

Numerical errors are treated separately from numerical uncertainties since errors can be eliminated and/or quantified by bug fixing and asymptotic parametric procedures (Iaccarino, 2009). In Table 3, the potential errors in the lift coefficients that occur in the 2D-AISE code for a pitching airfoil configuration are listed.

Table 3.

Approximate maximum error/uncertainty of the lift coefficient $c_{l}$ (%) for the deep-stall case determined by 2D-AISE and experiments.

Numerical	Error (%)	Uncertainty (epistemic, model based, %)
Panel resolution and time step	3%
Higher order temporal discretization	⩽1%
Vortex merging criterion	3%
Convergence of viscous-inviscid coupling	⩽1%
Convergence criterion of boundary layer	3%
Boundary-layer skin-friction model		18%
Vortex-core model		⩽1%
Experimental	Error (%)	Uncertainty (%)
Correction of $α$ and $c_{l},$ $c_{d},$ and $c_{m}$	10%
Convergence of averaging to mean flow	2%
Pressure probes	3%
Inconsistent inflow		3% (aleatoric)
Prescribed angle of attack via linear motor		2% (epistemic)
Water level		1%
Flow field parameters		⩽0.01% (aleatoric)

Number of panels and time steps

The spatial and temporal resolution, that is, the number of panels and time steps per pitching motion has to be chosen sufficiently high to keep the error low. A spatial discretization of 150 panels and 50 time steps per pitching motion resulted in a maximum error of about 3% in $c_{l}$ coefficients when simulating an inviscid 2D-AISE configuration using the intervals of $k$ and $α$ given above. The reference value to compute the error is provided by a simulation with 600 panels and 150 time steps. A higher number of panels and time steps do not alter the resulting force coefficient anymore. Furthermore, the results do not strongly depend on the reduced frequency since the number of time steps is set per pitching motion. In viscous flows, that is, when the inviscid-viscous coupling of the panel and the boundary-layer solution is applied, the choice of the time step has also an impact on the force coefficients if flow separation occurs. By comparing pitching motion with flow separation at different spatio-temporal discretization in terms of number of panels and time steps, the error remains similar to the inviscid one when the result of 150 panels and 50 time steps is again compared to a configuration with 600 panels and 150 time steps. That is, the number of panels ensures the smoothness of the $c_{p}$ distribution and the impact of the trailing vortices are rather affected by the number of time steps.

Higher order temporal discretization

The influence of a higher order temporal scheme, that is, second-order scheme, on the integral error compared to the first-order Euler time-stepping scheme is only marginal, that is, less than 1% in case of inviscid and viscous simulations on the force coefficients.

Merging curtain

The merging curtain is used to effectively reduce the number of vortex cores in the simulation and hence to reduce the computational effort. Simulations showed that a distance between the trailing edge and the curtain in the range of 6–1000 chord lengths causes a maximum error in $c_{l}$ of about 3%. Here, the simulation based on a distance of 1000 chord lengths is the reference result. The error occurs due to the phase shift caused by the merging vortices near the airfoil. However, the error quickly drops with increasing distance.

Viscous-inviscid coupling

To couple the viscous and the inviscid formulation, a semi-inverse coupling algorithm (Zanon et al., 2013) mentioned in section “Panel method” is applied. The displacement thickness is iteratively adjusted to ensure a converged solution between the boundary-layer code and the panel method. The convergence criterion is set such that a deviation of 0.1% between the integral lift force of the current and the last iteration step is allowed. The iteration limit is chosen such that the error can be neglected.

Convergence criterion of boundary layer

The boundary-layer code operates with an internal convergence criterion and a maximum number of iteration steps. The latter is modified to investigate the influence on a separated flow case by applying 10–100 internal iteration steps. It was observed that the error is small, that is, less than 1%, if the flow is globally attached even when only 10 iterations are applied. However, it may significantly increase up to 10% at such low iteration numbers if strong stall occurs. Using 50 steps and more, the error is approximately 3%. Note that this quantitative analysis is only valid upstream of the separation.

Numerical uncertainty

Boundary-layer skin-friction model

The skin-friction model (Drela, 1989) used in the boundary-layer implementation is based on semi-empirical correlations to meet the experimental values of numerous measurement campaigns. To investigate the impact on the integral solution of a pitching airfoil, the parameters of the model are modified such that the flow separation behavior is significantly altered. Therefore, the empirical skin-friction criterion in the boundary-layer formulation is pre-multiplied with a coefficient varying between [0.8, 1.2]; see equations (22), (23), and (27) in the study of Drela (1989). This leads to a modified skin-friction criterion that still provides valid flat-plate boundary-layer solutions, that is, momentum thickness over Reynolds number varies between ±10% compared to direct numerical simulation results (Schlatter and Örlü, 2010). Furthermore, the separation criterion introduced in section “Panel method” is varied between [0, 0.1]. These epistemic uncertainties have been computed using the parameters given in section “Numerical errors” and grow with increasing reduced frequencies and pitching time steps. The force coefficients in high-amplitude pitching cases include a peak uncertainty up to 20%.

Vortex core model

The vortices shed from the separation and trailing-edge panel are described via Lamb–Oseen vortex model with constant viscous-core radius. It has been observed that due to the missing direct interaction between the airfoil and the shedding vortex, the impact of reasonable changes in the viscous-core radii, that is, variations of up to 20%, is negligible.

Experimental error

Correction of lift, drag, moment, and angle of attack

Due to the geometric ratios of the investigated wing model and the operating water tunnel, a set of semi-empirical correction expressions for the integral force coefficients can be applied to estimate the error of the expectation value. The principally associated effects are water tunnel blockage and streamline curvature. The correction formulas were applied and the results have been treated as potential error sources, which are a function of local force coefficients and angles of attack (Barlow et al., 1999).

Convergence to mean force coefficients

A sufficient number of pitching motions is measured to analyze the convergence of the mean integral force coefficients. To estimate the error, a statistically meaningful number of pitching motions has to be averaged and the standard deviation determines the error due to the measurement procedure and facility. This source of error can be several percent of the expectation value.

Pressure probes

The integral force coefficients are computed from the pressure data obtained by pressure probes that are distributed over the airfoil. That is, there are three sources of uncertainty of the expectation value of the integral forces. First, the probe itself is only accurate up to 170 mbar which leads to an error of about 20%. In the measurement campaign, the relative error of the expectation value was estimated to be much smaller, that is, approximately 3%. Second, the discretization of the pressure probes affects the spline reconstruction which is used to obtain the integrated force coefficients. It was found that the number of probes together with the spline reconstruction led to an error which was less than 1%.

Experimental uncertainty

Temporal and spatial inflow variation

The analysis of the flow field upstream of the airfoil based on PIV data showed that a low-frequency fluctuation is present in the water tunnel at 0.05 Hz and an amplitude of about 0.02 U_∞. Furthermore, a spatial inconsistency occurs due to the water tunnel setup. The time-averaged cross-sectional velocity distribution varies up to 20% where uniform velocities were expected. That is, two distinct effects in the freestream configuration of the numerical simulation to mimic their influence are considered

u (z, t) = \bar{u} (z) + 0.02 \bar{u} (z) \sin (ω_{w t} t)

(4)

where $ω_{w t}$ is $2 π f_{w t}$ and $\bar{u} (z)$ denotes the mean velocity as a function of the channel height z. Although being genuinely aleatoric in the process of estimating the impact due the complexity of the actual flow phenomena, this uncertainty is treated as epistemic, that is, only the influence of the major velocity component at one frequency is modeled. The mean velocity $\bar{u (z)}$ is determined by the following generic spatial variation of the reference velocity and imposed on the inflow of 2D-AISE, that is

\bar{u} (z) = \bar{u_{r e f}} (z_{r e f}) + 0.8 \bar{u_{r e f}} (z)

(5)

where the factor 0.8 is used to fit the mean flow velocity distribution upstream of the airfoil measured by PIV. It is observed that the global impact on both attached and separated flows over a wide range of reduced frequencies and amplitudes are about 3% based on the maximum force coefficients.

Prescribed angle of attack

The epistemic uncertainty of the prescribed angle of attack of the pitching motion is taken into account by the upper limit of the error due to mechanical and control issues which is $Δ α = 0.01 α_{r e f} .$ That is, the angle of attack $α$ or the ratio of $u_{\infty} / v_{\infty}$ U_∞ / V_∞ in the numerical simulation is modified by

α = α_{r e f} + \sum_{k = 1}^{p} a_{k} \sin (ω_{k} t)

(6)

where $a_{k}$ is the normally distributed random amplitude which leads to $0.01 α_{r e f}$ when summed up over $p .$ Similar to the uncertainty connected to the inflow condition, the 2D-AISE simulations revealed that a global uncertainty of about 2% arises from a modified angle of attack computed by equation (6).

Water level

Due to the dynamics of the pitching motion in the water tunnel, the water level and the static pressure at the pressure probes may vary in time. In the numerical simulation, this is taken into account by modifying the potential as a function of the streamwise chord position. This is done by virtually adding static pressure according to a modified water level over a maximum chord length of 0.05 m. The impact on the integral force coefficients was about 1%.

Flow field parameters

Due to temperature fluctuations of about 0.2 K, the viscosity and subsequently the Reynolds number vary by approximately 0.4%. This is considered sufficiently small such that it can be neglected as input parameter variation in the numerical simulation.

Resulting numerical and experimental force coefficients due to errors and uncertainties

The combined value of the discussed errors and uncertainties is computed by executing a first-order sensitivity analysis. Assuming that the expectation value of the force coefficient is obtained using the temporal average, and the variation is computed by

Var (C) = \sum_{i} {(\frac{\partial C}{\partial P_{i}})}^{2} Var (P_{i})

(7)

where $C$ is the force coefficient and $i$ denotes the source of error or uncertainty with the parameter space $P_{i} .$ The $Var (P_{i})$ uses a standard deviation about the mean defining a 5–95 percentile. Using the expected value of the local force coefficient $μ = \bar{C}$ and $σ = \sqrt{Var},$ the relative coefficient of variation at each angle of attack is computed by $ξ_{r e l} = μ / σ .$

Figures 11 –13 show the error and uncertainty values superimposed on the time-averaged results of a steady flow and a harmonic as well as a non-harmonic pitching airfoil, that is, the cases B and C discussed in sections “Stationary airfoil” and “Pitching airfoil” are considered. It is observed that the errors and uncertainties of the experiment and the 2D-AISE simulation are of the same order of magnitude almost over the entire range of angle of attack. That is, although the differences of the origin and the related impact of the error sources in the measurements and computations are significant, the overall error and uncertainty of the mean result are comparable. Simulations and experiments of the flow cases show a maximum relative variation $ξ_{r e l}$ of the local force coefficient value of about 15%–20% at high angle of attack.

Figure 11.

(a) Normal force coefficient $c_{n}$ and (b) drag coefficient $c_{d}$ of the stationary airfoil case described in section “Stationary airfoil”; 2D-AISE and experimental data including errors and uncertainties given in Table 3.

Figure 12.

(a) Lift $c_{l},$ (b) drag $c_{d},$ and (c) moment coefficient $c_{m}$ of the harmonic pitching motion (case B); 2D-AISE and experimental data including errors and uncertainties given in Table 3.

Figure 13.

(a) Lift $c_{l},$ (b) drag $c_{d},$ and (c) moment coefficient $c_{m}$ of the non-harmonic pitching motion (case C); 2D-AISE and experimental data including errors and uncertainties given in Table 3.

Note that the numerical uncertainties in the investigated flow cases are locally even smaller than those of the experimental results. Due to the quantified errors and uncertainties, the numerical and the experimental results are considered, that is, to match in between the upper and lower bounds, over the entire range of angles of attack.

Figure 14 shows the maximum relative error and uncertainty $ξ_{r e l}$ based on the lift coefficient for the experimental and numerical results. By increasing the maximum angle of attack and the reduced frequency, the maximum value of $ξ_{r e l}$ increases. Note that the uncertainties assessed by Abdallah et al. (2015) for static airfoils are similar to the current maximum values given in Figure 14.

Figure 14.

Maximum relative error and uncertainty $ξ_{r e l}$ for reduced frequency $k$ and maximum angle of attack $α_{m a x}$ in the (a) measurements and (b) 2D-AISE results.

The thorough quantification of errors and uncertainties of numerical simulations or measurements is standard in neither industry nor science. However, this step is necessary to enhance the confidence and to efficiently reduce resources in the design process that depends on mean analyses of the flow fields.

Computational resources and confidence estimation

The computational resources of the 2D-AISE and the URANS analyses are summarized in Table 4. Based on the arguments given in section “Error estimation and uncertainty quantification,” the maximum normalized confidence of the force coefficients is $I = 1$ . From the results given in section “Pitching airfoil” and the reviews in the studies of Ekaterinaris and Platzer (1998) and Richter et al. (2011), the URANS simulations show a similar confidence value compared to the 2D-AISE simulations presented in section “Resulting numerical and experimental force coefficients due to errors and uncertainties,” that is, $I_{A I S E} = 1 - ξ_{r e l} |_{m a x} \approx I_{R A N S}$ . However, the computational time using URANS is still about $O (10^{3})$ core hours on small clusters with less than 500 cores, whereas the 2D-AISE code requires a standard desktop PC core hours in the range of $O (10^{0})$ to $O (10^{1}) .$ Note that unlike simulation hours, core hours directly quantify computational costs. That is, in case of a pitching airfoil problem, the 2D-AISE code is about two orders of magnitudes more efficient than a URANS simulation at similar accuracy.

Table 4.

Computational resources and confidence of the results for different numerical configurations.

Simulation type	Confidence $I$	Core hours	Simulation hours
URANS	$~ ~ 0.80$	$O (10^{3})$	$O (10^{2})$ to $O (10^{3})$
2D-AISE	$~ ~ 0.80$	$O (10^{1})$	$O (10^{0})$ to $O (10^{1})$

URANS: unsteady Reynolds-averaged Navier–Stokes similation; AISE: aeroelastic integrated simulation environment.

Conclusion

In this study, a two-dimensional numerical design tool for pitching airfoils is assessed. Harmonically and non-harmonically pitching airfoil configurations in a Reynolds number range representative for wind turbines are numerically and experimentally investigated. For a stationary airfoil configuration, an LES simulation at high angles of attack is considered as a reference solution of this problem and for the pitching airfoil problem, the aerodynamic data are compared to the solutions of a modified double-wake vortex-panel method with boundary-layer formulation and to URANS simulations.

The results of the 2D-AISE panel-boundary-layer code show satisfactory agreement with experimental and LES data for stationary airfoils at high angles of attack. Furthermore, the same quality as measurements and 2D-AISE comparisons is achieved for pitching airfoil configurations. The RANS/URANS results, however, are less convincing since stall onset and reattachment are poorly predicted. Note that this result is valid only for the current RANS/URANS setup and cannot be generalized to the RANS concept.

The rigorous analysis of the error and uncertainty sources in the numerical methods and the experiments proves that for design applications, the current low-cost double-wake model is well suited in the first stages of the design process. That is, the resulting errors and uncertainties of the first-order sensitivity analysis of the double-wake panel method 2D-AISE are in the same range as the measurements and the clearly more expensive URANS method. The unsteady panel-boundary-layer method provides a convincing confidence at very low computational costs in the investigated range of reduced frequencies and angles of attack for standard airfoils such that it is an appealing approach in the preliminary design stage of, for example, wind turbines.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Abdallah

Natarajan

Sørensen

(2015) Impact of uncertainty in airfoil characteristics on wind turbine extreme loads. Renewable Energy 75: 283–300.

Barlow

Rae

Pope

(1999) Low-Speed Wind Tunnel Testing. Toronto, ON, Canada: John Wiley & Sons.

Basu

Hancock

(1978) The unsteady motion of a two-dimensional aerofoil in incompressible inviscid flow. Journal of Fluid Mechanics 87(1): 159–178.

Celić

Hirschel

(2006) Comparison of eddy-viscosity turbulence models in flows with adverse pressure gradient. AIAA Journal 44(10): 2153–2169.

Corke

Thomas

(2015) Dynamic stall in pitching airfoils: Aerodynamic damping and compressibility effects. Annual Review of Fluid Mechanics 47: 479–505.

Danao

Edwards

Eboibi

, et al. (2014) A numerical investigation into the influence of unsteady wind on the performance and aerodynamics of a vertical axis wind turbine. Applied Energy 116: 111–124.

Deck

Laraufie

(2013) Numerical investigation of the flow dynamics past a three-element aerofoil. Journal of Fluid Mechanics 732: 401–444.

Drela

(1989) XFOIL: An analysis and design system for low Reynolds number airfoils. In: Mueller

(ed.) Low Reynolds Number Aerodynamics. London: Springer, pp. 1–12.

Drela

Giles

(1987) Viscous-inviscid analysis of transonic and low Reynolds number airfoils. AIAA Journal 25(10): 1347–1355.

10.

Ekaterinaris

Platzer

(1998) Computational prediction of airfoil dynamic stall. Progress in Aerospace Sciences 33(11–12): 759–846.

11.

Ferreira

Madsen

Barone

, et al. (2014) Comparison of aerodynamic models for vertical axis wind turbines. Journal of Physics: Conference Series 524: 012125.

12.

Galbraith

Gracey

Leith

(1992) Summary of Pressure Data for Thirteen Aerofoils on the University of Glasgow’s Aerofoil Database (Aero Report). Glasgow: University of Glasgow.

13.

Hedenström

Johansson

Wolf

, et al. (2007) Bat flight generates complex aerodynamic tracks. Science 316(5826): 894–897.

14.

Iaccarino

(2009) Quantification of uncertainty in flow simulations using probabilistic methods (technical report). Brussels: Von Karman Institute of Fluid Mechanics.

15.

Kármán

(1921) Über laminare und turbulente Reibung. ZAMM: Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 1(4): 233–252.

16.

Katz

Plotkin

(2001) Low-Speed Aerodynamics, Vol. 13. Cambridge: Cambridge University Press.

17.

Koumoutsakos

Leonard

(1995) High-resolution simulations of the flow around an impulsively started cylinder using vortex methods. Journal of Fluid Mechanics 296: 1–38.

18.

Lee

(2013) An experimental study of stall delay on the blade of a horizontal-axis wind turbine using tomographic particle image velocimetry. Journal of Wind Engineering and Industrial Aerodynamics 123: 56–68.

19.

Lee

Gerontakos

(2004) Investigation of flow over an oscillating airfoil. Journal of Fluid Mechanics 512: 313–341.

20.

Leishman

(2006) Principles of Helicopter Aerodynamics (2nd edn). Cambridge: Cambridge University Press.

21.

Leishman

Beddoes

(1986) A generalised model for airfoil unsteady aerodynamic behaviour and dynamic stall using the indicial method. In: Proceedings of the 42nd Annual Forum of the American Helicopter Society, Washington, DC, 2–4 June, pp. 243–265. Alexandria, VA: American Helicopter Soc.

22.

Leishman

Beddoes

(1989) A semi-empirical model for dynamic stall. Journal of the American Helicopter Society 34(3): 3–17.

23.

McCroskey

(1982) Unsteady airfoils. Annual Review of Fluid Mechanics 14(1): 285–311.

24.

Meinke

Schröder

Krause

, et al. (2002) A comparison of second- and sixth-order methods for large-eddy simulations. Computers & Fluids 31: 695–718.

25.

Menter

(1994) Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal 32(8): 1598–1605.

26.

Paraschivoiu

(2002) Wind Turbine Design: With Emphasis on Darrieus Concept. Montreal, QC, Canada: Presses Inter Polytechnique.

27.

Richter

Le Pape

Knopp

, et al. (2011) Improved two-dimensional dynamic stall prediction with structured and hybrid numerical methods. Journal of the American Helicopter Society 56(4): 1–12.

28.

Riziotis

Voutsinas

(2008) Dynamic stall modelling on airfoils based on strong viscous-inviscid interaction coupling. International Journal for Numerical Methods in Fluids 56: 185–208.

29.

Roidl

Meinke

Schröder

(2013) A reformulated synthetic turbulence generation method for a zonal RANS-LES method and its application to zero-pressure gradient boundary layers. International Journal of Heat and Fluid Flow 44: 28–40.

30.

Roidl

Meinke

Schröder

(2014) Boundary layers affected by different pressure gradients investigated computationally by a zonal RANS-LES method. International Journal of Heat and Fluid Flow 45: 1–13.

31.

Schlatter

Örlü

(2010) Assessment of direct numerical simulation data of turbulent boundary layers. Journal of Fluid Mechanics 659: 116–126.

32.

Voutsinas

(2006) Vortex methods in aeronautics: How to make things work. International Journal of Computational Fluid Dynamics 20(1): 3–18.

33.

Voutsinas

Riziotis

(1999) A viscous-inviscid interaction model for dynamic stall simulations on airfoils. In: 37th aerospace sciences meeting and exhibit, Reno, NV, 11–14 January.

34.

Whitfield

(1978) Analytical description of the complete turbulent boundary layer velocity profile. AIAA Paper 78: 1158.

35.

Zanon

(2011) A vortex panel method for VAWT in dynamic stall. PhD thesis, Universita Degli Studi Di Udine, Udine.

36.

Zanon

Giannattasio

Simão Ferreira

(2013) A vortex panel model for the simulation of the wake flow past a vertical axis wind turbine in dynamic stall. Wind Energy 16(5): 661–680.