4tanguilig-danao.pmd N.R.P. Tanguilig and L.A .M. Danao 5 SCIENCE DILIMAN (JULY-DECEMBER 2017) 29:2, 5-31 A Hybrid CFD-BEM Analysis of the Aerodynamic Performance of a Cut-Out Hollow Pipe Horizontal Axis W ind Turbine Blade Neil Richard P. Tanguil ig* University of the Philippines Diliman Louis Angelo M. Danao University of the Philippines Diliman ABSTRACT This study investigated the performance of a cut-out hollow pipe blade prof ile in small horizontal axis wind turbines (HAWT). Although this type of blade was expected to have losses in eff iciency, such blade prof ile can be easily manufactured locally, and could potentially have a lower cost compared to conventional blades with aerofoil prof iles. Numerical simulations using Unsteady Reynolds Averaged Navier Stokes (URANS) w e r e u s e d t o d e r i v e t h e a e r o d y n a m i c c h a r a c t e r i s t i c s o f t h e c u t - o u t hollow-pipe sections. The Blade Element Momentum (BEM) method was then used to investigate the performance of the HAWT’s rotor. Numerical results show that cut-out hollow-pipe sections have poor aerodynamic c h a r a c t e r i s t i c s d u e t o t h e i r s i m p l e g e o m e t r y a n d c r u d e d e s i g n . B E M demonstrate that rotor with cut-out hollow pipe blades can still extract the kinetic energy of the wind but only at low tip speed ratios. Parametric studies show that the performance can be improved by altering the pitch of the blades and by adding additional blades to the rotor. Keywords: Wind energy, CFD, BEM, blade analysis, wind turbine, pipe blade _______________ *Corresponding Author ISSN 0115-7809 Print / ISSN 2012-0818 Online A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 6 Figure 1. Distribution of installed units of small wind turbine generators worldwide (Gasenger and Pitteloud 2014). INTRODUCTION According to the Global Climate Risk Index of 2015, the Philippines is the f ifth most affected country by extreme weather events in the past two decades. In fact, in 2013, the Philippines topped the list because of Typhoon Haiyan, which is the strongest tropical cyclone on record to hit land. Because of the effect of climate change in the country, the Philippines has been actively participating in international mitigation programs for climate change. One of the remarkable steps of the Philippine government in addressing the issue of climate change is the enactment of RA 9513 or the Renewable Energy Act of 2008. The said law promotes the development, utilization, and commercialization of renewable energy (RE) resources in the Philippines (Congress of the Philippines 2008). RE resources are energy resources that are naturally replenished very rapidly, such that their availability is existent over an indef inite period of time, which includes, but are not limited to, biomass, solar, wind, ocean, and hydropower (Congress of the Philippines 2008). These resources are considered as “clean” energy because they do not emit harmful pollution during their operational stage. Since the implementation of the said law, 11242.6 MW of grid-connected RE project contracts have been awarded by the Department of Energy (DOE) to potential developers. This includes 1608.9 MW of wind energy potential and 1858.43 MW of solar energy potential. Unlike solar-photovoltaic, installation of small wind turbine generators is still in its infancy stage as most of these installations happen only in three countries, namely, China, USA, and United Kingdom (Gasenger and Pitteloud 2014). Figure 1 shows the distribution of small wind turbine installation worldwide. N.R.P. Tanguilig and L.A .M. Danao 7 Small wind turbine generators, like their large scale counterparts, can either be horizontal axis wind turbine (HAWT) or vertical axis wind turbine (VAWT). Currently, HAWTs are more popular than VAWTs in terms of manufacturing production and research mainly because HAWTs can be more eff icient and have simpler aerodynamic characteristics compared to VAWTs. HAWT is a lift-type machine which uses an aerofoil prof ile as cross-section of its rotor blades to help generate a torque enough to make the rotor turn. One of the main causes that affects the dissemination of small wind turbine generators is the cost of commercially available machines, which cost an average of US$6040/KW (Gasenger and Pitteloud 2014). In the Philippines, this value is still too high for the small communities or even the local government (baranggay and/or municipal) to shoulder. Another factor that affects the mobilization of small wind turbine generators is the diff iculty of manufacturing it locally as highly eff icient wind turbines require great craftsmanship to create an optimized rotor that can harness the energy of the wind eff iciently. Optimized rotor is composed of blade sections with aerofoil prof ile and twisted accordingly to maximize the aerodynamic characteristics of the aerofoil. This is usually made from molded f iber glass or carbon f ibers to ensure the accuracy of the design. Wood can also be used to signif icantly reduce the cost but it requires a highly skilled craftsman to do so. A good example of this is the Hugh Piggott machine which was named after its creator. This machine uses a flat bottom aerofoil prof ile without any twist to reduce the diff iculty of carving the wood; however, this makes it less eff icient than its commercially available counterparts. Although this machine has a great potential of entering the global market for small wind turbine generators, it has only been marketed worldwide as a do-it-yourself wind turbine as mass production seems impossible because of the carving-skill requirement for labor. Another rotor alternative for small wind turbine generator is blades created from hollow pipe section. Unlike the Hugh Piggott machine, this does not require great craftsmanship, as the blades with its desired twist can just be cut-out completely from the hollow- pipe. The downside of this is that the eff iciency of the rotor will signif icantly suffer as blade cross-section is not the usual rounded leading and pointed trailing edges aerofoil prof ile. A small wind turbine generator with rotor made from cut-out hollow pipe is not a new concept. In fact, the online resource Wind and Wet (www.windandwet.com) provides its visitors a blueprint of a blade made from hollow pipe section by requiring certain variables, such as length and diameter of hollow pipe. They also detailed the geometry computations of the cut-out hollow pipe blade. The blade A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 8 design is created by assuming a coeff icient of lift of 0.85 regardless of the geometry of the cut-out hollow pipe. Unfortunately, they did not mention any literature that explains how they arrived at such a value. This could pose a problem as the resulting rotor could perform poorly against its expected production if the actual lift coeff icient is below the assumed value. The technical contribution of this study is to present baseline knowledge of the aerodynamics and performance of cut-out hollow pipe rotor by developing a Computational Fluid Dynamics (CFD)-based numerical model of the aerofoil prof ile from cut-out hollow pipe, and to eventually design and evaluate the performance of a cut-out hollow pipe rotor through the well-established Blade Element Momentum (BEM) Theory method. MATERIALS AND METHODS Numerical Modell ing of Aerofoil HAWT A two-dimensional CFD model was used to represent the aerofoil and the free- stream condition. The basis for the 2-D model is that the spanwise velocity component of the wind turbine blade is much lower than that of the streamwise component (Hansen 2008). This is also supported by the review of relevant literature which show that the 2-D model is suff icient for this type of study (Hills and Sládek 2005; Menter and Langtry 2006; Castelli et al. 2012; Esfahanian et al. 2013; Lanzafame et al. 2013). Moreover, the lift, drag, and moment coeff icients for the 2-D models are def ined by the following equations: where U is the incoming wind velocity; c is the chord length of the aerofoil; ρ is the density of air; A is the projected aerofoil area (chord x span); and C l , C d , and C m are the lift, drag, and moment coeff icients, respectively. 2 l 1 L= ρU cC 2 2 d 1 D= ρU cC 2 2 m 1 M= ρU AcC 2 , (1) (2) (3) N.R.P. Tanguilig and L.A .M. Danao 9 Figure 2 shows the resulting blade geometry based on a blade radius (R) of two meters and tip speed ratio (λ) of six. It has a maximum chord length (c) of 308 mm at the base and a minimum chord length of 100 mm at the tip. Based on this geometry, the authors selected the chord lengths of 100, 150, 200, 250, and 308 mm for this study with maximum camber ratios of 0.01, 0.08, 0.11, 0.19, and 0.36, respectively. Aerofoils a200, a250, and a308 can be classif ied as high-cambered aerofoils because their respective maximum camber ratios exceed 0.10 (Milgram 1 9 7 1 ) . 2-D Numerical Flow Field For this study, the calculation of the Reynolds number (Re) was based on a wind speed of 12 m/s and tip-speed ratio (λ) of six. Moreover, the standard condition for the thermophysical properties of air was considered (air density ρ = 1.225 kg/m3, kinematic viscosity μ = 1.7894x10-5 m2/s). The calculated Re for the aerofoils were 4.99×106 (a100), 3.07×10 6 (a150), 2.70×106 (a200), 2.60×10 6 (a250), and 2.64×106 (a308). However, for the density based solver with far-f ield conditions, the Mach number (Ma) must be def ined instead of the Re number. Ma is the ratio of the relative velocity to the speed of sound (340.29 m/s). The calculated Ma for the above Re were 0.2145, 0.0878, 0.0579, 0.0447, and 0.0368, respectively. Structured, unstructured, and hybrid grids were used to create all the meshing requirements of this study. The near-wall region is basically divided into three parts, namely the viscous sublayer, the buffer zone, and the fully-turbulent layer, as shown in numerous experimental studies (ANSYS 2011a). In this study, modeling the near-wall region correctly will signif icantly increase the accuracy of the simulation, since it is where the solution variables have large gradients (ANSYS 2011b). In order to correctly model the near-wall region, f iner grids should be Figure 2. Aerofoil geometries from the hollow pipe. A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 10 used to resolve the solution up to the wall (including the viscous sublayer). Therefore, choosing the f irst cell height of the grid from the wall is very important because it can signif icantly impact the accuracy of the solution. An important parameter to consider in correctly modeling the computational domain is the node density around the aerofoil. The node density study avoids the over- prediction of the node spacing of the computational grid, which can result in a longer computational time of the solver, without sacrif icing the accuracy of the solution. Initial value of the number of nodes in each upper and lower surfaces of the aerofoil was spaced at 1% of the chord length of the aerofoil, while the number of nodes on both blunt leading and trailing edges was 20. This initial value in the upper and lower surfaces was then halved and doubled to consider the effect of decreasing and increasing the node spacing, respectively. Table 1 shows the resulting number of nodes for the a100 aerofoil used in the node density study. Table 1. Number of nodes around the a100 aerofoil 1% c 2% c 0.5% c Upper and lower surfaces nodes 200 100 400 Blunt leading and trailing edge surfaces nodes 40 40 40 Total 240 140 440 Figure 3 shows the results of the lift, drag, and moment coeff icients for different node densities along the surface of the aerofoil. It can be observed that increasing the node density had little effect on the drag and moment coeff icients, especially in low angles of attack (α). In terms of the lift coeff icient (c|), the f inest node density deviated to a smaller value at α = 80° and 120° with a maximum Δc| = 0.1 at α = 80° from the coarsest node density. However, the maximum deviation of the f inest node density to mid-level node density was Δc| = 0.03 at α = 80° which can be deemed as very small and negligible. Therefore, the chosen node density for the rest of the simulation runs was computed based on a spacing of 1% of the chord length of the aerofoil to reduce simulation time while still having adequate spatial resolution to accurately solve the problem. A structured O-type grid was adapted for a100 (Figure 4), a150, a200, and a250 (Figure 5), and was created by extruding the cells normally from the face of the aerofoils. The O-type grid is more suitable for aerofoils with blunt trailing edges (Cooperman et al. 2010) compared to the conventional C-type grid, and minimizes N.R.P. Tanguilig and L.A .M. Danao 11 skewness especially for surfaces with high curvature. The f irst cell height for the grid was estimated such that the y+ value does not exceed 1, which is a requirement to accurately model the viscous sublayer of the boundary layer. (a) (b) (c) Figure 4. a100 aerofoil grid: (a) computational domain, (b) leading edge mesh, (c) surface mesh. (a) (b) (c) Figure 3. Results of the node density study: (a) lift coeff icient C l , (b) drag coeff icient C d , (c) moment coeff icient C m . A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 12 The hybrid O-type mesh was adapted for the a308 aerofoil because of the diff iculty in creating a high quality structured O-type mesh, particularly with regards to the skewness of the cells. The main factor for this is the very high camber characteristics (max camber ratio = 0.36) of the a308 aerofoil. Basically, the hybrid mesh is composed of a structured mesh from the surface of the aerofoil up to a certain height in which the skewness is tolerable, and an unstructured mesh for the rest of the computational domain. Similar to the structured O-type mesh, the f irst cell height distance was predicted by using a value of y+ equal to 1. The succeeding cells were generated by using a growth factor of 1.1 up to a grid height of 0.025 which produces cells with skewness of already up to 0.65. The unstructured section Figure 5. a250 aerofoil grid: (a) computational domain, (b) leading edge mesh, (c) surface mesh. (a) (b) (c) Aside from the node density discussed above, the nodes were also distributed such that the f irst and last cells on the upper and lower surfaces of the leading and trailing edges have a maximum spacing of 0.0005 m. This is to ensure a higher spatial resolution on the regions where higher gradients of pressure and flow are expected. Moreover, the growth rate was set at 1.1 and extruded up to a total length of 20 times the chord length of the aerofoil for both the inlet and outlet boundaries. Table 2 shows the f irst cell height and the total number of cells generated for each aerofoil type. Aerofoil Type First Cell Height (m) Number of Cells a100 0.000004 26,904 a150 0.00001 36,624 a200 0.00002 48,832 a250 0.00002 57,352 Table 2. First cell height and number of cells for each aerofoil N.R.P. Tanguilig and L.A .M. Danao 13 i i u 0 x    (4) (5)jii i j j i j tu u p u , t x x x              (6) i j i jt 2 s , (7) ji i j j i uu1 s = + . 2 x x        (8)   i i j i j j i j u p + u u = + 2 s t x x x             of the mesh was created by connecting isotropic triangle cells to the outer surface of the structured mesh, continuing up to a total grid height (structured plus unstructured) of 20 times the chord length of the aerofoil. The length of the connecting triangle cells is equal to the length of the outer quadrilateral cells to preserve the smoothness of the growth rate of the cells towards the outermost region of the domain. The total number of structured cells was 38,532, whereas the total number of unstructured cells was 21,928. To ensure proper modeling, equiangle skewness, a mesh quality parameter, can be checked. Equiangle skewness is def ined as the maximum ratio of the cell's included angle to the angle of an equilateral element. Cells with highly skewed angles can decrease the accuracy and destabilize the solution. Skewness value can range from 0 (good) to 1 (bad), and for quadrilateral cells, the value is recommended to be limited to 0.8 to qualify as a good grid. For this study, the grid skewness was limited to 0.7 due to the high-camber characteristics of the aerofoils. 2-D Numerical Simulation In this study, the problem is well within the incompressible region. The equations for conservation of mass and momentum for incompressible flow are where u i is velocity, x i is position, t is time, p is pressure, ρ is density, and t ij is the viscous stress tensor def ined by where μ is the molecular viscosity and s ij is the strain-rate tensor, Rewriting and simplifying the previous equations yield the Navier-Stokes equation in conservation form. A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 14 T ime averaging equations (4) and (6) yields the Reynolds Averaged equations of motion in conservation form, Rewriting equation (8) in its reverse yields its most recognizable form. Equation (11) is usually referred to as the Reynolds-averaged Navier-Stokes equation, where the quantity is known as the Reynolds-stress tensor. The averaging process effectively introduces unknowns through the Reynolds-stress components without any additional equations. The closure problem of turbulence is essentially devising approximations for the unknown correlations in terms of flow properties that are known, so that a suff icient number of equations exist. Ansys Fluent 14.5 was the CFD package used in all simulations performed in this study. It implements 2-D Reynolds-Averaged Navier-Stokes equations by using a f inite-volume-based solver. A density-based solver with second order implicit formulation was selected for all flow computations, since most flows are assumed to be not in line with the mesh (ANSYS 2011b). Residuals were set to 10 -6 as a global convergence criterion for improved accuracy. Generally, the ideal time-step Δti is calculated as 1/10 of the time the flow crosses the aerofoil, as shown in the following equation: However, to prevent any loss of unsteady physics of the flow f ield, suff icient temporal resolution is necessary. Therefore, additional time-step sizes with values of half and twice of the calculated ideal time-step were also tested. The a100 aerofoil in this particular study was used with time-step sizes of 0.0001 (ideal time-step size Δti), 0.00005 (0.5Δti), and 0.0002 (2Δti). The simulation for each case was performed for up to 1500 time-steps or until a periodic convergence is achieved. Moreover, the number of iterations for each time-step (internal loop) was set at 50. All three cases were tested for angles of attack α from 0° to 20° in 2° intervals. i i u 0 x    (9) (10)   i j i j i ji j i j U P U U + u u 2 S t x x x                 (11) i ij ji j i j i j U U P U 2 S u u t x x x                   j i u u   i R c t 10U   (12) N.R.P. Tanguilig and L.A .M. Danao 15 Figure 7 shows the resulting lift, drag, and moment coeff icients simulated from the different time-step sizes. The values resulting from the largest time-step size Δt highly deviated from the two smaller Δt, especially in terms of the lift coeff icient. On the other hand, the results from the two f iner Δt are in good agreement, particularly in the lower ranges of angle of attack in which the design of the blade Figure 7. Results of the time-step size study: (a) lift coeff icient, (b) drag coeff icient, (c) moment coeff icient. (a) (b) (c) Figure 6. a308 aerofoil grid: (a) computational domain, (b) leading edge mesh, (c) surface mesh. A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 16 in this study is limited to. Based on these results, the time-step size adopted for the rest of the simulation runs was computed from the ideal time-step size Δti equation above in order to minimize simulation time without sacrif icing the accuracy of the solution. BEM Method The Blade Element Momentum (BEM) theory is an old but widely used method in predicting the performance of a HAWT rotor. It was developed by Glauert in 1935 by combining momentum theory and the events taking place at the actual blades (blade element theory) to analyze airplane propeller performance (Glauert 1935). Momentum theory uses the conservation of momentum to derive the forces and flow conditions acting on the HAWT rotor. Applying the conservation of linear and angular momentum will result to the expressions for the differential thrust and differential torque, respectively: where r is the local radius; dr is the annular section thickness; Ω is the blade rotational speed; and a and a' are the axial and tangential induction factors, respectively. The axial induction factor is def ined as the fractional decrease of wind velocity between the free stream and the rotor plane: where U 1 is the velocity at free stream, U 2 is the velocity right before the rotor plane, and U 4 is the downstream velocity. On the other hand, the blade element theory expresses the forces on the blade in terms of the lift and drag coeff icients and the angle of attack (Manwell et al. 2009). The procedure subdivides the span length of the blade into N elements which are assumed to act independently of each other. The total force and moment acting on the blade is then computed by summing up the resulting forces and moments acting on each individual N element. Based on these procedures and assumptions, the derived expressions for the normal force and differential torque are as follows:  2dT ρU 4a 1 a πrdr    3dQ 4a' 1 a ρUπr Ωdr,  (13) (14) 1 2 1 U U a U    4 1U U 1 2 a ,  (15) (16) N.R.P. Tanguilig and L.A .M. Danao 17 where ϕ is the angle of relative velocity of the wind, and σ' is the local solidity ratio. For a rotor with B number of blades, the solidity ratio is def ined by: By equating equations (13) and (17) and (14) and (18), the expressions for the axial and tangential induction factors can be derived: where λ r is the local tip speed ratio, and F is the Prandtl's tip loss factor given by: Moreover, the Power Coeff icient Cp can be derived by taking the ratio of the power produced by the rotor to the power derived from the wind: where λh is the tip-speed ratio at the hub. QBlade is an open source BEM software used for the horizontal and vertical axis wind turbine design and optimization. It is based on the same BEM discussed above and is integrated with XFOIL to enable the user to design and compute the performance of a custom aerofoil. However, in this study, the aerofoil performance was predicted using the numerical methods discussed earlier. In general, the pre-stall lift and drag coeff icients data (α ≈ ±15°) of aerofoils were used in the design of HAWT blade, since turbines nowadays are stall-regulated to maximize rotor eff iciency, and in some cases, limit the produced power. In the (17) (18)   2 2 N l d2 U (1 a) d F σ'πρ C cosφ C sinφ rdr sin φ      2 2 2 l d2 U (1 a) d Q σ 'πρ C sinφ C co sφ r d r, sin φ    B c σ . 2πr  (19)  l d 2 σ C s inφ C c o sφa 1 a 4 F co s φ    (20)  l d 2 r σ C c o sφ C sinφa' , 1 a 4 Fλ c o s φ    (21)   1 B r12 2 RF c o s e x p . rπ co sφ R                     (22) (23)    h λ 3 d 2P r r l λ C8C Fλ a 1 a 1 co tφ dλ Cλ ,            A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 18 QBlade iteration process, angles of attack in the post-stall range may occur temporarily, and thus, polar data for all possible 360° angles of attack should be made available to ensure continuation of the BEM algorithms (Marten and Wendler 2013). The polar data for all 360° angles of attack can actually be obtained in the CFD simulations but it will not be computationally economical because of the much longer time required to complete the simulation. For this study, an extrapolation algorithm embedded in QBlade was used to obtain the full 360° angle of attack range. Table 3 details the blade prof ile based on the guidelines from Wind and Wet (www.windandwet.com). Since polars were only created for aerofoils a250, a200, a150, and a100, the polars of the other aerofoil prof iles were adopted from the next smaller aerofoil with available polar, except for a143.75 and a137.5 polar, which were adopted from the a150 polar, since these aerofoils differ very slightly (in terms of geometry) with the latter. Table 4 lists the variables and their corresponding values used in the rotor simulation. Aerofoil Name Chord (mm) Twist (Degrees) Position (mm) a250 250 34.0181 0 a250 250 34.0181 260 a200 200 20.9124 440 a175 175 15.238 600 a150 150 9.92682 760 a143.75 143.75 8.64081 915 a137.5 137.5 7.37031 1070 a131.25 131.25 6.11335 1225 a125 125 4.8689 1380 a118.75 118.75 3.63603 1535 a112.5 112.5 2.41386 1690 a106.25 106.25 1.20158 1845 a100 100 0 2000 Table 3. Blade section characteristics Table 4. Variables used in rotor simulation Variables Values fluid density, ρ 1.225 kg/m3 kinematic viscosity, μ 1.7894×10-5 m2/s number of blade elements, N 50 maximum number of iterations 100000 maximum ε for convergence 0.00001 relaxation factor, ω relax 0.35 N.R.P. Tanguilig and L.A .M. Danao 19 (a) (b) Figure 8. Periodic convergence plots of lift coeff icient: (a) a100 at α = 0°, (b) a150 at α = 20°, (c) a200 at á = 2°, (d) a250 at α = 6°. (d)(c) RESULTS The aerodynamic coeff icients derived for different cut-out hollow pipe aerofoil prof iles were based on the procedure described earlier. The simulations were performed at an initial number of time-steps of 1500 and were increased accordingly until periodic convergence was achieved. The internal loop or number of iterations per time-step was set at 50 in order to achieve a global convergence criterion of 1×10 -6. All simulations were performed at an α range of -20° to 20° with Re values from 2.64×106 (for a250) up to 4.99×106 (for a100). In steady simulations, a flat line of convergence is expected as more iteration is consumed. However, in unsteady or transient simulations, fluctuating force coeff icients are expected, and periodic convergence is achieved when the peaks and troughs of the curve are starting to settle into a f inal value. This periodic convergence is indicative of highly transient flow behavior with time-dependent flow separation and reattachment on the blade surface. Figure 8 shows different unsteadiness behaviors of the lift coeff icient relative to time. Figure 8a shows that convergence was already achieved within the f irst 1500 time-steps as there was already a constant fluctuation of the values of the lift coeff icient. Initial run of the 1500 time-steps for a single aerofoil at 10 different α's took approximately f ive days to f inish in an Intel i7 machine. On the other hand, the remaining three plots needed additional time-steps before a periodic or cyclical behavior of the coeff icients was achieved. Shown f inite values of the lift coeff icients were computed by taking the average during the time range in which periodic convergence was already apparent. Values of other force coeff icients (i.e. , drag and moment) at different α's were derived using the same flow process. A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 20 (b)(a) Figure 9. Numerically predicted force coeff icients of the different aerofoil sections: (a) lift, (b) drag. Figure 9 shows the lift coeff icient Cl and drag coeff icient Cd of the four aerofoils at different α. For regularly shaped aerofoil prof iles, Cl usually increases linearly with increasing α until a maximum value is reached. Hereafter, Cl decreases and the aerofoil is said to have stalled. On the other hand, Cd is often predicted to be relatively constant at low α, and then increases abruptly as the aerofoil begins to stall. However, in the case of the subject aerofoils, Cl and C d behaved differently with increasing α, because, among other things, of the geometric characteristics of the aerofoils having blunt leading and trailing edges, as well as being symmetric with respect to the vertical. It can be seen from the graphs that both C l and Cd tended to increase when the aerofoils rotate clockwise (i.e. , increasing positive α). On the other hand, while they were rotating counter-clockwise (i.e. , increasing negative α), Cl increased f irst then started to drop at higher values, while Cd tended to increase from the beginning. It can also be seen that shorter aerofoils tended to produce lower lift but higher drag than their longer counterparts. Considering that the thickness of each aerofoil is the same, shorter aerofoils appear to be “thicker” overall and flatter. This results in a much drastic flow separation from the wall compared to longer and “thinner” aerofoils, as shown in Figure 10. Flow separation on the lower side of the aerofoil decreases the pressure on that side, thereby increasing the drag and lowering the lift. (a) (b) (c) (d) Figure 10. Streamtraces of the four different aerofoils at α =10°: (a) a100, (b) a150, (c) a200, (d) a250. N.R.P. Tanguilig and L.A .M. Danao 21 Figure 12. Blade geometry created from cut-out hollow pipe. Figure 11 shows the Cl/Cd plotted against α. It can be observed that the Cl to Cd ratio was very small, especially for aerofoils a100 and a150. For a100 aerofoil, in particular, the ratio was almost unity in positive angles of attacks, which means that while Cl was increasing, Cd was also increasing. This is very important to note, as aerofoil performance, in general, depends on its capacity to generate high lift while minimizing the effects of drag. For this case, the drag could overwhelm the effects of the lift, hence degrading the performance of the blade made from these aerofoils. This characteristic could be improved by either choosing a much smaller pipe thickness to make the shorter aerofoils appear “thin” and minimize separation in the flow regime, or by altering the geometry itself of the subject aerofoils by rounding off its leading edge and thinning its trailing edge. However, these methods are not within the scope of this study. The blade geometry was based on the Wind and Wet (www.windandwet.com) procedure but was slightly altered at the base due to the non-convergence of the a308 aerofoil in the static aerofoil studies. Figure 12 shows the geometry of the blade and the location of the different aerofoils along the span. The twist angle at the hub was 34.018° and was uniformly reduced to get a 0° twist at the tip. Figure 11. Lift/drag ratio at different angles of attack. A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 22 Figure 13. Output power at different tip speed ratios. The initial analysis was performed in the Rotor BEM Simulation module of QBlade using the simulation variables in Table 4. The simulation was run in a wind speed of 12 m/s and tip speed ratio range of 0 to 8. At each λ, the radial rotor experiences different angles of attack because of the different relative velocities and wake inductions. Hence, considering the initial pitch of the blades, there is a certain λ in which the rotor will perform optimally. In this case, the optimum λ was at 0.5 which produces a power of 24.12 Watts (Figure 13). It can also be observed that, at other λ, the power is negative which could only mean that the drag force acting on the blades was much higher than the lift force being experienced by the blades. A closer look at this phenomenon is presented in Figure 14 which shows the angle of the relative velocity of the wind or inflow angle and the tangential blade force coeff icient at different blade positions. The thick green line represents the plot for the tip speed ratio of 0.5. For high tip speed ratios, the inflow angle was almost parallel to the plane of rotation which makes the direction of the drag vector to point almost directly opposite the direction of rotation. Because of this phenomenon and the fact that the blade has aerofoils with low lift to drag ratio, the resulting tangential blade force coeff icients in almost two-thirds of the span of the blade were negative and served to counter the positive tangential force coeff icients acting on the remaining one-thirds of the blade. Parametric Stud ies Parametric studies were conducted to investigate how the variation of certain parameters will affect the performance of the rotor. The wind speed, blade pitch, and number of blades of the rotor were the parameters investigated. Similar to the analysis in the initial study, all simulations were conducted in a λ range of 0 to 8 at 0.5 interval. N.R.P. Tanguilig and L.A .M. Danao 23 (a) (b) Figure 14. Results for the simulation at 12 m/s wind at different tip speed ratios: (a) inflow angles, (b) tangential blade force coeff icients. Unlike in conventional blades with regular aerofoil sections, the twist angles of the aerofoil sections for the blade created from cut-out hollow-pipe section were f ixed. Because of this, optimum twist angles for the aerofoil sections to produce maximum lift cannot be obtained. The only option left for designing the blade in order to achieve optimum performance is to vary the pitch in which the resultant of the force coeff icients acting on each aerofoil section will produce maximum lift and minimum drag forces. For this study, the investigated pitch range is from 0° to 90° at 5° interval. Figure 15 shows the power produced by the rotor at different λ when the blades were pitched at different angles. The performance of the rotor can be seen as very dependent on the pitch angle of the blades. It can also be observed that the rotor is operating only at low λ range. The highest operating λ was only at 2.5 when the blades were pitched at 20°. Moreover, the highest power output was 202.65 W (CP is 0.013) when the blades were pitched at 25° and operating at a λ of 2.0. This behavior of the rotor can be explained by looking at the axial induction factors and tangential blade force coeff icients along the span of the blade. At higher values of A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 24 λ, the rotor is partially acting like a propeller, as can be observed in the negative values of the axial induction factors (Figure 16). Instead of decelerating the flow behind the rotor, the energy from the rotor's rotation is converted into kinetic energy which then accelerates the wake rotation behind the rotor, resulting in a much less energy extraction of the rotor. On the other hand, the blades are experiencing positive values of tangential blade force coeff icient along most of the span of the blade at the lower λ levels compared to levels at higher λ (Figure 17). The overall effect of the lift vector is higher than that of the drag vector at lower λ which makes the rotor harvest the kinetic energy of the flow, and thus, produce power. Figure 16. Values of axial induction factors at different blade positions pitched at 25° (heavy line is at ë = 2.0). Figure 15. Effects of varying the pitch angle on the performance of the rotor. N.R.P. Tanguilig and L.A .M. Danao 25 Figure 17. Tangential blade force coeff icient at different positions along the blade span pitched at 25° (heavy line is at ë = 2.0). Typically, horizontal axis wind turbines are limited to three-bladed rotors. While increasing the number of blades can also increase the power, the difference is minimal or of no practical importance for high-speed machines. However, since previous studies showed that the rotor is only operating at low λ or at slow rotations, adding additional blades to the rotor can help increase the torque, and thus, the power output signif icantly. Figure 18 shows the power output of the f ive-bladed rotor. Similar to the previous analysis with three blades, the maximum power output was achieved when the blades were pitched at 25°. However, in this case, the value for maximum power was 317.20 Watts (CP is 0.02), over a hundred watts more compared to the maximum output of the three-bladed rotor. Similar to the three-bladed rotor, when the blades were pitched at 30°, the maximum power output was almost similar to that when pitched at 25°, although at a much slower rotation (at λ of 1.5). One of the reasons for the increase in power output is the increase in the axial induction factor of the blades in the f ive-bladed rotor, as can be observed in Figure 19 for the blades pitched at 25°. This indicates that an increase in the induced velocity of the rotor contributes to the rotor's rotation. Another parameter that was affected by increasing the number of blades is the tip loss correction factor. As the rotor of the wind turbine rotates, multiple helical structures in the wake were created due to tip vortices, which can signif icantly affect the induced velocity distribution at the rotor. The tip loss correction factor addresses this influence of tip vortices shed into the wake on the induced velocity f ield. The spanwise distribution of the tip loss correction factor was much higher in A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 26 Figure 18. Power performance of the f ive-bladed rotor at different pitch angles. Figure 19. Difference in the axial induction factors for f ive- and three-bladed rotors with blades pitched at 25°. rotor with more blades (Figure 20). The effect of tip vortices shed is reduced by increasing the number of blades because the wake screw downstream are kept closer together, making it rigid and thus containing the flow in the cylindrical slipstream (Branlard 2011). Figure 21 shows the lift and drag coeff icients of the a250 aerofoil section compared to a similar cambered plate section with rounded leading and pointed trailing edges. The aerodynamic coeff icients of the cambered plate section were derived from XFoil which imposes a steady state analysis (time-independent). Although the lift values were comparable, the drag values for the blunt aerofoil were signif icantly larger than that of the rounded leading edge and pointed trailing edge aerofoil. Such low drag values can have signif icant improvements in overall N.R.P. Tanguilig and L.A .M. Danao 27 performance as the tangential component of drag will be lower, resulting in higher positive torque values. Figure 21. Force coeff icients of a250 aerofoil derived from CFD analysis and a similar cambered plate derived from XFoil: (a) lift coeff icient, (b) drag coeff icient. Figure 20. Tip loss factors for f ive- and three-bladed rotors with blades pitched at 2 5 ° . (a) (b) A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 28 CONCLUSIONS An investigation of the aerodynamic performance of a crude design of cut-out hollow-pipe blade for HAWT rotor was conducted. A combination of CFD and BEM analysis was used for the investigation of the performance. A 2-D unsteady RANS model was used to derive the aerodynamic coeff icients of the cut-out hollow pipe aerofoils. The derived aerodynamic coeff icients were then used as inputs to BEM to compute the overall performance of the rotor made from cut-out hollow pipe blades. The unsteady periodic convergence of the flow is indicative of highly transient flow behavior with time-dependent flow separation and reattachment on the blade surface. Due to the geometric characteristics of the subject aerofoils, the force coeff icients at different angles of attack behaved differently. While the Cl increased, the Cd also increased with positive angles of attack, which is in contrast with the aerodynamic characteristics of conventional aerofoils in which the Cl is maximized and the effect of C d is minimized. This suggests that these aerofoils have a poor aerodynamic performance if used in a wind turbine blade. The BEM simulations were conducted using the open-source software QBlade. The force coeff icients, which were derived through CFD, were used as polars to def ine the aerofoils' aerodynamic characteristics. As expected, the crude design of the cut- out hollow pipe blade performed poorly compared to its conventional counterparts, on account of the drag acting on the blade and overcoming the effects of lift. Moreover, for higher tip speed ratios, the values of the axial induction factors along the span of the blade were negative, indicating that the rotor is acting like a propeller, wherein the energy from the rotor's rotation is converted into kinetic energy which then accelerates the wake rotation behind the rotor. Parametric studies were also conducted to observe how certain variables affect the performance of the rotor. The f irst parameter investigated was the pitch angle of blade. Adjusting the pitch angle resulted in higher power output. The maximum output power was obtained when the blade was pitched at 25°. The other parameter that was investigated is the number of blades. Increasing the number of blades (from 3 to 5) also resulted in a much higher power output at low tip speed ratios. The methods presented in this study can serve as baseline knowledge in the investigation of the aerodynamics of cut-out hollow pipe turbine blades. Although N.R.P. Tanguilig and L.A .M. Danao 29 such methods have been validated and used before in the study of conventional HAWT blades, this is the f irst time, to the author's knowledge, that such methods have been extensively used to study cut-out hollow pipe blades. Thus, to further advance the study of the subject, the following are suggested: • An experimental study with the use of wind tunnel can be performed to further validate the numerical results of the 2-D unsteady RANS model. In literature, these RANS models were used to investigate conventional aerofoils with rounded leading edge and pointed trailing edge. However, in this study, the aerofoil had blunt leading and trailing edges which can greatly affect the flow physics of the solution. Hence, doing an experimental study can lessen the exposure of using RANS models. • The crude design of the cut-out hollow pipe blade can be further improved by rounding the leading edge and pointing the trailing edge of the aerofoil. As seen in the XFoil steady simulation, doing such design alteration can drastically improve the aerodynamic coeff icients of such aerofoils, and hence, the overall performance of the rotor. This could also result in a much lower cost per power produced since no additional materials for the blade will be needed. • Full three-dimensional CFD model of the rotor and flow regime will eliminate the simplifications brought about by the two-dimensional model and the assumptions in the BEM. However, it is computationally expensive because millions of cells are required for creating a good quality mesh and time-steps will be much smaller in order to not overlook the effects of unsteady physics. A supercomputer is usually used to do such simulations. ACKNOWLEDGMENTS The authors would like to thank the Department of Science and Technology (DOST)- Engineering Research and Development for Technology (ERDT) and the Off ice of the Vice Chancellor for Research and Development (OVCRD), University of the Philippines - Diliman for the funding provided for this work. A Hybrid CFD-BEM Analysis of the Pipe Horizontal Axis W ind Turbine Blade 30 REFERENCES Airfoil Tools [Internet]. 2015. Airfoil tools; [cited: 2015 November 15] Available from: http://airfoiltools.com/airfoil/details?airfoil=cp-180-050-gn. [ANSYS] ANSYS Inc. (USA). 2011a. ANSYS FLUENT theory guide. Pennsylvania (PA): ANSYS Inc. [ANSYS] ANSYS Inc. (USA). 2011b. ANSYS FLUENT users guide. Pennsylvania (PA): ANSYS Inc. Branlard E. 2011. Wind turbine tip-loss corrections [dissertation]. Denmark: Risø DTU National Laboratory for Sustainable Energy. Castelli MR, Grandi G, Benini E. 2012. Numerical analysis of laminar to turbulent transition on the DU91-W2-250 airfoil. World Academy of Science, Engineering and Technology. 6(3):914-921. [Congress] Congress of the Philippines (Phil). 2008. RA 9513: Renewable energy act of the Philippines. Metro Manila: Republic of the Philippines. Cooperman A, McLennan A, Baker J, Dam CV, Chow R. 2010. Aerodynamic performance of thick blunt trailing edge airfoils. In: 28th AIAA Applied Aerodynamics Conference; Illinois. Reston, Virginia (VA): AIAA . Esfahanian V, Salavati Pour A , Harsini I, Haghani A, Pasandeh R, Shahbazi A, Ahmadi G. 2013. Numerical analysis of flow f ield around NREL phase II wind turbine by a hybrid CFD/BEM method. Journal of Wind Engineering and Industrial Aerodynamics. 120:29- 36. Gasenger S, Pitteloud J. 2014. Small wind world report 2014. Bonn, Germany: World Wind Energy Association. Glauer t H. 1935. Airplane Propellers. In Durand WF, editor. Aerodynamic theory. Vol. IV Div. L. Berlin: Springer. p. 169-360. Hansen MO. 2008. Aerodynamics of wind turbines. London: Earthscan. Hills J, Sládek A. 2005. Numerical modelling of turbulent flow past an airfoil. Prague: Czech Technical University. Ingram G. 2005. Wind turbine blade analysis using the blade element momentum method. Durham: Durham University. Kreft S, Eckstein D, Junghans L, Kerestan C, Hagen U. 2014. Global climate risk index 2015. Berlin: Germanwatch e.V. L a n z a f a m e R , M a u r o S , M e s s i n a M . 2 0 1 3 . W i n d t u r b i n e C F D m o d e l i n g u s i n g a correlation-based transitional model. Renewable Energy. 52:31-39. N.R.P. Tanguilig and L.A .M. Danao 31 Le a r y J . 2 0 1 0 . Co m p u t a t i o n a l f l u i d d y n a m i c s a n a l y s i s of a l ow co s t w i n d t u r b i n e . Sheff ield: The University of Sheff ield. Manwell J, McGowan J, Rogers A . 2009. W ind energy explained—theory, design and application. 2nd ed. United Kingdom: John Wiley and Sons Ltd. Marten D, Wendler J. 2013. QBlade guidelines. v0.6. Berlin: TU Berlin. Menter FR, Langtry RB. 2006. A correlation-based transition model using local variables- par t II: Test cases and industrial applications. Journal of Turbomachinery. 128(3):423- 434. Milgram JH. 1971. Section data for thin, highly cambered airfoils in incompressible flow. Massachusetts (MA): Massachusetts Institute of Technology. Sørensen NN. 2009. CFD modelling of laminar-turbulent transition for airfoils and rotors using the model. Wind Energy. 12(8):715-733. Wind and Wet [Internet]. 2015. Wind and Wet; [cited: 15 November 2015]. Available from www.windandwet .com. _____________ Neil Richard P. Tanguil ig is presently Assistant Project Site Manager for Alternergy Wind One Corporation. He was involved in the pre- development, construction, and commissioning of the 54 MW Pililla Wind Farm in Pililla, Rizal. Currently, his responsibilities include wind resource assessment, wind and power forecasting, and project management. He obtained both his BS Civil Engineering and MS Energy Engineering degrees from the University of the Philippines Diliman. Louis Angelo M. Danao is presently an Associate Professor of the Department of Mechanical Engineering, University of the Philippines Diliman. He obtained his Ph.D. in 2012 from the University of Sheff ield, Sheff ield, UK, where he conducted experimental and numerical work related to the performance and aerodynamics of vertical axis wind turbines. Prior to this, he was involved in several solid mechanics research including the analysis of abdominal aortic aneurysms and the deformation of external f ixator for tibial fractures using f inite element analysis. He also developed mathematical models for torsion of solid and hollow rectangular sections using analytical and FEA approach. Presently, he is carrying out work on design and performance analysis of horizontal axis tidal turbines. He is a member of the American Society of Mechanical Engineers, the Philippine Society of Mechanical Engineers, and the International Association of Engineers.