Numerical Simulation of Unsteady Flow to Show Self-starting of Vertical Axis Wind Turbine Using Fluent

MATERIALS AND METHODS

Computational fluid dynamic, CFD tool solving the Reynolds-averaged navier-stroke equations based on moving mesh technique was conducted on October 8, 9 and 10, 2010 using high speed processor computer. The laboratory is located in the zoo campus of Harbin Institute of Technology, Department of Manufacturing and Automation. The analysis incorporates 2D unsteady flow of vertical axis wind turbine model with NACA0018 modified airfoil blade section at different tip speed ratios and steady flow conditions at three different orientations of the airfoils. The procedures followed in the analysis are as follows.

Computational analysis: The basic structure of the vertical axis wind turbine selected for the computational analysis is as shown in Fig. 4. It is a fixed pitch type with three straight blades of modified airfoil. The basic analysis system of the turbine is shown in Fig. 5.

As the blades turn about the central shaft, they encounter an incident wind that is composed of the ambient local wind velocity and the blade rotational velocity as indicated in Fig. 5.

Fig. 3:	Modified geometry of NACA0018 airfoil

Fig. 4:	Schematic representation of VAWT model created with solid works

Fig. 5:	Topographic view of the model with lift and drag component

This incident wind velocity generates lift and drag forces on the blades, which can be decomposed into a thrust force and a radial force on the turbine arms.

Fig. 6:	Boundary conditions

Fig. 7:	Mesh near airfoil

Computational domain: The modified NACA0018 airfoil was used for the analysis of 2D unsteady flow. The airfoil was set to 0.2 m chord length and the turbine radius was set to 2 m. Gambit 2.3.16 version modeling software was used to create 2D model and to generate mesh. The domain size was created with a rotating sub-domain surrounding the blades and stationary sub domains in the remaining region. Mesh for the rotating sub-domain, central stationary sub-domains and stationary sub domain located to the right side of rotating sub-domain were generated with square pave. Frontal stationary sub-domain and the last stationary sub-domain located at far end of the outlet are generated with structured type of grid. The rotating sub-domain was set to a major diameter of 5 m and minor 3 m. The inlet width and out let width were set to 3 times the major diameter of rotating sub-domain. The inlet was located 3 times the major diameter of rotating sub-domain upstream and outlet 6 times the major diameter of the rotating sub-domain downstream.

Computational method: The model and mesh generated in gambit modeling software were read into the commercial CFD code, fluent V.6.2.30 for numerical iterative solution. The RANS equations were solved using the green-gauss cell based gradient option and the sliding mesh method was used to rotate the sub-domain for the turbine blades. For pressure-velocity coupling, the simple algorithm was used.

Fig. 8:	Rotating sub-domain

Standard was set as pressure discretization and first order upwind was set for momentum. Time integration was done implicitly and the minimum convergence criteria were set to 1e-06. The RNG k-epsilon model was adapted for the turbulence closure.

For the moving mesh simulations, the computational sub-domain is split into a moving part around the turbine and a fixed part for the fixed environment. The rotational motion is simulated by allowing the mesh block around the wind turbine to rotate at constant angular velocity. The mesh movement is defined explicitly by specifying time-varying positions for all of the moving mesh block cell vertices. An interface boundary surrounding the moving mesh part within the model slides at the specified velocity. This represents the relative motion between the rotating wind turbine and the fixed environment.

Boundary conditions: The boundary conditions are shown in Fig. 6. The inlet was defined as a velocity inlet, which has constant inflow velocity while the outlet was set as a pressure out let, keeping the pressure constant. The velocity at the out let was determined by the extrapolation from inside. The no slip shear condition was applied on the turbine blades, which sets the relative velocity of blades to zero. There were four domains in the computational domain with 522 cells around each blade. The cells were concentrated near the blades as shown in Fig. 7 for better result.

The flow condition used for the analysis is shown in Table 1. Time step size was set corresponding to 2 degree for each rotational speed of the rotor (ω) given in Table 1 corresponding to each TSR.

Fig. 9:	Lift coefficient comparison, for Re=360,000

Table 1:	Flow conditions

One hundred and twenty iterations were used per time step.

The operating speed of the turbine, expressed as tip speed ration (TSR) was set between 0.1 and 1.

(1)

where, R is the turbine radius, ω is angular velocity and V∞ is the free stream velocity

Code verification: To verify the reliability of the computational method, hydrodynamic forces acting on the conventional symmetrical airfoil, NACA0018 was computed from different angles of attack. The cord length was set to 1m and the corresponding cord -based Reynolds number Re, was 360,000. The converged solution was obtained after 1500 iterations. Figure 9 shows the comparison of the computational solutions and experimental data (Sheldahl and Klimas, 1981). The computed lift forces are in good agreement with the experimental data for angles of attack between -10 and 10 degree which is considered to be the normal operating range of turbine blades. The drag coefficient comparison is also shown in Fig. 10.

Fig. 10:	Drag coefficient comparison, for Re=360,000

RESULTS

Figure 11 shows simulated torque values for the modeled NACA0018 modified airfoil at lower TSR. It shows the torque values at different azimuth angle in N-m for complete revolution of TSRs 0.1, 0.25, 0.5, 0.75 and 1. The torque values were obtained from coefficient moment (Cm) of the modeled airfoil, air density, turbine area, free stream velocity chosen and the radius of the turbine modeled. The graph shows that the average torque values at each of the TSR simulated are positive.

Figure 12 shows simulated torque values for NACA2415 cambered airfoil modeled in the same way as the modified NACA0018 using the same parameter for comparison. It shows the torque at different azimuth angles in N-m for complete revolution of TSRs 0.1, 0.25, 0.5, 0.75 and 1. The torque values were obtained using the same principle as mentioned above. As can be seen from the graph, the torque values are higher for the modified airfoil compared to the camber airfoil.

Figure 13 shows the coefficient of moment (Cm) of the simulated model for NACA0018 and NACA2415. This coefficient determines the average torque of wind turbines. The Cm values were obtained from the average moment of the three airfoils modeled through CFD computational analysis. As can be seen, Cm near zero is higher and seems to reduce up to TSR = 0.5 and then starts to rise. It is also clear that the Cm of modified NACA0018 is greater than that of NACA2415 camber airfoil.

Figure 14 shows the Coefficient of Power (Cp) for the modified NACA0018 airfoil. The graph is generated by combining the performance of turbine trailing edge inclined modeled for TSR 0.1 to 1 and without inclination of the trailing edge for TSR greater than 1.

Fig. 11:	Torque for modified NACA0018

Fig. 12:	Torque for NACA2415 camber airfoil

Fig. 13:	Comparison of Cm for NACA2415 and NACA0018 modified airfoil

The Cp was obtained from the ratio of the modeled turbine power to the available wind power in the air and used to determine the performance of the turbines.

Figure 15 shows the steady state torque values at TSR 0 for different wind speeds at three different orientations of the blades. The blade orientations were taken at three different azimuth angles of 0, 45 and 90° as shown in Fig. 16. This helps to show the performance of the turbine at its steady state. The simulation result shows that the torque values are positive at all the orientations and increases with increase of wind velocity.

Fig. 14:	Cp Curve of the modified airfoil

Fig. 15:	Torque versus velocity at three locations of blades

Fig. 16:	Blades orientation at different azimuth angle

With this result it is possible to conclude that the turbine can self start from steady state at any orientation of the blades.

In Fig. 16, Blade 1, 2 and 3 represent the steady state orientation of the turbine at azimuth angle of 0°. B1′, B2′, and B3′, represent the orientation of the turbine at azimuth angle of 45°. Lastly B1, B2 and B3′′ represent the third steady state orientation of the turbine at azimuth angle of 90°.

DISCUSSION

Blade aerodynamic forces and torque are computed from the solution of RANS equations through the integration of the pressure and shear stress over the blade surface. The total wind turbine force components and torque are obtained by adding the 3 blade force components and torque. The driving torque is obtained by calculating the average of the instantaneous values corresponding to the last revolution of the rotor. These values were used to derive the expected wind turbine power coefficient (Cp). The turbine was allowed to turn until stable torque was created and the minimum number of turns used for this model was 6 turns.

Power of the turbine is defined as:

(2)

where, p is expected output power of turbine, ω is angular velocity of turbine, T averaged torque.

Cp is then calculated as:

(3)

where, ρ is the air density, V∞ free stream velocity, A frontal area of the turbine

Symmetrical airfoils NACA0012, NACA0015 and NACA0018 are the conventional airfoil sections used in Darrieus type VAWTs. However, the main drawbacks with these types of sections are their minimum or negative torque generation at lower TSRs. For comparison of the modified airfoil result, NACA 2415 camber airfoil was numerically analyzed with computational fluid dynamic analysis taking the same parameter as the modified airfoil. The simulation result is in agreement with previous works (Lazauskas, 1992; Barker, 1983; Kirke and Lazauskas, 1991; Healy, 1978; Kotb, 1990; Habtamu and Yingxue, 2011) that camber airfoil can increase staring torque but reduction in peak efficiency and operating in narrower ranges are the draw backs.

The self- starting capability of the modified airfoil is then compared to the cambered airfoil as shown in Fig. 13. As can be seen from the simulation result, the modified airfoil has shown a better coefficient of moment for the simulated model than the camber airfoil. This indicates that the turbine can accelerate at lower TSR which cannot be possible using symmetrical airfoils.

The steady state result shown in Fig. 15 at TSR 0 also indicates that the modeled turbine can generate positive torques at all the selected three orientations. This indicates that the turbine can start to turn from stationary point.

Fig. 17:	Static pressure contour at 0° angle of attack in Pascal

Fig. 18:	Static pressure contour at 10° angle of attack in Pascal

Fig. 19:	Static pressure contour at 16° angle of attack in Pascal

The coefficient of moment and the steady result implies that the modified airfoil has shown good self- starting for the modeled turbine at low tip speed ratios.

Figure 17 to 19 shows the static pressure distribution around the airfoil in Pascal for angle of attacks 0, 10 and 16 degrees, respectively. Depending on the color of the static pressure contour, one can predict where the pressure is high and lower around the airfoil using the value of pressure displayed on the status bar for the modeled airfoil. Referring back the relationship between pressure and velocity from the Bernoulli equation one can also able to predict where the flow velocity of the wind high and low around the airfoil.

The present work has attempted to analyze the performance of Darrieus type VAWT with three fixed straight blades from different perspective through modification of existing airfoil at low tip speed ratios. The research is aimed to contribute to the current literature in renewable energy and self staring studies of vertical axis wind turbine which is a promising design for diversified applications.

CONCLUSIONS

2D unsteady flow of VAWT with NACA0018 modified geometry blade section based on fixed pitch three blades was analyzed using computational fluid dynamics. The model was analyzed at lower TSRs 0.1, 0.25, 0.5, 0.75 and 1. The model was also analyzed for its performance at its steady state or TSR = 0. The coefficients of moments were then compared with NACA2415 camber airfoil which is self- starting that is analyzed using the same parameter. Coefficients of moments for the higher torque were also analyzed without the inclination of the trailing edge for performance prediction. The simulation result shows that the modified airfoil has shown better performance in self- staring of the turbine at lower TSRs for the modeled turbine. The power coefficient of the simulated model obtained through combination is also in the normal range of turbine performance. Optimization of the airfoil through CFD is an extension of future work.