Influence of Turbulence Model for Wind Turbine Simulation in Low Reynolds Number

In designing a wind turbine, the validation of the mathematical model’s result is normally carried out by comparison with wind tunnel experiment data. However, the Reynolds number of the wind tunnel experiment is low, and the flow does not match fully developed turbulence on the leading edge of a wind turbine blade. Therefore, the transition area from laminar to turbulent flow becomes wide under these conditions, and the separation point is difficult to predict using turbulence models. The prediction precision decreases dramatically when working with tip speed ratios less than the maximum power point. This study carries out a steadiness calculation with turbulence model and an unsteadiness calculation with laminar model for a three-blade horizontal axis wind turbine. The validation of the calculations is performed by comparing with experimental results. The power coefficients calculated without turbulence models are in agreement with the experimental data for a tip speed ratio greater than 5.


INTRODUCTION
Computational fluid dynamics (CFD) modeling and experiments both have advantages and disadvantages.If both can be done in conjunction, we can expect more effective understanding of the phenomenon.Although CFD is more advantageous for predictions where experiments are difficult to carry out, e.g.free stream conditions, it is generally difficult to obtain reliable results compared with experimental data.However, it is possible to obtain useful CFD results with verification by experimental results.Moreover, experiments cannot deliver correct results for any arbitrary condition due to limitations of experimental equipment, measurement errors and problems with measurement systems.CFD is an efficient tool for turbo-machinery and can complement uncertain experimental results [1].However, CFD simulation generally takes a long time for practical designs.It is necessary to reduce the calculation time for many design conditions.
The Reynolds number of general wind tunnel tests is lower than the one available for turbulence models.Hence CFD with and without a turbulence model have been attempted to solve for wind turbine performance and provide data for various detailed characteristics.
In this paper, it is attempted to solve for more accurate characteristics of a wind turbine in a short time possible even on a personal computer, using coarse grid.The reliability of the experimental results and the CFD results are discussed.Furthermore flow visualization is carried out to obtain more detailed information concerning flow around the blade.The results analyzed by 3dimensional CFD simulation are investigated to determine the factors contributing to wind turbine characteristics in detail.
The CFD code is an in-house incompressible finite volume Navier-Stokes solver which was developed by the author.The solver is based on structured grids and the use of curve-linear boundary fitted coordinates.The SIMPLE algorithm is used for pressure-velocity coupling.The convection term is calculated using the QUICK scheme and the other terms in space are calculated using 2nd order difference schemes.It is well known that sophisticated turbulence models do not always produce better results than the very simple models.Therefore the proven and computationally efficient Launder-Sharma low-Reynolds-number k-ε turbulence model is used in this report.Detailed characteristics are solved for by 3-dimensional CFD.
The most important points in this research are to make clear the difference in results between the calculations and the experiments, and to make the calculation results contribute to the improvement of experimental method.

NUMERICAL METHODS
The in-house code used is a self-developed incompressible finite volume Navier-Stokes solver.The solver is based on structured grids and the use of curvelinear boundary fitted coordinates.The grid arrangement is collocated (Perić et al., 1988) [2] and the Rhie and Chow interpolation method (1983) [3] is used.The SIMPLE algorithm (Patankar, 1980) [4] is used for pressure-velocity coupling.The convection term is calculated using the QUICK scheme (Leonard, 1979) [5] and the other terms in space are calculated using 2nd order difference schemes.It is well known that sophisticated turbulence models do not always produce better results than the very simple models.For practical applications it is often wiser to use a simple approach than the computationally expensive.Therefore the proven and computationally efficient Launder-Sharma low-Reynolds-number k-ε turbulence model (1974) [6] is used in this report.In this paper the calculations with the turbulence model are only steady, and the calculations with the laminar flow model are only unsteady.

Experimental Condition
Figure 1 shows the experimental apparatus.The experiments of the wind turbine were carried out by Mie University [7].A three bladed wind turbine is situated in front of the wind turbine.The experimental wind turbine has a diameter of 2.5m and the blades consist of four airfoil profiles as shown Table 1.The pitch angle of the tip is θ tip =-2°.The internal diameter of the wind tunnel is 3.6m.The experiment was conducted at a wind velocity of 7m/s, and the measured data were the wind velocity, the number of rotation, the torque, and the thrust.
Figure 2 shows the relation of the fluid force acting on the blade element of the wind turbine at radius, r; the angle of pitch,θ; the angle of attack, a; the lift, L; the drag, D; the tangential force, F t ; the axial force, F a ; i.e. thrust force, the wind velocity, V a ; tip speed, U=rΩ; the rotational speed, Ω; and the relative velocity, W.
Relation among the lift, L, the drag, D, the tangential force, F t , and the axial force, F a , in Figure 2 are written by: sin( ) cos( ) cos( ) sin( ) Under no stall conditions, i.e. a small angle of attack, Eq. ( 1) is approximated as follows: The axial force, F a , i.e. thrust force, is predicted with the same accuracy as lift.On the other hand, since the tangential force, F t , i.e., torque, is strongly influenced by drag, D, and it serves as the difference of the force by the lift and the drag, the produced force becomes small.For this reason, the predicted accuracy of torque is less than that of the thrust force.
The Reynolds number Re=V a R/ν is expressed by the turbine radius, R, the wind velocity, V a , the kinematic viscosity of air, ν, and Re=6×10 5 .The characteristics of the wind turbine are expressed by the tip speed ratio, λ, the power coefficient, C P , and the thrust coefficient, C a .

Computational Grid
Figure 3 shows the computational grid around the wind turbine rotor.Periodic boundary conditions are applied to both side surfaces and calculation is performed at the region around one blade, that is, one third of a sphere domain.The radius of this sphere is twenty times the rotor radius.The O-O type grid enables a suitable grid arrangement, being able to arrange many grid points along the wing surface without distributing many points to unnecessary parts.The number of grid points is 130 around the configuration, 63 points spanwise, 52 points normal to the surface direction, and the 425,880 points in total.The grid is generated using an algebraic grid generation method (Eriksson, 1982) [8] based on the transfinite interpolation method which gives 5×10 -6 in a direction normal to the near-wall grid spacing to unit rotor radius, and y + values of less than 1.0.

GENERAL PERFORMANCE
In general, the thrust force is nearly equal to the lift, while the tangential force is strongly influenced by drag.
Because it becomes the difference of L sin(a+θ) and D cos(a+θ) from Eq. ( 1), it turns into a small force of less than 10 percent of the thrust force.For this reason, it is easy not only in numerical computation but also in an experiment, for large errors to occur in the tangential force result.For the experiment accurate measurement of the thrust is very difficult, because the experimental apparatus sets up the force transducer on the base under the tower of the wind turbine model.The force transducer takes the gravitational force and the moment becomes far larger than that of the thrust of the turbine, and the large capacity transducer is selected.The measurement accuracy is thus worsened.The experimental results show the following characteristics: the maximum power coefficient C P =0.426 appears at a tip speed ratio of λ=5.25; the stall region appears below λ=4.2; the power coefficient decreases above λ=5.5, because the angle of attack becomes smaller as the tip speed ratio increases.The big difference between the k-ε turbulence model and the experimental results is produced near just after stall angle where the tip speed ratios are λ=4∼2.This is because the present turbulence models cannot fully predict the transition from laminar to turbulent flow.In this respect, we anticipate the development of a more accurate turbulence model in the future.Therefore we attempted to obtain data using the laminar model.In the k-ε turbulence model, the leading edge separation occurs completely in the region where the tip speed ration is less than λ=3.For data obtained at a tip speed greater than λ=5 in which the angle of attack becomes small so as not to stall, the power coefficients can be fully predicted by the laminar model and are in good agreement with experimental results.
The computational results of the thrust coefficient agree well with the tendency of the experimental results for all regions, displaying an upward shift of about 0.1.
Figure 6 shows the convergence of the unsteady calculation using the laminar model.The results show the tendency for calculations to take a long time as tip speed ratio increases.

DETAIL CHARACTERISTICS
It is shown in the previous section that the performance of a wind turbine is predicted with sufficient accuracy in spite of having few grid points.Detailed characteristics, such as the spanwise distribution of characteristics, the flow visualization of blade, and the pressure distribution around the blade are shown in this section.

Characteristics of Spanwise Distribution
Figures 7 and 8 show the spanwise local power coefficients for λ=5 and λ=6.The integration of the local power coefficient along the radius presents the power coefficient of Eq. ( 3).The results of the laminar model are lower than the k-ε turbulence model over the entire span.The results of the laminar model for λ=5 in Fig. 7 show a sudden drop from r/R=0.86 to 0.90.The most different result between the laminar and the k-ε turbulence model appears at r/R=0.90.For λ=6 in Fig. 8 the most different result between the laminar and the k-ε turbulence model appears at r/R=0.80.The reasons for this are described in the later subsection.

Flow around the Blade
Flow visualization of the calculation results is carried out to grasp more detailed information in this section.Figures 13 and 14 show the flows around the blade analyzed using the k-ε turbulence model.The pressure side experiences a separation bubble from r/R=0.6 to 0.82 for λ=5.The suction sides for λ=5 and 6 experi-ence trailing edge separation near the hub.
Figures 15 and 16 show the flows around the blade analyzed using the laminar model.The pressure side experiences a separation bubble at r/R<0.82 for λ=5 and all over for λ=6.The suction side experiences leading edge separation at r/R>0.85 for λ=5, and at 0.2 chord length around r/R=0.65, leading edge separation at r/R<0.2.The suction side for λ=6 experiences separation at the leading edge for r/R<0.35and at 0.

Pressure Distribution on the Blade
This subsection discusses the connection between the pressure distribution on the blade and the separation that is affected by having and not having a turbulence model, while making it correspond to the local power and thrust coefficients, and the flow around the blade.The influences are investigated at a tip speed ratio higher than λ=5, where the results of the laminar model are in agreement with the experimental results.
Figures 17 and 18 show the pressure distribution on the blade for λ=5 and 6.The k-ε turbulence model shows the states in which the separation is not produced; therefore the pressure-drop near the leading edge and the pressure recovery near the trailing edge are sufficiently represented.In contrast, it is shown that the power and thrust coefficients for the laminar model become lower than the k-ε turbulence model by the separation, in spite of the thickness of the separation being thin.
The pressure coefficient on the pressure side be-tween the k-ε turbulence model and the laminar model do not appear remarkably different, although the separation exists.
A remarkable difference appears on the suction side.Figure 17(a) shows the pressure coefficients at r/R=0.80.The flow of the laminar model detached at 0.2 chords and the pressure drop at the leading edge and the pressure recovery near the trailing edge become less sufficient than the k-ε turbulence model, that is, without separation.Figure 17(b) shows the pressure coefficients of r/R=0.91.The flow of the laminar model causes the detachment at the leading edge, and there is little pressure drops in the neighborhood of leading edge.The high pressure is maintained from the leading edge to 0.5 chords, and the pressure recovery from 0.5 chords to the trailing edge is not done enough.The lift is obtained to integrate pressure along the blade surface.It can be understood from the pressure distribution that the difference of the thrust coefficient in Figs.9−10 becomes small regardless of whether or not the turbulence model is used.In contrast, a remarkable difference is seen in the power coefficient in Fig. 6, because the pressure drop at the leading edge and the pressure recovery at the trailing edge have a big influence on the rotational force.Since the flow detachment from the leading edge is not produced for λ=6, the same remarkable difference at ratios such as λ=5 are not seen.A similar tendency appears in the power coefficient and the thrust coefficient in Figs.7−10.
In the calculation of the laminar model, it has been suggested that it may be possible to acquire proper results at a tip speed ratio higher than λ=5.However, it takes many calculation times.

CONCLUSIONS
The calculation results have been obtained in agreement with the experiment in spite of using the coarse grid in this study.The calculation grid is very coarse, thus there is a question of whether the turbulence phenomenon can be expressed accurately enough in representing the physical shape.However, the utility value is high if one considers the practical use of calculations which solves the force acting on the blade with precision.The unsteady laminar model is useful for the performance prediction of low Reynolds number systems such as model tests and small wind turbines.The thrust coefficients have little influence whether tests are performed with or without a turbulence model.On the other hand the power coefficients are sensitive to the presence of a turbulence model, because the pressure at the leading edge and trailing edge are strongly influenced by the separation.

Figure 3 .Figures 4
Figure 3. 3D computational grid around 3 blade wind turbine.Figures 4 and 5 show the power coefficients and the thrust coefficients, which explain the influence of the k-ε turbulence model and the time step ∆t of laminar model.The time steps are 1/4000, 1/8000, 1/16000 and 1/32000 rotation.The inflow turbulence intensity is set to 1% for the k-ε turbulence model.The solutions have almost converged at the time step that is less than 1/16000 rotation.The power coefficients of the laminar model are in agreement with the experimental results at a tip speed ratio greater than five.The experimental results show the following characteristics: the maximum power coefficient C P =0.426

Figure 4 .
Figure 4. Power coefficient compared with experimental data.

Figure 5 .
Figure 5. Thrust coefficient compared with experimental data.

Figure 6 .
Figure 6.Time history of power coefficients.

Figures 9
Figures 9 and 10 show the spanwise local thrust coefficients for λ=5 and λ=6.The integration of the local thrust coefficient along the radius presents the thrust coefficient of Eq. (3).The local thrust coefficients are in agreement between the steady calculation of the k-ε turbulence model and the unsteady calculation of the laminar model.Small differences appear around the r/R=0.8and near the hub.

Figure 10 .
Figure 10.Local thrust coefficient for tip speed ratio λ=6.Figures 11 and 12 show the spanwise axial velocity which is circumferentially averaged at the rotor plane for λ=5 and λ=6.The steady calculation of the k-ε turbulence model and the unsteady calculation of the laminar model almost coincide, because the local thrust coefficients are similar in Figs.9-10.
2 chord length for r/R>0.35.(a) Pressure side (b) Suction side Figure 13.Limiting streamlines solved on blade surface by the steady k-ε turbulence model for tip speed ratio λ=5.

Figure 14 .
Figure 14.Limiting streamlines on blade surface solved by the steady k-ε turbulence model for λ=6.
Figure 2. Fluid force acting on a blade.