Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 1, pp. 203-212, Warsaw 2018 DOI: 10.15632/jtam-pl.56.1.203 STEADY AND UNSTEADY ANALYSIS OF NACA 0018 AIRFOIL IN VERTICAL-AXIS WIND TURBINE Krzysztof Rogowski Warsaw University of Technology, Institute of Aeronautics and Applied Mechanics, Warsaw, Poland e-mail: krogowski@meil.pw.edu.pl Martin O.L. Hansen Technical University of Denmark, DTU Wind Energy, Department of Wind Energy, Lyngby, Denmark e-mail: molh@dtu.dk Ryszard Maroński Warsaw University of Technology, Institute of Aeronautics and Applied Mechanics, Warsaw, Poland e-mail: maron@meil.pw.edu.pl Numerical results are presented for aerodynamic unsteady and steady airfoil characteristics of the NACA 0018 airfoil of a two-dimensional vertical-axis wind turbine. A geometrical model of theDarrieus-typewind turbine and the rotor operating parameters used for nume- rical simulation are taken from the literature. Airfoil characteristics are investigated using the same mesh distribution around the airfoil edges and two turbulence models: the RNG k-ε and the SST Transition. Computed results for the SST Transition model are in good agreement with the experiment, especially for static airfoil characteristics. Keywords: airfoil characteristics, vertical-axis wind turbine, computational fluid dynamics 1. Introduction Generally, with respect to the orientation of the rotor shaft, wind turbines can be divided into two main groups: horizontal-axis wind turbines (HAWTs, or axial flow turbines) and vertical- -axis wind turbines (VAWTs, or cross-flow turbines) (Maroński, 2016). Wind turbines can al- so be divided with respect to the principle of operation: lift-driven and drag-driven machines (Rogowski, 2014). Although, HAWTs are now widely used in the industry, large-scale VAWTs are designed as offshore units – floating wind turbines (Madsen et al., 2013; Borg et al., 2014). Aerodynamic efficiency (power coefficient) of drag-driven wind turbines is low, therefore, they are used relatively rarely (Rogowski and Maroński, 2015). In 1931, Georges J.M. Darrieus, a French aeronautical engineer, patented his invention – a new type of windmill designed for po- wer generation (Blackwell, 1974). TheDarrieuswind turbine is a lift-drivenwind turbine having two or more blades. The rotor of the Darrieus wind turbine can achieve relatively high aero- dynamic efficiency (Hau, 2006). Originally, the Darrieus wind turbine had curved blades with a symmetrical airfoil in their cross sections. The curved blade shape, so-called troposkien, was designed to avoid large bending stresses of the blades, especially when applied to large units (Paraschivoiu, 2009). Darrieus-type wind turbines are designed both as large- and small-size wind turbines with both curved and straight blades. The characteristics of the Darrieus-type wind turbines are: slightly lower power coefficient thanHAWTs (Amet et al., 2009); the gearbox and the power generator can be installed at the ground level; the yaw system is not needed because the rotor operates regardless of the wind direction. The main shortcomings of these wind turbines are: low starting torque and vibrations of the structure during rotor operation 204 K. Rogowski et al. (Paraschivoiu, 2009). The growing demand for decentralized electricity generation in urban and rural areas is the motivation for studying wind turbines in a small scale. Darrieus-type vertical-axis wind turbines are relatively simple devices. The movement of a single wind turbine blade is similar to the movement of the pitching blade. During rotation of the rotor, the blade angle of attack, the local Reynolds number and the relative wind velocity vary according to the rotor azimuthal angle. These variations also depend on the relationship between the tangential velocity of the wind turbine blade and wind velocity. For these reasons, many nonlinear phenomena occur in a single cycle of the blade (Laneville and Vittecoq, 1986). The rotor power coefficient defined as the ratio of the power absorbed by the rotor shaft divided by the power available from the air streamflowing through the rotor swept area (Hansen, 2008), depends on the tip speed ratio defined as the ratio of the tangential blade velocity to the wind speed. Typical Darrieus wind turbine achieves the maximum power coefficient of about 0.4 at the tip speed ratio of 5-6 (Hau, 2006). Dynamic effects associatedwith dynamic stall phenomena occur at low tip speed ratios (below4).Aerodynamic effects of the rotor elements such as blades, tower, struts, etc., play important role in reduction of the rotor power coefficient at high tip speed ratios (above 6) (Paraschivoiu, 2009). Although,measurement techniques havebeen improved in the recent years, only a fewexperi- mental tests of unsteadyaerodynamicblade loads have beenperformed.Measurement difficulties are particularly associated with the tangential blade load component (tangential to the rotor swept area) which is responsible for creation of the rotor torque. This is because the tangential blade load is very low compared with the normal blade load component (normal to the rotor swept area). Experiments referring to aerodynamic blade loads of theDarrieus-type vertical-axis wind turbines were performed in a water towing tank at Texas Tech University (Strickland et al., 1979, 1981). Laneville and Vittecoq (1986) conducted investigations of lift and drag airfoil characteristics of a small-size vertical-axis Darrieus-type wind turbine in a wind tunnel. Ferre- ira et al. (2011) showed that it was possible, though crudely, that the blade loading could be extracted from velocity flow fields using amethod that they had developed. Streamtube models and single-wake vortex models are often used in simulations of aero- dynamic blade loads of VAWTs (Paraschivoiu, 2009; Ferreira, 2009). Nowadays, computational methods of fluid dynamics (CFD) have becomepopular inmany areas of engineering as they can provide very accurate resultswhen referring to the experiments performedon a full-scaled object (Lichota, 2013; Lichota, 2016). The incorporated turbulence models are in numerical computa- tions a compromise between the available hardware capabilities and accuracy of computations. In order to resolve all scales of turbulence, it is necessary to apply an appropriate mesh with very small grid elements. Using a space-timemesh fine enough to compute all scales of turbulen- ce is still a very difficult task for modern supercomputers. However, the increase in computing power ofmodern computers has led to the development of computationally expensive turbulence models (Ferreira et al., 2007). Ponta and Jacovkis (2001) investigated the Darrieus-type wind turbine using a combined method consisting of a classic free vortex model and finite element techniques. Amet et al. (2009) performed CFD analysis of the two-bladed rotor basing on the experiment of Laneville and Vittecoq (1986) at tip speed ratios of 2 and 7. Many numerical simulations of two-dimensional Darrieus-type wind turbines using different turbulence models were made by Rogowski (2014). 3D simulations of a straight-bladed vertical axis tidal turbines were performed byMarsh et al. (2013) using the SST k-ω turbulencemodel. Generally, Darrieus wind turbines operate at lowReynolds numbers. The range of the blade angle of attack is very large. Streamtube models and vortex models require CL and CD airfoil characteristics in order to computeaerodynamicblade loads.Aerodynamic characteristics canbe computedusingCFDmethods (Rogowski, 2014; Sarlak et al., 2014) or performedexperimentally (Sheldahl and Klimas, 1981; Laneville and Vittecoq, 1986). The three main objectives of this work are as follows: Steady and unsteady analysis of NACA 0018 airfoil... 205 • Determination of aerodynamic coefficients for unsteadyflowaround thewind turbineusing a hybrid mesh consisting of a structured quadrilateral mesh close to airfoil edges and an unstructured triangle mesh elsewhere. • Investigation of steady characteristic of the NACA 0018 airfoil using the same mesh di- stribution as during the unsteady flow simulation of the VAWT. • Comparison of the aerodynamic characteristics for two turbulence models: the RNG k-ε and the SSTTransition. 2. Wind turbine parameters In this paper, the authors present computed airfoil characteristics of a rotating wind turbine blade and of a stationary airfoil. Computed airfoil characteristics are compared with the expe- riment of Laneville and Vittecoq (1986). The experiment was conducted in an open jet wind tunnel at the Universite de Sherbrooke. The main objective of Laneville and Vittecoq was to measure aerodynamic blade loads for a two-bladed rotor with zero offset pitch angle using strain gauges. Basic geometrical parameters of the investigated wind turbine are given in Table 1. In the central part of the rotor, a torsion-free steel shaft supported by two ball bearingswasmoun- ted. The rotor blades made of balsa wood were supported by horizontal arms at the lower part of the rotor and by two guitar wires stretched between the shaft and the blades at the upper part of the rotor. Measuring devices such as force transducers and amplifiers were placed in the lower horizontal arms. During the experiment, a special variable-speed electric motor was used tomaintain the correct rotational velocity. The effect of centrifugal forces on aerodynamic blade loads were removed from experimental data. The experimental measured data was not correc- ted for blockage effects. Themethod of measurement of aerodynamic blade load components is presented in Fig. 1. Table 1.Basic parameters of the investigated wind turbine Parameter Value Rotational speed n [rpm] 300 Rotor radiusR [m] 0.3 Chord c [m] 0.061 Airfoil NACA 0018 Number of blades N 2 Tower diameter d [m] 0.0381 Tip speed ratioTSR 5 Wind velocity V ∞ [m/s] 1.88 The static NACA 0018 airfoil characteristics CL and CD were measured in the experiment using the samewind turbine andusing the samemeasuring systemas described above (Laneville and Vittecoq, 1986). 3. Lift and drag coefficients In this paper, the angle of attack is an angle between the tangential velocity of the rotor bladeVT (VT = ωR, where ω is angular velocity of the rotor, R – rotor radius) and relative velocity VR which is a resultant of the wind speed V ∞ and the tangential velocity VT taken with theminus sign (Fig. 2) VR =V∞−VT (3.1) 206 K. Rogowski et al. Fig. 1. Silhouette of the turbine rotor and the method of measurement of aerodynamic blade loads (Laneville and Vittecoq, 1986) Fig. 2. Geometrical parameters of the rotor. Velocity vectors, angles and aerodynamic loads From geometrical considerations (Fig. 2), the vector components VR in the tangential and normal directions to the blade trajectory are respectively VRt =VT +V∞cosθ VRn =V∞ sinθ (3.2) where θ is the azimuth angle. The tangent angle of attack α is tanα= VRn VRt = V ∞ sinθ VT +V∞cosθ (3.3) Dividing the numerator and denominator of this equation by V ∞ we get tanα= sinθ VT V∞ +cosθ (3.4) TSR is the tip speed ratio defined as TSR= VT V ∞ = ωR V ∞ (3.5) Steady and unsteady analysis of NACA 0018 airfoil... 207 Taking into account the above formula in equation (3.4), the angle of attack is α=tan−1 ( sinθ cosθ+TSR ) (3.6) The relative velocity VR can be defined as VR = √ (ωR+V ∞ cosθ)2+(V ∞ sinθ)2 (3.7) Duringwind turbine operation, the blade angle of attack varies with the azimuth θ, whereas the relative wind velocity is associated with a variation in the angle of attack (Fig. 3). The lift and drag coefficients are given by CL = L 1 2 ρc(ωR)2 CD = D 1 2 ρc(ωR)2 (3.8) whereL is the lift force,D – drag, ρ – air density, c – chord. Fig. 3. Evolution of the angle of attack and the relative velocity vs azimuthal angle at the tip speed ratio of 5 The definitions of the lift and drag coefficients contain the tangential velocity of the blade VT = ωR instead of the relative velocity VR. In the case of VAWTs, the relative velocity is constantly changing in bothmagnitude and incidence. The use of the constant reference velocity in the dynamic pressure is desirable since it is possible to compare force coefficients for different airfoils (Danao et al., 2012; Amet et al., 2009). 4. Numerical model One of the main objectives of this study is to investigate unsteady aerodynamic characteristics of the wind turbine airfoil and steady aerodynamic characteristics of the same airfoil. The numerical two-dimensional model of the vertical-axis Darrieus-typewind turbine consists of two NACA 0018 airfoils and a tower which has beenmodeled as a circle (Fig. 2). Simulations of the steady airfoil characteristics have been performed using only a single NACA 0018 airfoil with the same chord. The model of the wind turbine rotor has been enclosed in a square area of a virtual wind tunnel. According to the previous investigations of the authors (Rogowski, 2014; Rogowski and Maroński, 2015), the length of the virtual wind tunnel should be at least equal to ten rotor diameters. The static characteristics of NACA 0018 have been obtained using one airfoil placed in a square area of the virtual wind turbine with the same length as in the case of the rotating rotor. 208 K. Rogowski et al. Themesh near airfoils has been created using structural quadrilateral elements. The height of the first layer of the structural grid is 7 ·10−7mgiving y+ ¬ 1. The growth rate of each layer of the structured mesh is 1.13. The airfoil edges are divided into small parts with lengths of 2·10−4m.The growth rate of the unstructuredmesh is 1.06. Themesh for unsteady simulations, presented in Fig. 4 contains of 133366 elements. In the case of the stationary airfoil, the same meshdistributionaround theNACA0018 is usedand thenumber ofmesh elements of the virtual wind tunnel is 74834. Fig. 4. Mesh distribution In this paper, two turbulence models are taken into account: the two-equation RNG k-ε and the four-equation SST Transition. The RNG k-ε turbulence model closes the average Navier-Stokes equations introducing two transport equations: one for turbulent kinetic energy and one for turbulent dissipation. The SST Transition turbulence model solves the transport equations for the turbulence kinetic energy, the specific dissipation rate, the intermittency and the transition onset criteria. More detailed description of these turbulencemodels can be found in the ANSYS, Inc.15.0 documentation. Turbulence parameters of the wind tunnel of the Universite de Sherbrooke are unknown. However, in the case of open jet wind tunnels, the turbulence intensity of the incoming flow is usually high. Therefore, in this simulation, the value of the turbulence intensity is assumed to be 5%. 5. Results 5.1. Unsteady airfoil characteristics of the wind turbine blade Figures 5a and 5b present drag and lift coefficients as functions of the angle of attack. The airfoil characteristics are computed using two turbulence models: the RNG k-ε and the SST Transition. The numerical results are compared with the experiment of Laneville and Vittecoq (1986). As it can be seen fromFigs. 5a and 5b, the computed results of the drag coefficients are more similar to the experimental results than in the case of the lift coefficient. The differences can be causedby the accuracy ofmeasuringdevices.According toLaneville andVittecoq (1986), with the increasing tip speed ratio, the precision of experimental data decreases, especially for Steady and unsteady analysis of NACA 0018 airfoil... 209 the lift coefficient. The precision of experimental data has been estimated as follows: CD ±5% andCL±12%. Even though the experiment was considered as a two-dimensional (large aspect ratio of the blades), the 3D aerodynamic effects such as tip vortices can reduce the efficiency of the device (Paraschivoiu, 2009; Scheurich et al., 2011). Analyzing the obtained numerical results of airfoil characteristics, hysteresis loops both of the lift and drag coefficients are visible (Figs. 5a and 5b). In the upwind part of the rotor, for the azimuth from 0deg to 90deg, a significant increase in the lift and drag coefficients can be observed.Moreover, at the azimuthal angle of zero, which corresponds to the zero angle of attack, the lift coefficients are 0.23 and 0.37 for the SST Transient and the RNG k-ε turbulence models, respectively. The non-zero value of the lift force may have several reasons. Firstly, the definition of the angle of attack assumed in this paper does not take into account effects associated with the slowdown of the flow close to the rotor. Secondly, during rotation of the rotor the virtual camber of the airfoil occurs at the zero angle of attack caused by curved flow around the rotor blade. This means that the symmetrical airfoil of the vertical-axis wind turbine behaves as a cambered airfoil (Akimoto et al., 2013).Moreover, in the case of the airfoil oscillating around the zero average angle of attack, CL cannot be equal to zero because of the momentum and the inertia of the fluid (Laneville and Vittecoq, 1986). With the increasing azimuth from 90deg to 180deg, the lift and the drag force coefficients decrease. In the downwind part of the rotor, the computed lift coefficients are still positive while experimental results are negative. The values of the lift force coefficients obtained by the SSTTransition turbulencemodel aremuch better compared with the RNG k-ε turbulence model. Fig. 5. (a) Drag and (b) lift force coefficient versus the angle of attack atTSR of 5 Figure 6 presents vorticity fields computed using the SST Transition turbulence model at three azimuth positions: 0deg, 60deg and 120deg. Analyzing these figures, it can be noticed that at the azimuthal angle of 0deg two interactions between the blades and the aerodynamic wake occur. At the azimuth of 120deg, the rotor blade located at the downwind part of the rotor interacts also with the aerodynamic wake from the rotor shaft. 5.2. Airfoil characteristics of NACA 0018 The second part of this paper concerns the analysis of static airfoil characteristics using the samemesh distribution close to the airfoil and the same turbulence models as during unsteady analysis. The obtained results of aerodynamic loads are compared with the experimental data taken from two independent sources (Laneville andVittecoq, 1986; Sheldahl andKlimas, 1981). The experimental data from the report of Sheldahl and Klimas (1981) are commonly used in simplified aerodynamic models for Darrieus vertical-axis wind turbines applications. The experiments were performed for the Reynolds numbers of 3.8 ·104 in the Laneville andVittecoq experiment (1986) and 4 ·104 in the experiment of Sheldahl and Klimas (1981). 210 K. Rogowski et al. Fig. 6. Evolution of the vorticity field, ω [1/s] – the SST Transitionmodel The results of aerodynamic force coefficients as a function of the angle of attack are presented in Figs. 7a and 7b. The aerodynamic derivatives ∂CL/∂α at the angle of attack range between 0deg and 5deg are: 6.0144 for the SST Transition model; 5.78 for the RNG k-ε turbulence model; 6.19 for the experiment of Laneville and Vittecoq (1986) and 4.72 for the experiment of Sheldahl and Klimas (1981). The maximum values of CL are: 0.71 at the angle of attack of 8.57deg for the experiment of Laneville andVittecoq; 0.473 at the angle of attack of 6deg for the experiment of Sheldahl andKlimas; 0.91 at the angle of attack of 12.5deg for the SSTTransition model and 1.22 at the angle of attack of 15deg for the RNG k-ε model. The minimum values of the drag coefficient at the zero angle of attack are: 0.034 for the experiment of Laneville and Vittecoq; 0.0214 for the experiment of Sheldahl andKlimas; 0.042 for the SSTTransitionmodel and 0.032 for the RNG k-ε model. The largest difference of CL data between all data series (Fig. 7b) is observed in the static-stall region. Lower drag coefficients by Sheldahl and Klimas (1981) in comparison with those by Laneville and Vittecoq (1986) can be caused by turbulence parameters of the wind tunnel. It is worth noting that the characteristics, both the lift and the drag coefficients, obtained using the SSTTransition turbulencemodel aremore comparable with the experimental results of Laneville and Vittecoq (1986) than the experimental results of Sheldahl andKlimas (1981). Fig. 7. (a) Drag and (b) lift force coefficients versus the angle of attack 6. Conclusions The main purposes of this paper are the investigation of steady and unsteady lift and drag coefficients of NACA 0018 airfoil using the same mesh distribution around the airfoil and two turbulence models: the RNG k-ε and the SSTTransition. The analysis shows that: Steady and unsteady analysis of NACA 0018 airfoil... 211 • For the samemesh around the airfoils, consisting of structural quadrilateral elements near the blades and triangle elements elsewhere, the steady airfoil characteristics are in good agreement with the experimental results. However, 3D effects may cause errors in the experiment. • The SST Transition turbulence model gives more realistic results of aerodynamic force coefficients than the RNG k-εmodel. • Comparison of the results of the lift force coefficient obtained during two independent experiments shows significant differences in the static stall zones. Thispaper gives somepreliminary results of steadyRANSmodelingof theflowpast a turbine rotor. The presented simulations are a part of a more extensive numerical study of vertical-axis wind turbines. The results of simulations presented in the paper can be a database for other investigations. Acknowledgments This work has been supported by the European Union in the framework of European Social Fund through the “Didactic Development Program of the Faculty of Power and Aeronautical Engineering of theWarsawUniversity of Technology”. References 1. AkimotoH.,HaraY.,KawamuraT.,NakamuraT., LeeY.-S., 2013,A conformalmapping technique to correlate the rotating flow around a wing section of vertical axis wind turbine and an equivalent linear flow around a static wing,Environmental Research Letters, 8, 044040 2. Amet E., Mâıtre T., Pellone C., Achard J.-L., 2009, 2D numerical simulations of blade- vortex interaction in a Darrieus turbine, Journal of Fluids Engineering, 131, 111103-1-15 3. Blackwell B.F., 1974, The vertical-axis wind turbine “how it works”, Report SLA-74-0160, Sandia Laboratories, USA 4. Borg M., Shires A., Collu M., 2014, Offshore floating vertical axis wind turbines, dynamics modelling state of the art. Part I: Aerodynamics,Renewable and Sustainable Energy Reviews, 39, 1214-1225 5. DanaoL.A.,QinN., HowellR., 2012,Anumerical study of blade thickness and camber effects on vertical axis wind turbines, Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 226, 7, 867-881 6. Ferreira C., 2009, The near wake of the VAWT: 2D and 3D views of the VAWT aerodynamics, Ph.D. Thesis, Delft University of Technology 7. Ferreira C., van Bussel G., van Kuik G., 2007, 2D CFD simulation of dynamic stall on a Vertical Axis Wind Turbine: verification and validation with PIV measurements, 45th AIAA Aerospace Sciences Meeting and Exhibit/ASME Wind Energy Symposium, Reno, 16191-16201 8. Ferreira C., van Bussel G.W., van Kuik G.M., Scarano F., 2011, On the use of velocity data for load estimation of a VAWT in dynamic stall, Journal of Solar Energy Engineering, 133, 1, 011006-011006-8 9. Hansen M.O.L., 2008,Aerodynamics of Wind Turbines, Second Edition, Earthscan 10. Hau E., 2006,Wind Turbines, Springer 11. Laneville A., Vittecoq P., 1986, Dynamic stall: the case of the vertical axis wind turbine, Journal of Solar Energy Engineering, 108, 141-145 12. LichotaP., 2013,MaximumLikelihood estimation: amethod for flight dynamics – angle of attack estimation, 14th International Carpathian Control Conference, IEEE, Rytro, Poland, 218-221 212 K. Rogowski et al. 13. Lichota P., 2016, Inclusion of the D-optimality in multisine manoeuvre design for aircraft para- meter estimation, Journal of Theoretical and Applied Mechanics, 54, 1, 87-98 14. Madsen H., Larsen T., Vita L., Paulsen U., 2013, Implementation of the actuator cylin- der flow model in HAWC2 for aeroelastic simulations on vertical axis wind turbines, 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Texas, USA 15. Maroński R., 2016,Wind Turbines (in Polish), OficynaWydawnicza PolitechnikiWarszawskiej, Warszawa 16. Marsh P., Ranmuthugala D., Penesis I., Thomas G., 2013, Performance predictions of a straight-bladed vertical axis turbine using double-multiple streamtube and computational fluid dynamics models, Journal of Ocean Technology, 8, 1, 87-103 17. Paraschivoiu I., 2009,Wind Turbine Design with Emphasis on Darrieus Concept, Presses Inter- nationales Polytechnique 18. Ponta F.L., Jacovkis P.M., 2001, A vortex model for Darrieus turbine using finite element techniques,Renewable Energy, 24, 1, 1-18 19. Rogowski K., 2014,Analysis of performance of theDarrieuswind turbine, Ph.D.Thesis,Warsaw University of Technology, Faculty of Power and Aeronautical Engineering,Warsaw 20. Rogowski K., Maroński R., 2015, CFD computation of the Savonius rotor, Journal of Theore- tical and Applied Mechanics, 53, 1, 37-45 21. Sarlak H., Mikkelsen R., Sarmast S., Sørensen J.N., 2014, Aerodynamic behaviour of NREL S826 airfoil at Re=100,000, Journal of Physics: Conference Series, 524, 1 22. Scheurich F., Fletcher T.M., Brown R.E., 2011. Simulating the aerodynamic performance and wake dynamics of a vertical-axis wind turbine,Wind Energy, 14, 159-177 23. Sheldahl R.E., Klimas P.C., 1981, Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines, Energy Report SAND80-2114, Sandia National Laboratories, Albuquerque, NewMexico 24. StricklandJ.H., SmithT., SunK., 1981,Avortexmodel of theDarrieus turbine: an analytical and experimental study, Sandia National Laboratories, Technical Report SAND 81-7017 25. Strickland J.H., Webster B.T., Nguyen T., 1979, A vortexmodel of the Darrieus turbine: an analytical and experimental study, Journal of Fluids Engineering, 101, 500-505 Manuscript received July 27, 2017; accepted for print August 31, 2017