Numerical Analysis for Dynamic Stall Control Elements on Vertical-Axis Wind Turbines

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Numerical Analysis for Dynamic Stall Control Elements on Vertical-Axis Wind Turbines Autor: Miguel Cuesto Rubio Tutor: Luis Balam Modesto López Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Numerical Analysis for Dynamic Stall Control Elements on Vertical-Axis Wind Turbines Autor: Miguel Cuesto Rubio Tutor: Luis Balam Modesto López Profesor Contratado Doctor Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Trabajo Fin de Grado: Numerical Analysis for Dynamic Stall Control Elements on Vertical-Axis Wind Turbines Autor: Miguel Cuesto Rubio Tutor: Luis Balam Modesto López El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos Es te trabajo fin de grado marca el final de una etapa comenzada hace 4 años, una etapa marcada por penas y glorias, pero de la que realizo un balance general muy positivo. Por supuesto, esto no habría sido posible sin la inestimable ayuda de mis padres, de mi hermana, y de mis amigos, de Sevilla, de Huelva, y de todas las partes del mundo. ¡El futuro es verde! I Resumen La energía eólica está demostrando ser un medio de generación de energía cada vez más fiable y eficiente, al ser posible una generación de energía en masa y libre de subsidios económicos. Esto está todo siendo conseguido por las turbinas eólicas de eje horizontal (HAWTs). Sin embargo, existe otro tipo de turbinas de eje vertical (VAWTs), que tras una pausa de una década, están volviendo a ser investigadas. Las VAWTs podrían posibilitar generación de energía silenciosa y eficiente en ambientes urbanos o muy cambiantes, donde la naturaleza inflexible de las HAWTs (por ejemplo, las HAWTs no pueden adaptarse a cambios grandes en la dirección del viento) las limita. Uno de los principales problemas a los que se enfrenta la implementación de VAWTs es la ocurrencia de entrada en pérdida dinámica a ciertos ángulos relativos a la dirección del viento. La entrada en pérdida afecta gravemente a la turbina, decreciendo el par dinámico y por tanto la energía generada, además de creando momentos flectores que pueden afectar la estructura de la turbina en general. Por tanto, el objetivo de este proyecto es estudiar la posibilidad de evitar la entrada en pérdida implementando en la pala elementos pasivos simples de control, y determinar el efecto de estos elementos pasivos en la generación de energía, comparado con el caso de turbina base. Al contrario de los dispositivos activos de control (por ejemplo, ”jets” sintéticos o control de cabeceo), los elementos pasivos tienen como objetivo la simplicidad estructural, lo que es un requisito para las aplicaciones domésticas. De esta manera, es posible ofrecer costes bajos de adquisición y mantenimiento. III 2Chapter 1. Introduction 1.2 Developments on urban Wind Energy In the last 20 years, focus in research has been shifted to wind energy production in urban environments. As per [4], three modes of integration in buildings have been mainly pursued as follows: • integration of traditional free-standing wind turbines in buildings. • retrofit of wind turbines onto already existing buildings. • integration of specifically-designed wind turbines for urban implementation purposes. Such wind turbines are the object of this study, the so-called vertical-axis wind turbines (VAWTs from now on). In contrast to horizontal-axis models (HAWTs), VAWTs do not need to yaw to adjust to wind direction, as energy production is independent of this factor. Moreover, their lower rotational speeds (and overall smaller sizes) hinder a noise pollution problem. As an example, in rooftop applications, the noise level produced by the turbines themselves is normally lower than the one produced by wind shear as the current tries to shape its way past the building [5]. 1.3 VAWT typology VAWTs can be divided into 3 main types: Savonius, Darrieus and hybrid wind turbines, as seen on Figure 1.2. Figure 1.2 Savonius, Egg beater (type 1 Darrieus), and H-rotor (type 2 Darrieus) turbines, compared to a HAWT [6]. 1.3.1 Savonius turbines Savonius wind turbines are based on drag forces which makes them spin, due to their specific shaping. This very shaping also implies a lower flow energy utilisation in contrast to Darrieus wind turbines, and thus a comparatively lower efficiency [ 7 ]. The main advantage to this type is its self-starting ability, due to good starting torques when compared to lift-based models. 1.3.2 Darrieus turbines These turbines feature high efficiencies but no self-starting abilities, as they are lift-based. However, their lift generation mechanisms can place a great level of strain on the blades and tower, even reaching breaking points for certain materials. One can classify Darrieus turbines in two larger groups. 1.3 VAWT typology 3 • Egg beater-shaped Darrieus VAWTs feature two curved blades which meet at a lower and upper hub at the tower. The geometry enables minimal bending moment on the blades, thus enabling larger wind turbines following this principle, and hence getting into megawatt-producing areas as well. This being said, this complex geometry also implies limited mass production and higher costs. • Straight-bladed wind turbines (also called H-rotors) are most often built up of two or three blades with constant airfoil features and rotating around a hub with a constant radius. Production costs are comparatively lower, while the experienced bending moments are highest. Said structural issues though don’t pose a problem in the kilowatt operational area, such as for urban, domestic applications , while severely limiting their prospects as a mass energy production medium. Furthermore, the blade geometry associated to H-rotors can induce problems associated with unsteady aerodynamics, which can decrease the energy yield strongly. The object of this project is thus to answer the question ”Is there a way to improve Straight-bladed VAWT performance, by means of mitigating the effects produced by dynamic stall?”. 1.3.3 Hybrid turbines Self-starting capabilities at low wind velocity are incorporated into a Darrieus turbine, this being enabled by an additional Savonius rotor in the setup. 2 Physical Background 2.1 The double-multiple Streamtube Theory (DMS) I. Paraschivoiu undertook in 1988 an analysis of Darrieus rotor aerodynamics through the double-multiple streamtube theory [ 8 ], taking the Blade Element Momentum (BEM) method as its base. In the former, each section of the overall rotor is made up of a pair of actuator disks in tandem (conversely, only one actuator disk is used for BEM). Henceforth, the study on the induced velocities is undertaken separately in the up- and downstream halves, and each half is further divided in multiple aerodynamically independent streamtubes. This means, the streamtubes are supposed not to interact with each other. A schematic view on the streamtubes and their induced velocities is as per Figure 2.1. Figure 2.1 Relative velocity diagram on the level of a Darrieus turbine [9]. 2.1.1 Upwind kinematic analysis The upwind velocity (V) is less than the local ambient air velocity ( V∞1 ), them being related through an interference factor u, less than 1: V0=uV∞1(2.1) Additionally, the equilibrium velocity in the plane between both upstream and downstream currents ( Ve ) is also less than V: Ve= (2u−1)V∞1(2.2) 5 6Chapter 2. Physical Background Therefore, the induced velocity (V’) decreases in the axial streamtube direction, so V’ <Ve . The upwind, induced and rotational velocities form a velocity triangle, which can be seen on Figure 2.1. For - π /2≤ θ ≤ π /2, the local velocity of the upstream is given by: Wu =qVu2[( λ −sin θ )2+cos2 θ ](2.3) With λ being the so-called local tip speed ratio (TSR), which links local blade speeds and air speeds: λ = ω R V0(2.4) The local angle of attack ( α ) is defined as follows: α =arcsin cos θ cos α 0−(X−sin θ )sin α 0 p[( λ −sin2 θ )2+cos2 θ ]!(2.5) In this case, α 0 is the residual pitch of the blade, alas, the angle between the tangent line pertaining to the rotation circumference and the blade chord. Furthermore, a local Reynolds number can be defined, which measures the viscous behaviour of the air: Re = ρ Wuc µ ,(2.6) where ρ is the density, Wu the local velocity, c the chord of the airfoil and µ the dynamic viscosity of the air. 2.1.2 Dynamic analysis Akin to lift and drag forces, we can define the tangential and normal forces on the turbine blade: FN( θ ) = 1 2 ρ cLW 2CN(2.7) FT( θ ) = 1 2 ρ cLW 2CT(2.8) In the latter expressions, ρ is the density of the fluid, c the blade chord, L its length (can also be interpreted as H, blade height), and CN and CT are the non-dimensional coefficients of the normal and tangential forces, respectively. CN=CLcos α +CDsin α (2.9) CT=CLsin α −CDcos α (2.10) CL and CD are the so-called lift and drag coefficients, respectively. They are obtained as dimensionless coefficients which characterise the lift and drag forces acting upon an airfoil. Their distribution on an airfoil with respect to the inflow velocity is illustrated by Figure 2.2. The lift coefficient is: CL=L 1 2 ρ V2S,(2.11) And the drag coefficient: CD=D 1 2 ρ V2S,(2.12) V in this case being the relative velocity between airfoil and inflow. The torque produced locally is calculated as per the following (R being the turbine radius, from the turbine tower to the rotation circumference): Q( θ ) = FT( θ )R(2.13) 2.2 Dynamic Stall 7 Figure 2.2 Aerodynamic forces on an airfoil [10]. And ultimately, the torque provided by the upstream half is as per: Qu=N 2 π Z π /2 − π /2 Q( θ )d θ (2.14) A non-dimensionalisation further provides CQu: CQu=Qu 1 2 ρ V02SR (2.15) 2.1.3 Power analysis The power coefficient for the upstream half is calculated through: CPu=CQu λ r(2.16) In this case, λ r is the rotor TSR, defined through λ r=R ω /V∞1 . Following the same procedure for the downstream half (for π /2≤ θ ≤3 π /2 ), CPd shall be obtained, and thus, the total power coefficient shall be CP=CPu+CPd . As an additional remark, the DMS theory does not assume unsteady flow and thus there is no consideration for phenomena such as dynamic stall. 2.2 Dynamic Stall Stall is a phenomenon which occurs when the angle of attack perceived by the airfoil increases past the angle at which the lift coefficient is maximum ( CLmax ). As can be seen on Figure 2.3, CL behaves linearly until CLmax , when the coefficient reaches a plateau and then plummets. This decrease is produced by a separation of the current away from the airfoil. A VAWT blade is constantly rotating at a speed ω , and therefore the perceived angle of attack by the blade is also constantly changing, under the assumption that the wind direction remains constant. This perceived angle of attack nonetheless changes with the turbine TSR, as per the formula in section 2.1. This is shown on Figure 2.4, which illustrates the perceived angles of attack for λ =2,3,4,5 . Taking α CLmax ≈13 as an example, it can be determined that the turbine will experience stall for λ <4. This constant change of angle of attack along the azimuth position induces large contrasts in lift and drag forces in the turbine operational area. It is so, that for certain instants the turbine may experience deceleration, as for high angles of attack drag forces will prevail. This disparity along the azimuth also puts a sizeable structural strain on the tower and the turbine overall, as certain blades might be stalled and others might not at any given point. Vibratory loads are just as important under dynamic stall conditions, leading to structural fatigue and noise. 8Chapter 2. Physical Background Figure 2.3 CLagainst α for a DU 06-W-200 airfoil. Figure 2.4 α against θ for a VAWT airfoil. Under dynamic stall, the turbine experiences shed vorticity, VAWTs being especially sensitive to such phenomenon. A vortex created and shedded on the upstream passage of the blade may affect the blade further on the downstream half. This vortex induces an unfavourable current situation on the blade, even thought it may not be under stall conditions due to its azimuth position. 2.3 The Lifting Line Theory (LLT) The double-multiple streamtube method is nevertheless just valid as a first approximation, as convergence problems arise in case of high blade loading, high solidities or high TSRs. Therefore, for further preliminary study the LLT can be used, as such convergence problems lessen, and computational costs are lower than those usual for RANS- or LES-based Computational Fluid Dynamics (CFD) simulations. This method was developed by Marten et al. [ 11 ], and incorporated into their open-source turbine design software suite, QBlade, which will be used in this study. It is worth noting that for VAWTs there is still no certification-grade software which is low on computational cost, like GH Bladed for HAWTs. That is why for model verification, CFD has to be used and ultimately wind tunnel experiments should also be undertaken. 2.4 Computational Fluid Dynamics (CFD) 9 2.3.1 Physical fundamentals For each time step, blade forces are calculated using the 1/4-chord line as the lifting line, which contains distributed vortices (hence the overall circulation in the blade being Γt , calculated through tabulated specific parameters for the airfoil, as well as the inflow angle for the timestep). Downstream from the airfoil, the wake is calculated through free-flowing vortex elements which separate from the trailing edge of the airfoil. For the case of HAWTs, wake vortex elements can be said to have a negligible influence after a certain downstream distance from the airfoil and are therefore truncated at a certain age. This is mainly due to the fact that for the wake of a HAWT, the trailing vorticity is dominant. Such an approximation cannot be made for a VAWT, as the angle of attack is changing constantly and therefore, trailing and shed vorticity effects are equally significant. The strategy for VAWTs is therefore rather a reduction on the wake, as it can be assumed that for every time step, there are vortices separating from the trailing edge with negligible circulation. As an example, for straight-bladed turbines (the case which shall be further studied), shed vorticity is negligible at low blade loading positions, and trailing vorticity from inner blade sections [ 11 ]. Thus, reducing shed vorticity in the model enables this method to have a lower computation cost, as a 50% reduction on vorticity decreases the computation time by 75%, while it can be ensured that the difference is not any greater than 1%, in comparison to there not being any vortex reduction. This method is nonetheless best suited for studies in high TSR cases, as there is no consideration for dynamic stall situations, and a successful calculation is dependent on good polar data. As an example, after setting up an airfoil, QBlade calculates the airfoil polar for a restricted azimuth through XFOIL and then extrapolates using the Viterna-Corrigan method for 0 ≤ θ ≤ 2 π . This yields an unrealistic result for low TSRs when compared to an experimentally-obtained aircraft polar, which takes phenomena such as dynamic stall into account. This study was undertaken in [ 12 ]. Figure 2.5, also taken from the same study, shows the effect of dynamic stall on the power coefficient (in red), and subsequently the disparity with the stall-absent case (in blue). Figure 2.5 Influence of dynamic stall on the power coefficient. 2.4 Computational Fluid Dynamics (CFD) ANSYS Fluent will be used as the CFD analysis tool for the project. CFD is needed for a potentially more exact turbine performance calculation, especially due to the built-in turbulence models. The reason is that turbulence follows highly chaotic and time-dependent behaviours which the forementioned, simplified methods fail to take into account. 10 Chapter 2. Physical Background 2.4.1 Fundamentals of CFD The main principle behind CFD is the discretisation of variables. Through this discretisation: • Continuous variables are approximated into their values in a finite amount of points (these points being called nodes). • Derivatives are therefore converted into algebraic differences. • Equations which make use of discretised variables, work out values in the nodes. This means, the system of partial differential equations (PDEs, namely Reynolds-Averaged Navier-Stokes equations) associated with CFD will be applied to the given nodes. Furthermore, there are several methods for discretising: • Finite Difference Method (FDM): continuous variables are supposed to vary linearly between two given points, and therefore derivatives are calculated as slopes. As an example, ∂ x ∂ y≈xi+1,j−xi−1,j 2∆y. • Finite Element Method (FEM): the nodes form a coarse grid, which in turn is composed by nodes being linked to their closest ones (each polygon being its own ”finite element”. This way, the PDEs can be solved simply for every finite element. Later on, these solutions are used to further assemble a system of equations which obtains the global solution of the given problem. • Finite Volume Method (FVM): a finite volume is assembled around each node and the variables are supposed to make up a conservative flux (the input and output flows are the same). This way, volume integrals part of a PDE can be converted into surface integrals by the theorem of divergence. These newly-converted terms are then evaluated as fluxes at the surface of each finite volume. This is the method behind ANSYS Fluent and CFD. After designing the geometry and the flow field, the flow field is discretised and a mesh is formed. The PDEs are then evaluated at the nodes, through which the PDEs are transformed into a system of algebraic equations. These equations are then solved numerically, thus obtaining values for every node in the mesh. 2.4.2 Reynolds-Averaged Navier-Stokes equations (RANS) Assuming incompressible flow, the Continuity equation is: ∂ρ ∂ t+∇·( ρ u) = 0,(2.17) where ρ is the density of the fluid (constant due to incompressible flow), ∇ the nabla operator, and u(x,t) the velocity. Similarly, the Momentum equation is: ∂ u ∂ t=u·∇u− ν ∆u=−1 ρ ∇p+g,(2.18) where ν is the kinematic viscosity ( ν = µ ρ ), p(x,t) is the pressure, and g are the volume forces which act upon the flow (i.e gravity). Averaging in time, the Continuity equation transforms into: ∇·U=0,(2.19) while the time-averaged Momentum equation is: ∂ U ∂ t+U·∇U− ν ∆U=−1 ρ ∇P+G−1 ρ ( ρ u0 0 0·∇u0 0 0),. (2.20) where U,P, and G are the time-averaged components of the velocity, pressure and volume forces respectively, and u0 0 0 is the fluctuating component of the velocity. These are the so-called Reynolds Averaged Navier-Stokes equations, also RANS equations. 2.4.3 Reynolds turbulence As seen in the time-averaged Momentum function, velocity can be decomposed into an averaged component ( U ), and a fluctuating component ( u0 ), thus being u=U+u0 . The fluctuating component of velocity can 2.4 Computational Fluid Dynamics (CFD) 11 be traced back to turbulence, which is a 3D, transient, fluctuating movement associated to a time-averaged main flow. The term −1 ρ ( ρ u0 0 0·∇u0 0 0) , also defines the turbulence in the time-averaged Momentum equation. In tensor notation, can be written as ∂ ∂ xj( ρ u0 iu0j) . The expression between the brackets is the Reynolds tensor, τ i j = ρ u0 iu0j , where the first index indicates the plane of action of the stress, and the second index indicates the direction upon which the stress acts [13]. Reynolds stresses cannot be known beforehand for turbulent current cases, and therefore they need to be modelled. One method to model the Reynolds stresses is through eddy viscosity. The Boussinesq eddy viscosity approximation is as follows: − ρ u0 iu0j= ν t( ∂ Ui ∂ xj + ∂ Uj ∂ xi )−2 3k δ i j,(2.21) where k is the turbulent kinetic energy and ν t the eddy viscosity. Both variables are still unknown, but they can be modelled through several methods. In our research, they will have a common modelling equation for k, and a more sophisticated system of equations for ν t. 2.4.4 The k−equation It is sensible to begin the turbulence equation definition by shaping up a kinetic energy model, and thus modeling the characteristic velocity component of a turbulent field as the square root of the kinetic energy pertaining to the field: k=1 2(u0 xu0 x+u0 yu0 y+u0 zu0 z)(2.22) It can be seen that k equals the sum of the Reynolds stresses multiplied by 0.5. Thus, using the Reynolds stress equation (the calculation process goes beyond the scope of this study): ρ∂ k ∂ t+ ρ∂ ∂ xj (Ujk) = τ i j ∂ Ui ∂ xj− µ∂ u0 i ∂ x0 k ∂ u0 i ∂ x0 k + ∂ ∂ xj µ∂ k ∂ xj−1 2 ρ u0 iu0 iu0j−p0u0j(2.23) The term µ∂ u0 i ∂ x0 k ∂ u0 i ∂ x0 k is of special interest for this study on fluid dynamics, as it represents the dissipation of kinetic energy in terms of heat, at a molecular level. This dissipation is denoted by ε : ε = ν∂ u0 i ∂ x0 k ∂ u0 i ∂ x0 k (2.24) As ∂ u0 i ∂ x0 k will always be positive, it can be determined that ε will imply a loss of energy away from the flow, by means of dissipation. 2.4.5 The k− ε Turbulence Model The k− ε turbulence model adds two additional equations to the RANS. The k−equation: ρ∂ k ∂ t+ ρ Uj ∂ k ∂ xj = ∂ ∂ xjh µ + µ t σ ki ∂ k ∂ xj+ τ i j ∂ Ui ∂ xj− ρε ,(2.25) and the ε -equation: ρ∂ε ∂ t+ ρ Uj ∂ε ∂ xj = ∂ ∂ xjh µ + µ t σε i ∂ε ∂ xj+C ε 1 ε k τ i j ∂ Ui ∂ xj−C ε 2 ρε 2 k(2.26) The eddy viscosity µ t is to be calculated through µ t= ρ C µ k2 ε . C ε 1 , C ε 1 , C µ , σε and σ k are empirically obtained closure coefficients. 18 Chapter 4. Airfoil Analysis Figure 4.6 LLT, operational point. Figure 4.7 LLT, output, algorithm, wake parameters. 4.2 Base Configuration: the NACA 0021 airfoil 19 4.2 Base Configuration: the NACA 0021 airfoil The NACA 0021 airfoil was chosen specifically due to its improved capability to operate under dynamic stall. For similarly-sized turbines as the case in this study, wind tunnel tests have been undertaken, comparing NACA0015 and NACA0021 airfoils, as well as changes in blade pitch and its influence on performance. Blade pitch changes are nonetheless out of the scope of this study. The NACA 0021 airfoil is part of the NACA 4-digit series, developed by the National Advisory Commitee for Aeronautics (NACA) in the 1940s. The 4 digits define key features for the airfoil: • The first digit defines relative camber, as a percent of the chord • The second digit defines the relative distance of the maximum camber point from the leading edge, as a percent of the chord • The last two digits define the maximum thickness of the airfoil, as a percent of the chord For a given position along the chord (as a percent of the chord), the local half thickness is given by: y=5t[0.2969√x−0.1260x−0.3516x2+0.2843x3−0.1015x4](4.1) Thus, points can be calculated in order to plot the airfoil as a first step. These values are illustrated on Appendix A.1. Figure 4.8 NACA0021, clean airfoil. Under null pitch conditions, [19] showed that the point of CPmax for the NACA0015 airfoil was given for λ =1.6 , while for the NACA0021 airfoil, CPmax would be obtained at λ =1.45 . It was speculated that a lower value for λ could therefore be linked with a more adequate dynamic stall performance. Furthermore, in the same study, it was determined that with a higher blade pitch (namely β = +7.8deg ), the turbine featuring a NACA0021 airfoil was still able to yield power in spite of dynamic stall. In the meantime, under the same conditions, such a blade pitch surpassed the operational conditions for the NACA0015 airfoil. This phenomenon could as well be justified by improved dyamic stall performance. As a state-of-the-art analysis, most present day studies on VAWTs and dynamic stall are being undertaken on several NACA airfoils, for the most part NACA0012 and NACA 0018. The NACA0015 airfoil would then be at an intermediate point, in terms of performance. After entering data for the NACA0021 airfoil on QBlade, the airfoil polars were also generated through this same tool. QBlade uses XFoil as its airfoil generator. Its using XFoil will pose problems later, as XFoil is unable to converge for complicated geometries, and another airfoil generator, namely JavaFoil, was found. XFoil Polars As generated by QBlade, the Cl plotted against Cd and Cl plotted against α can be seen on Figures 4.9 and 4.10, respectively. It can be seen that for null Cl , Cd equals 0.17 and Cd rockets for Cl=1.1 . This increase in drag coefficient can be attributed to the occurrence of stall on the airfoil, as is shown on Figure 4.10. Moreover, using these polars QBlade calculates the power and torque coefficients. These values are illustrated on Figures 4.11 and 4.12. The latter is of special interest, as it shows the least favourable azimuth, at which dynamic stall is most present, producing negative torque, and therefore power, on the turbine. On the former, it is seen that the maximum power coefficient value is 0.52. This is a highly irrealistic value which shows the overestimation expected due to the use of the DMS theory. 20 Chapter 4. Airfoil Analysis Figure 4.9 Clagainst Cd, clean configuration, XFoil. Figure 4.10 Clagainst α , clean configuration, XFoil. Figure 4.11 CPagainst λ , clean configuration, XFoil. JavaFoil Polars JavaFoil was also used for additional airfoil polar calculation. JavaFoil[ 20 ] first determines the velocity distribution through a panel method. This panel method divides the airfoil into multiple straight panels, with each inducing unknown velocities on themselves and others. These velocities are then calculated through a system of linear equations. Here, a boundary condition is that the velocities are always tangent to the 4.2 Base Configuration: the NACA 0021 airfoil 21 Figure 4.12 Ctagainst θ , clean configuration, XFoil. panels. Later on, a viscous boundary layer analysis is undertaken, so as to calculate the friction drag, as well as lift, drag and moment corrections in case of flow separation. JavaFoil obtained fairly similar polars to those obtained by XFoil. These JavaFoil polars were then also used by QBlade to calculate similar performance coefficients, when compared to the XFoil case. These polars and values can be seen on Figures 4.13, 4.14, 4.15 and 4.16. Irregularities in the polars can be seen in these Figures, and also a post-processed alternative. Peaks and buckets in the graphs were smoothed out using MATLAB (see List of Codes), using a 4th-degree polynomial best fit. These irregularities in the polars can be attributed to convergence problems during polar generation. Post-processed polars shall be used for further analysis. Figure 4.13 Clagainst Cd, clean configuration, JavaFoil. 22 Chapter 4. Airfoil Analysis Figure 4.14 Clagainst α , clean configuration, JavaFoil. Figure 4.15 CPagainst λ , clean configuration, JavaFoil. Figure 4.16 Ctagainst θ , clean configuration, JavaFoil. Lifting Line Parameter adaptation The Lifting Line Theory method was used to simulate a regular load case, with λ =2.5and constant wind speed Vw=8m/s, obtaining CP=0.38. 4.2 Base Configuration: the NACA 0021 airfoil 23 Problems were found while making use of the LLT method, as calculations are highly sensitive to its vortex parameters. Parametres have to be adapted to each user case and turbine. The programme developer validates the LLT in [21], analysing sensitivity to certain vortex modeling parametres. These tests (Figures 4.17 and 4.18) showed that a maximum wake age of 6 revolutions and a full wake age of 1 revolution is offers a good balance in terms of calculation time and accuracy. Figure 4.17 Effect of maximum wake age on CPaccuracy [21]. Figure 4.18 Effect of full wake age on CPaccuracy [21]. Wake reduction will also be analysed. Wake reduction intends to reduce computation time by truncating the vortices with the lowest circulation. Simulations were undertaken on QBlade and as can be seen on Figure 4.19, the best results can be seen for vortex reduction 0.2 and 0.15, in terms of stability and convergence. Even though the results for vortex reduction nearest to 0 are the most accurate, it can be seen that they fluctuate and do not converge in a meaningful manner within the simulation time (10 revolutions), so for a preliminary analysis higher wake reduction proportions are more appropriate. Even though the LLT method uses more sophisticated calculations, it depends heavily on a well-calculated polar, which is a problem for VAWT, as θ ranges from 0 to 2 π . These CP values calculated by the DMS and LLT methods should just taken as trend, serving as an orientation through this early stage of the study. The reason is that typical CP values for Darrieus turbines are usually in the 0.3-0.4 range. This sort of discrepancy between results yielded by LLT-DMS theories and CFD could also be seen on Figure 2.5. For the aforementioned parameters, LLT analysis yields CP=0.383 (Figures 4.20 and 4.21), which is closer than DMS to the expected value to be given by CFD, alas, lower. 24 Chapter 4. Airfoil Analysis Figure 4.19 CPagainst wake reduction, clean configuration, JavaFoil. Figure 4.20 CPfor λ CPmax at final status, clean configuration (LLT), JavaFoil. Figure 4.21 CPfor λ CPmax vs time, clean configuration (LLT), JavaFoil. 4.3 The Gurney flap 25 4.3 The Gurney flap The Gurney flap is a structurally simple device, made up of a short strip, which is placed perpendicularly to the chord line of an airfoil, near the trailing edge. This device was first developed by Gurney for motorsports applications, being later on used for aerospace operation. The Gurney flap is said to break the Kutta boundary flow condition, as it creates a pair of contrarotating vortices in the immediate vicinity of the flap, which makes for an altered trailing edge circulation [22]. Liebeck studied the device and found best results with a height comprising 1.25 % of the chord, as a compromise between lift and drag for the studied airfoil. Drag increases exponentially for heights greater than 0.02c. Multiple cases were studied: • 2%c height, placed 90%c from the front of the airfoil: 90200 • 2%c height, placed 95%c from the front of the airfoil: 95200 • 1.25%c height, placed at the rear of the airfoil (100%c): 100125 • 1.25%c height, placed 95%c from the front of the airfoil: 95125 Figure 4.22 Gurney flap-modified NACA0021 airfoil. The polars (obtained through JavaFoil, as the XFoil method cannot converge satisfactorily under the given geommetry) and the values calculated by QBlade are illustrated by Figures 4.23, 4.24, and 4.25. Figure 4.23 Clagainst Cd, Gurney flap-modified airfoil. 26 Chapter 4. Airfoil Analysis Figure 4.24 Clagainst α , Gurney flap-modified airfoil. Figure 4.25 CPagainst λ , Gurney flap-modified airfoil. Even though the calculated maximum CP remains similar to the clean configuration, this may be due to the nature of the calculation of the polar and its further computation on QBlade, as the DMS theory does not take dynamic stall into consideration, and the performance improvement that such devices would provide. Nonetheless, it can be seen that λ CPmax is lower than for the clean case, and that the curve is overall flatter. This might mean that the effects of dynamic stall are no longer as strong. Out of all the tested Gurney flap configurations, 2% heights performed best. A 0.95c placement can be speculated to be the most beneficial, as the CP decreases comparatively more slowly, while rising faster in the low λ range. Thus, the Gurney flap-modified 95200 airfoil will be used for further analysis. The plot details of such an airfoil are summarised in Appendix 4.2. 4.3 The Gurney flap 27 Meanwhile, Lifting Line Theory analysis yielded CP=0.461 , as per Figures 4.26 and 4.27. This would imply a performance improvement of 20.3% , relative to the already studied clean configuration case. This is the largest improvement out of all devices studied. Figure 4.26 CPfor λ CPmax , Gurney flap-modified airfoil (LLT). Figure 4.27 CPfor λ CPmax vs time, Gurney flap-modified airfoil (LLT). 5 CFD Simulation In order to verify the quality of the calculations already undertaken using the Lifting Line Theory, the next step is to use ANSYS Fluent to apply Computational Fluid Dynamics to the designed profiles. ANSYS Fluent was chosen for this project because it offers a comparatively shorter learning curve, as it provides the user with a simple graphical interface and configuration wizards. 5.1 Simulation Procedure 5.1.1 Design Modelling It is understood that to simulate the rotation of the turbine about a z-axis, the overall mesh needs to be decomposed into a static part and a sliding mesh, which rotates about the desired axis. Lanzafame et al. [ 14 ] showed that a sliding ring mesh containing the 3 blades, and of the least thickness possible would provide the most stable results. A sliding ring also simulates a VAWT closest to the experimental situation. Conversely, a larger domain does not increase the quality of the simulation or the validity of the result when compared to the domain size to be used in this project, but just increases computation time. In Figure 5.1 a detailed view of the simulation domain used in this project can be seen. Figure 5.1 Simulation domain, detailed view. 1. The first step was to import the airfoil data from MS Excel to CATIA V5 using built-in macros, and then turn the sketch into a solid part. This way, CATIA V5 can calculate the inertia parameters, of which the centre of gravity coordinates are most important. These data, as applied to the clean configuration, are as per Figure 5.2. 35 36 Chapter 5. CFD Simulation Figure 5.2 Centre of gravity coordinates, clean configuration. 2. The first blade was then put into place using these coordinates, aligning the CoG with the centre of coordinates and then displacing the blade 1500mm along the y-axis (the turbine radius). This process produces the result shown on Figure 5.3. Figure 5.3 Displaced single blade, clean configuration. 3. In order to model the two remaining blades, a rotational pattern was set up, offsetting the blades 120 degrees. Turbine dimensions are measured from the turbine z-axis to the blade center line. The ring (0.2*Diameter thick around the turbine radius circumference) was created next using the dimensions as per Figure 5.4, and using boolean operations it was filled as a solid, and the blades emptied. 4. The last steps were to sketch and fill the inner and outer static domains, with dimensions as shown on Figure 5.1. Subsequently, the model could be said to be ready for implementation, and the CATIA product file was imported into ANSYS DesignModeler. The three parts (ring, outer domain and inner domain) were imported separately, and then joined in DesignModeler using the Form New Part command. As a last step in the whole modelling phase, the Shared Topology Method was set as None within the individual parts and as Automatic for the general part, and the material as Fluid. The shared topology method is key for this study, as setting it to none allows for interfaces to be formed between the different meshing components. The model is seen as per Figure 5.4 on DesignModeler. 5.1 Simulation Procedure 37 Figure 5.4 Simulation domain, DesignModeler. 5.1.2 Meshing Meshing plays a crucial role in a CFD simulation, as a bad mesh can make the simulation not be representative. Thus, special care needs to be taken for meshing methods and their quality indicators. First of all, Named Selections have to be set up, so that boundary conditions are easier to deal with. The Named Selections set up for this study are: • The velocity inlet upstream of the turbine, named ”inlet”, • The pressure outlet downstream of the turbine, named ”outlet”, • The 3 blades of the turbine, named ”blade1”, ”blade2” and ”blade3” • The outer fluid domain, named ”exterior”, • The inner fluid domain, named ”interior”, • The sliding ring area, named ”corona”, • Within these 3 domain areas, two interfaces, ”interface1” between ”exterior” and ”corona”, and ”interface2” between ”corona” and ”interior”. Each of the interfaces is made up of two Named Selections, for each side of the interface belonging to each fluid domain. The following steps have been taken for meshing: • Edge Sizing for all 3 blades: ElementSize =10−3m . This allows for higher resolution along the blades, so that blade contours are correctly represented. • Meshing Method: All Triangles Method. Either square or triangular mesh is implemented in ANSYS Fluent. A triangular mesh needs less time to be produced, which is key when undertaking CFD simulations with a personal computer. • Inflation layers around all 3 blades: FirstLayerT hickness =3.82∗10−5m . This parameter adds parallel layers to the blades, which model the boundary layer next to them. The given first layer thickness was calculated using a yplus calculator, knowing that y+<1 (or in its range) for the turbulence models used for simulation. Furthermore, the inflation was set to 20 layers at a growth rate of 1.1, also used on [14]. A detail of said inflation layers can be seen on Figure 5.5. • Face sizing: changes the mesh coarseness within the chosen domain area. Element size was chosen 0.05 m for the ”corona” area and 0.1 m for ”exterior”. Otherwise, the model was meshed as per the general element size, which was set at 0.5 m. The overall mesh thus would look similar to that on Figure 5.6. A closer look into the turbine area is shown by Figure 5.7. 38 Chapter 5. CFD Simulation Figure 5.5 Inflation layers around a blade. Figure 5.6 General mesh. Figure 5.7 Turbine core mesh, detailed view. Once the meshing process is finished, it is then necessary to check the quality of the results. Out of many metrics, two parameters are most influential in the successful undertaking of the simulation: 5.1 Simulation Procedure 39 • Skewness: measures how far from the ideal (equilateral triangle or equiangular quad) each mesh element is. Thus, null skewness implies best quality, and Skewness =1 implies a degenerate mesh element. [26] • Orthogonal quality: measures the degree of alignment between the edge normal vectors and the vectors from the face centroid to the edge centroids, the value for the mesh element being the minimum of this product for all edges. Thus, OrthogonalQuality =1 implies best quality (vectors completely aligned in all directions). 5.1.3 Setup and Simulation Once the mesh has been verified, the next step is to set up the simulation. For the simulations in this study, a personal computer featuring an Intel i7-4510U processor and 16 GB RAM was used. Therefore, in order to use the full capabilities of the processor, the simulation was set up as Double Precision and using both cores. General Setup Due to the nature of the simulations in this study, the solver was chosen as pressure-based and transient in time. Viscous Model The choice of a turbulence model for a 2D VAWT simulation is still not straightforward. Early CFD simulations were undertaken using κ − ε turbulence models, which were not sucessful in approximating flow behaviour near the wall. Thus, currently studies are prominently using κ − ω SST and SST Transition models [ 13 ][ 14 ]. In this study an uncalibrated Transition SST model will be used, as it enables the simulation of laminar sub-layers within the laminar flow, being a progression from the κ − ω SST model. Calibration of the model is out of scope for this study, as it needs extensive experimental data. Nonetheless, the uncalibrated Transation SST model proved to be sufficiently accurate for this study. Cell Zone Conditions In Cell Zone Conditions is where the sliding mesh concept becomes important. As such, the ”corona” region was set to have a rotating mesh motion of 13.3 rad/s (with an inlet velocity of 8 m/s amounting to TSR =2.5 , approximately the optimal TSR for all turbines). This motion is relative to the static mesh zone, this static mesh zone being made up of both ”exterior” and ”interior”. Boundary Conditions The following boundary conditions were set up in this study: • Inlet Velocity: the wind velocity. Magnitude of 8 m/s, Turbulent Viscosity of 5% and Turbulent Viscosity Ratio of 10. These last two values are default values on ANSYS Fluent 19.2. • Blade Rotation: for ”blade1”, ”blade2” and ”blade3”, the walls must be modelled as Moving Walls, with a rotational motion of 0 rad/s with respect to the adjacent cell zone ”corona”. This step is very important, as otherwise the blades are not modelled appropriately. The effect of this mistake can be seen on 5.3.1. Mesh Interfaces In this step, the Named Selections ”interface1i” and ”interface1e” are turned into ”interface1”, and ”interface2i” and ”interface2e” into ”interface2”. Reference Values These Reference Values are used for adimensionalisation of the variables calculated by the Fluent solver, which for this study is most important for the calculation of CP . The Reference Values used are as per Figure 5.8. 40 Chapter 5. CFD Simulation Figure 5.8 Reference Values, ANSYS Fluent. Solution Method The solution method used, also used by [14] is as follows: • Scheme: SIMPLE • Gradient: Least Squares Cell Based • Pressure: Second Order • Momentum: Second Order Upwind • Turbulent Kinetic Energy: First Order Upwind • Specific Dissipation Rate: First Order Upwind • Intermittency: First Order Upwind • Momentum Thickness Re: First Order Upwind • Transient Formulation: Second Order Implicit Monitors and Reports As CP is the objective variable in this study and it is directly related with the Cm as CP= λ Cm , a monitor and a report were created for Cm , using the z-axis from (0,0,0) as moment source for ”blade1”, ”blade2” and ”blade3”. Animations In order to illustrate the results and also the convergence of the simulation and the different flow situations during a revolution, a velocity contour animation was created, which created one image every 5 timesteps. Timestepping, Convergence and Solution Initialisation As a result of [ 27 ], Rezaeiha et al. found that for convergent results, a simulation encompassing 20 turbine revolutions produced acceptable results (a 2.41% difference compared to simulating 100 turbine revolutions). Furthermore, azimuthal increments of 0.5◦ also provided sufficient accuracy (less than 0.5% difference compared to 0.1◦per timestep). Under these premises, there thus being 720 timesteps in a revolution, each simulation took around 5 days in this study. Convergence criteria were set as 10−5 and maximum 50 timesteps. A hybrid solution initialisation was undertaken. 5.2 Results: Clean Configuration 41 General Parameter Summary For further clarification, an overview of the most relevant parameters used on the simulations is show on the following Table. Table 5.1 Key simulation parameters. Turbine radius (R) 1.5 m Turbine height 1m (2D study, unitary height) Number of blades 3 Inflow velocity (V∞) 8 m/s Tip Speed Ratio (TSR- λ ) 2.5 (approximately the optimum) Turbulence model used Uncalibrated Transition SST Simulation length 20 revolutions Azimuth discretisation 0.5◦ Timestep size 6.5∗10−4ms 5.2 Results: Clean Configuration Performing a similar analysis to the one carried out on QBlade for the base configuration and both simulated passive devices, the data obtained by the moment coefficient monitor was processed on MATLAB. Figure 5.9 shows the moment coefficient over time, and Figure 5.10 the moment coefficient, averaged the whole timespan. Nonetheless, it is more representative to show the evolution of Cm over the last revolution, as the values will be nearly converged at this time range. Thi evolution, as well as the average Cm value over this last revolution, is shown in Figure 5.11. Figure 5.9 Cmover time, Clean configuration. 42 Chapter 5. CFD Simulation Figure 5.10 Averaged Cmover time, Clean configuration. Figure 5.11 Cmover the last full revolution and average, Clean configuration. As can be seen in Figure 5.11, the average power coefficient over the last full revolution is CP= λ Cm= 0.288. It is worth noting, especially upon analysis of the later two cases which will be studied, that there is no apparent regularity on blade momentum for this clean configuration. It has been observed that for this blade configuration, dynamic stall occurs in an uncontrolled, inconsistent manner in every rotation, which is the reason why the relative minima on the Cm chart are not regular. Thus, it is demonstrated why it is necessary to study and control the dynamic stall phenomenon further for a VAWT, as turbine performance is highly sensitive. 5.2 Results: Clean Configuration 43 Figure 5.12 Velocity contour for a NACA0021-based turbine at its least favourable azimuth. Figure 5.13 Velocity contour for a NACA0021-based turbine at its most favourable azimuth. In order to illustrate the inconsistencies caused by dynamic stall, it is then worth including the velocity contour details for two different moments at which a blade is at its least favourable azimuth. Note the different velocity situations on Figures 5.14 and 5.15. The velocity contour at its most favourable azimuth (Figure 5.16) remains constant for all 3 blades during the permanent regime. Figure 5.14 Velocity contour for a NACA0021-based turbine at its least favourable azimuth, blade detail (1). In these two figures it is clear that the blade experiences dynamic stall , as the blade is stalled only for a certain azimuth range. Even though the situations illustrated by both figures are similar, maximum velocity 50 Chapter 5. CFD Simulation Figure 5.30 Cmover the last full revolution and average, Stepped airfoil. As can be seen in Figure 5.30, the average power coefficient over the last full revolution is CP= λ Cm= 0.272 , which is 5.56% lower than the clean configuration, thus representing a worsening in the operational conditions. When compared to the clean configuration, it can be seen that there is no apparent improvement in flow separation or stall avoidance , which is what produces negative Cm values for the least favourable azimuth. Lower average Cm values relative to the clean configuration may also be caused because of the drag increase produced by the step at the most favourable azimuth, which reduces the amplitude between the minima and maxima in the Cmcurve. Figure 5.31 Velocity countour for the Stepped airfoil at its least favourable azimuth. 5.4 Results: Stepped Airfoil Modification 51 Figure 5.32 Velocity countour detail for the Stepped airfoil at its most favourable azimuth. Figure 5.33 Velocity countour for the Stepped airfoil at its most favourable azimuth. Figure 5.34 Velocity countour detail for the Stepped airfoil at its most favourable azimuth. Upon analysis of the velocity contours 5.33 and 5.34, it can be asserted that flow separation for the least favourable azimuth is not reduced by this step configuration, when compared to the base case. Moreover, at its most favourable azimuth, the step creates some additional drag, which is why the maxima on the Cm curve over its last revolution are lower than for the clean configuration. The test shows that a Stepped airfoil with the used placement and sizing can be asserted not to be of any benefit on the performance of the given turbine for the given operational point. Nonetheless, as was already discussed during the QBlade simulation phase, the step provides an improved performance behaviour within the low-TSR range, so such a configuration might be more beneficial for turbines where the optimal CPis obtained for lower TSRs. 52 Chapter 5. CFD Simulation 5.5 Performance Summary and Comparison with previous Studies For clarity, a summary of the resulting performance data obtained is as follows: Table 5.2 Performance summary. Turbine Configuration Power Coefficient (CP) value Improvement (%) Clean configuration 0.297 Base config. Gurney flap 0.43 44.78% Stepped airfoil 0.28 -5.56% In addition, for further illustration and verification, the results obtained in this project will be compared to those obtained through experimental and computational studies by other authors. As no study on the implementation of Gurney flaps and steps on the NACA 0021 airfoil is available, results on similar airfoils were investigated. 5.5.1 Clean configuration In 2015, Dossena et al. [ 28 ] undertook extensive experimental study on a similarly sized turbine, featuring three straight blades based on a NACA 0021 airfoil. The main sizes of this turbine are as follows: Table 5.3 Dossena et al. experimental turbine. Turbine radius (R) 1.028 m Turbine height (H) 1.5 m Blade chord (c) 0.086 m Optimal TSR approx. 2.5 As can be seen on Table 5.3, this turbine is operationally similar to the clean configuration in this study, both because of its sizing and its operational optimum (see Figure 5.31). Dossena et al. took aerodynamic measurements 0.75D and 1.5D downstream from the turbine, and 0.75D upstream from the turbine. The results of the measurements are shown in Figure 5.31. Figure 5.35 Clean configuration experimental results [28]. 5.5 Performance Summary and Comparison with previous Studies 53 The optimal CP obtained experimentally by Dossena et al. was λ =0.28 . This is very close to the result obtained numerically in this study ( 0.297 ). The higher result can be attributed to it being a 2D simulation, as typical 3D effects and tower interference are not taken into consideration. 5.5.2 Gurney flap The addition of a Gurney flap has not yet been studied extensively in a wind tunnel environment, in its application for vertical-axis wind turbines. Nonetheless, several computational results are available. One of such studies was undertaken by Ismail et al. [ 29 ], using a 1.25% chord length Gurney flap, with a dimple right before it (see Figure 5.32), applied to a NACA0015 airfoil. Ismail. et al used the dimple to reduce wake intensity, and thus reduce the drag component on the airfoil. Figure 5.36 Detail of a 1.25%c Gurney flap-equipped NACA0015 airfoil with a dimple [29]. Using the κ − ω SST turbulent model, Ismail et al. found increases in torque relative to the clean configuration were of 35% in the steady case and 40% in the transient case , approximately. The evolution of the tangential forces on the airfoil can be seen on Figure 5.33. Thus, this improvements found by Ismail et al. are similar to the improvements found in this project (44.78%). Figure 5.37 Tangential forces on a NACA0015 airfoil with a 1.25%c Gurney flap [29]. 5.5.3 Stepped airfoil For the stepped or notched airfoil, efforts have been concentrated in its self-starting capabilities, alas, improved behaviour for lower Tip Speed Ratios. Also based on a NACA0021 airfoil, Kumar et al. [ 30 ] undertook an experimental study on a notched airfoil, which had a 19%c deep, 50%c long step (see Figure 5.38). The characteristics of the general turbine tested were as follows. 54 Chapter 5. CFD Simulation Table 5.4 Kumar et al. experimental turbine. Turbine radius (R) 0.2 m Turbine height (H) 0.3 m Blade chord (c) 0.1 m Optimal TSR range approx. 1-1.5 Figure 5.38 Notched airfoil, Kumar et al. experimental study [30]. The main concept behind such an airfoil is the formation of a trapped vortex bubble inside of the airfoil, for a certain azimuth range. For example, Kumar et al. found that the bubble started forming for an azimuth of 30◦, and was strongest for an azimuth of 45◦. This pattern would repeat every 90◦. It was found that such a configuration does improve low-TSR performance and also improves the maximum CP , but only for low wind speeds (0-8m/s). This sensitivity to wind speed is shown on Figure 5.39. It can be seen that for higher Reynolds numbers, a notched airfoil is no longer beneficial for turbine operation. Figure 5.39 Performance variability with Reynolds number, Kumar et al. [30]. Nonetheless, this notched airfoil seems most interesting for low-cost, low-power domestic applications at sites which are not especially windy. Its design can be likened to a hybrid between a Darrieus and a Savonius turbine (due to the notch), and thus has better startup capabilities , as drag is more important within its operation. This importance of drag resides in the shear layer between the notch and the airfoil surface. 6 Inconsistencies Found During Simulation 6.1 Torque instability within the clean configuration As has already been discussed, during the clean configuration the torque coefficient was found to be highly unstable when compared to all other simulations, and even though a general trend could be deduced, maximum and minimum values would change for every simulation. It can be attributed to two main causes: • Insufficient simulation time. Even though for the case studied on [ 27 ] 20 revolutions attained a stable result which did not differ significantly to simulating 100 revolutions, a longer simulation size might have been needed for the studied turbine to converge to a higher level. • Dynamic stall control. It might be that because the clean configuration does not have any stall control device implemented, stability cannot be reached, as dynamic stall is a highly unstable process. Moreover, this has to be added to the fact that in a CFD simulation, convergence is limited. For instance, after 50 iterations in a timestep, the continuity relative error stood at 10−3 for the torque maxima and minima and 10−5for the intermediate values. 6.2 Asymmetrical Gurney flap simulation When the first simulation for the Gurney flap turbine was completed, it was observed that there was an asymmetry in the animation, and that one of the blades worked differently to the other two, especially when rotating past the most negative point. Later on it was detected that this was due to not setting up the Cell Zone Condition correctly for this single blade. The moving wall was set as translational instead of rotational, also relative to the ”corona” cell zone. 55 56 Chapter 6. Inconsistencies Found During Simulation Figure 6.1 Incongruencies found during Gurney flap turbine CFD simulation, Cmover time graph. Figure 6.2 Incongruencies found during Gurney flap turbine CFD simulation, averaged Cm over time graph. What is most noticeable about the figure below is that its best fit line would not be of the form C(1−e−kt ) like in the other cases, but it has a maximum and then its asymptote towards c (the ”final” Cm value as can be said), alas, it could be likened to C(1−e−at −e−bt ) . Furthermore, on Figure 6.3 it can be seen that even though all three maxima and minima should be the same or at least similar over the last revolution, there is a maximum and minimum which is different to all others. This can be attributed to the bad setup on one of the three blades, while the other two behave similarly as can be seen. As can be seen on Figure 6.3, the average power coefficient over the last full revolution is CP= λ Cm= 0.473 , which is an improvement of 64.2% relative to the clean configuration. Nonetheless, this result is very much invalid, due to the causes which have already been discussed. 6.2 Asymmetrical Gurney flap simulation 57 Figure 6.3 Cmover the last full revolution and average for the Gurney flap-equipped turbine. Figure 6.4 Velocity contour at the least favourable azimuth for the Gurney flap-equipped, wrongly-configured blade. Figure 6.5 Velocity contour at the least favourable azimuth for the two remaining Gurney flap-equipped blades. Upon analysis of Figures 6.4 and 6.6, which shows a detail on the velocity contour of the wrongly set up blade at its least favourable azimuth, it can be seen why it is not realistic, as there are no signs of flow separation, unlike it would happen for the remaining two blades. The velocity contours for the other two blades would be as per Figures 6.5 and 6.7. The situation at the most favourable azimuth would be similar for all 3 blades, see Figure 6.8. 58 Chapter 6. Inconsistencies Found During Simulation Figure 6.6 Detailed velocity contour at the least favourable azimuth for the Gurney flap-equipped, wronglyconfigured blade. Figure 6.7 Velocity contour at the least favourable azimuth for the two remaining Gurney flap-equipped blades. Figure 6.8 Velocity contour at the most favourable azimuth, as applicable for all Gurney-flap equipped blades. 7 Conclusion and Recommendations for Further Research In this project, it has been verified that for a turbine based on the NACA 0021 airfoil, the effects of dynamic stall on its performance can be most significantly mitigated using Gurney flaps, improving its performance over the clean configuration by 40%. On the other hand, Liebeck’s simple stepped airfoil format has proven to be highly ineffective and even negative for performance improvement. Thus, a way forward to work on this project is to develop innovative passive devices, based on Gurney flap principles. One such concept was the implementation of a dimple on the forward side of the flap. This concept was discussed by [29], an overview of which was given in Chapter 5.5.2 of this project. Other main points of this project have been: • Numerical methods for airfoil polar data generation are highly unstable for non-standard geometries, and results can vary significantly between different numerical methods (see XFoil and JavaFoil). • Using these numerically-generated polars, QBlade obtains results whose quality is uncertain. Nonetheless, even under these circumstances, QBlade can prove effective for preliminary analysis and for determining differences between different turbines and airfoils, due to its low computational costs and model complexity, when compared to DMS simulations. • Also using numerically-generated polars, DMS obtained unrealistic results, for instance obtaining very high CP values near the theoretical Betz’s limit. DMS was only useful for a very low computational cost analysis of the optimal operational point. • The uncalibrated Transition SST turbulent model used in this project has shown a behaviour close to experimental studies for the clean configuration. This turbulent model can therefore be verified for this configuration. Experimental studies are needed for the remaining two airfoils. 7.1 Recommendations for Further Research The main points related to this project which need to be addressed are further model verification tasks, as well as iterative experimental design: • An experimental wind tunnel study of the turbine, using particle image velocimetry (PIV), to obtain the torque behaviour in a steady state and verify the model. • Calibration of the Transition SST model using experimental data, and comparison with the results of this project. • Iterative experimental design, using the same base airfoil, and implementing different passive devices on each test, for instance changing Gurney flap placement, shape, sizing. 59 66 Chapter B. MATLAB Code Compilation B.2 ANSYS Fluent result processing Code B.2 CFD post-processing. clear all clc load('step.out') data=step; figure(1); hold on; plot(data(:,3)); title('Cm values over time'); xlabel('Timesteps'); ylabel('Cm [-]'); hold off; for i=1:length(data(:,3)) mean1(i)=mean(data(1:i,3)); end figure(2); hold on; plot(mean1); title('Averaged Cm over time'); xlabel('Timesteps'); ylabel('Cm [-]'); hold off; figure(3); hold on; plot(data(13100:13830,3)); m=mean(data(13100:13830,3)); m=m*ones(720); plot(m); title('Cm values over time'); xlabel('Timesteps'); ylabel('Cm [-]'); hold off; figure(4); hold on; plot(mean1*2.5); title('Averaged Cm over time'); xlabel('Timesteps'); ylabel('CP [-]'); hold off; figure(5); hold on; plot(data(13100:13830,3)); m=mean(data(13100:13830,3))*2.5; m=m*ones(720); plot(m); title('CP values over time'); xlabel('Timesteps'); ylabel('CP [-]'); hold off; List of Figures 1.1 Evolution of yearly wind energy capacity installation 1 1.2 Savonius, Egg beater (type 1 Darrieus), and H-rotor (type 2 Darrieus) turbines, compared to a HAWT [6] 2 2.1 Relative velocity diagram on the level of a Darrieus turbine [9] 5 2.2 Aerodynamic forces on an airfoil [10] 7 2.3 CLagainst α for a DU 06-W-200 airfoil 8 2.4 α against θ for a VAWT airfoil 8 2.5 Influence of dynamic stall on the power coefficient 9 2.6 Turbulence models available in Fluent [15] 12 3.1 Influence of the Tip Speed Ratio on the Reynolds number 13 3.2 Influence of the Tip Speed Ratio on the power coefficient 14 4.1 XFoil settings 16 4.2 QBlade polar 16 4.3 Viterna-Corrigan extrapolation 16 4.4 Rotor design, Base configuration 17 4.5 DMS, configuration parameters 17 4.6 LLT, operational point 18 4.7 LLT, output, algorithm, wake parameters 18 4.8 NACA0021, clean airfoil 19 4.9 Clagainst Cd, clean configuration, XFoil 20 4.10 Clagainst α , clean configuration, XFoil 20 4.11 CPagainst λ , clean configuration, XFoil 20 4.12 Ctagainst θ , clean configuration, XFoil 21 4.13 Clagainst Cd, clean configuration, JavaFoil 21 4.14 Clagainst α , clean configuration, JavaFoil 22 4.15 CPagainst λ , clean configuration, JavaFoil 22 4.16 Ctagainst θ , clean configuration, JavaFoil 22 4.17 Effect of maximum wake age on CPaccuracy [21] 23 4.18 Effect of full wake age on CPaccuracy [21] 23 4.19 CPagainst wake reduction, clean configuration, JavaFoil 24 4.20 CPfor λ CPmax at final status, clean configuration (LLT), JavaFoil 24 4.21 CPfor λ CPmax vs time, clean configuration (LLT), JavaFoil 24 4.22 Gurney flap-modified NACA0021 airfoil 25 4.23 Clagainst Cd, Gurney flap-modified airfoil 25 4.24 Clagainst α , Gurney flap-modified airfoil 26 4.25 CPagainst λ , Gurney flap-modified airfoil 26 4.26 CPfor λ CPmax , Gurney flap-modified airfoil (LLT) 27 4.27 CPfor λ CPmax vs time, Gurney flap-modified airfoil (LLT) 27 4.28 T-strip-modified NACA0021 airfoil 28 67 68 List of Figures 4.29 Detailed view, T-strip-modified airfoil 28 4.30 Clagainst Cd, T-strip-modified airfoil 29 4.31 Clagainst α , T-strip-modified airfoil 29 4.32 CPagainst λ , T-strip-modified airfoil 30 4.33 CPfor λ CPmax , T-strip-modified airfoil 30 4.34 CPfor λ CPmax vs time, T-strip-modified airfoil 30 4.35 Stepped NACA0021 airfoil 31 4.36 Clagainst Cd, Stepped airfoil 31 4.37 Clagainst α , Stepped airfoil 31 4.38 CPagainst λ , Stepped airfoil 32 4.39 CPfor λ CPmax , Stepped airfoil (LLT) 32 4.40 CPfor λ CPmax vs time, Stepped airfoil (LLT) 33 5.1 Simulation domain, detailed view 35 5.2 Centre of gravity coordinates, clean configuration 36 5.3 Displaced single blade, clean configuration 36 5.4 Simulation domain, DesignModeler 37 5.5 Inflation layers around a blade 38 5.6 General mesh 38 5.7 Turbine core mesh, detailed view 38 5.8 Reference Values, ANSYS Fluent 40 5.9 Cmover time, Clean configuration 41 5.10 Averaged Cmover time, Clean configuration 42 5.11 Cmover the last full revolution and average, Clean configuration 42 5.12 Velocity contour for a NACA0021-based turbine at its least favourable azimuth 43 5.13 Velocity contour for a NACA0021-based turbine at its most favourable azimuth 43 5.14 Velocity contour for a NACA0021-based turbine at its least favourable azimuth, blade detail (1) 43 5.15 Velocity contour for a NACA0021-based turbine at its least favourable azimuth, blade detail (2) 44 5.16 Clean configuration, most favourable azimuth, blade detail 44 5.17 Detail on the most and least favourable azimuths 44 5.18 Cmover time, Gurney flap, second simulation 45 5.19 Averaged Cmover time, Gurney flap, second simulation 45 5.20 Cmover the last full revolution and average, Gurney flap, second simulation 46 5.21 Velocity countour for the Gurney flap-equipped airfoil at its least favourable azimuth 46 5.22 Velocity countour detail for the Gurney flap-equipped airfoil at its least favourable azimuth 47 5.23 Velocity countour for the Gurney flap-equipped airfoil at its most favourable azimuth 47 5.24 Static pressure countour for the Gurney flap-equipped airfoil at its most favourable azimuth 47 5.25 Static pressure contour detail for the Gurney flap-equipped airfoil at its most favourable azimuth 48 5.26 Static pressure countour for the Gurney flap-equipped airfoil at its least favourable azimuth 48 5.27 Static pressure countour detail for the Gurney flap-equipped airfoil at its least favourable azimuth 48 5.28 Cmover time, Stepped airfoil 49 5.29 Averaged Cmover time, Stepped airfoil 49 5.30 Cmover the last full revolution and average, Stepped airfoil 50 5.31 Velocity countour for the Stepped airfoil at its least favourable azimuth 50 5.32 Velocity countour detail for the Stepped airfoil at its most favourable azimuth 51 5.33 Velocity countour for the Stepped airfoil at its most favourable azimuth 51 5.34 Velocity countour detail for the Stepped airfoil at its most favourable azimuth 51 5.35 Clean configuration experimental results [28] 52 5.36 Detail of a 1.25%c Gurney flap-equipped NACA0015 airfoil with a dimple [29] 53 5.37 Tangential forces on a NACA0015 airfoil with a 1.25%c Gurney flap [29] 53 5.38 Notched airfoil, Kumar et al. experimental study [30] 54 5.39 Performance variability with Reynolds number, Kumar et al. [30] 54 6.1 Incongruencies found during Gurney flap turbine CFD simulation, Cmover time graph 56 6.2 Incongruencies found during Gurney flap turbine CFD simulation, averaged Cmover time graph 56 6.3 Cmover the last full revolution and average for the Gurney flap-equipped turbine 57 6.4 Velocity contour at the least favourable azimuth for the Gurney flap-equipped, wrongly-configured blade 57 List of Figures 69 6.5 Velocity contour at the least favourable azimuth for the two remaining Gurney flap-equipped blades 57 6.6 Detailed velocity contour at the least favourable azimuth for the Gurney flap-equipped, wronglyconfigured blade 58 6.7 Velocity contour at the least favourable azimuth for the two remaining Gurney flap-equipped blades 58 6.8 Velocity contour at the most favourable azimuth, as applicable for all Gurney-flap equipped blades 58 List of Tables 4.1 Turbine characteristics 15 5.1 Key simulation parameters 41 5.2 Performance summary 52 5.3 Dossena et al. experimental turbine 52 5.4 Kumar et al. experimental turbine 54 A.1 NACA 0021 airfoil 61 A.2 NACA 0021 airfoil, Gurney flap h=0.02c, x=0.95c 62 A.3 NACA 0021 airfoil, T-strip 63 A.4 NACA 0021 airfoil, Stepped 64 71 List of Codes B.1 Polar polynomial best fit 65 B.2 CFD post-processing 66 73 Bibliography [1] European Wind Energy Association, ‘‘Wind energy installed capacity.’’ http://www.windeurope.org, 2017. [2] Asociación Española para la Energía Eólica, ‘‘La eólica en españa.’’ https://www.aeeolica.org/es/ sobre-la-eolica/la-eolica-espana, 2017. [3] Agencia Andaluza de la Energía, ‘‘Energía eólica marina.’’ https://www.agenciaandaluzadelaenergia.es, 2016. [4] A. Grant, C. Johnstone, and N. Kelly, ‘‘Urban wind energy conversion: The potential of ducted turbines,’’ Renewable Energy, vol. 33, no. 6, pp. 1157–1163, 2008. cited By 61. [5] F. Toja-Silva, A. Colmenar-Santos, and M. Castro-Gil, ‘‘Urban wind energy exploitation systems: Behaviour under multidirectional flow conditions—opportunities and challenges,’’ Renewable and Sustainable Energy Reviews, vol. 24, pp. 364 – 378, 2013. [6] P. Mehrpooya, Improvement of Vertical-axis Wind Turbine Performance Via Turbine Coupling. PhD thesis, Illinois Institute of Technology, 2014. [7] M. M. A. Bhutta, N. Hayat, A. U. Farooq, Z. Ali, S. R. Jamil, and Z. Hussain, ‘‘Vertical axis wind turbine – a review of various configurations and design techniques,’’ Renewable and Sustainable Energy Reviews, vol. 16, no. 4, pp. 1926 – 1939, 2012. [8] I. O. N. PARASCHIVOIU, ‘‘Double-multiple streamtube model for studying vertical-axis wind turbines,’’ Journal of Propulsion and Power, vol. 4, pp. 370–377, Nov. 2018. [9] P. Tchakoua, R. Wamkeue, M. Ouhrouche, T. A. Tameghe, and G. Ekemb, ‘‘A new approach for modeling darrieus-type vertical axis wind turbine rotors using electrical equivalent circuit analogy: Basis of theoretical formulations and model development,’’ Energies, vol. 8, no. 10, pp. 10684–10717, 2015. [10] Wikipedia, ‘‘Airfoil lift and drag.’’ https://upload.wikimedia.org/wikipedia/commons/a/a2/Airfoil_ lift_and_drag.jpg. [11] D. Marten, G. Pechlivanoglou, C. N. Nayeri, and C. O. Paschereit, ‘‘Nonlinear lifting line theory applied to vertical axis wind turbines: Development of a practical design tool,’’ Journal of Fluids Engineering, vol. 140, no. 2, p. 021107, 2018. [12] D. Marten, A. Bianchini, G. Pechlivanoglou, F. Balduzzi, C. N. Nayeri, G. Ferrara, C. O. Paschereit, and L. Ferrari, ‘‘Effects of airfoil’s polar data in the stall region on the estimation of darrieus wind turbine performance,’’ Journal of Engineering for Gas Turbines and Power, vol. 139, no. 2, p. 022606, 2017. [13] M. Kortleven, ‘‘Simulation verication and optimization of a vertical axis wind turbine using cfd,’’ Master’s thesis, TU Delft, 2016. 75