Sensitivity analysis of airfoil aerodynamics during pitching motion at a Reynolds number of 1 . 35 × 105

Numerically reproducing wind tunnel experimental tests of flow behavior around a pitching airfoil is a challenge, especially under the occurrence of dynamic stall at relatively high angles of attack. This not only requires the application of advanced turbulence models, but also asks for examining the influence of the most relevant model parameters in detail. This contribution presents the results of an extensive computational study of the unsteady flow around a pitching NACA0012 airfoil at a Reynolds number of 1.35 105, as obtained by first performing 2D Unsteady Reynolds-Averaged Navier-Stokes (URANS) simulations whereby the flow characteristics are simulated by the Transition SST turbulence model. The influence of a large number of parameters on the numerical results is investigated, namely the blockage ratio, computational grid resolution and yþ value, time step size, inlet freestream turbulence properties, airfoil trailing edge detailing and turbulence model. Integral aerodynamic force coefficients and detailed flow patterns are analyzed and compared with measurements from wind tunnel experiments presented in the literature. For the best-performing parameters, an adequate agreement with the experimental tests is obtained for the upstroke phase, while for the downstroke phase some differences appear. These differences are investigated in more detail by considering the results from a 2.5D Large Eddy Simulation (LES), which provides deep and complementary insights into the flow behavior during dynamic stall at the selected Reynolds number.


Introduction
Airfoils undergoing pitching motion experience a stall behavior that is distinctly different from the behavior at fixed angles of attack (AOA). This phenomenon is therefore referred to as "dynamic stall", for which the stall angle and integral forces and moments (lift coefficient C L , drag coefficient C D and momentum coefficient C M ) can all be far in excess of their static stall counterparts (McCroskey, 1981). The forces and moments also exhibit a hysteresis effect as a function of the actual AOA, which does not appear for static airfoils. The dynamic stall phenomenon is particularly relevant for applications involving rotating machinery, such as turbines, compressors and helicopter rotors. In the case of airfoils subjected to high-amplitude pitching motion, the obstacle body is not streamlined anymore as the AOA increases. Consequently, for moderate-thickness airfoils, flow at the leading edge starts to roll up thereby forming an energy-containing vortex (leading edge vortex, LEV), which provides an additional suction across the upper airfoil surface (Leishman, 1990). The added suction increases lift and postpones stall to a higher AOA. This LEV grows and affects the boundary layer thickness as it propagates downstream (McCroskey, 1981;Lee and Gerontakos, 2004). Next, a second vortex, which counter-rotates with respect to the LEV, develops at the suction side close to the leading edge (McCroskey, 1981;Gharali and Johnson, 2013). This secondary vortex interacts with the LEV and both are shed into the wake. In the meantime a trailing edge vortex (TEV) forms. While the LEV is advected downstream, it interacts with this TEV, which in turn is also shed into the wake. During the downstroke phase of the pitching motion this process is repeatedly followed by the generation of another LEV and TEV and their subsequent shedding into the wake, until flow reattaches at a lower AOA. Accordingly, the wake flow during the downstroke is characterized by an alternating vortex pattern originating from shed LEVs and TEVs. The integral forces and moments experience substantial variations and large hysteresis, resulting in a significant increase in stresses that make the blades more susceptible to fatigue damage. The main physical features involved in dynamic stall are discussed by McCroskey (1981) and more recently by Lee and Gerontakos (2004).
Dynamic stall occurs regularly in vertical axis wind turbines (VAWTs), due to the wide range of AOAs experienced by turbine blades during operation, especially at low tip speed ratios (Rezaeiha et al., 2018). This has a clear effect on the fatigue resistance, noise emission, and output efficiency of VAWTs. In urban applications an H-type Darrieus micro-VAWT normally operates at a chord-based Reynolds number (Re) on the order of 10 5 . The characteristics of the boundary layer on airfoils, such as the transition onset location, cause the aerodynamic performance of the VAWT as a whole to be very sensitive to a wide range of parameters, especially at Reynolds numbers in between 10 4 < Re < 10 6 (Kim and Change, 2014;Lissaman, 1983). These parameters are related to pitching kinematics (mean pitching angle α m , pitching amplitude α p and reduced frequency κ), blade shape (airfoil profile and chord length), and operating conditions (Reynolds number, Mach number and turbulence characteristics of the freestream flow). It is widely recognized that an accurate and reliable simulation of these parameters by computational fluid dynamics (CFD) analyses is challenging, as the CFD results may be very susceptible to the computational settings used. For this purpose, in the present communication an elaborate numerical study is performed for a range of parameters governing the aerodynamics of deep dynamic stall of an airfoil subjected to pitching motion. For achieving high-fidelity computational results, the numerical results are systematically validated with experimental data presented in the literature.
In the last decades extensive experimental work has been performed on the physics of dynamic stall and the aerodynamic performance of pitching airfoils (McCroskey et al., 1976(McCroskey et al., , 1982Carr et al., 1977;Lee and Basu, 1998;Lee and Gerontakos, 2004;Lind and Jones, 2016). Of the existing dynamic stall predictive methods, which were initiated in the helicopter industry and recently shifted for use in wind turbine studies (Leishman, 2002;Larsen et al., 2007), semi-empirical methods are typically used to predict the essential physics of airfoils under dynamic stall conditions, by correlating force and moment data obtained from wind tunnel tests (McCroskey, 1981;Paraschivoiu and Allet, 1988;Holierhoek et al., 2013). Although semi-empirical methods generally provide acceptable approximations for the unsteady aerodynamics, calculation of the aerodynamic coefficients for a specific airfoil shape is less accurate, especially for airfoils in the deep dynamic stall regime (Holierhoek et al., 2013). With the vast progress in computational methods and the tremendous increase in computational power, a large part of recent research on pitching airfoils has focused on numerical work (Osswald et al., 1991;Martinat et al., 2008;Wang et al., 2012;Gharali and Johnson, 2013;Kim and Xie, 2016;Rahromostaqim et al., 2016 (LES)), the turbulence model adopted, the dimension of the computational geometry and domain (2D, 3D) and the actual parameter for which a sensitivity analysis is performed. Although a lot of effort has been put into predicting dynamic stall, the large sensitivity of the flow characteristics to the airfoil geometry and conditions mentioned above has so far seriously hampered validation exercises, especially for airfoils pitching into deep dynamic stall with a large oscillation rate and operating at a low Reynolds number (Tank et al., 2017). As an example of the above, CFD model validations of the wind tunnel test data obtained by Lee and Gerontakos (2004) by different authors resulted into a large spread the numerical results (Wang et al., 2010(Wang et al., , 2012Gharali and Johnson, 2013;Kim and Xie, 2016).
In CFD simulations the trailing edge profile shape can have a significant impact on the flow field across airfoils, thereby affecting their overall aerodynamic performance (Standish and Van Dam, 2003;Almohammadi et al., 2015). However, the studies were performed with large blunt edges (e.g., a thickness to chord ratio at the trailing edge of 0.1 in Almohammadi et al., 2015). Given the importance of TEV in dynamic stall, one can ask how to geometrically reproduce the trailing edge in the pre-processing step of a CFD simulation: this is typically done by either extending or truncating the trailing edge (see Fig. 6). The uncertainty introduced by the chosen method on the dynamic stall characteristics (usually with a small thickness-to-chord-ratio, such as the 0.0025 adopted in this study) is not known. Tank et al. (2017) simplified the investigated NACA0012 geometry by truncating the trailing edge at a thickness of 0.01c in order to simplify grid generation for the investigation of the integral aerodynamic force coefficients at a low Reynolds number. In contrast, Standish and Van Dam (2003) argued that it is reasonable to extend the chord by 0.002c, and did so by fitting the trailing edge shape with a quadratic polynomial in a baseline simulation. However, the above simulations are based on static airfoils, and were aiming at quantifying the uncertainty arising from these simplifications rather than providing a comparison of both methods.
The turbulence modeling approach adopted generally plays a crucial role in the accurate simulation of fluid flow. Gharali and Johnson (2013) performed a 2D study around a NACA 0012 pitching airfoil using the Shear-stress Transport (SST) k-ω model, and the Transition SST model was employed in Wang et al. (2012) for the same calculation. The simulation results showed that the lift coefficient benefited from the Transition SST model's more accurate transition modeling approach. Zanotti et al. (2014) conducted a comparison of 2D and 3D URANS simulations. They concluded that the 3D results were in better agreement with measurements, leading to an improved accuracy of the lift coefficient curve slope prediction; the 3D simulations also properly captured the movement of vortices along the airfoil upper surface during the downstroke. A study by Martinat et al. (2008) further showed that 3D models were better able to predict the decrease in the coherent structure amplitude during the downstroke than their 2D counterparts.
Although it can be expected that the URANS approach will continue to play an important role in many fields of application, there is general consensus that scale resolving approaches such as DES, LES and Wallmodeled LES (WMLES) have the capability to provide improved flow predictions (Menter, 2012;Blocken, 2014). Wang et al. (2012) evaluated the behavior of the lift coefficient of a pitching airfoil by using the 2D URANS approach with the Transition SST turbulence model and the SST k-ω based Delayed DES (DDES) model, respectively. The DDES failed to show advantages compared to the 2D URANS results in terms of the prediction of stall angle and the peak value of the lift coefficient. Recently, Kim and Xie (2016) studied the effect of freestream turbulence on the aerodynamic characteristics of a pitching airfoil using the LES modeling approach. Three reduced frequencies of 0.025, 0.05 and 0.1 were used. They found that results for the cases with small κ corresponded well with experimental data in terms of dynamic forces, while notable deviations were observed for LES results for the case κ ¼ 0.1. As remarked in Kim (2013), at the chosen Reynolds number (Re ¼ 1.35 Â 10 5 ) the flow characteristics are exceedingly sensitive to the problem details, as a result of which experimental results could vary significantly between different flow facilities and experimental setups. Another reason relates to the observation that the leading edge vortex bursting in the case of κ ¼ 0.1 is much more violent than in cases with lower κ, which most likely requires a higher computational resolution to properly resolve this aspect. In this research we specifically focus on the challenging case with κ ¼ 0.1, for which the numerical results can be compared in detail with the wind tunnel test data obtained by Lee and Gerontakos (2004). The goal of the present paper is to perform a high-fidelity validation study with quantification of the discretization error, input uncertainty and physical modeling uncertainty (Casey and Wintergerste, 2000;Oberkampf and Trucano, 2002;Franke et al., 2007;Tominaga et al., 2008;Blocken, 2015). Computational and physical parameters investigated in this paper include blockage ratio, computational grid resolution, y þ value, temporal resolution, freestream turbulence properties, airfoil trailing edge detailing and turbulence model.
The paper is organized as follows. Section 2 starts with outlining the experimental work of Lee and Gerontakos (2004) used for validating the numerical simulations, including the experimental setup, incoming flow properties and measurement methods. Subsequently, as a preliminary validation study, the experimental results are numerically reproduced with a 2D URANS simulation; the impact of computational parameters on force coefficients is detailed through a sensitivity analysis and quantified by four performance metrics. Section 3 presents a 2.5D LES model, for which the computational results are compared to the experimental results and the results computed by the 2D URANS model. The aerodynamic force coefficients and flow visualizations, such as the LEV formation and the vortex shedding characteristics, as computed from the 2D URANS and 2.5D LES models, are compared at various time instances. Final conclusions are given in Section 4.

Quality assessment of 2D URANS simulations
2.1. Review of wind tunnel experiments Lee and Gerontakos (2004) conducted wind tunnel measurements to investigate the unsteady boundary layer and stall events on a pitching NACA 0012 airfoil with a 0.15 m chord length (c) and a 2.5c span mounted horizontally in the 0.9 m (height) Â 1.2 m (width) Â 2.7 m (length) suction-type wind tunnel at McGill University, Canada. In order to obtain a 2-dimensional flow across the pitching airfoil and to eliminate the effect of highly unsteady swirls at the airfoil edges, two splitter plates with a diameter of 2c were mounted at the ends. The blockage ratio of the wind tunnel experiment, which is defined as the ratio of the projected frontal area of tested obstacles to the cross section of the wind tunnel, was 2.2% when the airfoil was at its maximum AOA (α ¼ 25 ). Pre-stall, during stall and post-stall flow physics and mechanisms responsible for stall, were measured and visualized for both a sinusoidal pitching motion  with an angle of attack αðtÞ ¼ α m þ α p sinðΩtÞ and for a static airfoil. Five sets of pitching modes were reported by prescribing different kinematic parameters: αðtÞ ¼ 10 þ 15 sinðΩtÞ with κ ¼ 0.025, 0.05 and 0.1; αðtÞ ¼ 0 þ 7:5 sinðΩtÞ with κ ¼ 0.05 and αðtÞ ¼ 0 þ 15 sinðΩtÞ with κ ¼ 0.05, is the reduced frequency, with Ω the angular velocity and U ∞ the freestream velocity. In this paper we focus on the mode αðtÞ ¼ 10 þ 15 sinðΩtÞ (corresponding with κ ¼ 0.1) in which the airfoil moves through a deep dynamic stall regime with large flow unsteadiness. The pitching axis is 0.25c downstream of the leading edge on the chord line as depicted in Fig. 1. The freestream velocity is 14 m/s, which results in a Reynolds number of 1.35 Â 10 5 based on the chord length. A turbulence intensity (TI) of 0.08% was measured at a free stream velocity of 35 m/s but no value was provided at the measured velocity. It was also not indicated at which position in the wind tunnel this TI was measured.
Multiple hot-film sensor arrays, surface pressure transducers, hot-wire probes placed in the wake, and smoke visualization were used to investigate critical boundary-layer events (laminar separation, laminar-turbulent transition, flow reversal, reattachment, etc.), aerodynamic forces evolution and flow pattern development during a pitching stroke. Both the instantaneous wake velocities and measured pressure signals were ensemble-averaged over more than 100 oscillations to obtain the mean and fluctuating velocities, and the unsteady aerodynamic forces.

Computational domain and grid
The rectangular computational domain contains a circular subdomain with a radius of 2c as shown in Fig. 2. The width of the computational domain is 6c, which is identical to the height of the wind tunnel test section in the experiments of Lee and Gerontakos (2004). Based on the size of the computational domain adopted in a previous study by Wang et al. (2010Wang et al. ( , 2012, the length is chosen to be 28c, which is larger than the length characterizing the wind tunnel test section, in order to allow for the full development of the wake flow. The distance from the domain inlet to the airfoil is 8c to avoid an unphysical influence of the inlet boundary on the induction field upstream of the airfoil. The strategy employed to construct a high-quality grid for the URANS simulations is as follows: (1) a dense grid is used in the boundary layer, which is gradually coarsened away from the airfoil surface. The first cell height is taken as 6.00 μm (4 Â 10 À5 c) to capture the laminar and transitional boundary layer regions. This results in a maximum y þ value along the airfoil during one complete pitching cycle of 0.001 < y þ < 1. Values of y þ falling outside this range might affect the transition onset location, as observed in a study on the effect of y þ on the location of transition for a flat plate (ANSYS Fluent, 2015); (2) grid points are clustered towards the leading and trailing edges in order to capture small vortices and large flow gradients as shown in Fig. 3(c and d); (3) a C-type grid is used in areas near the airfoil with the distance from the C-type grid boundaries to the airfoil of 1c, as depicted in Fig. 3(b) to facilitate the process of generating a structured grid; (4) a maximum edge length expansion ratio of 1.07 is used in the wall normal direction, which nowhere exceeds 1.2 in the freestream direction. The above modeling strategy is executed in the pre-processor Gambit 2.4.6, resulting in a grid with about 1:18 Â 10 5 quadrilateral cells, shown in Fig. 3(a).
In addition to the grid used in the above reference case, two more grids with the same grid topology were generated in order to perform a grid-sensitivity analysis, see Section 2.4. Accordingly, the number of cells  in both the streamwise and wall normal directions were subsequently changed by approximately a factor of ffiffiffi 2 p and 1= ffiffiffi 2 p in every direction, resulting in a fine grid with 236,460 cells and a coarse grid with 64,398 cells, respectively, as depicted in Fig. 4. The properties of the coarse, medium, and fine grids are listed in Table 2.
In order to understand the effect of an increasing y þ value on the calculated aerodynamic forces acting on the pitching airfoil, the first layer height of 6.00 μm in the medium grid (reference case), corresponding with a maximum y þ of 0.9, is further increased to 30.00 μm and 90.00 μm without changing the growth ratio or the grid structure outside the boundary layer. This generates another two grids with maximum y þ values of 4.1 and 10.1, respectively. The grid details at the near wall area are shown in Fig. 5, together with the reference mesh (also shown in Fig. 4 (b)). The grid properties of the lower y þ , medium y þ , and higher y þ grids, which apply to the y þ sensitivity analysis, are listed in Table 3.
Based on a description of the airfoil profile shape with 200 node coordinates, the generated NACA0012 profile considered has a trailing edge thickness of 0.0025c at x/c ¼ 1, as shown in Fig. 6. Two trailing edge configurations are generated to evaluate their impact on the dynamic force coefficients (see Section 2.4). The sharp trailing edge utilized in the reference case is obtained by extending the upper (suction) surface and lower (pressure) surface using straight lines tangential to the profile surface at both end points, resulting in an extension of the chord by 1.0%    Table 3 Some grid properties of the lower y þ , medium y þ , and higher y þ grids for the y þ sensitivity analysis. as also shown in Fig. 6(left). For the blunt trailing edge, the two points at the suction surface and pressure surface are connected by a vertical line, as shown in Fig. 6(right).

Boundary conditions
At the inlet, a freestream velocity of 14 m/s is imposed with a turbulence intensity of 0.080%. The incident turbulence intensity, which is computed at the location just upstream of the airfoil in a domain with the airfoil absent, is found to be 0.079%, showing that the decay is minimal due to the absence of turbulence generation along this distance. The turbulence length scale l is set in accordance with the physical dimensions of the flow problem as l ¼ 0:07D h =C 3=4 μ ¼ 0:438 m, in which D h denotes the hydraulic diameter of the test section and C μ is a model constant equal to 0.09. The airfoil surface is modeled as a no-slip wall with zero roughness, and its kinematics are imposed on the circular subdomain through a UDF (user-defined function), which prescribes its motion relative to the rest of the computational domain. To achieve a parallel flow at the top and bottom of the domain, a symmetry boundary condition is applied, assuming zero normal velocity and zero normal gradients of the flow quantities. At the outlet of the domain a zero static gauge pressure is imposed.

Turbulence model
In order to limit the computational demand of the analyses in the sensitivity study, an incompressible 2D URANS modeling approach is employed. Typical turbulence models in this approach are the k-ε or k-ω models in their different forms. Flow fields predicted by k-ε models are generally inaccurate under adverse pressure gradients and boundary layer separations (Wilcox, 1993), and these models are therefore not recommended for application to pitching airfoils where both effects are of critical importance. Although this downside can be remedied by employing an equation for the specific dissipation rate ω instead of the turbulent dissipation ε, the drawback of the standard k-ω model is that the solution is too sensitive to freestream values of k and ω outside the shear layer (Menter, 1992). The SST k-ω model is designed to avoid such freestream oversensitivity. The Transition SST model (also known as the γ-Re θt model) is based on the combination of the SST k-ω model with two other correlation-based transport equations, one for the intermittency γ and one for the transition momentum thickness Reynolds number Re θt (Menter et al., 2006). For wall-bounded flows this model is better able to predict the laminar-turbulent transition. Hence, in this study the Transition SST model is used for the 2D URANS calculations.

Numerical procedure
The simulations were performed with the commercial CFD solver ANSYS Fluent version 16.1 on an HP DL360P Generation8 server node equipped with 16 Intel Xeon E5-2690 2.90 GHz cores and 347 GB fully buffered DDR2 memory. In order to achieve results with minimal numerical effort and to circumvent convergence difficulties, the simulation is performed by two subsequent calculations, namely (1) a steady-state solver is used to initialize the flow field at an AOA of 10 , during which the motion of the pitching airfoil is mimicked by the Moving Reference Frame (MRF) approach. In this approach the flow around the moving airfoil is modeled in steady state by incorporating an additional source term in the equations of motion; (2) using the steady state flow pattern as the starting point, a transient moving grid simulation is carried out, in which the sliding mesh method is employed to couple the rotational subdomain to the rest of the domain. The pressure-velocity coupling is warranted by the SIMPLE algorithm, the pressure interpolation is performed with a second-order scheme, and second-order discretization schemes are used for both the convection terms and the viscous terms of the governing equations. Iterations are terminated when all scaled residuals are below 1 Â 10 À5 . The following guidelines should be taken into account to determine a suitable time step size: (1) To limit the interpolation error at the sliding interface, the marching length for each time step is chosen below the minimum cell length along the interface. The maximum rotation speed at the interface Ω max ¼ dα dt j t¼0 ¼ 4:89 rad=s : The minimum cell length along the interface in the aforementioned grid is Δx min ¼ 2:45 Â 10 À3 m : The maximum time step size Δt max is determined by the condition Δt max Ω max 2c ¼ Δx min ; giving Δt max ¼ 1:67 Â 10 À3 s : A much more restrictive criterion based on the Courant-Friedrichs-Lewy (CFL) number at the conservative interface was introduced by Trivellato et al. (2014): Δα is the angular marching step, R g the radius of the rotating grid, which equals 2c as depicted in Fig. 2, Δx ¼ 0:0115 m is the cell length at the upwind interface and Λ is the tip speed ratio of the rotating grid, defined as Λ ¼ ΩRg U∞ . According to this criterion, the time step size should meet the condition Δt < 1:1 Â 10 À4 s, which corresponds to a maximum angular marching step of 1/32 (0.031 ) during one complete pitching cycle.
(2) The time step size is expected be small enough to guarantee convergence for each time step within a prescribed number of iterations. With a maximum of 20 iterations per time step, this requirement could not be met at higher AOAs when specifying a time step as large as Δt ¼ 1:1 Â 10 À4 s. A time step size of 3:0 Â 10 À5 s was found to be necessary to achieve residuals below 1 Â 10 À5 over a complete pitching cycle. (3) After obtaining a preliminary time step size from satisfying the above criteria, the effect of the time step size on the results is determined in a sensitivity analysis, which is presented in Section 2.4.
Hence, a time step size of 3:0 Â 10 À5 s is used in the reference case, which corresponds to 11,218 time steps per pitching cycle. The calculation is first run until a statistically converged state is reached, after which statistic measures are sampled for post-processing.
The parameters considered in the sensitivity analysis in Section 2.4 include the size of the computational domain, the turbulence properties of the incoming flow, the temporal and spatial resolutions, the y þ value, the turbulence model, and the trailing edge detailing. An overview of the studied computational parameters is presented in Table 4. In the sensitivity analyses the specific parameters considered are varied individually with respect to the reference case.

Convergence of simulations and results of the reference case
The time history curve of the lift coefficient is depicted in Fig. 7 for the reference case, showing a sudden stall after the lift coefficient reaches TI ¼ turbulence intensity, l ¼ turbulence length scale, gridres ¼ grid resolution, timres ¼ temporal resolution, maximum y þ denotes the maximum y þ during a complete pitching cycle, T.E. ¼ trailing edge.
its peak value, followed by a generally unsteady behavior during the downstroke phase. This is qualitatively in agreement with the general lift force features for dynamic stall (Gharali and Johnson, 2013). In order to quantify the statistical convergence of this simulation, the standard deviations (σ) of the differences in force coefficients between two successive pitching cycles are calculated and listed in Table 5. Together with Fig. 8, where the instantaneous force coefficients are plotted against the pitching cycle at six AOAs (15 ↑, 23 ↑, 22 ↓, 15 ↓, À5 ↓, À3 ↑, where ↑ denotes upstroke and ↓ designates downstroke), this illustrates that force coefficients in the second pitching cycle are different from those in the first pitching cycle, and that the differences between the force coefficients in the second cycle and the rest of the cycles are negligible. Hence, sampling of results is started after obtaining a statistically steady state at the second cycle. The integral aerodynamic force coefficients are sampled for a total flow time of 4T with T ¼ 2π/Ω the oscillation period, during which a maximum CFL number of 5.2 and a maximum y þ of 0.9 along the airfoil surface are observed.
The streamwise velocities at sampled AOAs on the vertical line m located one chord length downstream of the trailing edge are monitored and depicted in Fig. 9. This figure indicates that the largest variation of velocity occurs in the downstroke (22 ↓ and 15 ↓ in the figure). The velocity variation along this line is also quite strong at incipient stall (23 ↑). Less significant deviations are found at small AOAs during both upstroke and downstroke phases (15 ↑, À5 ↓, À3 ↑). These findings agree qualitatively with observations in the experimental work of Lee and Gerontakos (2004).
The calculated force coefficients are compared with wind tunnel measurements (Lee and Gerontakos, 2004) and other URANS numerical validation studies on this experimental dataset published by Wang et al. (2012) and Gharali and Johnson (2013), see Fig. 10. It should be noted that the experimentally determined aerodynamic forces were obtained by ensemble-averaging the measured pressure signals over 100 oscillation cycles. The current CFD results, however, are ensemble-averaged over only 4 oscillation cycles, since, without any perturbation in the incoming flow, the solver predicts a repetitive force pattern (as discussed above), in contrast with experimental observations. Thus, an improved result would not be expected by ensemble-averaging the data over more cycles.
In terms of the lift coefficient the current CFD results show an improved agreement with the experiments compared to previous numerical works, in that (1) the predicted peak value of the lift coefficient is closer to the measurements, and (2) the lift coefficient in the downstroke phase oscillates around the mean value obtained from the measurements and does not show the significant underprediction shortly after stall (from α ¼ 23 ↓ to 17 ↓) observed by Wang et al. (2012) and Gharali and Johnson (2013). However, the maximum stall angle is underpredicted and a sharp secondary peak is observed at the maximum AOA, which is in contrast with the experimental observations. Interestingly, the results of Wang et al. (2012) correspond somewhat better with the experiment at this stage while the results of Gharali and Johnson (2013) closely mimic those of the current work, despite that the turbulence model adopted here is the same as that in Wang et al. (2012), and differs from that in Gharali and Johnson (2013). Furthermore, the behavior near the minimal angle of attack (α < 0 ) is qualitatively similar to that in the experimental observations, but is quantitatively different for all numerical studies, although the current results virtually coincide with those of Wang et al. (2012).
For the drag force there is generally a less satisfactory agreement between the numerical simulations and the experiment, particularly when α > 13 . This is due to persistent flow separations at high angles of attack, resulting in difficulties to accurately model the viscous effects adjacent to the airfoil surface when using RANS models. Additionally, current CFD results predict a secondary LEV that helps the recovery of lift and drag coefficients around the maximum angle of attack (α ¼ 25 ), while in the experiment the secondary LEV is observed at a smaller AOA (α ¼ 21.8 ↓), which is less dominant in the recovery of the lift coefficient, as indicated in Fig. 10.  The dynamic stall mechanisms can be illustrated through the visualization of flow structures at different airfoil AOAs. Fig. 11 depicts the instantaneous vorticity field at eight characteristic AOAs in which major flow features representative of the pitching airfoil are observed. At α ¼ 10 ↑ (Fig. 11(a)) the boundary layer is still attached with a thickened shear layer at the trailing edge compared to the leading edge, and the force coefficients vary linearly with the AOA. As the airfoil oscillates further upwards, an area of reversed flow spreads from the trailing edge towards the leading edge, a symptom that is shown in Fig. 11(b). Afterwards, a strong energy containing LEV forms at the leading edge and moves downstream over the airfoil surface. This results in an increase in integral forces, as represented by a larger slope of the force coefficient curves plotted in Fig. 10(a and b). At α ¼ 22.8 ↑ (Fig. 11(c)) the LEV spans roughly the whole chord, which is associated with the force coefficients reaching their peak values. This is followed by abrupt stall caused by the breakdown of the LEV, as shown in Fig. 11(d) for α ¼ 24.2 ↑. Simultaneously, a TEV is generated rapidly during the stall. This vortex is soon shed downstream, and replaced by a secondary LEV that spans most of the chord, as shown in Fig. 11(e). This leads to the recovery of the force coefficients observed in Fig. 10(a and b). The secondary LEV is subsequently shed into the wake again and a secondary TEV has formed at 22.8 ↓ (Fig. 11(f)), the flow pattern of which is very different from that at the same AOA upstroke (Fig. 11(c)). The alternating LEV and TEV formation and shedding continues throughout the downstroke, leading to strong fluctuations in the force coefficients. The local peak values are characterized by the maximum span of LEVs across the airfoil surface, and local minima corresponding to the breakdown of LEVs and the formation of TEVs. Around α ¼ 0 ↓ (Fig. 11(g)) the separated boundary layer reattaches to the airfoil, after which the force coefficients curves become less unsteady. However, small oscillations in the force coefficient curves are observed even at low AOAs; a representative vorticity field for α ¼ À2.5 ↓ is shown in Fig. 11(h), where small vortices are observed in the wake area. These are responsible for the small fluctuations in the lift coefficient around this AOA, which indeed are visible in Fig. 11(a).
In order to investigate the impact of the pitching motion on the aerodynamic performance of the airfoil, seven static pitch angle calculations were performed (at À5 , 0 , 5 , 10 , 15 , 20 and 25 ), using the same inflow conditions and computational domain as for the reference case. Likewise, a transient solver was employed in these seven cases. The obtained lift coefficient and drag coefficient are time-averaged over 2T ft after reaching statistical convergence, where T ft is the flow-through time (T ft ¼ L=U ∞ with L the length of the computational domain). The predicted values of the force coefficients are in good agreement with the measurements by Lee and Gerontakos (2004), as shown in Fig. 12. The force coefficients are also compared with their dynamic counterparts, which illustrates that the pitching motion postpones stall occurrence to a higher AOA.

Size of the computational domain
The blockage ratio of the 2D computational model, BR, is defined as the width of the airfoil at the maximum AOA divided by the width of the computational domain (BR ¼ sin25 c/W); its impact on the predicted aerodynamic forces is tested by varying the domain width W. A width of W ¼ 6c is used in the reference case, which coincides with the wind tunnel test section height and corresponds to a BR ¼ 7%. Additionally, widths of 12c (BR ¼ 3.5%) and 24c (BR ¼ 1.75%) are considered, whereby the upstream and downstream lengths of the computational domain remain fixed at 8c and 20c, respectively. Fig. 13 shows that the     Lee and Gerontakos (2004). The CFD results in pitching and static state are indicated as 'Pitching CFD' and 'Static CFD', respectively. blockage ratio influences the shape of the hysteresis loops of the force coefficients. For BR ¼ 7% the peak values of the force coefficients are slightly higher compared to those for BR ¼ 3.5% and BR ¼ 1.75%, resulting in a better agreement with the experimental peak value of the lift coefficient. In general, however, during the upstroke phase of the airfoil the lift coefficients for the different blockage ratios are quite similar. During the downstroke motion the lift coefficient following from the CFD calculations oscillates around the experimental data, whereby the oscillation appears to be stronger for a lower blockage ratio. For example, at α ¼ 12 ↓ a relative difference with respect to the experimental result is found of 45% for BR ¼ 1.75% and of 16% for BR ¼ 3.5%. This indicates that the agreement between CFD predictions and wind tunnel measurements decreases with an increasing width of the computational domain. It also shows that the aerodynamic forces are most sensitive to the domain width when the airfoil pitches downward, during which the flow is highly unsteady.

Grid resolution
The coarse, medium (reference case) and fine computational grids depicted in Fig. 4 with the grid properties listed in Table 2 are used for a grid-sensitivity analysis. It is found in Fig. 14 that the force coefficients during the upstroke phase are similar, including the lead peak values and the α values at which they occur. During the downstroke phase substantial differences are observed, especially in the range 10 ↓ α 20 ↓.
Here, the difference between the lift coefficient measured experimentally and that computed by the CFD simulations is the smallest for the medium mesh used in the reference case, and the largest for the coarse mesh. In this context, it may be concluded that the medium mesh results in an acceptable prediction of the experimental results.

y þ value
Apart from the overall grid resolution of the whole computational domain, the grid distribution near the airfoil wall is crucial for gaining an accurate representation of laminar and transitional boundary layers, and hence for the correct prediction of the flow transition and flow separation dynamics. The non-dimensional quantity y þ is directly linked to the grid resolution close to the wall. In the case of a flat plate, a maximum y þ value above 8 results in a transition onset upstream of the correct location (ANSYS Fluent, 2015). Due to a more complex flow pattern, the required maximum y þ value around an airfoil to ensure a high-fidelity prediction is expected to be lower than 8. This can be investigated by considering three different cases, with corresponding maximum y þ  values of 0.9, 4.1 and 10.1, see Table 3. Fig. 15 shows the effect of an increasing y þ value on the magnitudes of the aerodynamic forces. Surprisingly, the predicted force coefficients appear to be identical for the three cases, even for a maximum y þ value as high as 10.1. This indicates that the grid with a first layer height of 90 μm is still below the limit at which simulation results become inaccurate.

Turbulence model
In addition to the transition SST model used in the reference case, the ability of the SST k-ω turbulence model to accurately predict the force coefficients of the pitching airfoil is evaluated. The computational results obtained by the transition SST and SST k-ω turbulence models are shown in Fig. 16 together with the experimental results from the wind tunnel tests. The SST k-ω model underestimates the peak value of the lift coefficient by about 4% compared to both the experimental and transition SST model. Otherwise, a similar trend is obtained for both models in the upstroke phase. Conversely, during the downstroke phase the SST k-ω model predicts a significant drop in the lift and drag coefficients at 21.4 ↓, whereas a much smaller drop at 22.8 ↓ is observed in the transition SST results; this indicates a more intense shedding of the secondary LEV with the former turbulence model. Lastly, the SST k-ω model fails to predict the experimentally observed lift coefficient hysteresis for α < 0 , which indeed is captured by the Transition SST model. In general, it may be concluded that the results computed by the Transition SST model are significantly closer to the experimental results than those calculated by the SST k-ω model.

Trailing edge detailing
The effect of the trailing edge detailing on the dynamic force coefficients is depicted in Fig. 17. It is clear that this aspect has little effect on the lift and drag coefficients during the upstroke phase. During the downstroke phase, however, noticeable differences are observed starting at 20.8 ↓. The discrepancy remains along a large part of the downstroke phase, during which the flow is highly unsteady, and disappears for both the lift and drag coefficients when pitching down to about 2 ↓. The overall agreement with the experimental data is nevertheless similar for the sharp and blunt trailing edges, so that both models may be considered as acceptable representations.

Freestream turbulence properties and time step size
In the wind tunnel experiments of Lee and Gerontakos (2004), the measurement position of the turbulence intensity, the turbulence  . It turned out that these variations did not significantly change the time evolutions of the lift and drag coefficients, and therefore the graphical representation of these results is omitted here. In order to understand the influence of the temporal resolution on the dynamic force coefficients and to identify an appropriate time step size Δt, apart from Δt ¼ 3:0 Â 10 À5 s employed in the reference case, two additional values of 1:0 Â 10 À5 s and 1:2 Â 10 À4 s were investigated. These values correspond to maximum CFL numbers of 1.8 and 21, respectively, in comparison with a value of 5.2 for the reference case. As for the variation of the turbulent properties, the computational results obtained with this time step variation were virtually identical, so that their graphical representation is omitted here. Furthermore, the time step size used in the reference case (Δt ¼ 3:0 Â 10 À5 s) is also adopted for the remaining URANS simulations presented in this manuscript.

Performance evaluation
Four performance metrics are used to quantify the impact of the tested parameters presented in Section 2.4 and to determine the optimal computational parameters. These performance measures are the normalized mean square error (NMSE), the fractional bias (FB), the correlation coefficient (R), and the fraction of CFD results that are within a factor of 1.3 of the corresponding experimental data (FAC1.3), see Schatzmann et al. (2010). The metrics are defined as follows: Here, O i and P i denote observed (experimental) and predicted (CFD) values of the quantity of interest for the data point i, respectively, and angular brackets indicate the arithmetic mean over a set of data points, n is total number of points and σ indicates the standard deviation over a set of data points.
NMSE and FB cannot be used for the quantification of parameters that have both positive and negative values (Schatzmann et al., 2010;Gousseau et al., 2013), such as the negative lift coefficients observed when the airfoil oscillates in the range α < 0 . To circumvent this issue, the lift coefficients as well as the drag coefficients corresponding to α > 10 are used to generate these metrics, as these are strictly positive. Since most complex flow features such as large scale separation, dynamic stall and severe turbulence occur in this range of AOAs, and quantities of interest are generally most sensitive to the tested parameters in this range, this is expected to provide useful information. The ideal values indicating a perfect match between CFD results and experimental measurements would be NMSE ¼ 0, FB ¼ 0, R ¼ 1, and FAC1.3 ¼ 1. Table 6 and Table 7 respectively list the validation metrics for the lift coefficient and drag coefficient for each of the parameters studied in Section 2.4, with bold text indicating the statistics from the reference case. These metrics confirm that both dynamic force coefficients generally show a similar behavior for many of the studied parameters. The most sensitive parameter is the turbulence modeling strategy, where the transition SST model predicts a much better agreement with the experiment for the evaluated metrics compared with the SST k-ω model. For example, for the lift coefficient a difference of 0.079 and 0.115 is found in terms of R and FAC1.3; for NMSE and FB the differences are 0.022 and 0.034, respectively. In terms of the domain size, the FAC1.3 value for the lift coefficient of the reference case, with W ¼ 6c, is 0.101, and 0.026 higher than for the cases with W values of 12c and 24c (case 3), respectively, while no obvious differences are observed in NMSE, FB and R. However, the FAC1.3 value for the drag coefficient of the reference case is lower than for the other two domain sizes, making it impossible to objectively decide which is best. A clear difference between case 4 with TI ¼ 0.03% and the cases with TI ¼ 0.08% and 0.24% is observed in the FAC1.3 value for the lift coefficient, which is the lowest (0.694) for TI ¼ 0.03% and has a value of 0.824 and 0.823 for TI ¼ 0.08% and 0.24% respectively. However, again the opposite is observed for the drag coefficient, preventing an objective determination of the best parameter value. Regarding the grid resolution, a better agreement between medium grid results and experiments is observed for the lift coefficient than for both other resolutions, but once more the drag coefficient results show a different picture, favoring the fine grid instead. With respect to the maximum y þ and the time resolution, only minor differences between

. Computational domain and grid
The dimensions of the ground plane are taken the same as in the 2D URANS simulation (reference case) presented in Section 2. The geometry of the 2.5D computational domain results from constructing a ground plane and uniformly extending it in the third direction. The grid topology of the ground plane is somewhat different from that used in the 2D URANS case to ensure approximately square cells (Δx/Δy % 1) in the entire ground plane. The first cell height is 1:2 Â 10 À4 m (8 Â 10 À4 c) along the airfoil surface, resulting in an average y þ < 4 during the complete pitching cycle. Given that the force coefficients obtained from the 2D URANS calculations are not very sensitive to the trailing edge detailing on the scale studied here, see Section 2, and in order to simplify the creation of low aspect ratio cells at the trailing edge, the blade is modeled with a round trailing edge in the 2.5D LES simulation.
A structured grid is used in areas close to the airfoil while the overall grid in the rotational subdomain is non-conformal with 7 grid zones with different cell sizes. At least 5 cells are used for each zone in both the wall normal and streamwise directions; the cell growth ratio between two adjacent zones is limited to 2 within the rotational subdomain, and the cell growth ratio at the sliding interface, which connects the rotational subdomain and the rectangular domain, is set equal to 1 in order to minimize interpolation errors at the interface. In this way, the total number of cells is kept limited while maintaining a good accuracy, resulting in a total of 51,070 cells in the ground plane. The ground plane grid is then extended in the third dimension over a length of 0.5c, during which a maximum aspect ratio Δz/Δx 5 in each grid zone is guaranteed. The resulting 2.5D LES computational grid comprises 4.2 million hexahedral cells. Fig. 18 shows the grid topology in the near-wall area at two different levels of detail.

Boundary conditions
The boundary conditions used in the 2.5D LES calculations differ from those in the 2D URANS simulations in two ways: (1) periodic boundary conditions are imposed on both lateral faces, as depicted in Fig. 19; (2) at the inlet plane the Vortex Method (VM) is employed to generate a timedependent velocity profile that provides turbulent fluctuations on the mean velocity (14 m/s) with desired statistical properties in a relatively inexpensive yet accurate way (Mathey et al., 2006). The intensity and size of generated vortices depend on the local value of turbulent kinetic energy and the turbulence dissipation rate. The number of vortices used to generate turbulent fluctuations at the inlet plane is set to the default value of 190.

Numerical procedure
Due to the high computational cost the 2.5D LES calculations are performed on the Dutch national computing cluster (Lisa). A URANS calculation result with added random turbulent fluctuations is used as initialization. The standard Smagorinsky-Lilly SGS model with C S ¼ 0.1 is used, whereby the sub-grid eddy viscosity close to the wall is refined using a mixing length -related to the unresolved eddies -that is proportional to the wall normal distance. The discretization of the convection term in the filtered momentum equation is performed with a Table 6 Statistical performance of the computational parameters tested in Section 2.4 on the lift coefficient, quantified by the normalized mean square error (NMSE), fractional bias (FB), correlation coefficient (R), and fraction of data within a factor of 1.3 of observations (FAC1.3). Bold text indicates the statistics from the reference case.  bounded central difference scheme. In comparison to the central difference scheme, this spatial discretization scheme can prevent damping of small scale eddies by reduced numerical dissipation and is simultaneously more accurate for lower quality grids (Menter, 2012). A least squares cell-based gradient spatial discretization scheme is chosen, and the pressure interpolation is of the second order. The bounded second-order implicit scheme is used for time discretization, which provides a higher accuracy than a first-order scheme and is much more robust than the non-bounded version (Ferziger and Peri c, 2002). The non-iterative time advancement scheme is used to reduce calculation time. The PISO scheme is used for pressure-velocity coupling. The time step size used is 8 Â 10 À6 s. The resulting maximum CFL is well below 2 at most AOAs, while it exceeds 2 for the angle of attack from α ¼ 22 ↑ to 20 ↓ with a maximum value of around 15.28 at α ¼ 24.50 ↑ and then drops back below 2 again as the airfoil pitches down. The total sampled flow time is 13T, chosen as a compromise between computational cost and accuracy. Fig. 20 shows the dynamic force coefficients at several AOAs of interest as a function of the pitching cycle for the 2.5D LES; the 2D URANS results (see Fig. 8) are also included for comparison. It is found that (1) 2.5D LES results are consistent with 2D URANS results at low AOAs, (2) 2D URANS tends to predict higher lift and drag coefficients at high AOAs, e.g. α ¼ 23 ↑, α ¼ 22 ↓, and (3) the 2.5D LES results during the downstroke phase show a strong cycle-to-cycle variation around the rather smooth 2D URANS results. The only cycle-to-cycle variation observed in the URANS results is between the first and second cycles (in the initialization phase). In addition, the starting vortex generated at the onset of pitching motion near the leading edge also contributes to the variation at the fist cycle in both the 2D URANS and 2.5D LES simulations, due to the impulsive rotation of the airfoil (Panda and Zaman, 1994;Ol et al., 2009). For these reasons, the initialization phase is essential for both 2D URANS and 2.5D LES simulations.    The force coefficients are ensemble-averaged over the last 9 pitching cycles computed by the 2.5D LES model, and depicted in Fig. 21 together with the 2D URANS results, LES results taken from Kim and Xie (2016), and experimental results (Lee and Gerontakos, 2004). It can be seen that the present LES model leads to a better prediction of the lift coefficient than that computed from the LES model of Kim and Xie (2016). For the drag coefficient the 2D URANS prediction and the LES predictions are in good correspondence at low AOAs, but show significant differences at high AOAs. Some further observations are:

Results of 2.5D LES
(1) CFD results (both 2D URANS and LES) for both lift and drag coefficients are well in line with experimental observations at low AOAs, especially for the upstroke phase during which no significant boundary layer separations and vortex shedding are present.
(2) A clear increase of the slope of the lift coefficient beyond α ¼ 19 ↑ predicted by both the current LES simulation and the URANS simulation, which corresponds to the vortex detachment on the suction surface, is in agreement with the trend observed experimentally, although it there occurs at a higher AOA. Note that this slope increment is not clearly observed in the LES model of Kim and Xie (2016). (3) Both the 2D URANS (reference case) and the current LES predict a comparable peak lift coefficient as observed in the experiment, though at a lower dynamic stall angle. This lower dynamic stall angle is the reason for the difference in force coefficient recovery between CFD results and experimental data, since the second LEV responsible for the force coefficient recovery is generated upon further increment of the AOA beyond the dynamic stall angle. (4) The LES results (both the current results and those obtained by Kim and Xie, 2016) show a seemingly less unsteady behavior during the downstroke phase than the 2D URANS results. This is likely to be due to the ensemble-averaging of the LES data, considering that the unsteadiness predicted by URANS simulations does not palliate after ensemble-averaging (as observed in Fig. 8). In fact, the LES results during the downstroke phase are very sensitive to small physical and/or numerical perturbations, resulting in cycle-to-cycle variations of the force coefficients shown in Fig. 20. Note that in the downstroke phase the current LES results show a closer agreement with the experimental data than the results of Kim and Xie (2016), both in terms of the lift and drag coefficients. In addition, the hysteresis loop of the lift coefficient for α < 0 , which is not present in the Kim and Xie (2016) results, is predicted in the current 2.D LES simulation, although quantitative differences with the experiments can be observed. (5) The current work shows that 2.5D LES is able to predict force coefficients in deep dynamic stall with reasonable accuracy, although a significant improvement over the 2D URANS simulation results in terms of the agreement with experimental data is lacking. Possible reasons for the discrepancy between CFD simulations and the experimental observations are discussed in Section 3.3.
It is noted that the maximum CFL number substantially exceeds 2 at high AOAs. In order to test if this high maximum CFL number has an effect on the values of the force coefficients, a second simulation was performed of a single cycle using a reduced time step of 8 Â 10 À7 s in the range from α ¼ 22 ↑ to 20 ↓. This time step adaptation warrants a maximum CFL value strictly below 2 throughout the calculation; however, this did not improve the results compared to the case with Δt ¼ 8 Â 10 À6 s. Fig. 22 shows the comparison of wall pressure coefficients (C p ) along three sampling lines in the spanwise direction, taken at a mutual distance of 0.25c (line 1: z ¼ À0.25c, line 2: z ¼ 0; line 3: z ¼ 0.25c). The results are depicted for four AOAs, namely two angles in the upstroke phase ((a) α ¼ 15 ↑, (b) α ¼ 23 ↑) and two angles in the downstroke phase ((c) α ¼ 20 ↓, (d) α ¼ À5 ↓). For each AOA the numerical data along each sampling line are ensemble-averaged over four pitching cycles. It can be seen that the wall pressure coefficients on the pressure side generally coincide for the three lines, whereas some differences appear on the suction side. The differences are quite small during the upstroke phase of the pitching motion, see Fig. 22(a) and (b), and become significant during the downstroke phase, see Fig. 22(c) and (d). For example, for α ¼ 20 ↓ the wall pressure coefficient at z ¼ 0.9c downstream of the leading edge at line 3 is 60.4% higher than that at the same streamwise location at line 1. This shows again that there is a much larger flow unsteadiness during the downstroke phase than during the upstroke phase. It is possible that significant 3D effects due to massive flow separations are also present during the downstroke phase, as was found by Martinat et al. (2008); averaging the response over significantly more pitching cycles would be necessary to test this hypothesis.
In Fig. 23, the obtained wall pressure coefficients from the 2.5D LES are compared to those from the 2D URANS (reference case). Both the 2.5D LES predictions and 2D URANS predictions are ensemble-averaged over four pitching cycles. Furthermore, the 2.5D LES results are also averaged over the three lines in the spanwise direction indicated in Fig. 21. Comparison of force coefficients obtained from the current calculations (2D URANS and 2.5D LES), the LES simulations performed by Kim and Xie (2016), and the wind tunnel tests of Lee and Gerontakos (2004): (a) lift coefficient; (b) drag coefficient. The 2D URANS results correspond to the reference case presented in Section 2. The results from the current LES simulation are ensemble-averaged over 9 pitching cycles. Fig. 22. In general, similar wall pressure distributions are observed at α ¼ 15 ↑ (Fig. 23(a)), α ¼ 23 ↑ (Fig. 23(b)) and α ¼ À5 ↓ (Fig. 23(d)). In contrast, a remarkable difference is found at α ¼ 20 ↓, see Fig. 23(c), at which the wall pressure coefficient predicted by the URANS approach is significantly more negative, with a maximal difference of 0.51 at z ¼ 0.38c. This discrepancy is especially noticeable in the interval from α ¼ 24 ↑ to 18 ↓ through a significant variation of the lift coefficient, see To visually inspect the turbulence structures, two iso-surfaces with the Q- where Ω v is the vorticity and S is the strain rate as indicated above, at α ¼ 10 ↑ corresponding to Q ¼ 5 Â 10 3 s À2 and α ¼ 20 ↓ corresponding to Q ¼ 5 Â 10 5 s À2 , are shown in Fig. 24.
The colors in the plots indicate the instantaneous streamwise velocity. It can be seen that at α ¼ 10 ↑ the flow is attached close to the leading edge and starts to separate at about 0.3c from the leading edge, whereby the transition into a fully turbulent shear flow occurs at around 0.6c.
Conversely, at α ¼ 20 ↓ the flow is only attached at the leading edge, followed by substantial separations over the suction surface.

Discussion
The CFD modeling results for an airfoil during pitching motion were compared to the experimental data of Lee and Gerontakos (2004). Although an improved agreement is obtained compared with previous modeling work by other researchers (Wang et al., 2012;Gharali and Johnson, 2013;Kim and Xie, 2016), there are still some notable discrepancies. The results of the extensive parameter variation study performed in the present study imply that these discrepancies are not related to any of the parameters investigated. A potential source of discrepancy are the simplifications made to the 2.5D LES domain and grid in order to make the calculation computationally tractable, such as the presence of non-conformal grid interfaces, the use of cells with an aspect ratio up to 5, and the limited length in the spanwise direction. In fact, these drawbacks can only be resolved at the expense of a significant increase in computational time. The use of other scale-resolving simulation models, such as DES or Stress-Blended Eddy Simulation (SBES), in combination with a limited computational demand, would be an interesting topic for future research.
Additionally, flow behavior at a moderate Reynolds number, such as the value of 1.35 Â 10 5 adopted in this study, also creates large uncertainties in terms of flow aerodynamics during deep dynamic stall. At the chosen Reynolds number several flow phenomena are critically dependent on the exact flow conditions, such as the location of flow separation, the laminar to turbulent flow transition point and the height of the boundary layer. Hence, it is very challenging to capture surface pressure data with high fidelity both experimentally and numerically, as was recently demonstrated by Tank et al. (2017) in a validation exercise for the dynamic force coefficients on a static NACA0012 airfoil at Reynolds numbers ranging from 1 Â 10 4 to 1.5 Â 10 5 . In this respect, the following experimental uncertainties were noted, which might affect the comparison with CFD results presented in this communication: As discussed in the experimental study (Lee and Gerontakos, 2004) and a recent numerical validation study by Kim (2013), the length of the Tygon tubing, which connects orifices distributed in the airfoil to the external pressure transducers, has an effect on the unsteady pressure signals when measuring the surface pressure of the pitching airfoil. Such an effect results in a constant time lag in all pressure signals with an oscillation frequency κ above 2.95 Hz (Lee and Gerontakos, 2004), which sets a limit on the reduced frequency of 0.0993, and makes that for κ > 0.0993 the curves of lift, drag and moment coefficients can only be considered qualitatively. In the simulated case κ equals 0.1, which is near the above limit. However, in response to this concern, as argued by Kim (2013), a reduction of the tubing length did not show any improvement. The turbulence properties of the incoming flow are to some extent uncertain in the experiments of Lee and Gerontakos (2004). As mentioned in Section 2, it is unclear whether the reported value of the turbulence intensity of 0.08% represents the total turbulence intensity or only the streamwise component. Further, the location at which this quantity was measured is unknown. In the case of LES, additional information is required on turbulence length scales and the spectral distributions, which is also not available. It should be noted that this aspect may further be affected by the chosen LES inflow method. The pressure orifices for measuring pressure distributions across the airfoil in the experiment of Lee and Gerontakos (2004) were dispersed from the leading edge to 0.78c in the streamwise direction, by which the pressure data close to the trailing edge is excluded. This may influence the agreement between LES results and experimental results in terms of the peak values of lift and drag coefficients, as also argued by Kim and Xie (2016).

Conclusions
The present study focuses on the CFD modeling of the aerodynamic behavior of a NACA0012 airfoil subjected to a large pitching motion, which at a Reynolds number of 1.35 Â 10 5 is characterized by a deep dynamic stall. A 2D URANS formulation equipped with the Transition SST turbulence model was firstly used to investigate the sensitivity of the response to several computational parameters. It is observed that for an airfoil subjected to a pitching motion, during the downstroke phase the aerodynamic forces are found to be much more sensitive to changes in the computational parameters than during the upstroke phase. It is shown that the type of turbulence model chosen has a profound effect on the dynamic force coefficients, whereby the best agreement with the experimental results for the two turbulence models tested here is obtained for the Transition SST turbulence model. Furthermore, in the downstroke phase the predictions for the force coefficients are found to be rather sensitive to the variation of the grid resolution. Conversely, the range of studied values for the first cell height and the time step size do not result in a significant variation in the dynamic force coefficients. In addition, comparable CFD results are obtained for geometries with sharp and blunt trailing edges.
The 2.5D LES simulations show significant flow unsteadiness during the downstroke phase, but no significant improvement over the URANS simulations. This was possibly due to an insufficient grid resolution, which unfortunately could not be tested due to limitations in computational capacity. Despite this, with the ongoing increase in computational power 2.5D LES simulations have the potential to be applied for accurate modeling of deep dynamic stall. The results of the present study may serve as a guideline to further improve wind tunnel experiments and CFD simulations for pitching airfoils subjected to deep dynamic stall, a complicated, challenging and rich problem that deserves further research.