Numerical investigation of the flow evolution of an oscillating foil with the different reduced frequency

The power extraction systems based on oscillating foil are played more and more attention in alternative energy extraction systems. The dynamic stall problem on oscillating foil is related to general performance of energy extraction systems, where a varying flow filed is introduced by a rapid change in the effective angle of attack. The objective of this paper is to investigate the flow evolution of an NACA0012 foil undergoing sinusoidal pitching motion with the different reduced frequency. The k-ω SST turbulence model, coupled with a two-equation γ-Reθ transition model, is used for the turbulence closure. The results showed that the numerical method can well capture the dynamic stall process. With the increasing of the reduced frequency, the dynamic stall is delayed during the development of the Leading Edge Vortex stage. During post-stall vortex shedding stage, the mechanism of post-stall vortex shedding present different style at different reduced frequency K. When K=0.1, the secondary vortex forms and covers the whole suction face again along with the shedding of LEV. Compared to the case of K=0.1, no secondary vortex was induced for K=0.05. However, When K=0.2, the secondary vortex covers the whole suction surface again because the oscillating velocity is too fast to develop the secondary vortex.


Introduction
With the high cost of non-renewable fossil fuels, the development of renewable energy sources have been promoted [1][2] [3] . Among the various renewable energy, wind energy has attracted a lot of attention because it is resource-rich, widely distributed, clean. Additionally, the wind turbines that can convert the wind energy into electricity have also achieved great development [4][5] [6] . Among the wind turbine system, the wind turbine blades often play an important role. However, the dynamic stall problem on the wind turbine blades is related to general performance of the wind turbine, where a varying flow filed is introduced by a rapid change in the effective angle of attack [7][8] [9] . Hence, it is of significant importance to understand the dynamic interactions between the transient wing motion and unsteady flow conditions. Moreover, as one of the primary flow energy conversion technologies, an improved understanding will help the oscillating wing system to enhance the energy harvest performance for wind or hydro energy.
Many experimental and numerical investigations on aerodynamic performance of an oscillating foil have been conducted. Yang [10] utilized a videometric technique to measure the deformation of large-scale wind turbine blade in the full-scale tests and during operation. The results indicate the deformation of composite wind turbine blades would affect the performance of the wind turbine.
McKinney [11] applied a pitching wing to investigate a windmill to extract wind energy. The results presented the windmill can extract energy compared with rotary designs. LEE and Gerontakos [12] investigated the transient behaviour of the blade surface unsteady boundary layer and the characteristics of the dynamic stall events associated with an oscillating NACA0012 airfoil with multiple hot-film sensor arrays and smoke flow visualizations measurement. The results showed the unsteady flow mechanisms and the influence of the reduced frequency, oscillation amplitude on dynamic characteristics.
With the increasing of the computing resources, much attention has been paid to the computational fluid dynamics technique. Liu [13] studied the dynamic stall characteristics of an S809 airfoil undergoing translational motion as well as pitching motion with the Computational Fluid Dynamics (CFD) techniques. Simulation results indicated that both the out-of-plane and in-plane translational motions of the airfoil affect the unsteady aerodynamic forces significantly. Danao [14] investigated the effects of various rotor blade thickness and camber on the performance of a vertical axis wind turbine, where the NACA0012, NACA0022, NACA5522 and LS0421 are employed. The results showed slightly cambered foils can improve the overall performance of the vertical axis wind turbine because of generating higher values of torque in both the upwind and the downwind regions. Huang [7][15] [16] numerically investigated a hydrofoil undergoing transient up-and-down pitching motion at various pitching velocities. It revealed that the flow experienced complex boundary layer events such as laminar-to-turbulence transition and the pitching velocity had a significant effect on the transition and hydrodynamic load.
Although much work has been done to investigate the dynamic performance of oscillating wings [17][18] [19] , the significant uncertainty of the impact of the unsteady flow features on the power generation and also energy harvest system still exists. The objective of this study is to analyze the influence of the reduced frequency on the transient vortex structures around an oscillating foil, and to further improve the understanding of the interplay between unsteady flow, oscillatory motion of the foil and the dynamic performance. The numerical models are presented in Section 2, followed by the summary of numerical set-up in Section 3. In Section 4, detailed analysis of the influence on the flow evolution is presented. Finally, the conclusions are summarized in Section 5.

Basic governing equations
The flow field of oscillating foil is simulated by solving the 2D unsteady Reynolds-averaged Navier-Stokes equations, with the continuity and momentum equation listed below: where, ρ is the fluid density, t is time, u is the velocity, x is the coordinate, p is the pressure,  is the fluid viscosity, and subscripts i and j denote the directions of the Cartesian coordinates.

Turbulence model
The current simulation solves the URANS equations using the revised k-ω SST turbulence model, which couples the k-ω SST turbulence model proposed by Menter [20] and the -Re transition model [21][22] [23] .
where, k is the turbulent kinetic energy, ω is the specific turbulent dissipation, , Pk and Pω are production terms, Dk is the destruction term, F1 is blending functions.

Numerical setup
As shown in figure 1(a), the computational domain and boundary conditions are set according to the experimental setup [12], containing two domains connected with a sliding grid interface. In which, the static domain has a length of 18c and a height of 6c, and the dynamic domain has a diameter of 4c. The interaction between the static domain and the dynamic domain is controlled by CFX expression language (CEL) subroutine. The velocity of the inlet and the pressure of outlet, and non-slip wall were adopted as boundary conditions. The inlet is sufficiently far (7 times chord length) from the airfoil to simulate an effectively unbounded flow field. Figure 1(b) show the boundary layer mesh and the refined grids around the foil, respectively. About 450 nodes were placed around the airfoil boundary layer, which is selected to meet y + =yuτ/ν≈1,where y is the thickness of the first cell from the foil surface, and uτ is the wall frictional velocity [20] . The 2D fluid mesh is composed of 2.2*10 5 elements. The mean free stream velocity U∞ is 14m/s with turbulence intensity of 0.08%, which corresponds to a Reynolds number of Re=1.35*10 5 based on chord length of the airfoil. The evolution of the lift coefficients Cl (Cl =L/(0.5ρU∞ 2 sc)) and drag coefficients Cd (Cd =D/(0.5ρU∞ 2 sc)) with the geometric angle of attack α for the sinusoidal pitching motion are presented in figure 3, where L and D are, respectively, the lift and drag. In Figure 3, the dynamic experimental results are obtained from Ref. [12]. It is shown that in the upstroke from t1(α + =-5°) to t7(α + =25°), the numerical results agree well with the experimental data, and the maximum of predicted lift coefficient Cl is very close to the measured value. In the down stroke, the predicted lift coefficients present a small amplitude and lowfrequency oscillating behaviour from t7(α + =25°) to t13(α -=18.1°). This may due to the interaction between the induced secondary vortex on the leading and trailing edge, and the detailed analysis of the flow features is given in the following paragraphs.
To sum up, although the numerical results exist slightly discrepancies, two-dimensional studies by k-w SST model and the -Re transition model are able to capture the mean features and qualitative tendencies of the dynamic stall phenomenon. It is prominent that the difference of the mechanism of the LEV development at different reduced frequency. Figure 5 shows comparisons of the predicted pressure coefficients at same geometric angles of attack for the different reduced frequency.

Comparisons between
As K increased, the formation and rearward convection of the LEV was delayed to higher angles of attack, as shown in figure 4(a-c). At K=0.1, the formation and rearward convection of the LEV occurs at α + = 17.5°. The LEV has covered the whole suction surface at α + =23.5°. As the development of the LEV, the max value of pressure coefficient moves the trailing edge in figure 5(a-c). When K=0.05, the LEV has covered two thirds of the suction surface at α + = 17.5°. As the angle of attack reaches 20.4°, the LEV has shed away. Compared to the case of K=0.1, the pressure coefficient is larger. At K=0.2, the formation of the LEV occurs at α + = 20.4°. Compared to the case of K=0.1, the LEV only has covered the half suction surface at α + =23.5°. As shown in figure 5, the pressure coefficient is larger than the case of K=0.1.
Although the formation and rearward convection of the LEV was delayed with the increasing reduced frequency, the corresponding lift and drag coefficients all increased almost linearly with the increasing angle of attack α. Because the pressure coefficient become more and more larger with the development of the LEV.

4.2.2.
Post-stall vortex shedding stage. As shown in figure 6, the contours of z-vorticity superimposed on the instantaneous streamlines in post-stall vortex shedding stage at the same angle of attack for different reduced frequency are presented. It is find that the vortex-shedding phenomenon exists discrepancy at different reduced frequency. Comparisons of the predicted pressure coefficients at same geometric angles of attack for the different reduced frequency are shown in figure 7. When K=0.1, it should be noted that at α + = 24.8° the LEV has shed away because of the interaction between the LEV and the TEV. Meanwhile the leading edge of the suction surface has induced a secondary vortex. In the downstroke phase at α -=24.7°, the secondary vortex has covered the whole suction surface, as shown in figure 4. With the development of the secondary vortex, the pressure coefficients become larger, as shown in figure 7(a-b). It should be noted that at α -= 22.8° a new TEV rolls up, grows larger and squeezes the secondary vortex, while the secondary vortex begins to shed away, which is responsible for the reduction in the slope of the lift coefficient curve. As α increased, the multiple interaction between the LEV and the TEV results in small amplitude oscillating behavior in the lift and drag coefficient curves. When the angle of attack reaches to α -=18.1°, the flow begins to transit from turbulent to laminar.
When K=0.05, the vortex-shedding phenomenon was different from as K=0.1. At α + = 24.8° in upstroke phase, the LEV has shed away. As α increases, it is obvious observe that no secondary vortex was induced. Compared to the case of K=0.1, the pressure and lift coefficients are lower, as shown in figure 3 and figure 5. Increasingly, when K=0.2, the vortex-shedding phenomenon was different of the cases of K=0.05, 0.1. At α -= 24.7° in downstroke phase, the LEV has shed away and the secondary vortex has spanned over approximately one third of the suction surface, while the trailing edge vortex is induced. Compared to the case of K=0.1, the secondary vortex begins to shed away at α -= 22.8°.