Aerodynamic Simulation of the MEXICO Rotor

CFD (Computational Fluid Dynamics) simulations are a very promising method for predicting the aerodynamic behavior of wind turbines in an inexpensive and accurate way. One of the major drawbacks of this method is the lack of validated models. As a consequence, the reliability of numerical results is often difficult to assess. The MEXICO project aimed at solving this problem by providing the project partners with high quality measurements of a 4.5 meters rotor diameter wind turbine operating under controlled conditions. The large measurement data-set allows the validation of all kind of aerodynamic models. This work summarizes our efforts for validating a CFD model based on the open source software OpenFoam. Both steady- state and time-accurate simulations have been performed with the Spalart-Allmaras turbulence model for several operating conditions. In this paper we will concentrate on axisymmetric inflow for 3 different wind speeds. The numerical results are compared with pressure distributions from several blade sections and PIV-flow data from the near wake region. In general, a reasonable agreement between measurements the and our simulations exists. Some discrepancies, which require further research, are also discussed.


Introduction
The unsteady nature of wind makes the flow over wind turbine blades extremely difficult to predict. Indeed, even with stationary inflow conditions, complex unsteady phenomena will occur, including flow separation, secondary flows and stall delay due to rotational augmentation [1]. CFD simulations are a very promising method for predicting accurately the flow over the blades. Once the flow characteristics are well predicted and understood, improved and efficient engineering models can be developed for avoiding the high computational cost of CFD. However, up to now only few CFD models of wind turbines have been validated in depth.
Detailed measurements of wind turbines operating under controlled conditions are required for performing extensive validations. The aim of the MEXICO project [2] was to provide the project partners with the afore mentioned data for validation purposes. The available measurements, which cover different operating conditions, have been used within the project MexNext IEA Task 29 by 20 research institutions from 11 different countries for validating numerical models and for the analysis of aerodynamic effects [3] [4].
In this work we summarize our efforts for validating a CFD model of the Mexico-turbine based on the open source toolbox OpenFoam [5]. We focus on axisymmetric inflow conditions under different wind speeds, keeping the rotational speed and pitch angle constant. The simulations correspond to a case with fully attached flow, a case at design conditions and a case in the post-stall region. The numerical results are compared with PIV flow data and pressure data from pressure transducers placed on several blade radial stations.

The MEXICO Experiment
The MEXICO project stands for Model Experiments in Controlled Conditions.
The measurement campaign took place in December 2006 and involved the extensive measuring of load, pressure and flow data from a 3 bladed wind turbine placed in the LLF (Large Low-Speed Facility) of the German-Dutch Wind tunnel DNW, which has an open section of 9.5x9.5 m 2 [6] [7]. The blockage ratio of the wind tunnel is 18 % and the use of breathing slots behind the collector helps to reduce tunnel effects. Prior studies by other project-partners show that tunnel effects do not influence dramatically the rotor flow [8].
The wind turbine has a rotor diameter of 4.5 m. The blades rotate in clockwise direction and are twisted and tapered. Their design is based on 3 different aerodynamic profiles: DU91-W2-250 from 0.20 to 0.45 span, RISO-A1-21 from 0.55 to 0.65 span and NACA 64-418 from 0.7 span to the tip region. A zig-zag tape placed at 5 % of the chord was used both on the suction and the pressure sides of the blades for triggering the laminar to turbulent flow transition. The material used for the blades is Aluminium 7075-T651 Alloy. The tower center is located 2.1 m downwind from the rotor, and its influence on the rotor flow is believed to be minimal.
The measurements were carried out under several flow conditions, including 2 different rotational speeds (324 rpm and 424 rpm), several wind speeds (ranging from 10 to 30 m/s), several pitch angles (from -5.3 deg to 1.7 deg) and several yaw-inflow angles (from 0 deg to 45 deg).
The pressure measurements were performed by means of Kulite transducers placed at 5 different blade sections corresponding to 25% , 35%, 60%, 82% and 92% blade span positions. Some pressure transducers located at the 0.25 blade span position were malfunctioning during most of the measurements, so they are disregarded in this work.
The PIV measurements were performed along multiple PIV windows upstream and downstream of the turbine. Two radial and two axial traverses were performed. Unfortunately, only outer radial positions can be analyzed with the available PIV data, since the innermost PIV-data are at the 45 % radial position. This is an important draw-back for the study of 3D effects, which play a major role in the inner region of the blades.
In this work three different wind speeds have been considered (10 m/s, 15 m/s and 24 m/s), while the rotational speed has been kept constant at 424.4 rpm. These conditions correspond to pre-stall, design conditions and post-stall respectively. The Reynolds number varies between 3.5×10 5 at the blade root to 7×10 5 at the blade tip. The tip speed is 100 m/s for all cases, resulting in a tip speed ratio of 6.7 at design conditions (15 m/s). The Mach number is well below 0.3 and compressible effects can be disregarded. The rotation of the blades can be simulated in different ways. In the present simulations the "frozen rotor" method was chosen for the steady state simulations. For the time-accurate simulations the method of the sliding meshes is used. In the frozen rotor approach the rotating flow field is solved in a non-inertial reference system, whereas the non-rotating parts are solved in an inertial reference system. The Coriolis and centrifugal forces are added to the momentum equation in the regions subjected to rotation. In this way, it is possible to simulate a rotating system without a moving mesh. This allows a considerable speed-up of the simulations. The relative position of the rotor and the stator remains constant for the whole simulation. This is an obvious disadvantage when a strong rotorstator interaction is expected. In the case of the MEXICO turbine, the interaction between the rotor (blades) and the stator (tower) is believed to be negligible. Consequently the tower is not modeled. One further limitation of this method is the fact that only steady-state (timeindependent) solutions are possible with it. This implies that the flow unsteadiness is not resolved in time. The use of unsteady simulations with sliding meshes is useful when the above mentioned limitations can influence negatively the solution. The solver used for the frozen rotor simulations is able of dealing with multiple reference frames (MRF) and uses the SIMPLE algorithm for pressure-velocity coupling. Hence its name MRFSimpleFoam. It is an incompressible, isotherm, turbulent and stationary flow solver.

Simulation model
The transient simulations make use of the so called General Grid Interface (GGI) between the rotor and the stator. This approach allows for a relative mesh motion between the rotor and the stator, avoiding at the same time mesh topology modifications. This is done by means of a special algorithm that allows the interpolation of the flow variables between the GGI interfaces. Both meshes (the moving and the static meshes) can be conformal or non-conformal. Since no tower was modelled for the MEXICO turbine, the interest of performing simulations with the sliding mesh approach was obviously not the study of any rotor-stator interaction. The reason for performing these simulations was rather to obtain a time-dependent solution. This should in principle allow a better prediction of unstationary effects. The use of GGI implies the calculation of weighting factors for every facet, originating from the geometric intersection of the faces of the rotor and the stator interfaces. This is done in order to ensure flux conservation between both parts. The price that has to be paid for the enhanced accuracy and more realistic flow prediction of the unsteady simulations with moving meshes is a much higher computational cost. The solver used for the transient simulations, named pimpleDyMFoam, is an incompressible, isotherm, turbulent solver which can deal with moving meshes. It makes use of the PIMPLE algorithm for the pressure-velocity coupling. The PIMPLE algorithm is a hybrid of the standard SIMPLE and PISO algorithms that allows the use large time steps.
The 2nd order linear-upwind discretization scheme has been used in this work for the convective terms of both the steady-state and the transient simulations. The time is discretized with a 2nd-order backward formulation in the case of the transient simulations.
The simulations, which are run fully turbulent, make use of the Spalart-Allmaras turbulence model [9].
The transient computations use the steady-state results as a pre-solution. For the case with an inflow wind speed of 10 m/s, no transient simulations have been run, since no significant flow separation or unsteady effects are expected. All the simulations have been run at the Flow-Cluster [10] of the University of Oldenburg using 240 processors. The steady-state computations reached convergence in less than 5 hours. For the transient simulations, however, more than 1 week was required for solving 4 rotor revolutions.
The wind speed at the inlet and the pressure at the outlet were set to Dirichlet boundary conditions. Three different inflow wind speeds were simulated (u= 10 m/s, 15 m/s and 24 m/s). The pressure at the outlet was always set to atmospheric pressure. No-slip boundary conditions were set for the blades and the nacelle. For the wind speed at the outlet and the pressure at the inlet Neumann conditions (zero gradient) were used. At the interfaces between the moving and the static part of the mesh, a special GGI boundary condition was used.

Mesh description
The same mesh was used for all the simulations. As explained above, only the blades and the nacelle were simulated.The high-computational cost of the transient simulations made it necessary to keep the number of grid cells as low as possible. On the other hand, the high requirements on the solution quality made it advisable to use a fine mesh. After several computations with different meshes, it was decided that the best compromise between computational cost and solution accuracy was offered by a mesh with 11×10 6 cells. The mesh was made with the tool snappyHexMesh, which is included in the OpenFoam package. In fact, the mesh consist of a non-moving mesh merged to a moving mesh. The domain is cylindrical, with a radius of 16 m and a length of 52 m. The rotating region has a radius of 6 m and a length of 14.5 m. Most of the cells are hexaedrons, although split hexaedrons are also present in the proximity of the boundary layer. The boundary layer of the blades was resolved with 10 prism-layers. For the nacelle, 3 prism-layers were used. Since wall functions were employed in the computations, an extremely fine resolution of the cells in the boundary layer could be avoided. The y+ value on the blades oscillated between 50 and 200. The wake and the regions where high gradients were expected were accordingly refined. In figure 1, showing a slice of the mesh at the inboard region of one blade, a refinement region in the vecinity of the blade can be seen.  However, for the wind speed 24 m/s, the computed pressure distributions (not shown in this paper) are not completely consistent with the measurements. The angle of attack rises with increasing wind speed owing to the fact that the rotational speed and pitch angle are kept constant for all cases (as it would occur in a stall-regulated wind turbine). At 15 m/s, corresponding to the design conditions, the flow is still attached, as it can be seen from the pressure distributions. Indeed, at 24 m/s the blade is massively stalled, with the suction side of the blades presenting no significant underpressure in large blade regions from the root to the tip. The Spalart Allmaras turbulence model is known to perform well with attached flows. However, when the flow is separated, as it is the case of the Mexico turbine at 24 m/s, the prediction capabilities of the mentioned turbulence model are limited.

Flow traverses
The flow upwind and downwind of the wind turbine has been analyzed with PIV measurements. An axial traverse and a radial traverse are compared with measurements. Both traverses were taken when blade 1 was pointing upwards (12 o'clock position). The PIV windows were placed in a horizontal plane at the 9 o'clock position.   Figure 3 shows the axial wind speed for the traverse corresponding to the 0.61 radial position (1.37 m away from the rotor axis). As it can be seen, a general good agreement between measurements and simulations exists for all wind speeds.
At 15 and 24 m/s the measurements show ripples in the axial induction. This is attributed to the presence of increased vorticity in that region as a result from flow separation on the blades. The simulation at 24 m/s predicts satisfactorily this effect. At 15 m/s, however, the numerical results do not present the measured oscillations in the axial induction. Furthermore, at 15 m/s the axial induction is underpredicted. This is certainly related to the fact that no flow separation was predicted by the simulations at design conditions (15 m/s). Apparently, at that wind speed the stall-onset begins in the experiment, whereas the simulations still predicts attached flow at those conditions. At 24 m/s, however, the simulation predicts massive stall, and the agreement with the experimental data is much better.
The radial traverse situated 1.53 m behind the rotor shows an interesting behavior. As it can be seen in figure 4 for the 24 m/s case, the measured axial speed oscillates strongly along the radial direction. At the innermost measured region, corresponding approximately to the radial position 1.2, the axial speed suddenly drops. The simulation results offer a satisfactory qualitative prediction of the flow pattern, although the measured velocity drop is less drastic in the numerical than in the experimental results. Unfortunately, the evolution of the axial speed at the inboard regions of the blade can not be analyzed from the experimental data, since only the central and outer regions of the blade were measured. According to the simulations, the  velocity recovers rapidly just after the sudden drop. This pattern is attributed to the presence of a strong vortex, which under normal conditions would not be expected at that radial position (approx. 60 % radius). Interestingly, this region corresponds to the transition zone between the RISO-A1-21 and the NACA 64-418 profiles. It seems obvious that both profiles present different aerodynamic characteristics with a poor compatibility between them.

Conclusions
The open source CFD toolbox OpenFoam has been validated against the MEXICO data-set. In general, a good agreement between numerical results and experimental data exists. The pressure distributions along 5 blade sections reveal more details of the rotor flow: the attached flow found at low wind speeds evolves into massively separated flow at high wind speeds due to the fact that the rotational speed of the rotor is kept constant. The numerical results are here in good agreement with the measurements for the cases with attached flow (10 and 15 m/s). When the blades are stalled (24 m/s) the limitations of the Spalart Allmaras turbulence model make the prediction of the pressure distributions less reliable. The PIV data provide an excellent insight into the rotor wake. The studied traverses suggest the presence of a strong vortex at the 0.6 r/R region of the blade. Its existence was confirmed by the numerical results. The release of the mentioned vortex, taking place at the transition between the RISO and the NACA profiles is attributed to the poor aerodynamic compatibility between both profiles. In a next stage of this work, new simulations with other turbulence models, including a laminar-turbulent transition model, are planned. Furthermore, a detailed study of 3D effects is also aimed.