Numerical Simulation of a Pitching Airfoil Under Dynamic Stall of Low Reynolds Number Flow

1.Ferdowsi University of Mashhad – Faculty of Engineering – Department of Mechanical Engineering – Mashhad – Iran. Correspondence author: javareshkian@um.ac.ir Received: Jan. 6, 2018 | Accepted: Jan. 23, 2019 Section Editor: Panagiotis Tsoutsanis ABSTRACT: In this research, viscous, unsteady and turbulent fl uid fl ow is simulated numerically around a pitching NACA0012 airfoil in the dynamic stall area. The Navier-Stokes equations are discretized based on the fi nite volume method and are solved by the PIMPLE algorithm in the open source software, namely OpenFOAM. The SST k – ω model is used as the turbulence model for Low Reynolds Number fl ows in the order of 105. A homogenous dynamic mesh is used to reduce cell skewness of mesh to prevent non-physical oscillations in aerodynamic forces unlike previous studies. In this paper, the effects of Reynolds number, reduced frequency, oscillation amplitude and airfoil thickness on aerodynamic force coeffi cients and dynamic stall delay are investigated. These parameters have a signifi cant impact on the maximum lift, drag, the ratio of aerodynamic forces and the location of dynamic stall. The most important parameters that affect the maximum lift to drag coeffi cient ratio and cause dynamic stall delaying are airfoil thickness and reduced frequency, respectively.


INTRODUCTION
By 2030, energy consumption will be increased more than two-third of present condition worldwide (Dorian et al. 2006). According to the announcement of the International Energy Agency, oil and gas sources will be depleted in the next 41.8 and 60.3 years, respectively (Hosseini et al. 2013). Moreover, excessive utilization of fossil fuel has signifi cant undesired impacts on global climate change (Farhad et al. 2008). Renewable energy is the best solution to overcome these problems (Mahmoudimehr et al. 2016). Wind energy is a renewable resource which has been improving considerably in recent decades (Dai et al. 2015).
Dynamic stall can be considered as a delay in separated fl ow over wings and airfoils due to rapid variation of the angle of attack in the unsteady motion (Carr 1988). Dynamic Stall Phenomenon (DS) can be simulated by various models. Many researchers (Merz et al. 2012;Holierhoek et al. 2013;Liu et al. 2014;Howison and Ekici 2014;Dyachuk et al. 2014 Dyachuk andGoude 2015a;2015b;Antonini et al. 2015;Zanon et al. 2015;Wang and Zhao 2015;Elgammi and Sant 2016) presented new models to simulate DS, some of them comparing their results with the famous Leishman-Beddoes model. Based on the turbulent nature of the fl ow, it is necessary to couple turbulence equations with basic ones. Simão Ferreira et al. (2010), Wang et al. (2012) Buchner et al. (2015) and Almohammadi et al. (2015) compared different types of turbulence models with each other. Wang et al. (2010) analyzed a pitching NACA0012 airfoil at Reynolds number (Re) of 1.35 × 10 5 . Their results demonstrated significant oscillations in instantaneous forces because they divided the solution domain into two zones, the dynamic and the static one. However, CFD results had an acceptable agreement with experimental results, except for conditions with a maximum angle of attack. The effects of accelerated flows were shown by Choudhry et al. (2013), Gharali and Johnson (2015) and Karbasian et al. (2016a). Some researchers (Mulleners and Raffel 2013;Melius et al. 2016) experimentally analyzed DS. Many researchers experimentally represented leading edge vortex (LEV) shedding. They realized that fluid behavior strongly depends on some parameters such as airfoil shape, mean angle of attack (α mean ), oscillation amplitude (d), reduced frequency (k) and specially Reynolds number or Mach number (Dickinson and Goetz 1993;Lehmann 2004). The effects of unsteady parameters like k, α mean , pitch axis location and Re on a NACA0012 airfoil were investigated by Akbari and Price (2003) for 3 × 10 3 ≤ Re ≤ 10 4 . Based on their results, the increment of reduced frequency causes a delay in flow separation. Also, moving the location of the pitch axis backwards, from the quarter-chord to the mid-chord, has a minor effect on the flow field in terms of the flow separation. Amiralaei et al. (2010) concluded that d, k and Re can change the maximum value of aerodynamic forces for 555 ≤ Re ≤ 5 × 10 3 . Moreover, it was shown that k and d do not have a significant impact on lift coefficient slope. Yu et al. (2010) studied unsteady aerodynamic characteristics of a S809 airfoil under sinusoidal pitch oscillation for a various range of k, α mean and d at Re = 10 6 . Their results indicated that by increasing the reduced frequency, maximum lift coefficient (C l,max ) occurs at a higher angle of attack. The effects of pitching location on the aerodynamic performance of a NACA0012 airfoil were described by Ferrari (2012). He showed that lift coefficient (C l ) enhances by an increment of pitching location. However, this increment in pitching location results in drag increase. Danao et al. (2012) conducted a 2D numerical simulation to analyze the influence of rotor blade thickness and its camber on the performance of a 5 KW vertical axis wind turbine (VAWT). They proved that for symmetric airfoils the thinner one has higher maximum pressure coefficient. Also, the airfoil with a small camber, like LS0421, generally has better performance, whereas NACA5522 with a 5% camber is not appropriate. Gharali and Johnson (2012) investigated the effect of an oscillating free stream on a stationary S809 airfoil. The simulation was performed for Re = 10 4 , 10 5 and 10 6 and a vast range of reduced frequencies from 0.026 to 18. They pointed out that for k = 0.002 and Re = 10 6 , all force coefficients of the stationary airfoil in an oscillating free stream have approximately the same trend as an oscillating airfoil in a constant free stream. Furthermore, they revealed that the behavior of aerodynamic forces, vortices and velocity field are very sensitive to the reduced frequency. By performing simulation on a NACA0012 airfoil in Re = 1.35 × 10 4 , Lu et al. (2013) revealed that the instantaneous lift coefficient, thrust generation and flow structure are affected by reduced frequency, pitching amplitude and airfoil curvature. It was shown by Raeisi and Alighanbari (2014) that trailing edge (TE) curvature and airfoil thickness have a major role in early flow separation of a pitching NACA23012 airfoil. Karbasian et al. (2016b) investigated the dynamic stall of a S809 airfoil for different reduced frequencies. The results showed that each increase in the reduced frequency rapidly enhances aerodynamic characteristics of the oscillating airfoil and also decreases the size of vortices formed around the airfoil. Using CFD methods and with help of Menter-SST turbulence model,  investigated the flow structure of a single blade VAWT under dynamic stall. They highlighted that the blade vortex seriously affects the flow behavior near the blade and leads to power production loss.  also simulated a pitching NACA0012 airfoil with help of k -ε URANS turbulence model to damp the turbulence production near the wall. Furthermore, they visualized vortices at several chord lengths downstream of the airfoil in order to discuss its flow characteristics. By performing a research study on a pitching NACA0018 airfoil, Hand et al. (2017) stated that increment of Re causes dynamic stall to be delayed at higher angles of attack. Geissler and van der Wall (2017) numerically focused on the control of the dynamic stall by airfoil deformation. They showed that not only strong vortex can be avoided by proper airfoil deformation, but also propulsion efficiency can be improved.
To summarize, the effects of unsteady parameters have been discussed by many researchers for different ranges of Re under DS. However, the authors could not address the parameters that have the most influence on the maximum lift to drag coefficient ratio and the dynamic stall delay. The main objective of this paper is to investigate the effects of k, d, Re and airfoil thickness on aerodynamic force coefficients and dynamic stall delay to specify the most effective ones on maximum lift to drag coefficient (C l/d,max ) and delay of dynamic stall for 1 × 10 5 ≤ Re ≤ 2 × 10 5 . Unsteady flow fields around airfoils in low Reynolds number are xx/xx 03/14 dominant with unsteady vortex shedding and viscous flow, which can lead to dynamic stall, leading edge vortices, wake capturing mechanism and clap and flying mechanism. This range of Reynolds number is applicable in the micro aerial vehicles and the wind turbines (Dickinson and Goetz 1993;Osborne 1951). To the authors' knowledge, not enough attention has been paid to modify non-physical oscillations in aerodynamic forces in the previous studies. Therefore, in this paper, a homogenous dynamic mesh is used to reduce cell skew factor (CSF) angle of mesh to prevent mentioned oscillations, which is validated with experimental data. The fluid flow is assumed turbulent, incompressible, and viscous. The Navier-Stokes equations are solved by using CFD techniques, and discretization is done based on the finite volume method (FVM) in the OpenFOAM software. The SST k -ω model is considered to simulate the turbulent flow.

DYNAMIC STALL PHENOMENON (DS)
This phenomenon results in an intensive increment of lift and delays the stall. DS is started with turning back of flow at TE. Then large annular flow appears in the boundary layer and reversed flow propagates up to LE. Suddenly, the reversed flow will be converted to the vortex at LE. By forming the vortical flow at LE, lift coefficient gradient increases rapidly due to the increment of the angle of attack until deep stall. Then abrupt drop is observed in Cl. Subsequently, the flow is attached to the airfoil LE and this process is repeated again in the next cycle (Ferrari 2012). Various steps of dynamic stall are demonstrated in Fig. 1 in more detail.

DYNAMIC MESH AND GOVERNING EQUATION
One procedure to simulate the oscillatory motion of foils is dynamic mesh. In this method, the mesh will be corrected by airfoil motion in each moment. The linear tension and bending spring technique is implemented in this approach to model motion alteration at the boundary. In other words, each line between two nodes in the mesh is converted to a spring, which stiffness is conversely proportional to the distance of nodes. Consequently, longer lines have less stiffness.
Mass, momentum and scalar equations in vector type and as independent of the coordinate system are as follows (Eq. 1): where ρ = density (kgm -3 ); φ = scalar; U g = mesh velocity; Γ φ = diffusion coefficient; Sg = source term. Equation 1 is discretized by FVM and is solved implicitly by a pressure-based algorithm. In this research study, discretization of the diffusion and convection terms is performed by central difference and second order linear upwind methods, respectively. Discretization of the unsteady term is conducted by the backward implicit scheme. The velocity and pressure equations are solved by PIMPLE method, which is a large time-step transient solver for incompressible flow and is a combination of SIMPLE and Piso algorithms. The SST k -w model is considered as turbulence model because, as Wang et al. (2010) revealed, this model showed an improvement compared to standard k -w model to follow experimental data of Wernert et al. (1996). Furthermore, it could properly predict the flow characteristics at very high angles of attack where the flow is detached. Additionally, Martinat et al. (2008) performed two and three dimensional simulations on a pitching NACA0012 airfoil under DS in Re = 10 5 and 10 6 in order to compare different turbulence models. The results indicated that the SST k -w model is the most accurate one to predict stall angle as well as lift and drag coefficients in upstroke motion. In this paper, pitching motion is applied by Eq. 2: where α = angle of attack. According to Eq. 2, the instant variation of the airfoil is obtained by angular velocity and oscillation amplitude. Airfoil movement is prescribed by a sinusoidal motion. Angular velocity and reduced frequency are achieved by Eqs. 3 and 4, respectively:

MESH AND BOUNDARY CONDITIONS
A C-type mesh is used to model the flow around the airfoil due to better performance than H and O-types for simulation of pitching motion (Esmaeilifar et al. 2017). In other words, the quality of other mesh types is lower than that of C-type (Ferrari 2012). The C-type mesh utilized in this study is illustrated in Fig. 2. (2) Constant velocity and constant pressure are considered as boundary conditions at the entrance and outlet, respectively (Samiee et al. 2018). In order to prevent any disturbances at the boundaries, the computational domain extended 20C from the airfoil in all directions (Bos et al. 2008;Lentink and Gerritsma 2003). The boundary conditions around the airfoil are illustrated in Fig. 3. Turbulent kinetic energy at the entrance boundary and initial condition is computed by Eqs. 5 and 6, respectively: e = (3(Ti × U ∞ ) 2 )/2 ω = √e/l where e = turbulent kinetic energy (m 2 /s 2 ); Ti = turbulent intensity (%); l = turbulence length scale. Turbulence intensity (Ti) and free stream velocity in Eq. 5 can be seen in Table 1. The obtained turbulence kinematic energy is used as the entrance boundary xx/xx 05/14 condition. In Eq. 6, omega will be quantified by turbulence kinematic energy (e) and the turbulence length scale (l), which is usually 0.01 times of basic length in external flows. Basic length is considered equal to the chord length (Versteeg and Malalasekera 2007). The height of the first cell on the airfoil surface is chosen in order of 10 -5 to make the y + criterion less than 1. Therefore, first cells on the airfoil are located in viscous sub-layer. y + is a normalized number defined as follows (Eqs. 7 and 8): where u* = friction velocity (ms -1 ); ϑ = kinematic viscosity; τ = shearing stress (N/m 2 ); w = wall. It should be mentioned that the shear stress in Eq. 8 is computed the same as laminar flow.

RESULTS AND DISCUSSIONS
Initial accurate value of pressure, velocity, kinetic energy and omega are needed at each cell of mesh to solve unsteady equations of continuity, momentum and turbulence. Steady results are used as the initial condition for the unsteady problem. Grid study is conducted to ensure the independency of the results from the mesh. Therefore, grids that consist of 1.0 × 10 5 , 1.5 × 10 5 , 2.0 × 10 5 and 2.5 × 10 5 cells are generated at Re = 1.35 × 10 5 and 10 degree angle of attack. As illustrated in Fig. 4, by raising the number of cells from 2 × 105 to 2.5 × 105, although computational cost increases, computations accuracy is almost constant. Therefore, the mesh with 2 × 10 5 cells is chosen for the rest of simulations. C l,max is selected to analyze the effects of different time steps for the mesh with 2.0 × 10 5 cells. According to Table 2  To validate the present simulation, obtained results are compared with the numerical results of Wang et al. (2010), Karbasian and Kim (2016) and the empirical data of Lee and Gerontakos (2004). The flow and geometric conditions of these simulations tabulated in Table 1. As demonstrated in Fig. 6, there are some oscillations in force coefficients in the results of Wang et al. (2010) and Karbasian and Kim (2016) in comparison with the empirical data. However, these oscillations are reduced significantly in the present paper. Cell skew factor is an important parameter in mesh quality and it remained 0.2% for points at LE and TE,   while it is almost zero for other points. According to Fig. 6, obtained results have an acceptable agreement with experimental data in upstroke motions. In downstroke motions, some oscillations around the experimental data exist due to the flow separation, reattachment and transition from laminar to turbulent. Mean content of lift coefficient is 0.7 and 0.65 for empirical results and present simulation, respectively. Also, drag coefficient (C d ) is equal to 0.2 and 0.19 in this research and empirical data, respectively. Based on the present results, the stall angle is 13° in the static condition. However, in the dynamic condition and according to Figs. 6a and 6b, even up to 16° there is not any reversed flow. Indeed, the first vortex will be created in α = 16.5° at LE. In other words, an increase in momentum of fluid, obtained from airfoil motion, helps the flow to undergo the intensive pressure gradient in a higher angle of attack without any flow separation.
In the lift coefficient hysteresis diagram (Figs. 6a and 6b) there is a sharp increase in Cl at α = 21° up to stall moment in upstroke motion. According to the comparison between α = 21° and 22° in upstroke motions, shown in Fig. 7, regions with rapid growth of vortex cause pressure reduction and finally lift coefficient enhancement. According to Figs. 6a and 6b, there is a tremendous drop in Cl at α = 22.6° in upstroke motion. The reason for this is that the flow escapes from the lower surface to upper surface at TE of the airfoil. This can be seen in Fig. 7 for AOA = 24 in the upstroke. Subsequently, there would be a velocity increment at lower surface of the airfoil, resulting in a sudden pressure reduction and also an intensive lift coefficient decrement. According to Fig. 8, at α = 24.7° in the downstroke motion, low-pressure region which is created at the lower region of the airfoil causes lift coefficient to be reduced. Conversely, a high-pressure region will be created at the lower surface of the airfoil after that.
Two types of flow rotation can be seen in each part of Fig. 7 in the upstroke motion. According to this figure, the generation of a clockwise vortex begins from LE. This vortex gets bigger by increasing the angle of attack until the airfoil approaching the end of the pitching cycle. In this point, the motion velocity of the airfoil reduces and the secondary flow is created. After that, pressure decreases at the lower surface of the airfoil. On the other hand, growth and separation of counter-clockwise vortex begins at TE. At this time, dynamic stall occurs. Then some oscillations are observed in the lift coefficient diagram because of flow separations until the airfoil approaches the maximum velocity. According to Fig. 5, oscillations of aerodynamic coefficients repeated similarly as the first cycle.

EFFECTS OF REYNOLDS NUMBER
Herein, the effects of various Re of 1 × 105, 1.35 × 105 and 2 × 105 on Cl are investigated. When the reduced frequency is fixed, the physics of flow is not changed. Therefore, Reynolds number is altered by dynamic viscosity. Simulation is done for k = 0.1, α mean = 10° and = 15° around the pivot point, which is located 0.25 times of NACA0012 airfoil chord from leading edge. C l,max and dynamic stall angle are increased from 2.15 to 2.40 and from 21.75° to 22.40°, respectively by an increment of Re from 1 × 105 to 1.35 × 105 (Fig. 9a). However, C d is increased from 0.8 to 1. Therefore, Cl/d,max is equal to 2.69 and 2.4 for AOA (deg) 3.2 3.txt Numerical (Wang et al., 2010) Experimental (Lee and Gerontakos, 2004) 2 Present study Numerical (Karbasian and Kim, 2016) Experimental (Lee and Gerontakos, 2004) Present study Experimental (Lee and Gerontakos, 2004)  Re of 1 × 105 and 1.35 × 10 5 , respectively. Moreover, the lift coefficient behavior at Re of 1 × 10 5 is not the same as the behavior of two other Re at downstroke motion from the angle of 20°. This is because of Re decrement, which is a factor for the flow separation.
There is no change in Cl by an increase in Re from 1.35 × 105 to 2 × 10 5 , but dynamic stall angle is increased from 22.40° to 22.80°. According to Fig. 9, there exists a discrepancy between aerodynamic coefficients in the same angle of attack. This difference can be justified by strong diffusion of vortices during circulation around the airfoil, which is called aerodynamic phase delay. Gordon (2006) explained this unsteady impact in more details.

EFFECTS OF REDUCED FREQUENCY
The effect of reduced frequency on the unsteady pitching motion is investigated by lift and drag coefficient diagrams for k = 0.1, 0.15, 0.2, Re = 1.35 × 10 5 , α mean = 10° and d = 15° around the pivot point of 0.25 chord length (C). As demonstrated in Fig. 10, C l,max is not changed by reduced frequency increment from 0.1 to 0.15, but there is a 15% increment, from 0.98 to 1.14, in C d . Reduced frequency enhancement from 0.1 to 0.15 causes 2.1° delay in DS. C l/d,max is equal to 1.98 in this condition, indicating a 17% reduction compared to k = 0.1. By comparing the two mentioned conditions, it is revealed that in the downstroke motion and from the angle of attack of 25° to 18°, there are more reductions in instantaneous Cl in k = 0.15 because in this reduced frequency the flow becomes more turbulent and, as a result, the vortex will be stronger and finally Cl is reduced more. By comparing two reduced frequencies of 0.15 and 0.2, there is no significant difference in the instantaneous lift coefficient when the airfoil is in the upstroke motion unless in the stall point. It is worth mentioning that C d in k = 0.2 has a significant reduction (22%) in comparison with k = 0.15 in the downstroke motion (Fig. 10b). According to the diagrams of aerodynamic coefficients, in a low angle of attack (before deep stall), aerodynamic coefficients have a similar behavior. In addition, by increasing the reduced frequency from 0.1 to 0.2 C l,max remains constant. However, the dynamic stall delayed as a result of increasing the reduced frequency. xx/xx 10/14

EFFECTS OF AIRFOIL THICKNESS
Three airfoils of NACA family, 0018, 0015 and 0012, are selected to examine the effects of airfoil thickness on the aerodynamic coefficients. Simulations are accomplished for α mea n = 10°, k = 0.1 and Re = 1.35 × 10 5 . Pitching location is assumed 0.25 times of chord length. As demonstrated in Fig. 11a, by increasing the airfoil thickness in the upstroke motion, there is no significant difference in the behavior of instantaneous lift coefficient before deep stall area. In addition, the flow pattern around mentioned airfoils is similar with each other. By entering to deep stall region, the trend of C l is changed significantly so that by an increase in thickness from 0012 to 0015 there is a reduction in C l,max from 2.4 to 1.86, while stall point is preserved. However, by increasing thickness from 0015 to 0018, C l,max and stall angle are increased from 1.86 to 1.96 and from 22° to 23.5°, respectively. According to Fig. 11a, in downstroke motion, there is a reduction in instantaneous lift coefficient by an increase in thickness. This reduction is significant in the stall area and while the flow is separated but it will disappear by the transition of flow to the fully attached condition. By analysis of instantaneous drag coefficient in Fig. 11b, although C l,max has been decreased by thickness increment, the maximum drag coefficient is also reduced. Maximum drag coefficient is equal to 0.99, 0.56 and 0.59 for NACA0012, 0015 and 0018 airfoils, respectively. C l/d,ma x is equal to 2.42, 3.32 and 3.32 for NACA 0012, 0015 and 0018 airfoils, respectively. Consequently, thickness increment in deep dynamic stall (DDS) region results in enhancement of C l/d,ma x and stall angle. That is to say, increased airfoil thickness has favorable effects on aerodynamics of pitching airfoils. This is because drag coefficient reduces much more than lift coefficient, resulting in an increase in maximum lift to drag coefficient.

EFFECTS OF OSCILLATION AMPLITUDE
The effects of airfoil oscillation amplitude are studied for d = 10°, 13° and 15°, k = 0.1, α mean = 10° and Re = 1.35 × 10 5 on a NACA0012 airfoil. As represented in Fig. 12a, results show that C l,max is increased with the increment of oscillation amplitude. Although there is no significant enhancement in C l,max among oscillation amplitudes of 13° and 15°, dynamic stall angle is increased from 21° to 22.4°. A huge drop in Cl is reported from α = 20° onwards, which is called the Deep Dynamic Stall (DDS). DDS is a flow regime in the dynamic stall. Vortex shedding phenomenon is the predominant characteristic of this flow regime and by increasing the angle of attack, lift coefficient drops intensively after it reached to C l,max (McCroskey 1981). This process can be changed by variation of Re, thickness and reduced frequency. However, in a turbulent flow with the range of Re in this paper, stall is called the Light Dynamic Stall (LDS) between α = 15° and 20°. C l/d,max is equal to 2.53, 2.82 and 2.69 for d = 10°, 13° and 15°, respectively. Moreover, augmentation of oscillation amplitude causes widening of hysteresis loops. It should be noted that lift coefficient enhancement in a LDS till DDS is significant at d = 13°. This result is represented for C d in Fig. 12b.

CONCLUSION
In this research, two dimensional, viscous, turbulent and unsteady flow around a pitching NACA0012 airfoil was investigated under dynamic stall. Fluid flow was analyzed for 1 × 10 5 ≤ Re ≤ 2 × 10 5 . Effective unsteady parameters such as reduced frequency, oscillation amplitude, Reynolds number and airfoil thickness and their impacts on the dynamic stall delay, flow pattern and instantaneous lift coefficient diagram were studied. The main findings of the present study can be summarized as follows: Modification of mesh quality by CSF can prevent non-physical oscillations in the aerodynamic coefficients for the pitching motions.
Aerodynamic force coefficients have a slightly different behavior at Re of 1 × 10 5 rather than two other Re at downstroke motion from the angle of 20°. This is because of Re decrement, which is a factor for the flow separation.
In low angles of attack aerodynamic coefficients have a similar behavior. By increasing the reduced frequency from 0.1 to 0.2 C l,max remains constant. However, the dynamic stall delayed as a result of increasing the reduced frequency. Consequently, reduced frequency is the most effective parameter on delaying of dynamic stall.
Airfoil thickness is the most effective parameter on the maximum lift to drag coefficient ratio. By increasing the thickness from 0.12 chord to 0.18 chord, C l/d,max increases by 37%.