Numerical simulation of the flow around oscillating wind turbine airfoils Part 1 : Forced oscillating airfoil

This is the first part of a two part paper on flow around vibrating wind turbine airfoils. In this part 1, the unsteady, incompressible, viscous and laminar flow over a forced oscillating airfoil is computed using a method based on a commercial Computational Fluid Dynamics (CFD) code. Beforehand, the Navier-Stokes equations are solved for the unsteady flow around a NACA 0012 airfoil at a fixed 20° incidence and the low Reynolds numbers of 103 and 104 to check the reliability of the CFD computations. Then the flow around a pitching airfoil is simulated for prescribed values of the reduced frequency. The Navier-Stokes equations are expressed in ALE formulation and solved with moving mesh. The effects of the discretization scheme and the moving mesh technique are investigated. The hysteresis loops of the dynamic stall phenomenon are highlighted.


INTRODUCTION
Wind turbines are subjected to a hard environment as the atmospheric turbulence, the ground boundary layer, the rapid variations in wind speed and direction and the tower shadow for downwind turbine.With the architecture of the rotor, the time delay of the yawing system and the resulting yawing error, the blades are experiencing cyclic variations of the attack angle which lead to 3D unsteady aerodynamics and dynamic stall [1,2].The dynamic stall results in fluctuating blade force which leads to blade pitch oscillations, flap-wise or lead-lag oscillations, known as aeroelastic phenomena.Then problems of aeroelastic stability can be encountered in particular on the new large wind turbine blades.Indeed, the new blade lengths are more and more important and as they are fixed at only one extremity, they are subject to high bending and torsion stress.Aeroelastic stability problems are also encountered on parked wind turbine at high wind speeds [3,4].
Most of the wind turbine aeroelastic analyses were performed with the use of engineering methods where the blade forces were often computed by the BEM theory and the dynamic stall was modelled with empirical models such as the ONERA or the Beddoes-Leishman.The dynamic response of the wind turbine was then determined using aeroelastic tools with a reference rotor speed which was assumed constant [5][6][7].The problems of viscous fluid flow and of elastic body deformation were then studied separately.These weak coupled approaches are limited to small deformations.For the large displacements, the interaction between those two media and the modelling of the unsteady aerodynamics and dynamic stall are of great importance for the aeroelastic stability study.Then there is a need to perform aeroelastic computations by means of a strong coupled approach with Computational Fluid Dynamics (CFD) simulations of the viscous fluid flow.
To the authors' knowledge, the first coupled analysis of the aeroelastic stability of a wind turbine was accomplished in the frame of the European project KNOW-BLADE at RISOE: Ellipsys3D, an in-house 3D Navier-Stokes equations solver, was coupled to a 3D structural model based on the beam element theory.Then the problem of lead-lag oscillation of a wind turbine blade was investigated [8].Others Fluid Structure Interaction studies can be found in the literature.However, most of them are intended for high-speed turbo-machinery and the applied CFD codes solve the compressible Navier-Stokes equations.
In this paper, the problem of classical Flutter phenomenon is dealt with a strong coupled method where the dynamic response of the wind turbine blade is determined in time accurate sequences.The numerical simulations are performed using a commercial CFD code which is coupled to a computational structural code for the solution of the dynamic equations of the blade.The aim of this study is to propose a straightforward technique to be used with general computer tools.
All computations are performed for a 2D airfoil with moving mesh.Beforehand, to check the reliability of the CFD computations, the cases of the unsteady flow field around (i) a NACA 0012 airfoil at a fixed large angle of attack and (ii) around this airfoil in forced pitch oscillations are considered.The effects of the discretization schemes and the moving mesh technique are investigated.The hystereris loops of the dynamic stall are highlighted with the moving mesh computations.This is the purpose of this part 1 of the present paper.The problems of flow induced vibrations of a symmetrical NACA 0012 airfoil in pitch oscillations and that of a cambered NACA 632 415 airfoil in flutter are considered in the part 2 of the present paper.

NUMERICAL APPROACH
It is generally assumed that the flow around wind turbines is turbulent and incompressible.The Mach number is lower than 0.2 and the Reynolds number lies between 10 5 and 2 × 10 6 .However, according to Leishman [1], although the fluid velocities are relatively low compared to a helicopter rotor, the justification of an incompressible flow require also that the frequency of the sources of the unsteady effects must be small compared to the sonic velocity.In this paper, the computations are performed for low values of the reduced frequency to justify the assumption of an incompressible flow.In addition, the computations are implemented with the low Reynolds numbers 10 3 and 10 4 , and it is assumed that the flow is laminar.These simulations of a laminar flow allow us to perform direct numerical simulation (DNS): the equations are not averaged as with the RANS methods.The computations of the fluid flow around the oscillating airfoil are performed with a 2D moving mesh.

FLUID EQUATIONS
The governing equations are described in the Arbitrary Lagrangian Eulerian (ALE) coordinates where the grid is treated as a referential frame moving with an arbitrary velocity.A detailed description of the derivation of the ALE formulation can be found in Refs.[9,10] among others.
Let Ω f ∈ R 2 the spatial domain occupied by the fluid.The incompressible Navier-Stokes equations defined in ALE coordinates are written as, in tensor notation: (1) where t, is the time, x i the Cartesian coordinate of a point of Ω f , u i , the absolute velocity component in the direction i, u ci , the velocity of the moving grid, p, the pressure, ρ, the fluid density, the determinant of the metric tensor and τ ij , the stress tensor components defined as, in the case of laminar flows, (2) with µ, the molecular dynamic fluid viscosity.An additional equation called the Space -Conservation Law (SCL) is enforced for the moving coordinate velocity components.
The boundary conditions are prescribed as follows.On the West boundary, the inflow condition is applied having defined the free-stream velocity.The outflow condition is set on the East boundary.A symmetry condition is applied on the South and the North boundaries.At the airfoil surface, the fluid velocity is equal to the airfoil speed: where u is the fluid velocity adjacent to the airfoil surface and v b is the velocity of the airfoil.
Block topology applied for the oscillating airfoil.

COMPUTATION OF THE AERODYNAMIC FORCES
The resultant of the aerodynamic forces is obtained from the solution of the Navier-Stokes equations by integrating the pressure p and the shear stress τ w over the blade surface S: where: • is the outward-pointing vector normal to the airfoil surface • N is the number of cells on the airfoil surface • p i , τ wi and δS i are the wall cell pressure, shear stress and face area respectively • is the outward-pointing vector normal to δS i F X , the drag force and F Y , the lift force are then computed as: ( and are the unit vectors in the parallel and vertical directions of the flow, respectively.The lift and drag coefficients C L and C D are derived from the computed values of the aerodynamic force components and with: (6) A ref is the reference area and V ∞ is the reference velocity set to the free stream velocity.

DYNAMIC EQUATIONS OF THE AIRFOIL
The instantaneous angle of attack of the NACA 0012 airfoil in forced pitch oscillations is given by the relation: where α 0 is the mean angle of attack of the oscillating airfoil, α m is the oscillation amplitude and ω = 2π f with f, the frequency of oscillation.

MOVING MESH TECHNIQUE
An overview of the computational domain is shown in Fig. 1.The axis of rotation of the airfoil is located at 25% chord length from the leading edge.The grid is an O-H block structured mesh, allowing a meshing technique easy to implement.Moreover, it was shown in [11] that solutions obtained with structured meshes are of higher accuracy.To check the meshing technique implemented for moving the grid, results are compared to that obtained with the Arbitrary Sliding Interface (ASI) method which is built-in the CFD code.Both moving grid techniques are described in the following: (i) With the ASI technique, the computational grid is split into a moving part around the airfoil and a fixed part.The cells of the circular sub-domain of radius R 1 (Fig. 1) move at the airfoil speed.The mesh movement is defined explicitly by specifying time-varying position for all the moving mesh block vertices.An interface boundary surrounding the moving mesh part slides at the specified velocity.(ii) In the meshing technique, the coordinates of the moving cells are expressed in cylindrical coordinates (R i , θ i ).As with the ASI technique, the grid of the circular sub-domain of radius R 1 moves with the airfoil as it oscillates.An algebraic interpolation is applied for moving the vertices of the annular subdomain delimited by the radius R 2 and R 1 .The grid of the outer domain of radius R i ≥ R 2 remains stationary.At each time step t n , the angular position of the vertices is then given by the relation: (8) where: (9a) (9c) ∆t is the time step and is the airfoil angular velocity: (10) With this moving mesh algorithm, the mesh distortions are small and the original mesh quality is preserved.Given an appropriate choice of the radius R 1 and R 2 , this technique can be applied even for large displacements of the airfoil.The method is implemented in a User subroutine called by the CFD code at each time step.All simulations are carried out with StarCD, a commercial computer code based on the solution of the Navier-Stokes equations by the finite volume method.The PISO algorithm is applied for the solution of the coupled pressure-velocity equations.This algorithm is similar to SIMPLE but has more than one pressure corrector step.In the original model, there are two correction levels [12].In StarCD, the maximal number of PISO correctors is 20 but when the solution converges, only 3 or 4 corrector steps are applied.The UPWIND scheme, a first order accurate differencing scheme, is applied for the interpolation of the convective terms.The QUICK scheme, a high order differencing scheme would be more accurate.However, stability problems have been encountered as it is shown on Figs. 2 and 3 which compare the time histories of the lift and drag coefficients computed with both differencing schemes for a stationary airfoil at 20°angle of incidence.The lift and drag coefficients calculated with the QUICK scheme vary with large amplitudes and the amplitude variation of the force coefficients increases as the Reynolds number increases.These results are due to the instability of the QUICK scheme when the Peclet number is too high.As for the simulations of the flow around the oscillating airfoil it is important to distinguish the numerical instabilities of the solution from the unsteady flow associated with the airfoil movements, the computations are performed with the first order UPWIND scheme, less accurate, but more stable.For the temporal discretization of the equations, the implicit θscheme is applied with a factor equal to 0.8.

APPLICATION
These simulations are performed for a NACA 0012 airfoil which is used on vertical axis wind turbine.First computations are carried out for the simulation of the unsteady flow around a stationary NACA 0012 airfoil at the large angle of attack α = 20°.Then the case of this airfoil in forced pitch oscillations is considered.All computations are performed with an airfoil of 1m chord length.The properties of the fluid are those of air with a density ρ = 1.2 kg/m 3 and a dynamic viscosity coefficient µ = 1.8 Pas.

FLOW AROUND A FIXED NACA 0012 AIRFOIL AT α = 20°T
hese calculations are carried out with a constant integration time step ∆t =10 -2 s for the chord Reynolds numbers Re = U ∞ C/ν = 10 3 and 10 4 where U ∞ is the free-stream velocity, C, the airfoil chord and υ = µրρ, the kinematic viscosity.The computational domain extends on a distance equal to -8.1 × C at upstream and the outlet boundary is located at a distance equal to 15 × C. The North and South boundaries are located at ± 10 × C respectively.The grid is a C-H block structured mesh and consists of 66 600 cells with 320 cells around the airfoil.The height of the first cell row around the airfoil is y 0 /C ≈ 3.5 × 10 -3 (Fig. 4).It is assumed that the solution has converged when the temporal variations of the force coefficients are periodic.The force coefficients are then obtained by calculating the arithmetic average of the instantaneous values corresponding to the last iterations of the periodic solution.
Fig. 5 shows the instantaneous velocity contours obtained with Re = 10 3 at the dimensionless times τ * = U ∞ t/C = 30 and 45.It can be seen that the flow on the upper surface of the airfoil separates near x/C = 0.5 and that the shape of the separation zone changes with time.A regular vortex shedding appears which is attenuated in the wake.When Re = 10 4 , the flow around the airfoil and in the wake is more complex as it is shown on Fig. 6   represents the instantaneous contours of the velocity magnitude at the dimensionless times τ * = U ∞ t/C = 30 and 67.5.The shape and the changes of the separation zone are different to those obtained with Re = 10 3 and the wake vortices are not dampened.This reflects that the oscillations of the flow are of larger scale than with Re = 10 3 .
The computed lift and drag coefficients are summarized in Table 1 compared to the values calculated by Hoarau and al [13] and to data measured by Sunada et al [14] for Re = 4 × 10 3 , Laitone [15] for Re = 2.07 × 10 4 and Critzos et al [16] for Re = 5.0 × 10 5 : (i) A mismatch is found for the lift coefficient computed with Re = 10 3 and the computed values of Ref. [13].Note that a closely related value was obtained in our computations with the QUICK scheme.However, as the QUICK scheme is unstable, this result is uncertain.(ii) The force coefficients calculated here with Re = 10 4 are closer to the values found in [13].(iii) The computed C L and C D coefficients are larger to experimental data of Refs.[14], [15] and [16].
The low agreement of the computational results with measurements is attributed to uncertainties in the experimental data.Indeed, although the NACA 0012 airfoil was the subject of numerous experiments [15][16][17][18][19] and numerical simulations [13,19,20], there are very few data for low Reynolds and Mach numbers.Laitone [15] explained this gap by the insufficient precision (or sensitivity) of the strain gauges.Moreover, as shown in Table 1, data found in the literature are somewhat different.In addition, various values of the static stall angle of attack and the corresponding C L, max are given in Refs.[15-17, 19, 20].Nevertheless, it is generally found that the lift coefficient increases (and the drag coefficient decreases) as the Reynolds number increases even at the Reynolds numbers in the range 10 3 ≤ Re ≤ 10 5 .These suggest that our results are more reliable than those of Ref. [13].

NACA 0012 IN FORCED PITCH OSCILLATIONS
The computational domain extends are defined as previously.An overview of the block topology is shown on Fig. 1.The O-H mesh used for these computations consist of 111 000 cells with 500 cells around the airfoil and 720 cells at the periphery of the sub-domain of radius R 1 (Fig. 7).
The calculations are performed for the chord Reynolds numbers Re = 10 4 with a constant integration time step ∆t = 10 -4 s.The instantaneous value of the angle of attack is given by the equation (7) with α 0 = 0°, α m = 10°and ω corresponding to k* = 0.19 and 0.45 where k* = ω C/2U∞ is the dimensionless reduced frequency based on half chord and the free-stream velocity.At time t = 0 s, the airfoil incidence is zero.Preliminary computations are carried out for a fixed airfoil with the initial angle of attack to get the initial flow field.The simulations for the airfoil in forced oscillations are initiated when the solution becomes stable.The

Comparison of the two moving mesh techniques
Figs. 8 and 9 show that the contours of velocity magnitude obtained with both moving mesh techniques, the ASI method and the meshing method are similar, as well as over the airfoil (Fig. 8) than in the wake (Fig. 9).The time-wise variations of the calculated lift coefficient with both moving mesh techniques are compared on Fig. 10.This figure shows again that the obtained results with both moving mesh techniques are similar and confirms the suitability of the meshing technique.

Instantaneous velocity magnitude
On Fig. 9, it is shown that the wake vortices follow the airfoil oscillations.The velocity contours around the airfoil at selected times within the oscillating cycle, are depicted on Figs.11 and 12 for the reduced frequency k* = 0.19 and 0.45 respectively.Similar flow patterns are found for same positions of the airfoil when it is upstroke (or down stroke) cycle (Figs.11a and 11f or Figs.12a and 12e).At opposite positions, the flow patterns are similar but reversed (Figs.11b and 11d, 11c and 11e or Figs.12b and 12d).With the larger reduced frequency, the near wake vortex shedding is more important.
On Fig. 13 are depicted the time histories of the velocity magnitude and pressure at a point located near the airfoil surface, in the vicinity of the leading edge.The amplitude of the velocity magnitude variation is lower when the reduced frequency decreases.On the opposite, the time-wise variations of the pressure are similar however a higher peak is get with the lower reduced frequency.

Force and moment coefficients
The time evolutions of the force and moment coefficients are depicted on Fig. 14, compared to the oscillations of the airfoil.The lift and moment coefficients have the same frequency (T CL ) as that of the airfoil oscillation (T): T CL = T = 111.70s when k* = 0.19 and T CL = 46.54s with k* = 0.45.The variation of the lift and moment are in phase with the airfoil oscillations when k* = 0.45.The drag's frequency (T CD ) is one half the lift's frequency (T CD = T CL /2).Similar behaviour was found in [21] with Euler computations of the flow over an oscillating airfoil.In Brevins [22, p.47], it is reported that this has been shown in experiments and that it is a consequence of the geometry of the vortex street.
The hysteresis loops C L (α) and C D (α) are depicted on Fig. 15.When the airfoil is oscillating, both lift and drag coefficients increase compared to the 2D steady state: (i) The mean values of the drag coefficient corresponding to α = 0°are C D = 0.0775 and 0.0684 with k*=0.19 and 0.45 respectively.In previous computations of the drag coefficient for the fixed airfoil at zero incidence, it was found that C D = 0.0374.(ii) The lift coefficients vary roughly between ± 0.37 and ± 0.12 respectively.The lift hysteresis loops are then larger when the reduced frequency is lower.However, the predicted maximum lift and drag coefficients increase with the reduced frequency.

SUMMARY AND CONCLUSION
Simulations of the fluid flow around an airfoil in forced pitch oscillations have been performed with a commercial CFD code based on the solution of the Navier-Stokes equations by the finite volume method.The governing fluid equations are described in ALE coordinates and solved with a moving mesh technique that is implemented in a user subroutine called by the CFD code at each time step.The meshing technique is straightforward and suited to structured grids.The hysteresis loops of the dynamic stall phenomenon are highlighted with the moving mesh computations.It is found that the period of the drag coefficient C D is twice the period of the lift coefficient C L .Both coefficients increase with the unsteady motion in comparison to the 2D steady state.As expected, the hysteresis loops and time evolutions of the lift and drag coefficients vary with the reduced frequency k*.

370 Numerical simulation of the flow around oscillating wind turbine airfoils 1 Figure 2
Figure 2 Time histories of the lift and drag coefficients computed for a NACA 0012 airfoil at a fixed incidence of 20°and Re = 10 3 .

Figure 3
Figure 3 Time histories of the lift and drag coefficients computed for a NACA 0012 airfoil at a fixed incidence of 20°and Re = 10 4 computed with both differencing schemes.

Figure 4
Figure 4 Mesh near the airfoil at fixed incidence of 20°.

Figure 5
Figure 5 Instantaneous velocity contours around the airfoil at a fixed incidence of 20°with Re = 10 3 at the dimensionless times t * = 20 and t * = 30.

5 Figure 6
Figure 6 Instantaneous velocity contours around the airfoil at a fixed incidence of 20°with Re = 10 4 .

Figure 7
Figure 7 Mesh near the oscillating airfoil.

Figure 8
Figure 8 Contours of velocity magnitudes over the oscillating airfoil computed with both moving mesh techniques and k* = 0.45.

igure 9
Contours of velocity magnitudes in the whole computational domain computed with both moving mesh techniques and k* = 0.45.

378 Numerical simulation of the flow around oscillating wind turbine airfoils 1 Fig 10
Fig 10 Time evolution of the lift coefficient computed with both moving mesh method and k* = 0.45.

Figure 13
Figure13 Time histories of the velocity magnitude and pressure at a point located near the airfoil surface, in the vicinity of the leading edge.

382 Numerical simulation of the flow around oscillating wind turbine airfoils 1 Figure 14
Figure 14 Time evolution of the force and moment coefficients of the oscillating airfoil with a) k* = 0.19 and b) k* = 0.45.

Figure 15
Figure15 The hysteresis loops for the lift, drag and moment coefficients for the oscillating airfoil with k* = 0.19 and 0.45.

Table 1 .
Lift and drag coefficients computed for the NACA 0012 airfoil at 20°fixed incidence convergence of the solution for the flow field around the oscillating airfoil is checked both with the number of PISO correctors and the maximal value of the CUNO parameter, a number analogous to the CFL number but whose calculation is based on the face fluxes.It is recommended to check that the peak values of the CUNO do not exceed about 300 (http://www.adapco-online.com/).However, this criterion proved to be insufficient and values lower than 20 have been set.All results are depicted as a function of the dimensionless time t* = t/T where T = 2π/ω is the period of oscillation.