This paper presents large eddy simulations (LESs) of symmetric and asymmetric (cambered) airfoils forced to undergo deep dynamic stall due to a prescribed pitching motion. Experimental data in terms of lift, drag, and moment coefficients are available for the symmetric NACA 0012 airfoil and these are used to validate the LESs. Good agreement between computed and experimentally observed coefficients is found confirming the accuracy of the method. The influence of foil asymmetry on the aerodynamic coefficients is analyzed by subjecting a NACA 4412 airfoil to the same flow and pitching motion conditions. Flow visualizations and analysis of aerodynamic forces allow an understanding and quantification of dynamic stall on both straight and cambered foils. The results confirm that cambered airfoils provide an increased lift-to-drag ratio and a decreased force hysteresis cycle in comparison to their symmetric counterparts. This may translate into increased performance and lower fatigue loads when using cambered airfoils in the design of vertical axis turbines (VATs) operating at low tip-speed ratios.

## Introduction

Wind and tidal turbines work under highly turbulent conditions and are subjected to a large range of velocities and turbulence intensities compromising their survivability, specifically from a structural integrity and a service life point of view. Material fatigue in these turbines as a result of repetitive dynamic loads can lead to eventual failure [1]. The loads are mainly due to dynamic stall of the blades characterized by large flow separation from the blade and shedding of energetic vortices. There is an obvious need to understand dynamic stall in vertical axis turbine (VAT) rotors with the goal to damp or reduce its impact on the structure [2]. Over the past decades, the study of dynamic stall has been focused mainly on the design of helicopters [3], micro-aerial vehicles, as well as wind and tidal turbines [4] for which further research is still needed in order to improve current designs.

The complexity of dynamic stall on moving airfoils was investigated experimentally throughout different motion patterns, such as heaving, plunging, and pitching, or combinations of them, respectively. During extensive experimental studies [5–9], it was observed that the aerodynamic behavior of an airfoil undergoing pitching motion is dominated by dynamic stall, which modifies its aerodynamic characteristics compared to the steady-state ones. Dynamic stall is related to the generation of a leading edge vortex (LEV) at high angles of attack, greater than the static stall angle, which overshoots the lift generation capabilities of the airfoil. Under deep dynamic stall conditions, the shedding of the LEV provokes large flow separation over the entire suction side of the airfoil, and this is subjected to successive generation and shedding of a series of leading and trailing edge vortices (TEVs). In this poststall regime, the airfoil loses its aerodynamic capabilities and experiences large force fluctuations until front-to-rear flow reattachment is achieved during the downstroke motion recovering its ability to generate lift.

Despite all the research undertaken to date, it has not been possible to extract a unique conclusion about the nature of dynamic stall. Choudhry et al. [10] remarked that dynamic stall in moving air/hydrofoils depends on many factors such as blade geometry (e.g., thickness or cambering), Reynolds and Mach numbers, oscillation frequency, or movement pattern (pitching, heaving, plunging or ramp-type motion), among others. Concerning the pitching motion, the three key factors are the pitching oscillation frequency, preset angle of attack, and pitch amplitude. These determine whether deep or light stall occurs, i.e., the flow separation region is extended over the entire suction side or just over a smaller section of the foil closer to the trailing edge [8], respectively.

Experimental work on dynamic stall started in the 1970s. However, it was not until the 2000s when computational fluid dynamics models were employed as a complementary tool to experiments aiding in the analysis of dynamic stall. Since then, many authors have studied dynamic stall using numerical simulations as they are able to provide more detailed information than experiments, such as detailed velocity and pressure fields around airfoils as well as visualizations of turbulent flow structures. Akbari and Price [11] analyzed the dynamic stall of a NACA 0012 using a vortex method for a range of Reynolds numbers (Re* _{c}* = 3 × 10

^{3}to 1 × 10

^{4}). Martinat et al. [12] simulated a NACA 0012 and compared different turbulence models in two-dimensional (2D) and three-dimensional (3D) Reynolds-averaged Navier–Stokes (RANS) simulations, which appeared to have a noticeable influence on the predicted aerodynamic performance of the foils. Wang et al. [13] performed a similar analysis using the experimental work from Lee and Gerontakos [7] to validate their model. Their results showed that the computed aerodynamic coefficients are sensible to the chosen turbulence model especially during the shedding of the LEV. Their extension from 2D RANS to 3D RANS and then 3D detached eddy simulation (DES) provided different predictions in terms of the airfoil's aerodynamic behavior with the largest differences at high angles of attack, i.e., when flow separation is more predominant, highlighting that 3D models are required to accurately resolve dynamic stall due to its 3D nature.

Eddy-resolving approaches, such as large eddy simulation (LES), are able to compute the instantaneous flow field, which is required to accurately represent the time dependence of turbulence structures [14], such as the LEV present in the flow around pitching airfoils. Kim and Xie [15] obtained a good match with the experimental data from Lee and Gerontakos [7] using LES for various pitching frequencies. Their results evidenced that a remarkable improvement in the prediction of the LEV or laminar-to-turbulent shear layer transition is obtained when using LES instead of RANS. Visbal [16] studied the behavior of a heaving airfoil using implicit LES (ILES). The computed flow field matched closely the vortical structures that were observed during laboratory experiments. They revealed important flow phenomena such as Kelvin–Helmholtz instabilities during the downstroke motion developing on the upper side of the airfoil, as well as other instabilities within the LEVs or TEVs. The use of ILES on moving airfoils was extended by Visbal et al. [17] with an emphasis on the generated flow structures, which again agreed remarkably well with those observed in the experiments. The eddy-resolving nature of LES (and ILES) is particularly important in the simulation of complex flows dominated by energetic large-scale structures, such as those found during dynamic stall, which RANS is not capable of resolving due to its time-averaging nature. Finally, the recent work of Rosti et al. [18] with the direct numerical simulation (DNS) of an airfoil subjected to a ramp-up motion has provided detailed insights into the vortex generation and development during the different phases of the airfoil motion.

The numerical simulation of a moving airfoil is a challenging task for any high-accurate numerical method and requires careful treatment. Body-fitted models are often adopted for the simulation of moving bodies but they are sometimes limited to relatively small body displacements due to numerical stability reasons. Additionally, variable reallocation is required at each time-step [19], which notably increases the computational load especially when moving bodies are simulated using LES or DNS. In this work, an immersed boundary (IB) method [20] is adopted to represent the moving airfoil geometry, which treats the solid body as a detached Lagrangian grid that communicates with the fluid mesh using interpolation functions. The fact that the solid mesh is not embedded in the fluid mesh avoids additional computations, which lead to a smaller computational cost than that of body-fitted models. Additionally, the use of the IB method in rectangular Cartesian meshes together with fast Poisson equation solvers, e.g., multigrid methods, allows to perform expensive high-fidelity simulations with a reasonable amount of computational resources.

In the IB method, the fluid mesh does not conform to the immersed solid body, which compromises the accurate representation of the boundary layer in the flow around airfoils. Recent publications have demonstrated that the IB method is able to accurately reproduce the flow around airfoils even at medium-to-large Reynolds numbers whenever adequately fine fluid meshes are adopted. Castiglioni et al. [21] performed a LES of a NACA 0012 with a fixed angle of attack at Re* _{c}* = 5 × 10

^{4}and achieved good agreement of aerodynamic coefficients computed with the IB method compared with those obtained from body-fitted model simulations. Tay et al. [22] studied different IB methods for the simulation of flapping wings and showed, when comparing the results to experimental measurements, that the accuracy of their IB model is similar to that of a body-fitted model. Zhang and Schluter [23] focused on the influence of the LEV on the lift generation of a flat plate undergoing a sinusoidal motion for a range of Reynolds numbers between 440 and 21,000. Their results indicated that the IB method performs well in the prediction of flow separation and vortex shedding mechanisms. Ouro and Stoesser [24] proved the accuracy of the IB method in the LES of vertical axis tidal turbines (VATTs) whose blades experience a loading cycle similar to that of pitching airfoils. They highlighted how the moving turbine blades undergo dynamic stall along most of their rotation cycle, which affects its performance in particular at relatively low rotational speeds.

This study aims at providing an appreciation of dynamic stall in both symmetric and asymmetric (cambered) pitching airfoils, and quantifies the aerodynamic properties of a cambered airfoil in comparison to those of the symmetric equivalent when undergoing deep dynamic stall. The data and flow visualizations are then interpreted in the context of lift-driven Darrieus-type VATs for which the pitching airfoil is an ideal surrogate system. The performance of VATs is compromised by the occurrence of light and deep stall, which is the result of the constantly-changing angle of attack between the oncoming fluid flow and the rotating rotor blades [24]. This work is motivated by the fact that past research on Darrieus-type VATs [25,26] has demonstrated that the turbine's performance improves when adopting cambered blade profiles. Additionally, Choudhry et al. [10] stated that asymmetric blade shapes tend to experience smaller force hysteresis cycles, i.e., lower difference between maximum and minimum load magnitudes, an additional benefit for VATs because it reduces load amplitudes and thus material fatigue.

The paper is organized as follows: Sec. 2 describes the dynamic stall phenomenon on pitching airfoils and how it dominates the driving physics of VATs. Section 3 presents the numerical framework together with the computational setup of various pitching airfoil simulations. In Sec. 4, the sensitivity to spatial and temporal resolution on the simulation results is presented together with the validation of the code. The effect of blade cambering on the aerodynamic properties of a pitching airfoil is analyzed by comparing lift, drag, and moment coefficients of NACA 0012 and NACA 4412 airfoils in Sec. 5. Conclusions and design criteria to be followed in the design of VATs are summarized in Sec. 6.

## Dynamic Stall in Vertical Axis Turbines

*N*), chord length (

_{b}*c*), and its radius (

*R*) determines its solidity

*σ*=

*N*/2

_{b}c*πR*, which is the proportion of the turbine's swept circumference length covered by the blades. According to Amet et al. [27], the range of rotational speeds at which the turbine operates depends on its solidity: the greater

*σ*, the lower the tip speed ratio at which power is extracted most efficiently and vice versa. During the rotation of a VAT, its blades face a constantly varying effective angle of attack,

*α*, relative to the oncoming fluid flow, and this angle is defined [24] as

*θ*denotes the rotated angle of the turbine rotor, and

*λ*is the tip speed ratio (= Ω

*R*/

*U*

_{0}, where Ω is the rotor's rotational speed, and

*U*

_{0}is the freestream velocity). The maximum angle of attack,

*α*

_{max}, a VAT blade attains is a function of

*λ*and is calculated from

Figure 1(a) presents the variation of the angle of attack *α* of a turbine blade over the upstream half of its revolution, i.e., 0 deg < *θ* < 180 deg. At all rotational speeds, the blade overcomes the static stall angle, *α _{ss}*, which means, in terms of physics, that the flow separates and forms a small recirculation or flow reversal zone on the suction side of the blade. At tip speed ratios lower than approx. 3.0, the dynamic stall angle,

*α*, is surpassed, and physically this means that the separated flow does not reattach on the blade leading to a dramatic reduction in lift. As it can be seen as well from this figure is that the lower the tip speed ratio, the earlier dynamic stall occurs, i.e.,

_{ds}*α*>

*α*, and hence the longer the rotor blade undergoes deep dynamic stall during its rotation. Increasing the value of

_{ds}*α*, for instance by selecting an airfoil shape that is less prone to flow separation reduces the extent of deep dynamic stall at given tip speed ratios, and this is beneficial to VATs. Note that if

_{ds}*α*is smaller than

_{ds}*α*

_{max}, then the blade does not undergo deep dynamic stall as it is the case for

*λ*> 3.0 [24].

The value of *α _{ss}* is intrinsic to the airfoil's geometry and flow regime, whereas

*α*is a function of the parameters defining the motion of the airfoil, e.g., pitching frequency and/or amplitude of pitch angle, respectively [10]. Cambering an airfoil profile aims at providing a larger value of

_{ds}*α*which reduces or delays full flow separation, i.e., postponing

_{ds}*α*>

*α*. Considering a VAT rotates at tip speed ratio of 2.0, Fig. 1(b) illustrates those regions over a revolution the turbine blades undergo deep dynamic stall (

_{ds}*α*>

*α*, gray region) and light- or no-dynamic stall (

_{ds}*α*<

_{ss}*α*<

*α*or

_{ds}*α*<

*α*, respectively).

_{ss}Past research suggested that tidal versions of VATs, also known as VATTs, perform at their peak efficiency when *λ* = 2.0 and the maximum torque is generated at *θ* ≈ 90 deg [24,27,28]. Figure 1(b) evidences that during maximum power generation, VATT blades are under deep dynamic stall, i.e., *α*(*θ* = 90 deg) > *α _{ds}* according to Fig. 1(a). This is supported by the visualization of the flow separation from VATT blades as depicted in Fig. 1(c), which is based on the work of Ouro and Stoesser [24]. In the early stages of the revolution, the blade operates at an effective angle of attack

*α*smaller than the dynamic stall angle and thus there is no flow separation. Once the blade surpasses

_{I}*α*, a leading edge vortex forms, here at

_{ds}*α*

_{II}. This large-scale vortex controls lift and drag forces of the blade and it grows in size until the maximum angle of attack

*α*

_{max}is attained. At

*α*

_{III}, the leading edge vortex separates and the turbine drops dramatically in efficiency [24] due to the sudden loss of lift and rapid increase in drag of the blade. What follows is that the accurate prediction of dynamic stall is critical in order to accurately simulate the torque generated by a VATT during optimal operational conditions [29].

*λ*, and the reduced pitching frequency,

*κ*. Laneville and Vittecoq [30] introduced an equivalent reduced rotational frequency to characterize VATs, here denoted as

*κ** and defined as

The parameter *κ** depends on the geometrical properties of the VAT (e.g., solidity *σ*) and its motion in terms of tip speed ratio. Figure 2 demonstrates how *κ** varies with the curvature parameter, *c*/*R*, and also with the rotor's solidity. According to Amet et al. [27], a value of *κ** = 0.05 can be adopted as threshold above which a pitching airfoil experiences strong dynamic stall. VATTs are often designed with high solidities and operate at lower tip speed ratios (commonly at *λ* < 3) than their wind counterparts vertical axis wind turbines, which are mostly designed with low solidities allowing to avoid deep dynamic stall [31]. As a result, VATTs suffer more from deep dynamic stall conditions and hence increased fatigue loads [26].

## Numerical Framework and Computational Setup

where *u _{i}* and

*x*(

_{i}*i*or

*j*= 1, 2, 3) are the fluid velocity and coordinates in the three coordinates of space, respectively,

*p*denotes pressure, and Re

*is the Reynolds number set as Re*

_{c}*=*

_{c}*cU*

_{0}/

*ν*, where

*U*

_{0}is the inlet velocity,

*ν*is the fluid kinematic viscosity, and

*c*is the airfoil's chord length. The freestream velocity and chord length are the velocity and length scales used for the normalization, i.e., equal to 1, while the Reynolds number in Eq. (5) is set according to the value from the experiments [7]. The subgrid scale stress tensor,

*τ*, is approximated using wall-adapting local eddy viscosity subgrid scale model from Nicoud and Ducros [34]. The source term

_{ij}*f*is the forcing term of the IB method, which is employed to resolve the moving boundaries in a fixed Eulerian field.

_{i}Hydro3D has been validated for various complex hydrodynamic flows such as in compound channels [35], or around hydraulic structures [36,37]. Recent implementations include a Lagrangian forcing IB method and a Lagrangian particle tracking algorithm to accomplish fluid–structure interaction [24,38] and bubbly flow simulations [39,40], respectively. Hydro3D is a finite differences-based Navier Stokes solver operating on locally refined, staggered Cartesian grids [41]. The fractional-step method [42] is used with low-storage three-steps Runge–Kutta predictor for time advancement [43]. A fifth-order weighted essentially nonoscillatory scheme is used to approximate convective terms and diffusive terms are approximated with central differences. The solution of a Poisson pressure-correction equation using a multigrid technique is adopted in the final step as a corrector. A refined version [44] of the direct forcing IB method of Uhlmann [20] is adopted to represent the airfoil geometry [24,38].

### Pitching Airfoil Kinematics.

The geometrical parameters and forces considered during the simulation of a pitching airfoil are sketched in Fig. 3. The pitch angle at the time *t* is calculated as $\alpha (t)=\alpha 0+\Delta \alpha \xb7\u2009sin((2\kappa U0/c)t)$, where *α*_{0} is the preset pitching angle, Δ*α* is the angle amplitude, and *κ* is the reduced pitching frequency. Note that upstroke movement, i.e., increase of pitching angle, is denoted with ↑ while ↓ indicates downstroke movement.

where *A* = *Hc* corresponds to the projected area considering the spanwise length (*H*) used in the computational domain, *ρ* is the fluid density, and *U*_{0} is the freestream velocity.

### Setup of Simulated Cases.

Two different cases of an airfoil describing a sinusoidal pitching motion undergoing deep dynamic stall conditions are simulated. Table 1 provides details for these two cases concerning: airfoil shape, Reynolds number (Re* _{c}*), pitching motion parameters, maximum (

*α*

_{max}) and minimum (

*α*

_{min}) pitch angle, and reference of the experiments used for the validations. The experimental work undertaken by Lee and Gerontakos [7] is selected for the baseline simulation, as it was reproduced using RANS [13,45,46], DES [13], and more recently with LES [15].

In the baseline case, a NACA 0012 is simulated and this is adopted for the mesh and time-step sensitivity study. The effect of blade cambering on dynamic stall due to airfoil pitching is analyzed using a NACA 4412 under the same flow and kinematic conditions. The simulations are run over four pitching cycles with the first cycle discarded from phase averaging of aerodynamic coefficients (similar to Ref. [15]), although during the experiments [7], averaging was performed for more than 100 cycles.

These experiments were carried out in a suction-type wind tunnel of dimensions *L* = 2.7 m (18*c*) × *H* = 1.2 m (8*c*) × *B* = 0.9 m (6*c*), and the numerical domain is identical to these dimensions as presented in Fig. 4. During previous numerical studies of this case [13,15,45,46], the domain length in the *y*-direction (*H*) was set to 20*c*, which reduced flow blockage in comparison to the experiment. The pitching center is at 0.25*c* from the leading edge of the airfoil and is situated 4*c* downstream of the inlet and centered in the transversal direction (see Fig. 4). In the experiment, the airfoil had a spanwise length of 2.5*c* and was equipped with end plates. The numerical domain extends 0.2*c* in the spanwise direction, which is adequate to reproduce accurately the three-dimensionality of the flow, in fact Lee and Gerontakos [7] stated that the flow can be assumed quasi-two-dimensional. This was also supported by the numerical results of Kim and Xie [15], who studied the influence of the spanwise length of the airfoil using values of 0.5*c* or 1.0*c* predicting very similar aerodynamic coefficients. Using finite spanwise domains together with periodic boundary conditions is a common practice in the study of airfoils using LES [47] and DNS [48].

Concerning boundary conditions, a uniform freestream velocity is set at the inlet while a convective boundary condition is set at the outlet. Periodic boundary conditions are set in the spanwise direction, and no-slip conditions are imposed on the upper and lower bounds of the numerical domain representing the wind tunnel walls. The computational domain is decomposed into 264 subdomains, as depicted in Fig. 4. Four levels of local mesh refinement (LMR) are used to achieve a very fine mesh resolution close to the airfoil while coarser grids are employed far away from the airfoil, as shown in Fig. 5. The simulations run on 76 central processing units (CPUs), using Supercomputing Wales facilities, and each simulation using the finest mesh resolution require approx. 38,000 CPU hours.

## Code Validation

The present computational approach is first validated with a mesh resolution and time-step sensitivity analysis using as reference the experimentally obtained aerodynamic coefficients of the NACA 0012 case. Table 2 presents the details of the three Eulerian fluid mesh resolutions of the finest LMR level, namely Δ*x*_{1}, Δ*x*_{2}, and Δ*x*_{3}, and number of Lagrangian IB markers along the airfoil's boundary, *N _{L}*. The mesh resolution is uniform in

*x*- and

*y*-directions so the number of points distributed along upper and lower surfaces of the airfoil is identical, whereas in the spanwise direction, the mesh resolution is Δ

*z*= 2Δ

*x*. Note that there is one Lagrangian marker per Eulerian cell in order to accomplish that the total force exchanged between solid and fluid grids is constant, as required by the direct forcing IB method using delta functions [20].

The mesh resolutions can be compared with simulations that employed body-fitted meshes in terms of the number of cells or solid markers *N _{L}* along the airfoil's surface. In RANS simulations, the finest resolution was used by Gharali and Johnson [46] with 500 cells covering the entire airfoil surface. Kim and Xie [15] performed most of their LESs using a mesh with 579 division (386 on the upper side), although they tested a finer mesh with 893 divisions, which did not provide noticeable improvements of the prediction of aerodynamic coefficients. Hence, meshes Δ

*x*

_{2}and Δ

*x*

_{3}with 642 and 802 divisions, respectively, are deemed to be fine enough to ensure good resolution of the flow over the airfoil's surface. The time-step sensitivity analysis is performed using the finest mesh resolution (Δ

*x*

_{3}) with three different values: $\Delta t1*=8\xd710\u22124,\u2009\Delta t2*=4\xd710\u22124$, and $\Delta t3*=2\xd710\u22124$, where $\Delta t*=\Delta tU0/c$ is the normalized time-step.

Figures 6 and 7 present the phase-averaged lift and drag coefficients computed from the simulations on the three different meshes and time steps, experimental data [7] and 3D DES results from Wang et al. [13]. The results obtained with the coarser mesh fail to predict *C _{L}* and

*C*along most of the pitching cycle, featuring a premature drop in the lift forces. From flow visualizations (not shown here), it is appreciated that this mesh is not fine enough to accurately resolve the velocity gradients on the airfoil's surface, which are fundamental to correctly reproduce flow separation from and reattachment on the airfoil. An increase in mesh resolution provides an appreciable improvement in the prediction of both aerodynamic coefficients. The medium mesh, Δ

_{D}*x*

_{2}, gives accurate

*C*predictions for

_{L}*α*< 20 deg with a very good agreement during the downstroke motion. Further improvements are achieved using the finest mesh, Δ

*x*

_{3}, as the

*C*prediction gets closer to the experimental data also during the upstroke motion. The simulation with Δ

_{L}*x*

_{2}predicts a premature shedding of the LEV, at

*α*≈ 18 deg ↑, while the simulation with Δ

*x*

_{3}predicts shedding of the LEV to occur at

*α*≈ 21 deg ↑, achieving a closer match to the experiment.

The present results deviate from the experimental data for *α* > 15 deg during the upstroke and the downstroke. This is probably because of slightly different approach flow conditions and/or flow blockage. Also data from only three cycles are used to average the aerodynamic coefficients; however, their variation during each cycle (not shown here) is quite small except during the part of the cycle between LEV shedding and flow reattachment. In fact, simulating a larger number of pitching cycles would only smoothen the phase-averaged curve rather than improve it, which is similar to what Kim and Xie [15] have found. The disagreement between experimentally obtained and computed coefficients may also be explained by the uncertainty in the experimental measurements, such as constant time delay on the pressure signals for the selected pitching frequency [7], pressure calculations during LEV separation as highlighted in Ref. [49], or the fact that the experimental coefficients do not account for skin friction although Kim and Xie [15] argued that this has little influence.

Nonetheless, the predicted aerodynamic coefficients agree quite well not only with the DES [13] for *α* < 20 deg, but also with the ones from other numerical works using 2D RANS [45], 3D RANS [13], and LES [15]. In addition, the LES predicts well the lift overshoot due to the LEV and poststall conditions, as explained later in Sec. 5, albeit in the LES this occurs at *α* = 20 deg while in the experiment, the maximum lift value takes place at *α* = 24 deg. Noteworthy is that the present and other numerical studies coincide in the fact that the shedding of the LEV occurs at *α* ≈ 20 deg ↑ (large drop in lift coefficient), all deviating from the experimental results where the LEV is shed at the end of the upstroke motion. The early part of the downstroke motion features large flow separation, consequently larger differences in the aerodynamic forces are observed in comparison to the upstroke before stall occurs. On the finest mesh, the coefficient of drag is predicted in good agreement with the experiments until *α* ≈ 15 deg. The overprediction of the drag coefficient for *α* > 15 deg compared to the experiment is similar to the LES results of Kim and Xie [15]. Overall, the predicted *C _{L}* and

*C*distributions of the present LES using the IB method agree well with those predicted by Kim and Xie [15] who used a body-fitted mesh.

_{D}The sensitivity of the simulation results to the time-step is studied on the finest mesh only (Δ*x*_{3}), and results of *C _{L}* and

*C*are presented in Fig. 7. The simulations using the largest time-step, $\Delta t1*$, predict an early shedding of the LEV and consequently the poststall flow is not accurately predicted. Decreasing the time-step improves the prediction of the formation of the LEV. Some differences in the predictions using $\Delta t2*$ and $\Delta t3*$ are observed during the downstroke motion, with the latter providing a better match with experiments at 5 deg ↓ <

_{D}*α*< 18 deg ↓. However, the differences between the simulations using $\Delta t2*$ and $\Delta t3*$ are very small. These results evidence that small time steps together with fine mesh resolutions are required to achieve a good representation of this complex flow. Considering that the total number of time steps to obtain one pitching motion cycle using $\Delta t3*$ is double than that of $\Delta t2*$, i.e., double the computational cost, mesh Δ

*x*

_{3}and fixed time-step $\Delta t2*$ are adopted as the best configuration for the following simulations.

## Results and Discussion

The effect of blade cambering on the development of dynamic stall is analyzed comparing the aerodynamic behavior of the previously validated NACA 0012 with that of the cambered NACA 4412 under the same flow and kinematic conditions as provided in Table 1.

### Flow Visualization.

The flow field generated over the NACA 0012 and NACA 4412 is visualized allowing qualitative comparisons of the aerodynamics of the two airfoils over the entire pitching motion cycle. This motion is divided into three stages: upstroke pitching from the minimum angle of attack until shedding of the LEV (Fig. 8), upstroke pitching under deep stall conditions (Fig. 9), and downstroke motion (Fig. 10). The flow is visualized using contours of instantaneous normalized *z*-vorticity (*ω _{z}c*/

*U*

_{0}) in a plane at half the spanwise extent of the domain, i.e.,

*z*/

*c*= 0.1.

#### Upstroke Pitching Previous to Deep Stall.

During the first instances of the upstroke motion, *α* > *α*_{min} = –5 deg (shown later in Fig. 10(e)), the entire suction side of the NACA 0012 features a laminar shear layer without flow separation, whereas a short laminar-to-turbulent transition along the trailing edge develops in the NACA 4412 due to its cambered shape. At pitch angle, 6.23 deg ↑ (Fig. 8(a)), the shear layer remains attached over most of the NACA 0012 airfoil with evidence of the onset of laminar-to-turbulent transition toward the trailing edge and some subsequent laminar vortex shedding close to the tail of the airfoil. The flow over the NACA 4412 starts to develop rear-to-front flow reversal progressively shortening the laminar shear layer on the suction side with generation of coherent vortices. As Fig. 8(b) shows, at 9.99 deg ↑, the flow reversal extends over more than half of the suction side of the NACA 4412 while the NACA 0012 also experiences flow reversal that breaks down the laminar shear layer due to an adverse pressure gradient induced by the pitching motion. This provokes the onset of Kelvin–Helmholtz instabilities with generation of shear layer vortices. Similar results have been reported based on experimental work by Carr et al. [6] and McAlister and Carr [50], and numerical simulations from Visbal [16].

Figure 8(c) visualizes the flow at 11.54 deg↑ just before the airfoil reaches the static stall angle of *α _{ss}* = 13 deg [7]. The flow over the suction side of either airfoil is very similar with a train of coherent vortices convected along the upper surface, and small recirculation bubbles forming at the leading edge of the airfoils and extending over

*x*/

*c*< 0.10. This turbulent bubble is the premature formation of the LEV that is characteristic of dynamic stall. Figure 8(d) exhibits the flow field at 16.22 deg ↑ where the LEV starts to develop and already extends over half of the NACA 0012's length and approximately over the first quarter of the NACA 4412. Noticeable differences in the flow field of the two airfoils are observed in Fig. 8(e) at 19.45 deg ↑ especially in terms of the LEV. For the symmetric airfoil, this large-scale flow structure is detached from the upper surface and there is a recirculating bubble that covers the first quarter of the airfoil. Meanwhile, the asymmetric airfoil permits the LEV to stay attached to the upper surface albeit it appears less coherent.

#### Upstroke Pitching at Onset of Stall and Poststall Conditions.

The increase of pitch angle induces the LEV to grow in size and to extend over almost the entire upper surface of both airfoils, as depicted in Fig. 9(a). At 21.17 deg ↑, the LEV is not completely detached from either airfoil and is accompanied by a recirculating area (enclosed bubble) extending over the first half of the suction side. At this stage, the straight airfoil overcomes the dynamic stall angle of $\alpha dss\u224820.8deg\u2191$ and hence the flow is considered to be in the poststall regime. In contrast, the NACA 4412 dynamic stall angle is $\alpha dsc\u224821.3\u2009deg\u2191$ so at *α* = 21.17 deg ↑, the LEV is still attached but close to be shed. The LEV formation and detachment process is critical in terms of forces on the airfoil and is responsible for the lift overshoot characteristic of dynamic stall and its shedding results in a dramatic drop of the aerodynamic coefficients, as shown later in Sec. 5.2.

During the remaining upstroke motion, both airfoils are under deep dynamic stall condition and thus a lack of lift generation capability. At 22.74 deg ↑, Fig. 9(b), the LEV moves away from both airfoils, and its clockwise rotation (negative vorticity) induces the generation of a TEV that features a counter clockwise rotation with positive vorticity. This primary TEV increases gradually and becomes the dominating large-scale flow structure at 23.83 deg ↑, Fig. 9(c). Before the maximum pitch angle is reached, the TEV is eventually shed allowing the enclosed recirculating bubble to extend over the upper surface of either airfoil, as observed in Fig. 9(d). Near to the completion of the upstroke motion, at *α* = 24.95 deg ↑, a secondary TEV is formed at the NACA 4412 trailing edge, while in the NACA 0012, this is not appreciated and the enclosed recirculating area dominates the airfoil suction side.

#### Downstroke Pitching.

Figure 10(a) shows the development of a secondary TEV on the NACA 0012 shortly after the downstroke movement starts at 24.52 deg ↓. This flow feature evidences the rapid formation and shedding of large-scale structures during poststall condition. In the cambered airfoil, the TEV is still attached. Once the airfoil continues to pitch down, the TEV is shed and the flow field becomes again dominated by the recirculating bubble followed by front-to-rear flow reattachment, and irregular shedding of LEVs and TEVs. At 13.90 deg ↓, Fig. 10(b), the shear layer developed from the leading edge extends similarly over the suction side of both airfoils. Figure 10(c) shows that reducing the pitch angle further results in flow reattachment until *x*/*c* ≈ 0.1. At 0.00 deg ↓, the NACA 4412 exhibits a laminar shear layer until *x*/*c* ≈ 0.85 while for the straight airfoil, it extends over the first half of the suction side. In the flow field at the minimum angle of attack (*α* = –5 deg ↓, Fig. 10(e)), turbulent flow phenomena are absent on the upper surface of either airfoil, while on their pressure side, the shear layer breaks down generating roll-up vortices that are eventually shed. These are more predominant in the NACA 4412 due to the convex shape of its pressure side, which shortens the laminar shear layer compared to the flow over the NACA 0012.

#### Flow Three-Dimensionality.

Turbulence structures and three-dimensionality of the flow over the pitching NACA 4412 airfoil are visualized in Figs. 11 and 12. Three-dimensional views of isosurfaces of the spanwise vorticity and contours of streamwise velocities at half spanwise domain width (*z*/*c* = 0.1) are depicted in Fig. 11. At 11.54 deg ↑, the free shear layer at the leading edge of the airfoil exhibits two-dimensionality and laminar separation until it becomes unstable and transitions to three-dimensionality in the form of rollers due to Kelvin–Helmholtz instability. Such vortices are fairly coherent in the spanwise direction initially due to the absence of turbulent instabilities. Nonetheless, spanwise instabilities emerge as the shear layer rollers exhibit some undulation in this direction. Figure 12 presents two isosurfaces of streamwise vorticity with opposite sign and identifies the onset of coherent periodic instabilities in the rollers that are close to the leading edge. A total of four instabilities are depicted along the spanwise domain (*H* = 0.2*c*) whose wavelength is *d _{w}* = 0.2

*c*/4 = 0.05

*c*being constant in the first three rollers depicted here. It is noteworthy that the spanwise wavelength remains constant irrespective of their size. Similar pattern of these perturbations was observed in the DNS of a static NACA 0012 by Jones et al. [51], and Visbal [16] identified an analogous onset of spanwise perturbations in their ILES of a SD7003 plunging airfoil with the size of the developed spanwise instabilities

*d*≈ 0.04

_{w}*c*, which agrees well with the present observations.

At 20.22 deg ↑, the large-scale LEV vortex dominates the flow over the airfoil's suction side as streamlines and isosurfaces plotted in Fig. 11(b) show. The *z*-vorticity isosurfaces suggest an almost instant transition from the 2D shear layer to 3D structures, which are predominant around the LEV. At this pitch angle, the flow over the upper surface is dominated by flow reversal whose interaction with the freestream velocity above results in strong velocity shear at the interface and causing strong turbulence.

During downward pitching, at 8.04 deg ↓ (Fig. 11(c)), the laminar shear layer extends further downstream along the upper surface. Full flow separation over the second half of the upper surface is observed during the entire pitch-down cycle. Shear layer turbulence in the form of fairly incoherent small scale structures is found toward the tail of the airfoil. Such complex front-to-rear reattachment is notably different from the smooth laminar-to-turbulent transition experienced during the pitch-up motion as shown in Fig. 11(a). These visualizations suggest that after the airfoil undergoes deep stall, with the separated flow over the upper surface being fully 3D, the process of flow relaminarization and hence the airfoil's ability to generate lift is delayed in comparison to the pitch-up process.

Flow phenomena such as shear layer transition, flow separation, or reattachment are key in the aerodynamics around pitching airfoils. Figure 13 shows various of these events developed at different stages of the pitching cycle using isosurfaces of *Q*-criterion = 300 [52] colored with streamwise velocities with top views of the airfoil from (*a*) to (*g*), whereas (h) shows a bottom view at *α* = –5 deg ↓. During upstroke motion previously to the generation of the LEV (*α* < *α _{ss}*), the turbulent structures above the upper surface of the airfoil feature a roller-like shape as shown in Figs. 13(a)–13(c). The onset of spanwise instability is again observed from the top view of Fig. 13(c). During the development of the LEV at 21.17 deg ↑, coherent isosurfaces of

*Q*-criterion are found until

*x*/

*c*= –0.1 corresponding to the stable free shear layer. Following a quick turbulent transition, the flow becomes unstable shortly after the shear layer breakdown and noncoherent small-scale structures are distributed over the LEV influence area.

The complexity of the front-to-rear flow reattachment during the downstroke cycle is depicted in Figs. 13(e) and 13(f) for which after *x*/*c* = 0.0 the shear layer breaks into 3D smaller scale structures that evidence the chaotic turbulent structures distribution in the area of full flow separation. Close to completion of the downstroke motion, at 0.00 deg ↓ the flow over the upper surface of the airfoil still features some separation due to the convex shape of the NACA 4412 s upper surface. Finally, the bottom view of Fig. 13(h) presents the flow development along the pressure side of the airfoil at *α*_{min}. Due to the cambered shape of the NACA 4412 there is a prompt flow separation on this side compared to that exhibited during pitch-up motion at a similar angle of attack.

### Aerodynamic Coefficients.

The aerodynamic loads for both airfoils are analyzed with the goal to link the instantaneous flow field to the generated forces. The effect of blade cambering on the magnitude and distribution of the aerodynamic coefficients is presented in Fig. 14. The plots in the left column present the coefficients as a function of the angle of attack over the entire pitching motion cycle and the plots in the right column plot the coefficients between 19 deg < *α* < *α*_{max}. In these figures, $\alpha dss$ and $\alpha dsc$ are the dynamic stall pitch angle for the straight and cambered airfoils, respectively.

Figure 14(a) quantifies the increase in lift when using a cambered airfoil and this is for both upstroke and downstroke motions. This is appreciated over the entire pitching cycle except when the airfoil stalls, i.e., *α* > *α _{ds}*. The NACA 4412 features a tighter hysteresis loop, i.e., the difference of

*C*values between upstroke and downstroke at the same angle of attack is smaller compared to that of the NACA 0012. These findings agree with Choudhry et al. [10] who stated that the hysteresis loops of cambered airfoils are smaller.

_{L}The onset of the LEV provokes a lift overshoot for both airfoils starting at *α* ≈ 17 deg until its shedding at *α* = *α _{ds}*. During this stage of high lift generation, the NACA 4412 generates greater lift compared to that of the NACA 0012. The difference in the generation of maximum lift between the two airfoils is due to the LEV remaining closer to the suction side of the NACA 4412 in comparison with a quick and significant detachment of the LEV from the NACA 0012 (as depicted in Fig. 8(e) with

*α*= 19.45 deg ↑). The fact that the LEV is further away from the airfoil's upper surface makes it more vulnerable to the freestream flow and hence is easier to be shed. As a result, the dynamic stall angle of the straight airfoil $\alpha dss$ is approx. 20.7 deg ↑, whereas for the cambered, it is $\alpha dsc\u224821.3\u2009deg\u2191$. Cambering the airfoil shape appears to delay the shedding of the LEV resulting in an extra lift overshoot. Maximum

_{ds}*C*values can be quantified with Fig. 14(b) showing that at

_{L}*α*the NACA 0012 has its peak at

_{ds,}*C*= 2.1 while the peak of the NACA 4412 is

_{L}*C*= 2.45, i.e., the cambered airfoil generates approx. 15% more maximum lift than the straight airfoil.

_{L}The poststall flow (*α* > *α _{ds}*) is characterized by the shedding of the LEV causing a dramatic drop in

*C*,

_{L}*C*, and

_{D}*C*. The coefficients exhibit very similar patterns for both airfoils during poststall conditions. This is in line with the findings from McCroskey et al. [53], who stated that under poststall conditions, the blades lose their aerodynamic capabilities, i.e., they behave like simple bluff bodies.

_{M}During the downstroke motion, the flow starts to develop front-to-rear reattachment from *α* ≈ 20 deg ↓ onward during which the coefficients stabilize. Flow recovery is characterized by the shear layer forming at the leading edge and expanding along the upper surface as the airfoil's pitch angle is decreasing. The cambered airfoil improves this boundary layer reattachment process, and hence increasing its capability to generate lift. This is related to the flow field observed in Figs. 10(b)–10(d), where the flow relaminarization on the upper surface of the NACA 4412 is taking place earlier than for the NACA 0012. Consequently, from *α* ≈ 10 deg ↓ onward, the NACA 4412 generates *C _{L}* ≈ 0.5, in contrast to the NACA 0012 that generates almost no lift. At negative pitch angles, the latter airfoil always produces negative lift with the minimum

*C*≈ –0.5, while the former crosses over from positive to negative

_{L}*C*at

_{L}*α*≈ –2.5 deg, and the minimum is at

*C*= –0.25.

_{L}The hysteresis loops of the coefficient of drag are plotted in Figs. 14(c) and 14(d). The NACA 4412 produces slightly greater drag for *α* < 10 deg, i.e., before overcoming the static stall angle, and the values of *C _{D}* are similar for both airfoils until the shedding of the LEV, i.e., when

*α*≈

*α*. The overshoot of lift forces is accompanied by a drag increase, and similar to the lift the cambered airfoil features a greater maximum drag. Under poststall conditions, both airfoils generate approximately the same amount of drag. During the pitch-down motion and after

_{ds}*α*< 20 deg ↓, the quicker boundary layer reattachment on the NACA 4412 leads to higher

*C*values, which are maintained until

_{D}*α*≈ 12 deg ↓ when both bodies generate similar values of

*C*.

_{D}The *C _{M}* hysteresis loops are depicted in Figs. 14(e) and 14(f). The cambered shape of the NACA 4412 leads to a larger pitching moment until poststall conditions. Figure 14(e) shows that once the angle of deep stall is attained, i.e.,

*α*=

*α*, the

_{ds}*C*slope increases notably due to the LEV action. Carr et al. [6] observed an analogous situation during their experimental work and denoted this stage as moment stall. During the downstroke motion, the NACA 0012 experiences a variation from negative to positive

_{M}*C*values, whereas the values of the NACA 4412 are always negative. The cambered airfoil always generates larger

_{M}*C*than the straight airfoil, i.e., the hysteresis cycle of moment (and also lift) are reduced with the cambered airfoil in comparison to the NACA 0012.

_{M}Results show that the NACA 4412 not only exhibits a larger peak value of *C _{L}* but also produces greater drag in comparison with the NACA 0012. The greater lift force of the NACA 4412 on the upstroke is an important finding in terms of designing Darrieus-type VATs due to the fact that they are driven by the lift force to generate torque and in particular VATs generate the majority of the power during their upstroke motion. However, the lift increase is accompanied by an increase in drag and this is detrimental to the turbine's performance, because it adds drag torque to the rotor shaft and reduces the turbine's ability to self-start. Therefore, it is important to quantify the lift-to-drag ratio (

*C*/

_{L}*C*) of the two airfoils, in particular at large angles of attack, i.e.,

_{D}*α*> 19 deg. The

*C*/

_{L}*C*curve is presented in Fig. 15. The effect of significant flow unsteadiness during LEV separation is noticeable for 19 deg <

_{D}*α*< 21 deg, which is where the curve is not perfectly smooth. It is apparent that the cambered NACA 4412 achieves a larger

*C*/

_{L}*C*ratio compared to that of the NACA 0012 over the selected angles of attack. The benefits of increasing the lift generation are greater than the drawbacks associated with the drag increase. Hence, the improved aerodynamics of cambered airfoils, such as a NACA 4412, when subjected to severe pitching motion, including dynamic stall, suggest that VATs designed with cambered airfoils perform better than VATs equipped with symmetric airfoils. This, of course, only applies to VATs that are operating at relatively low TSR, such as vertical axis water turbines, and for which dynamic stall is unavoidable.

_{D}## Conclusions

The results of large eddy simulations of the flow over two pitching airfoils undergoing deep dynamic stall have been presented. The computational approach to represent moving bodies in the flow is based on a refined immersed boundary method whose accuracy has been validated first in this study. The sensitivity to spatial and temporal resolution has been assessed in terms of the prediction of aerodynamic coefficients aided by comparisons with experimental data and numerical results of other studies. The main objective was to study the effect of airfoil cambering on the flow over and aerodynamic performance of a symmetric NACA 0012 airfoil and an asymmetric NACA 4412 airfoil under the same flow and kinematic conditions.

The simulation data provided quantitative (flow visualization) and qualitative (aerodynamic coefficients) evidence that the aerodynamic performance of the NACA 4412 is superior to the NACA 0012 during prestall conditions, generating more lift and featuring a smaller force hysteresis cycle than the straight airfoil. The dominating large-scale leading edge vortex produced a larger lift overshoot for the cambered airfoil accompanied by a slight delay in its shedding. Under poststall conditions, both airfoils behaved almost identically with successive generation and shedding of LEVs and TEVs that led to a rapidly varying distribution of aerodynamic forces at high angles of attack. The LES-predicted and visualized flow under deep stall for both airfoils agrees well with other experimental findings, outlining that during full flow separation the airfoils lose their aerodynamic capabilities, and thus the blade cambering effect of increased lift generation is neglected. During the front-to-rear flow reattachment of the flow during their downstroke motion, the NACA 4412 featured a quicker boundary layer reattachment together with constantly greater lift in comparison with the NACA 0012.

The results demonstrate that a cambered airfoil shape provides pitching airfoils a higher lift-to-drag ratio and a short delay in the shedding of the dynamic stall vortex, which is responsible for the lift overshoot during dynamic stall. These findings confirm experimental findings that cambered airfoils improve the performance of Darrieus-type vertical axis turbines that operate at low tip speed ratio due to: (a) favorable lift-to-drag ratio during the maximum power generation phase on the upstroke of the turbine and (b) a reduced force hysteresis loop, which diminishes cyclic fatigue loads on the turbine rotor thus improving their survivability.

## Acknowledgment

The authors acknowledge the support of the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via Welsh Government, where the simulations were performed. Information on the data underpinning the results presented here can found at Cardiff University data catalogue.^{2}

## Funding Data

Engineering and Physical Sciences Research Council (Grant No. EP/K502819/1).

## Nomenclature

*A*=airfoil's projected area

*c*=airfoil's chord length

*C*=_{D}drag coefficient $(=2D/\rho U02A)$

*C*=_{L}lift coefficient $(=2L/\rho U02A)$

*C*=_{M}moment coefficient $(=2M/\rho U02Ac)$

*D*=drag force

*H*=airfoil's spanwise length

*L*=lift force

*M*=moment around pitching center

*N*=_{b}vertical axis turbine's number of blades

*p*=pressure

*R*=vertical axis turbine's radius

*Re*=_{c}Reynolds number based on chord length (=

*cU*_{0}/*ν*)*u*,*v*,*w*=streamwise, transverse and spanwise velocity components

*U*_{0}=freestream velocity

*x*,*y*,*z*=coordinate system

*α*=effective angle of attack

*α*(*t*) =airfoil's pitch angle

*α*=_{ds}dynamic stall angle

*α*=_{ss}static stall angle

*α*_{max}=maximum effective angle of attack

*α*_{min}=minimum effective angle of attack

*κ*=reduced pitching frequency (= Ω

*c*/2*U*_{0})*κ** =reduced rotational frequency

*θ*=angle rotated by turbine blades

*λ*=tip speed ratio

*ν*=fluid kinematic viscosity

*σ*=vertical axis turbine's solidity, (=

*N*/2_{b}c*πR*)- Ω =
turbine rotor rotational speed