DUCT OPTIMIZATION FOR A DUCTED VERTICAL AXIS HYDRO CURRENT TURBINE by MAHMOUD ALIDADI B.Sc., The Sharif University of Technology, Iran, 1999 M.Sc., The Ferdowsi University of Mashad, Iran, 2001 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) June 2009 © Mahmoud Alidadi, 2009 ABSTRACT Turbines developed for the power production from the energy of tidal and river currents are commonly referred to as hydro current turbines. This research focuses on duct optimization for a ducted vertical axis hydro current turbine. As a first step, a discrete vortex method, referred to as “DVM”, is developed for the performance analysis of unducted vertical axis turbines. It is shown that the loads for an unstalled blade of a vertical axis turbine may be adequately approximated by a quasi-steady component and a pitching component. A new approach is developed for the calculation of nascent vortex position relative to the blade. Besides, a new approach is developed to take into account the effects of connection point between a blade and its supporting arm. Subsequently, a new numerical model, referred to as “BEM-DVM”, is developed for the performance calculation of ducted vertical axis turbines by combining a boundary element method referred to as “BEM” with the DVM model. Comparison of the results for unducted and ducted turbines showed that ducting significantly improves the power coefficient of a vertical axis turbine. In addition, the results showed that ducting decreases the torque fluctuations significantly, especially at high tip speed ratios. Based on the method developed for ducted turbines, a new approach is developed to calculate the effects of towing tank wall on experimental results. The results showed that the power coefficient of a turbine increases when it operates in a tank. This is especially true for a ducted turbine which has a higher blockage ratio than the unducted turbine in the tank. The results also showed that the amount of increase in power coefficient grows when the tank width decreases. An optimization study is conducted to find the duct shape that maximizes the output power of a vertical axis turbine. The optimized ducts results showed a higher power coefficient for the turbine relative to the results of experimentally tested ducts. It is also shown that there is an upper limit to the amount of increase in power coefficient due to ducting. ii TABLE OF CONTENTS ABSTRACT ........................................................................................................................ ii TABLE OF CONTENTS ................................................................................................... iii LIST OF TABLES ............................................................................................................. vi LIST OF FIGURES ........................................................................................................... vii LIST OF SYMBOLS, ABBREVIATIONS AND NOMENCLATURE ........................... xi ACKNOWLEDGMENTS................................................................................................ xiii CHAPTER ONE: INTRODUCTION ................................................................................. 1 1.1: Research objectives .................................................................................................. 1 1.2: Research methodology ............................................................................................. 2 1.3: Thesis outline ........................................................................................................... 3 1.4: Scope of research ..................................................................................................... 4 CHAPTER TWO: BACKGROUND .................................................................................. 5 2.1: Renewable energy technologies ............................................................................... 5 2.2: Vertical axis turbines ............................................................................................... 6 2.3: Flow dynamics in a vertical axis turbine ................................................................. 7 2.4: Advantages of the vertical axis turbines .................................................................. 7 2.5: Disadvantages of vertical axis turbines .................................................................... 8 2.6: Ducted vertical axis turbines .................................................................................... 8 2.7: Literature review .................................................................................................... 10 2.7.1: Numerical study of unducted vertical axis turbines ...................................... 10 2.7.1.1: Potential flow models........................................................................... 10 2.7.1.2: RANS modeling ................................................................................... 13 2.7.2: Numerical study of ducted vertical axis turbines .......................................... 14 2.7.3: Experimental study of unducted and ducted vertical axis turbines ............... 14 CHAPTER THREE: UNDUCTED VERTICAL AXIS TURBINE ................................. 17 3.1: Introduction ............................................................................................................ 17 3.2: Loads on a blade of vertical axis turbine ............................................................... 20 3.3: Development of DVM for vertical axis turbines ................................................... 25 3.3.1: Wake representation in DVM ........................................................................ 26 3.3.2: Calculation of nascent vortex position relative to foil ................................... 27 3.3.3: The velocity induced by a discrete vortex ..................................................... 29 iii 3.3.4: Angle of attack calculation for a blade of VAT ............................................ 31 3.3.5: Effects of connection point between supporting arm and blade .................... 32 3.3.6: Calculation of turbine power coefficient ....................................................... 33 3.3.7: Computational procedure for DVM............................................................... 34 3.4: RANS simulation for the calculation of foil characteristics .................................. 35 3.4.1: Grid generation for the RANS simulations ................................................... 35 3.4.2: Solver specifications for the RANS simulations ........................................... 36 3.4.3: Results of the RANS simulations .................................................................. 37 CHAPTER FOUR: DUCTED VERTICAL AXIS TURBINE ......................................... 38 4.1: Introduction ............................................................................................................ 38 4.2: Boundary element methods for internal-external flows ......................................... 39 4.2.1: General solution for potential flow problem ................................................. 40 4.2.2: Boundary conditions ...................................................................................... 45 4.2.2.1: Neumann problem ................................................................................ 45 4.2.2.2: Dirichlet problem ................................................................................. 45 4.2.3: Kutta condition .............................................................................................. 47 4.2.4: Modeling of the vortices in the flow field ..................................................... 48 4.2.5: Calculation of influence coefficients ............................................................. 51 4.2.5.1: Constant-strength source distribution .................................................. 52 4.2.5.2: Constant-strength doublet distribution ................................................. 53 4.3: Development of the BEM-DVM for ducted vertical axis turbines ........................ 54 4.4: A method for the calculation of tank wall effects on the turbine performance ..... 56 CHAPTER FIVE: DUCT OPTIMIZATION .................................................................... 57 5.1: Introduction ............................................................................................................ 57 5.1.1: Search-based methods ................................................................................... 58 5.1.2: Gradient-based methods ................................................................................ 58 5.2: Duct optimization for ducted vertical axis hydro current turbines ........................ 59 5.2.1: Duct representation ........................................................................................ 59 5.2.2: Objective function ......................................................................................... 60 5.2.3: Optimization constraints ................................................................................ 60 5.2.4: Optimization methodology ............................................................................ 61 CHAPTER SIX: RESULTS AND DISCUSSION............................................................ 62 6.1: Introduction ............................................................................................................ 62 iv 6.2: Experimental testing of a vertical axis turbine at UBC towing tank ..................... 62 6.3: Comparison of the numerical results with the UBC experimental results ............. 66 6.4: Experimental testing of a vertical axis turbine at IOT ........................................... 72 6.5: Comparison of the numerical results with the IOT experimental results .............. 75 6.6: A discussion on the sources of discrepancies ........................................................ 79 6.7: A parametric study on the performance of a VAT................................................. 80 6.7.1: Effects of TSR on the torque and the angle of attack (AOA) of a blade ....... 80 6.7.2: Ducting effects on the torque fluctuations of a turbine ................................. 81 6.7.3: Tank wall effects on the power coefficient of a turbine ................................ 82 6.8: Results of duct optimization .................................................................................. 84 6.8.1: Symmetric duct optimization results ............................................................. 84 6.8.2: Asymmetric duct optimization results ........................................................... 89 CHAPTER SEVEN: CONCLUSIONS AND RECOMMENDATIONS ......................... 93 7.1: Conclusions ............................................................................................................ 94 7.2: Original contributions ............................................................................................ 96 7.3: Recommendations for the future work................................................................... 96 7.3.1: Foil shape study ............................................................................................. 96 7.3.2: Inclusion of a dynamic stall model into the numerical models ..................... 97 7.3.3: Inclusion of a boundary layer method for duct calculations.......................... 97 7.3.4: Inclusion of three-dimensional effects........................................................... 97 7.3.5: Using a high-order boundary element method .............................................. 97 References ......................................................................................................................... 98 Appendix A: Gradient-based optimization methods ....................................................... 105 A.1: Sequential quadratic programming (SQP) ..................................................... 106 Appendix B: Introduction to spline ................................................................................. 107 B.1: Bezier curves .................................................................................................. 107 B.2: B-spline curves ............................................................................................... 108 v LIST OF TABLES Table 2.1: Available Davis et al. reports ................................................................................ 15 Table 6.1: The characteristics of the UBC-Turbine. ............................................................... 63 Table 6.2: CTF values for an unducted turbine and a ducted turbine at U∞=1.5 m/s. ............. 82 Table 6.3: The design variables and their bounds in symmetric duct optimization. ............... 85 Table 6.4: The design variables and their bounds at asymmetric duct optimization. ............. 89 vi LIST OF FIGURES Figure 2.1: A VAT with helical blades [9]................................................................................ 6 Figure 2.2: A VAT with straight blades. ................................................................................... 6 Figure 2.3: Schematic of circular movement for a blade element of VAT. .............................. 7 Figure 2.4: Streamtubes in an unducted turbine and a ducted turbine. ..................................... 8 Figure 2.5: Schematic of a ducted turbine................................................................................. 9 Figure 2.6: Single streamtube model [1]. ................................................................................ 11 Figure 2.7: Multiple streamtube model [1]. ............................................................................ 11 Figure 2.8: Schematic of double-multiple streamtube model [1]............................................ 12 Figure 2.9: Modeling a blade by a vortex filament [1]. .......................................................... 13 Figure 3.1: The local and global coordinate system................................................................ 20 Figure 3.2: Mapping of a flat plate to a circle. ........................................................................ 21 Figure 3.3: The loads on a foil are placed in its quarter-chord. .............................................. 24 Figure 3.4: A schematic of vertical axis turbine vs. DVM modeling of this turbine. ............. 25 Figure 3.5: Circulation around any contour remains constant in a potential flow field.......... 26 Figure 3.6: The wake is modeled by discrete vortices in DVM. ............................................. 26 Figure 3.7: The vortex sheet from θ1 to θ 2 . ............................................................................ 27 Figure 3.8: The equivalent discrete vortex in DVM. .............................................................. 27 Figure 3.9: The vortex sheet and its equivalent discrete vortex are along the chord. ............. 27 Figure 3.10: Velocity variation of a vortex. ............................................................................ 29 Figure 3.11: Vortex shedding from a foil. ............................................................................... 30 Figure 3.12: Velocity components for a blade. ....................................................................... 31 Figure 3.13: The blade connected to the arm at a point other than the mid-chord. ................ 32 Figure 3.14: The loads on a blade of VAT at the quarter-chord of blade. .............................. 33 Figure 3.15: The flowchart of computational procedure for DVM......................................... 34 vii Figure 3.16: Profile of turbulent boundary layer [1]. .............................................................. 35 Figure 3.17: Domain of solution in simulations...................................................................... 36 Figure 3.18: Fine grid around the foil in simulations. ............................................................. 36 Figure 3.19: The characteristics of NACA 63(4)-021............................................................. 37 Figure 4.1: Schematic of a ducted turbine............................................................................... 38 Figure 4.2: A body submerged in a potential flow.................................................................. 40 Figure 4.3: The potential and its gradient inside and outside of body. ................................... 44 Figure 4.4: Flow leaves the T.E. along the bisector line based on kutta condition. ................ 47 Figure 4.5: The body surface is represented by some panels. ................................................. 50 Figure 4.6: Panel i in the local coordinates of panel j . ........................................................... 51 Figure 4.7: Schematic of modeling for a ducted vertical axis turbine. ................................... 54 Figure 4.8: The flowchart of BEM-DVM. .............................................................................. 55 Figure 4.9: Schematic of a ducted turbine operating in a tank. ............................................... 56 Figure 5.1: Schematic of a symmetric duct. ............................................................................ 59 Figure 5.2: Schematic of an asymmetric duct. ........................................................................ 60 Figure 5.3: Flowchart of the optimization process.................................................................. 61 Figure 6.1: The towing tank facility at UBC [3]. .................................................................... 63 Figure 6.2: The turbine assembly (all dimensions are in inches) [3]. ..................................... 64 Figure 6.3: Gearbox drive-train configuration [3]. ................................................................ 64 Figure 6.4: Symmetric duct dimensions (all dimensions are in inches) [3]. ........................... 65 Figure 6.5: Schematic of ducted turbine in the tank (dimensions are in inches) [3]............... 66 Figure 6.6: A schematic drawing of the ducted turbine operating in a tank. .......................... 66 Figure 6.7: The profile and dimensions of the arms (all dimensions are in inches) [3]. ......... 67 Figure 6.8: Experimental Cp vs. TSR values for the supporting arms [3]. ............................. 67 Figure 6.9: Results of Cp vs TSR for unducted turbine at U∞=1.5 m/s. ................................. 68 viii Figure 6.10: The torque of unducted turbine at TSR=3 and U∞=1.5 m/s. .............................. 69 Figure 6.11: The torque of one blade of the unducted turbine at TSR=3 and U∞=1.5 m/s. .... 69 Figure 6.12: Results of Cp vs TSR for ducted turbine at U∞=1.5 m/s. ................................... 70 Figure 6.13: The torque of the ducted turbine at TSR=3 and U∞=1.5 m/s.............................. 71 Figure 6.14: The torque of one blade of the ducted turbine at TSR=3 and U∞=1.5 m/s. ........ 71 Figure 6.15: The asymmetric duct and its positions relative to the turbine. ........................... 72 Figure 6.16: The tank and the carriage in IOT. ....................................................................... 73 Figure 6.17: Main beams of carriage. ..................................................................................... 73 Figure 6.18: The frame designed for experiments [4]............................................................. 73 Figure 6.19: The turbine in its initial position relative to the duct. ......................................... 74 Figure 6.20: The turbine in its secondary position relative to the duct. .................................. 74 Figure 6.21: Cp vs TSR comparison for the turbine in its initial position at U∞=1.5 m/s. ..... 75 Figure 6.22: Cp vs TSR comparison for the turbine in its initial position at U∞=2 m/s. ........ 75 Figure 6.23: Torque vs. θ for the turbine in its initial position at TSR=3, U∞=1.5 m/s. ......... 76 Figure 6.24: Torque vs. θ for the turbine in its initial position at TSR=3 and U∞=2 m/s. ...... 76 Figure 6.25: Cp vs TSR for the turbine in its secondary position at U∞=1.5 m/s. .................. 77 Figure 6.26: Cp vs TSR for the turbine in its secondary position at U∞=2 m/s. ..................... 77 Figure 6.27: Torque variations for turbine in its secondary position at TSR=3,U∞=1.5 m/s .. 78 Figure 6.28: Torque variations for turbine in its secondary position at TSR=3,U∞=2 m/s ..... 78 Figure 6.29: Variations of AOA and torque for a blade of turbine at U∞=1.5 m/s. ............... 80 Figure 6.30: Torque vs. theta (θ) for the unducted turbine at U∞=1.5 m/s.............................. 81 Figure 6.31: Torque vs. theta (θ) for the ducted turbine at U∞=1.5 m/s.................................. 81 Figure 6.32: Cp vs. TSR for the unducted turbine at different tank size at U∞=1.5 m/s......... 83 Figure 6.33: Cp vs. TSR for ducted turbine at different tank size at U∞=1.5 m/s................... 83 Figure 6.34: Symmetric duct and its control points. ............................................................... 84 ix Figure 6.35: The results of symmetric duct optimization. ...................................................... 86 Figure 6.36: The optimized symmetric duct for WD=1.4D. .................................................... 87 Figure 6.37: The optimized symmetric duct for WD=1.6D. .................................................... 87 Figure 6.38: The optimized symmetric duct for WD=1.8D. .................................................... 88 Figure 6.39: The optimized symmetric duct for WD=2.0D. .................................................... 88 Figure 6.40: Asymmetric duct and the control points. ............................................................ 89 Figure 6.41: The results of asymmetric duct optimization...................................................... 91 Figure 6.42: The optimized asymmetric duct at WD=1.4D. .................................................... 91 Figure 6.43: The optimized asymmetric duct at WD=1.6D. .................................................... 92 Figure 6.44: The optimized asymmetric duct at WD=1.8D. .................................................... 92 x LIST OF SYMBOLS, ABBREVIATIONS AND NOMENCLATURE AOA : Angle of attack AR : Aspect ratio AR = CFD : Computational fluid dynamics HAT : Horizontal axis turbine IOT : Institute of ocean technology RANS : Reynolds average Navier-Stokes equations VAT : Vertical axis turbine UBC : University of British Columbia WD : Width of the duct LD : Length of the duct TSR : Tip speed ratio TSR = ρ : Density of the fluid [kg/m3] Ω : Angular velocity of the blade [rad/sec] U∞ : Free stream velocity[m/sec] R : Radius of the turbine [m] D : Diameter of the turbine [m] C : Chord length of the foil [m] N : Number of blades CP : Power coefficient of the turbine C P = T : Average torque of a turbine[N.m] Pmax : A : Frontal area of the turbine[m2] , A = 2 R ⋅ 1 Ft : Tangential force applied on a blade[N] Fn : Normal force applied on a blade[N] MC/4 : Moment at quarter-chord on a blade[N.m] α : Angle of attack of a blade(rad) bladespan C RΩ U∞ TΩ Pmax Maximum attainable power from the free stream flow in the frontal area of the turbine Pmax = 0.5 ρU ∞3 A xi α& : Rate of angle of attack change [1/sec] k : Non-dimensional pitch rate parameter, k = α&C 2U ∞ Cl : Lift coefficient Cl = Cd : Drag coefficient C d = L : Lift force [N] D : Drag force [N] W : Relative velocity of a blade [m/sec] Wt : Tangential velocity of a blade [m/sec] Wn : Normal velocity of a blade[m/sec] ∆t : Time step [sec], ∆t = ∆θ L 1 2 ρW 2 ⋅ C ⋅ 1 D 1 2 ρW 2 ⋅ C ⋅ 1 ∆θ (rad ) Ω The angle which blade rotates at each time step[rad] θ : Phase angle of the blade (rad) φ 2 : Potential velocity function [m /sec] w : Complex function [m2/sec] G : Boundary function [m2/sec] z : Complex variable [m], z=x+iy ΓC 2 : Circulation [m /sec] ΓQ.S . 2 : Quasi steady circulation [m /sec] Γp 2 : Pitching circulation [m /sec] γ : Strength of a discrete vortex [m2/sec] rC : Vortex core [m] σ µ : Source strength per length[1/m] Λ : Blockage ratio, Λ = At .c.a 2 : Cross section area of the tank[m ] Ad . p.a 2 : Net projected area of the device[m ] : Doublet strength per length [1/m] xii Ad , p ,a At .c.a ACKNOWLEDGMENTS First and foremost, I would like to thank Professor Sander Calisal for providing me this great opportunity to work under his supervision during my PhD studies. Second, I would like to thank Jon Mikkelsen for all his help during these years. Moreover, I would like to thank all of my colleagues at the Naval Architecture and Maritime Laboratory who have provided the excellent advice and support that has led to the completion of this work. Special thanks must go to Voytek Klaptocz, Bill Rawlings, James McRoberts and Kevin Gould for sharing their expertise and ideas with me. I would also like to thank the Institute of Ocean Technology in St John's, Newfoundland for hosting us during two weeks of testing. Furthermore, thank you to Yasser Nabavi and Ali Vakil who have been very good friends for me during the past few years. I would also like to thank Blue Energy Canada, especially Jon Ellison, for financial support of the project and generous help during the experiments. Last but not least, a very special thank you to my family for their endless support and encouragement in whatever I choose to do. xiii Chapter One Introduction 1.1: Research objectives The first objective of this research is to develop a tool that can be used for designing unducted vertical axis turbines. This method should be fast and able to quantify the effects of different design parameters, such as turbine geometric parameters and free stream velocity, on the performance of turbines. The second objective of this research is to develop a method to predict the performance of ducted vertical axis turbines at different operating conditions. This method should have the ability to determine the ducting effects on the output power of a vertical axis turbines without need to run a RANS (Reynolds Average Navier-Stokes) simulation that is very computationally expensive (each 2-D RANS simulation of a ducted turbine takes 2-3 weeks as reported by Nabavi [1]), or need to conduct physical experiments that are time consuming and expensive. The third objective of this research is to develop an optimization procedure to find the duct shape that maximizes the output power (i.e. power coefficient) of a vertical axis turbine. This procedure should be fast and have the ability to include different kinds of constraints such as the length and width of the duct. 1 1.2: Research methodology A new discrete vortex method referred to as “DVM” is developed for unducted vertical axis turbines. This method is two-dimensional and the flow is assumed to be potential. In this method, a new approach is developed to approximate the unsteady loads on the blades by mapping a flat plate to a moving cylinder and using the equations for loads on a moving cylinder as given by Milne-Thompson [2]. The unsteady wake of turbine is also modeled by discrete vortices shed from the blades. The interactions between the blades and discrete vortices are taken into account by representing each blade with a vortex filament. A new approach is also developed for the calculation of nascent vortex position relative to the blade; in addition, a new approach is presented to take into account the effects of connection point between a blade and its supporting arm. A numerical model referred to as “BEM-DVM” is developed for the performance calculation of ducted vertical axis turbines based on DVM model and a boundary element method referred to as “BEM”. This boundary element method is used to calculate the twodimensional potential flow field around the duct. Contrary to most boundary element methods, BEM model is fast and shows good mass flow continuity in the duct. In this model, the duct is represented by a series of panels which have constant distributions of sources and doublets. To quantify the towing tank wall effects on experimental results, a new approach is also developed based on BEM-DVM model. Understanding of these effects is important for the design of a full scale turbine operating in the ocean. The numerical models for unducted and ducted vertical axis turbines are validated by comparing the numerical results obtained from the FORTRAN codes developed for these models with the results obtained from the experiments and RANS simulations. The experimental results are from the tests conducted in the towing tank at UBC [3] and the towing tank at IOT [4]; the author formed part of the team who conducted these experiments. The RANS simulation results are from studies conducted by Nabavi [1]. To find the duct shape that maximizes the power coefficient of a vertical axis turbine, an optimization procedure is developed by coupling BEM-DVM model with the MATLAB optimization toolbox. The coordinates of the control points representing the duct shape are used as design variables. At each stage of the optimization process, the duct is created by fitting splines to these control points using the MATLAB spline toolbox. 2 1.3: Thesis outline A background on renewable energies is given in the first part of Chapter 2. Then, two devices designed for harnessing energy of hydro currents, horizontal axis turbines and vertical axis turbines, are introduced. Next, a discussion on vertical axis hydro current turbines studied in this research is provided. Finally, a literature review is provided to present the state of the art of vertical axis turbines. In the first part of Chapter 3, the flow dynamics in a vertical axis turbine are briefly discussed. Then a new method for calculating the loads and circulation on an unsteady blade of vertical axis turbine is developed by mapping a flat plate to a moving cylinder. Next, the approach adopted for modeling the vortex shedding from an unsteady blade is explained. The details of a new model developed for analysis of unducted vertical axis turbines, DVM, are then discussed. Since this model requires the characteristics of the foil used as blade profile, the approach for the calculation of these characteristics is explained in the last part of this chapter. Chapter 4 is focused on development of a new numerical model for the analysis of ducted vertical axis turbines. In the first part of this chapter, the boundary element method used for the calculation of potential flow field around the duct, BEM, is discussed extensively. The model developed for ducted vertical axis turbines, BEM-DVM, is subsequently explained. In the last part of this chapter, the procedure developed to calculate the effects of towing tank walls on experimental results is disscussed. In Chapter 5, an introduction to the numerical optimization methods is first provided, then a procedure is developed for finding the duct shape that maximizes the output power of a vertical axis turbine. In chapter 6, the new numerical methods for unducted and ducted vertical axis turbines are validated by the comparison of numerical and experimental results, followed by a discussion on the sources of discrepancies between the results. The results of a parametric study on the turbine performance are then provided. In the last part of this chapter, the results of duct optimization for a vertical axis turbine are illustrated. Lastly, in Chapter 7, conclusions of the research and the recommendations for the future work are presented. 3 1.4: Scope of research The numerical methods in this research are fast and reliable tools for the analysis and design of unducted and ducted vertical axis turbines for most operating conditions; however, the turbulence of the incoming flow is not taken into account in these methods. These methods do not work well at tip speed ratio (TSR) values lower than 2.5, as dynamic stall phenomenon which is dominant at these TSR values is not taken into consideration. Since these methods are two-dimensional, they do not provide good performance prediction for turbine configurations in which the blades have a low aspect ratio (AR), especially for AR<10, where the three-dimensional tip losses are significant. The model developed for ducted vertical axis turbines does not work well for ducts that have large separation zones, as the boundary layer of the ducts is not considered during the calculation of the flow field around the duct. Also, the interactions from the shaft and supporting arms of the turbine are not included in these models, though their experimental parasitic drag values are used in the post processing of the data to relate the numerical and experimental results. Since the optimization approach developed in this research uses a gradient-based optimizer, one should try different initial conditions to find the optimum duct for a vertical axis turbine, as the gradient-based optimization methods find only the local optimum and not the global optima. 4 Chapter Two Background 2.1: Renewable energy technologies Renewable energy technologies use natural resources such as sunlight, wind, rain, tides and geothermal heat, all of which are naturally replenished. Climate change concerns coupled with high oil prices are driving research and development on renewable energies [5]. The main types of renewable energies are as follows [6]: • Biomass energy • Bio-fuels energy • Geothermal energy • Wind energy • Solar energy • Hydro energy o Hydroelectric energy o Wave energy o Tidal energy o River current energy 5 Harnessing the energy of tidal and river currents, referred to as “hydro currents”, is among the new means of power production from renewable energies, compared to harnessing the wind energy. The interest for extraction of energy from tidal and river currents has been increased in the past few years, as they are highly reliable and predictable. The extraction of these energies using low-head turbines is expected to be environmentally benign. In this research, the turbines developed for capturing the energy of hydro currents are referred to as “hydro current turbines”. 2.2: Vertical axis turbines The modern types of these turbines, which rotate around a vertical axis, employ foil-section blades to generate lift. They convert lift into positive torque when their blades are traveling sufficiently fast relative to the free stream. These types of turbines originally developed by G. J. Darrieus for generation of electricity from the wind energy in 1931[7]. The first known wind tunnel testing of vertical axis turbine were carried out by P. South and R.S. Rangi at the National Research of Canada in 1971[8]. Over the past three decades, vertical axis turbines have undergone considerable research, development and large scale field trials. The blades of vertical axis turbines have a wide range of shapes. Figure 2.1 shows a turbine with helical blades, while Figure 2.2 shows a straight-bladed turbine. In this research, the vertical axis turbines with straight blades are focused. Figure 2.2: A VAT with straight blades. Figure 2.1: A VAT with helical blades [9]. 6 2.3: Flow dynamics in a vertical axis turbine Each blade element of VAT moves in a circular path, as shown in Figure 2.3. The blade experiences changes in its angle of attack as it rotates, leading to highly unsteady flow dynamics over the blade. Low torque or even negative torque is produced by the blade at phase angles near 0 and 180 degrees as it encounters small angles of attack. On the other hand, high torque is produced around phase angles of 90 and 270 where the blade sees high angle of attack. Figure 2.3: Schematic of circular movement for a blade element of VAT. 2.4: Advantages of the vertical axis turbines The main advantages of vertical axis turbines are discussed briefly below: • Its operation is independent of free stream direction The direction of incoming flow continuously changes as there is a lot of turbulence in the free stream. The omni-directional operation of vertical axis turbines eliminates the need for a system to adjust the orientation of the turbine blades to the flow direction, thus reducing the turbine cost. • Mounting of gearbox and generator is simple The vertical axis turbines are well suited to drive vertical transmission shafts; therefore it allows the gearbox and the generator to be mounted above water. This unique feature makes the system maintenance easier and therefore decreases the energy conversion cost using vertical axis turbines. In addition, this feature allows stacking of turbines under bridges or other suitable structures. 7 2.5: Disadvantages of vertical axis turbines There are some challenges that need to be addressed to make vertical axis turbines a reliable choice for harnessing the energy of hydro currents. These issues are explained below: • Not self-starting Vertical axis turbines may fail to self-start even under no-load conditions, depending on the blades position relative to the free stream direction, the foil characteristics, free stream velocity, etc. • High torque fluctuations As mentioned in Section 2.3, vertical axis turbines generate fluctuating torque. The magnitude of fluctuations is especially high at low tip speed ratios (TSR) [3]. This means that turbine components such as bearings, arms, blades, gearbox and generator need to be designed stronger and possibly bigger to withstand the peak fluctuating loads. This leads to a more expensive and possibly heavier turbine with higher parasitic drag. 2.6: Ducted vertical axis turbines Unducted turbines extract the flow energy by reducing the flow velocity with little or no pressure reduction as the fluid passes through the turbine rotor. Therefore the streamlines must expand to maintain continuity as shown in Figure 2.4. However, they cannot expand indefinitely; hence there is a theoretical limit to the percentage of kinetic energy that can be extracted from the flow. This limit has been shown by Betz [10] to be 59.3% using single actuator disk theory. Newman [11] showed that the corresponding limit is 64% for a double actuator disk such as a vertical axis turbine. Figure 2.4: Streamtubes in an unducted turbine and a ducted turbine. 8 However, if the turbine is confined in a duct as shown in Figure 2.4, the flow boundaries are defined and the streamline expansion is limited by the duct geometry. Contrary to unducted turbines, the flow energy in ducted turbines is extracted primarily by a pressure drop across the rotor [12]. The pressure drop in a ducted turbine depends on the duct shape and the flow field around it. Riegler [13] showed that the theoretical maximum power coefficient for a ducted turbine based on the turbine frontal area is 1.96, i.e. 3.3 times of the Betz limit. This research is mainly focused on development a method for analysis and designing efficient ducted vertical axis hydro current turbines. Schematic of a ducted turbine is shown in Figure 2.5. The advantages associated with ducting of VAT are outlined below: • Higher power can be extracted from a ducted turbine by concentrating the energy of a large area into a smaller area. In this way a smaller, lower cost turbine can be used for a given output power. • Ducting the turbine decreases the torque fluctuations, as shown in UBC experiments [3]. • Low speed sites can be considered for energy conversion, as the flow accelerates in duct. • Ducts can be used as a float for the structural support of the turbine. • In areas where there is a danger of divers or floating debris being drawn into the turbine, a grid supported by the ducting could be placed on the upstream opening, thus reducing risk to marine life as well as damage to the turbine from debris. • Higher velocity in a duct decreases the amount of biofouling on the turbine. Figure 2.5: Schematic of a ducted turbine. 9 2.7: Literature review Development of an efficient and reliable vertical axis turbine requires good understanding and insight into the performance of vertical axis turbines in different operating conditions including varying flow speeds, angular velocity, and turbine sizing. In order to obtain this knowledge, a number of numerical and experimental studies were conducted over the past few decades. These studies are reviewed briefly as follows: 2.7.1: Numerical study of unducted vertical axis turbines The numerical models developed for vertical axis turbines range from simple and fast models to complicated and time-consuming ones. They can be classified into two categories: potential flow models and RANS models. 2.7.1.1: Potential flow models Potential flow models have been used for decades as the primary design tool, as they are simple and computationally fast. The potential flow models developed for vertical axis turbines can be categorized as momentum models and vortex models. The following sections cover a short review of these models; a comprehensive review of these methods is provided by Paraschivoiu [7] and Nabavi [1]. • Momentum models Momentum models are based on Glauert’s Actuator Disc Theory and Blade Element Theory [7]. The concept used in these models is that the total change in the axial momentum across the actuator disc is equal to the aerodynamic forces exerted on the blades in the axial direction. This force is also equal to the pressure difference across the disk. Bernoulli’s equation is then used in each stream tube to obtain the velocity in the wake. The main drawback of these models is that they become invalid for high tip speed ratios and also for high rotor solidities, as the momentum equation become invalid in these particular cases [7]. The main momentum models developed are the Single Streamtube model, the Multiple Streamtube model and the Double-Multiple Streamtube model. In 1974, Templin developed the Single Streamtube model; the first and the simplest model for the performance calculation of a vertical axis turbine [14]. This model assumes that the rotor is enclosed in a single streamtube as shown in Figure 2.6. The flow velocity within the streamtube is assumed to be uniform. While this model is elegant in its simplicity and 10 predicts overall performance rather well for lightly loaded blades, it is incapable of adequately predicting the loads on the blades as it requires a more precise knowledge of the flow velocity variations across the rotor. These variations become increasingly large as blade solidities and blade tip speed increases. To account for this non-uniform flow, Wilson and Lissaman used a sinusoidal variation in the inflow velocity across the width of turbine [15]. Figure 2.6: Single streamtube model [1]. Figure 2.7: Multiple streamtube model [1]. A further improved model is the Multiple Streamtube Model which was developed in 1975 [16]. In this model the rotor is enclosed in a series of streamtubes as shown in Figure 2.7. The streamtubes are aerodynamically independent in this model. The momentum balance is carried out separately for each streamtube resulting in different velocities in each streamtube. Despite the fact that this model produces better results than the single streamtube model, the results are still only valid for lightly loaded blades. In 1981, Paraschivoiu [17] introduced the Double-Multiple Streamtube model for the vertical axis turbine. In this model the differences between the upwind and downwind passes of each blade are taken into account by dividing each streamtube into an upwind half and a downwind half, as shown in Figure 2.8. The momentum balance is then carried out separately for each half of the streamtube. This model gives a better correlation between the calculated and the experimental results compared to the results of Multiple Streamtube model, especially for the local aerodynamic blade forces. On the other hand, this model 11 appears to have convergence problems, especially on the downstream side and at higher tip speed ratios [18]. Figure 2.8: Schematic of double-multiple streamtube model [1]. • Vortex models The vortex models developed for vertical axis turbines calculate the velocity field about the turbine by including the effects of vorticity in the wake of turbine. One of the first vortex models was developed by Larson in 1975 [19]. He used this model to analyze a cyclogiro windmill, in which the blades periodically flipped from a set positive pitch angle to a set negative angle and back at each revolution. For the performance calculation, he used a simplified wake with only two vortices that shed into the wake at each revolution at the points at which the blades flipped, and calculated an average velocity by which the vortices convect downstream. In 1976, Holme presented a two-dimensional vortex model for a fast running vertical-axis wind turbine, having a great number of straight and very narrow blades [20]. His analysis is valid only for very lightly loaded turbines in an inviscid fluid flow, but approximations for heavier blade loading and corrections for viscous effects are also provided. A two-dimensional analysis for “Giromill” performance is presented by Wilson in 1980 [21]. Giromill is a vertical-axis wind turbine with straight blades that are articulated to produce maximum energy extraction from the wind. It is found that the power coefficient and windwise force coefficient for the Giromill have the same limit as obtained for the horizontal axis wind turbines. 12 Wilson and Walker proposed a model for VAT called as Fixed Wake model [22]. To determine the differences in induced flow between the upwind and downwind blades, a vortex sheet wake is used. Forces are determined from the Kutta-Joukowski theorem and the principle of circulation conservation. The computational cost of this model is of the same order as the streamtube models. Free Wake models are the most complex and accurate vortex models that have been developed for vertical axis turbines. In all of these models, the wake is modeled by discrete, force-free vortices convecting downstream with local flow velocity. However, the blades are modeled in different ways. Fanucci and Walters presented the first Free Wake model in 1976 [23] for a straight-bladed rotor. They modeled the blades by a distribution of vortices along the blade camber line. In the method proposed by Strickland [24], each blade is replaced by a vortex filament as shown in Figure 2.9. In 2008, Li from the University of British Columbia developed a method to predict the output power from a tidal current turbine farm, by modeling each turbine similar to that presented by Strickland [25]. Figure 2.9: Modeling a blade by a vortex filament [1]. 2.7.1.2: RANS modeling The potential flow models are computationally fast, but they suffer from different simplifying assumptions. With the exponential increase in the computational speed, the RANS simulations are becoming more popular, as they can provide valuable insight into the flow field and can facilitate the optimization process. However, still most of the RANS simulations for vertical axis turbines are two-dimensional, as they are very expensive and time consuming. 13 Rajagopalan in 1986 did a two-dimensional study of viscous flow field around a vertical axis wind turbine [26]. In this study, he solved the 2-D steady, laminar, Navier-Stokes equations. In 1995, Ishimatsu and his colleagues reported a RANS simulation study for the flow field around a two-bladed Darrieus wind turbine using NACA 0018 as the blade profile [27]. In 1997, Allet from École Polytechnique in Montreal used a RANS solver to investigate the flow field around a Darrieus turbine [28]. In 2007, Guerri et al. [29] and Jiang et al. [30] from the Zhejiang University separately investigated the flow field around a vertical axis turbine by solving the RANS equations. Also in 2007, Nabavi [1] from the University of British Columbia studied the performance of a vertical axis turbine at different operating conditions by solving the RANS equations using the commercial package FLUENT. 2.7.2: Numerical study of ducted vertical axis turbines Not many numerical studies for ducted vertical axis turbine are found in the literature. In most of the researches on ducting, the effects of wind tunnel walls on the model performance are examined. In 2000, Wang and Coton [31] studied the wall effects of wind tunnel on a horizontal axis wind turbine, using a Boundary Element method. In 2003, Lawn from the University of London analyzed the performance of a shrouded hydro current turbine (horizontal or vertical), using a one-dimensional analysis, treating the ducts upstream and downstream of the turbine as contractions or expansions having specified diffusion efficiencies [32]. In 2007, the effects of ducting on the performance of a vertical axis turbine were studied extensively by Nabavi at UBC [1]. He examined a wide range of duct geometries using the commercial package FLUENT. In 2008, Werle et al. developed an analytical model for ducted/shrouded propeller-based propulsion and power systems using the momentum theory [33]. 2.7.3: Experimental study of unducted and ducted vertical axis turbines Most of the researches on ducting have been done on diffuser-augmented horizontal axis wind turbines. In 1983, Gilbert and Foreman tested many diffuser configurations in a wind tunnel [34]. They could increase the output power of a wind turbine model up to 4.25 times when a diffuser-type duct was added. They were able to achieve this result with a short, wide angle diffuser; however one should note that the turbine they used had untwisted blades with low aspect ratio which its performance in free stream flow was well below the Betz limit. 14 The first experimental study on vertical axis turbines for hydro applications was conducted by Barry Davis and his co-workers from NRC-Ottawa in 1978. Their extensive research during the 1980’s, completed by Nova Energy in collaboration with NRC, resulted in a number of pilot projects and the publication of multiple reports as given in Table 2.1. They also conducted several independent assessments validating the technology, as well as experimental studies on ducting effects on the performance of vertical axis turbines. However, due to the low cost of fossil fuels and the lack of funding for further development of tidal energy at the time, neither Nova Energy nor its successor, Blue Energy, was able to establish any further major projects through the 1990s. Table 2.1: Available Davis et al. reports Report Title Synopsis Flume tank tests of vertical and NEL-002: Water Turbine Model Trials [35] NEL-021: horizontal axis water turbines. Ultra Low Head Hydroelectric Power Vertical axis water turbine flume Generation Using Ducted Vertical Axis Water tank tests with caissons, walls, and Turbines [36] NEL-022: vane duct configurations. Ultra Low Head Hydroelectric Power Generation Using Ducted Vertical Axis Water Turbines [37] NEL-038: Research and Development of a 50kW to 100kW Vertical Axis Hydro Turbine for a Restricted Flow Installation [38] NEL-070: The Ducted Vertical Axis Hydro Turbine for Large Scale Tidal Energy Applications [39] Continuation of NEL-021 with a more robust model. Installation of 70 kW turbine within a dam in Nova Scotia. Investigates application of vertical axis turbine in a 474 turbine tidal fence. NEL-081: Commissioning and Testing of a 100kW Examines repaired and enhanced Vertical Axis Hydraulic Turbine [40] version of model in NEL-038. 15 Gorlov patented a vertical axis turbine using helical blades to distribute the torque loading in 1994 (U.S. Patent 5451127), and continues development work in Korea [41]. Given the commercial nature of this venture, efficiency data and torque curve data are closely guarded. In 2000, the researchers at the Nihon University in Japan tested the effects of number of blades, turbine solidity and the starting torque (the torque required by the turbine at start) on the performance of a Darrieus turbine in a water tank [42]. In 2004, the Department of Trade and Industry in UK sponsored three reports on vertical-axis tidal turbines. Only one report attempted experimental trials for numerical model validation and provided no quantitative data due to a number of factors, including excess friction in the gearbox and a less than ideal experimental flume facility [43]. To develop a new type of two-way diffuser suitable for a fluid flow energy conversion system, an experimental study was carried out by the researchers at the Saga University in 2004 [44]. In this study, a new type of two-way diffuser was developed for tidal currents. In 2005, a research was undertaken in Italy by Ponte di Archimede S.P.A. Company and the University of Napoli on a turbine with a patented passive angle of attack adjustment mechanism [45, 46]. No publicly available torque curve data has been published to date. From 1999 through 2005, researchers at the University of South Australia conducted several experiments to investigate the ducting effects on the output power of hydro current turbines [12]. Although the overall performance of the unducted turbine was low, they could increase the performance of a turbine model by a factor of 3 when a diffuser-type duct was added. In 2007, a group of researchers from the University of Buenos Aires tested the ducting effects on the performance of vertical axis turbines. In these experiments, they tested a series of scale models with different combination of deflectors, nozzles and diffusers [47]. Researchers at UBC conducted a wide range of experimental studies from 2006 to 2008 [3]. In these studies, the effects of parasitic drag, tip losses, angle of attack, cambered blades, shaft fairing and ducting on the performance of a vertical axis turbine were investigated. Their experimental results showed significant increase in the output power of the turbine due to ducting. These experiments also showed that ducting can significantly decrease the torque fluctuations of the vertical axis turbine. 16 Chapter Three Unducted Vertical Axis Turbine 3.1: Introduction The unsteady flow dynamics in a vertical axis turbine (VAT) is one of the main challenges for researchers and developers looking to increase the efficiency of this turbine. The changes in the angle of attack of blades during their rotation cause a difference between the static and dynamic loading on the blades. Noll and Ham [48] stated that the maximum loads on the blades of a vertical axis turbine may be as much as three times greater under the dynamic conditions than under the static conditions, while the peak pitching moment may be up to five times greater. Other authors (Klimas [49]; Paraschivoiu and Allet [50]; Jiang et al. [51]) have emphasized the importance of accounting the unsteady loading from the perspective of structural design, fatigue life and drive-chain and generator sizing. The effect of unsteadiness on overall turbine performance may also be significant, especially at low tip speed ratios. In general, the increased peak loads and delayed stall that are associated with the dynamic conditions lead to an improvement in the performance of turbine. Unsteady loading on a foil can arise due to several effects. In absence of stall, these effects can be loosely categorized as “added mass effects” and “circulatory effects”. But when stall happens, the unsteady loading is mainly associated to what is referred to as “dynamic stall”. 17 • Added mass effects Added mass is the inertia added to a system because an accelerating or decelerating body must move some volume of the surrounding fluid as it moves through it, since the object and fluid cannot occupy the same physical space simultaneously. For simplicity this can be modeled as some volume of fluid moving with the object, though in reality all the fluid will be accelerated to various degrees [2]. • Circulatory effects The pitching motion of a foil gives rise to changes in its circulation as a result of required adjustments to satisfy the Kutta condition. On the other hand because of the total vorticity conservation constraint [52], the change of foil circulation due to its movement leads to the formation of an unsteady wake in the form of a continuous vortex sheet sheds from the trailing edge of the foil. This unsteady wake itself also changes the circulation of foil. Various researchers have attempted to approximate the added mass effects and circulatory effects. Fung [53] used a classical thin airfoil theory to calculate lift and drag of a pitching and plunging foil. Using a quasi-steady analysis and ignoring the unsteady wake, he showed that the lift is proportional to the angle of attack at 3/4 chord measured from the leading edge. Strickland et al. [54] used the analysis by Milne-Thomson [2] for forces on a moving and rotating cylinder to estimate the circulation and loads on a blade of a vertical axis turbine. His analysis yields expressions for the additional circulation around the foil due to its rotation, as well as the added mass forces that are independent of the circulation. The tangential force is found to be proportional to the angle of attack at the mid-chord position. The normal force consists of a term that is proportional to the angle of attack at the threequarter chord position and an added mass term proportional to the normal acceleration of foil. Strickland et al. stated that an order of magnitude analysis justifies the omission of the added mass term and pitching moment. These results give good insight into the effects of unsteady dynamics on the circulation and forces of blades. However, calculating the normal force based on the angle of attack at three-quarter chord might not be a good approximation of the normal force on a blade of a vertical axis turbine, as these results are obtained from the analysis of a moving cylinder which does not much resemble a foil. 18 Wilson et al. [55] presented another method for the calculation of the loads on a blade of vertical axis turbine by mapping a flat plate to a cylinder, and then using the results of analysis for a moving and rotating cylinder presented by Milne-Thomson. Despite the fact that, Wilson and Strickland gave different expressions for the forces and moments, both of them concluded that the pitching moment and added mass term does not have much effect on the net power extraction; in other words, the quasi-steady assumption is sufficiently accurate for the calculation of the net power from a turbine. • Dynamic stall Experiments with different foils show that under steady conditions, the lift increases in an approximately linear fashion with the angle of attack up to a critical angle called the static stall angle of attack. At this angle of attack, the flow begins to separate from the upper surface of the foil (suction side) which leads to the sudden loss of lift and an increase in drag. Under dynamic conditions, the angle at which stall occurs depends on the rate of angle of attack change. This phenomenon at which under the dynamic conditions stall occurs at a higher angle of attack than under the static conditions, with an associated increase in the maximum attained lift is called “dynamic stall” [56]. Considerable experimental attention has been devoted to the topic of dynamic stall by measuring forces and pressure distributions on the pitching foils in wind tunnels. These experiments show that the magnitude of the stall delay and the peak loads attained is a function of the non-dimensional pitch rate parameter k = α&C 2U ∞ [56]. Different dynamic stall models are proposed because of disagreements over the magnitude of the stall delay, the maximum lift, and the non-dimensional pitch rate parameter, k . Among different dynamic stall models which originally developed for dynamic stall study of helicopter blades, Gormont model [57], MIT model [58] and more recently BeddosLeishman model [59] are the most popular models. The dynamic stall phenomenon in a vertical axis turbine is more complex than the dynamic stall in a helicopter rotor because of higher thickness ratio of the blades and higher angles of attack seen by the blades in a vertical axis turbine. Dynamic stall is not considered in this research because of lack of a reliable model for different pitch rates, Reynolds numbers and maximum angle of attack. 19 3.2: Loads on a blade of vertical axis turbine Here, a new method is developed for the calculation of the loads on a blade of vertical axis turbine. The main goal for developing this method is to address the shortcomings of the previous methods. While this method is easy to implement, it contains the main features of unsteady flow dynamics on a blade of vertical axis turbine. The conjugate force and moment on a moving cylinder is obtained by Milne-Thompson [2]. By specializing those results to the case of a flat plate rotating with angular velocity Ω on a circular path with radius R (Fig. 3.2), the conjugate force components, Fx and Fy , as well as the moment at mid-chord, M C / 2 , are as follows [55]: 2 ∂ iρ dw Fx − iFy = ∫ dz − ρΩ ∫ z dw + iρ ∫ w dz + 2πik C ρQ 2 B dz ∂t B B (3.1) ρ dw 2 ∂ M C/2 = Real Part − ∫ z dz + ρQ ∫ zdw − ρ ∫ zwdz ∂t B 2 B dz B (3.2) Here z is the complex variable and z is its conjugate; where ( x, y ) is the local coordinate system as shown in Figure 3.1. B is the bounding curve of flat plate and 2πk C represents the circulation around this curve. w , w are complex potential and its conjugate respectively. The r derivative dw dz yields conjugate velocity Q = Q x − iQ y , where its components at each phase angle θ are: Q x = − RΩ − U ∞ cos θ (3.3) Q y = −U ∞ sin θ x y Ω θ U∞ R Y X Figure 3.1: The local and global coordinate system. 20 The integrals above can be evaluated by transforming the flat plate from the z plane to the ξ plane. In the ξ plane the flat plate is transformed into a circle, using the following transformation: z = ξ + a2 ξ (3.4) Where a is the radius of circle as shown in Figure 3.2. ξ plane ϕ z plane y C a x ψ T .E. Figure 3.2: Mapping of a flat plate to a circle. From the transformation of equation 3.4, the relation between the chord length of flat plate and the radius of circle is: C = 4a From the articles 9.63 and 9.66 of reference 2, the complex potential, w , for a moving cylinder with circulation 2πκ C is: w = G1 (ξ ) + ik C ln ξ (3.5) Where G1 (ξ ) is equal to [2]: G1 (ξ ) = 2iI (Q ) a2 ξ − iΩ a4 (3.6) ξ2 The term I (Q ) is the imaginary part of conjugate velocity Q . Equation 3.5 gives the complex velocity potential for a moving circle; it depends on the circulation around the cylinder, which is always constant. However, based on the Circle Theorem [2], the complex potential of an unsteady foil (i.e. flat plate) depends on its unsteady wake as well as on its circulation. As it is mentioned earlier, the unsteady wake forms as a result of total vorticity constraint [52]. Here, the dependency of complex potential on unsteady wake is ignored to simplify derivation of circulation and loads on an unsteady foil. In next sections, it is discussed how the unsteady wake is taken into consideration. 21 By ignoring the unsteady wake, the complex velocity potential for a flat plate becomes as: w = ikC ln ξ + 2iI (Q ) a2 ξ − iΩ a4 (3.7) ξ2 Using the complex velocity potential of equation 3.7, the results of integrals in equation 3.1 are as follows: 2 dw ∫A dz dz = 0 ∫ zdw = 4πI (Q )a 2 A ∫ w dz = − 4πrk C − 4πI (Q )a 2 A Therefore the conjugate force on the foil is: Fx − iFy = π 4 1 2 ρΩQ y C 2 + i ρC π ∂ − 2πk C + Q y C + iρ (2πk C )Q 2 ∂t (3.8) Similarly using the complex velocity potential of equation 3.7, the results of integrals in equation 3.2 are as follows: 2 dw 2 ∫A z dz dz = −i2πk C ∫ zdw = 4πI (Q )a 2 A ∫ zwdz = −2πk C a 2 + 2πΩa 4 A Therefore the moment about the mid-chord is as: ∂ M C / 2 = Real Part iπρk C2 + 4πρQI (Q )a 2 − ρ − 2πk C a 2 + 2πΩa 4 ∂t ( ) From the above equation, the moment for a blade with constant angular velocity Ω is: ( ) M C / 2 = ρQ x 4πI (Q )a 2 − ρ ∂ − 2πk C a 2 ∂t ( ) 22 By replacing the time derivative, MC/2 = − π 4 ρQ x Q y C 2 + ∂ ∂ , with Ω , the moment at mid-chord becomes: ∂t ∂θ 1 ∂ ρC 2 Ω (2πk C ) 16 ∂θ (3.9) The trailing edge is a singular point of w , this implies that the velocity at the trailing edge must be set to zero in the circle plain , then: dw =0 dξ ξ =a This gives the circulation around the cylinder as: ΓC = 2πk C = −πCQ y − π 4 C 2Ω (3.10) On the other hand, the Kutta-Joukowski theorem [2] states that the circulation about a stationary foil with relative velocity Vrel = (−Q x ,−Q y ) is: ΓQ.S . = 1 CC l Vrel 2 (3.11) Also the lift coefficient of the foil, C l , at small angles of attack α is: C l = 2πα = 2π sin α Then the circulation of a stationary airfoil ΓQ.S . is: ΓQ.S . = −πCQ y (3.12) Comparing equations 3.10 and 3.12, one can conclude that the circulation of a blade of vertical axis turbine is equal to the sum of a quasi-steady circulation and a pitching circulation. This means that: ΓC = ΓQ.S . + ΓP (3.13) Where the pitching circulation is: Γp = − π 4 C 2Ω (3.14) 23 From equations 3.8 and 3.10, the conjugate force on the blade is: Fx − iFy = π 4 ρC 2 ΩQ y + i 3π ρC 2 ΩU ∞ sin θ + iρΓC (Q x − iQ y ) 4 Replacing ΓC in the last term above from equation 3.10, gives the force in x direction as: Fx = ρΓQ.S .Q y (3.15) It is interesting to observe that the force in the x direction ( Fx ) for a pitching thin foil depends only on the quasi-steady circulation. The force in y direction is: Fy = − ρΓQ.S .Q x − π 4 ρC 2 ΩU ∞ (4 sin θ + TSR ) (3.16) From equation 3.9, the moment about the mid-chord is: MC/2 = − π 4 ρC 2 Q x Q y − π 16 ρC 3 ΩU ∞ sin θ (3.17) It is more appropriate to work with force and moment at the quarter-chord of foil as the center of pressure of a stationary thin foil is close to its quarter-chord position measured from leading edge (Fig. 3.3). Then, the moment at quarter-chord is: C M C / 4 = M C / 2 + Fy 4 This gives the moment at quarter-chord which is: MC/4 = − π 16 ρC 3 ΩU ∞ (5 sin θ + TSR ) (3.18) Fy Fy M C /2 M C /4 Fx L.E. Fx L.E. α α Vrel Vrel Figure 3.3: Loads on a foil are placed in its quarter-chord. 24 3.3: Development of DVM for vertical axis turbines As shown in section 3.2, the circulation and loads for an unsteady foil (without considering its unsteady wake) is composed of two components: a quasi-steady component and a pitching component. This decomposition of loads and circulation enables us to develop a new model for vertical axis turbine referred to as “DVM” which has the following assumptions: • The model is two-dimensional and the flow is potential. • The pitching components of loads and circulation are obtained using the result of mapping a rotating flat plate onto a rotating cylinder, while the quasi-steady components are obtained using the stationary characteristics of foil. • The unsteady wake of blades is modeled by discrete vortices shed from the blades. It is assumed that these discrete vortices influence only the quasi-steady components of loads and circulation of blades. • To calculate the interaction between the foils and the unsteady wake, each blade is represented by a vortex filament. • The arms and shaft are not considered in this modeling. Figure 3.4 shows a schematic of vertical axis turbine versus the DVM model of this turbine. In DVM, an unsteady blade of vertical axis turbine which is rotating with angular velocity Ω is represented by a vortex filament (stationary momentarily) with relative flow velocity W . Note that the velocity components U W , VW are the induced velocity components of unsteady wake on the blade. W θ α VW RΩ R θ UW U ∞ U∞ Ω Vertical Axis Turbine DVM Modeling of Vertical Axis Turbine Figure 3.4: A schematic of vertical axis turbine vs. DVM modeling of this turbine. 25 3.3.1: Wake representation in DVM As it is mentioned in section 3.1, the total vorticity conservation constraint [2] gives rise to an unsteady wake in the form of a continuous vortex sheet sheds from the trailing edge of unsteady foil. However, in DVM modeling, the vortex shedding from unsteady blade is modeled by a discrete process, which at each time step a discrete vortex sheds from the blade. Based on the Kelvin theorem [52], when the blade moves from its initial position with circulation ΓCn to a new position with circulation ΓCn +1 as shown in Figure 3.5, a vortex referred to as “nascent vortex”, sheds from the foil with strength γ equal to: γ = −(ΓCn +1 − ΓCn ) (3.19) From the Kutta-Joukowski theorem, the circulation of a quasi-steady blade with relative flow velocity W is: ΓC = 1 C l CW 2 ΓCn +1 γ ΓCn Figure 3.5: Circulation around any contour remains constant in a potential flow field. Figure 3.6 shows a schematic of the wake of turbine which is represented by discrete vortices. The new position of discrete vortices which move with the local flow velocity after a time step ∆t is calculated as: X vn +1 = X vn + (U ∞ + ∑ U vi )∆t ∑U , ∑V i v i v Yvn+1 = Yvn + (∑Vvi )∆t are the total velocity components induced by all vortices and the blades. Figure 3.6: The wake is modeled by discrete vortices in DVM. 26 (3.20) 3.3.2: Calculation of nascent vortex position relative to foil The vortex shedding is a continuous process and appears as a vortex sheet. Figure 3.7 shows a schematic of vortex sheet that sheds from the trailing edge of a blade when it moves from the phase angle θ1 to the phase angle θ 2 . But in DVM, the vortex shedding at the same period of time is modeled with a nascent vortex of strength γ as shown in Figure 3.8. However, the position of nascent vortex is not known and needs to be specified. Vortex sheet γ U∞ θ2 Figure 3.7: The vortex sheet from U∞ θ1 θ1 to θ 2 . θ2 θ1 Figure 3.8: The equivalent discrete vortex in DVM. To find an explicit expression for the position of nascent vortex, it is assumed that the vortex sheet remains along the chord when the blade moves a small angle ∆θ = θ 2 − θ1 as shown in Figure 3.9. Using this assumption implies that the movement of vortex sheet under the influence of wake is not taken into consideration. Note that the ratio of chord to radius and the angle ∆θ = θ 2 − θ1 is exaggerated in this figure in order to show a better figure. R·∆θ Vortex sheet θ2 S S´ γ θ Figure 3.9: The vortex sheet and its equivalent discrete vortex are along the chord. 27 Since the nascent vortex in DVM should have the same effect as the vortex sheet, then their induced velocities (based on Biot-Savart law [52]) at the mid-chord of blade is assumed to be θ2 equal. This implies that: dγ ′ γ ∫θ 2πS ′ = 2πS 1 Where γ is the strength of nascent vortex and dγ ′ is the strength of a vortex sheet element. S , S ′ are the distances of nascent vortex and the vortex sheet element from the mid-chord of blade respectively. From equations 3.10 and 3.19, the nascent vortex strength is: γ = πCU ∞ sin θ 2 − πCU ∞ sin θ1 Also, the strength of an element of vortex sheet with strength ω is: dγ ′ = ωdθ The vortex sheet strength calculated from equation 3.10 is as follows: ω= dΓC = πCU ∞ cos θ dθ Assuming the vortex sheet strength remains constant between angles θ1 and θ 2 , and is equal to the strength of vortex sheet at its middle; the vortex sheet strength becomes as: ω = ω m = πCU ∞ cos θ m Where θ m = θ1 + θ 2 2 Using the vortex sheet strength and the nascent vortex strength, the distance S becomes: S= sin θ 2 − sin θ1 θ2 dθ ′ θ1 S cos θ m ∫ After calculation of integral above, the position of nascent vortex relative to mid-chord is: χ S = C / 2 ln (1 + χ ) Where χ = (3.21) 2 R ⋅ ∆θ , ∆θ is in radian. C 28 3.3.3: The velocity induced by a discrete vortex The experimental studies conducted by Ciffone and Orloff [60] showed that a vortex rotates as a solid body within a distance rC from the vortex center referred to as “vortex core”. The tangential velocity increases linearly within the vortex core and reaches its maximum value at vortex core edge as shown in Figure 3.10. For points out of the vortex core, the tangential velocity changes inversely with distance r from the vortex center. However, the quantitative description of vortex core is not straightforward. Here the method presented by Strickland is used to calculate the radius of the vortex core [24]. Wvi WC Wi ∝ 1 r rC r Figure 3.10: Velocity variation of a vortex. For a point p located out of the vortex core, the tangential velocity of a counter clockwise vortex is [52]: Wvi = − γ 2π r for r ≥ rC (3.22) Where r is the distance of point p from the vortex center and γ is the vortex strength. Assuming the velocity WC is the velocity associated with the vortex sheet originating from the trailing edge of foil as shown in Figure 3.11, from the definition of circulation, we obtain: ∂Γ = 2WC ∂S → WC = 1 ∂Γ 2 ∂S Where S is distance along the vortex sheet as shown in Figure 3.11. But for nascent vortices, ∂Γ γ ≈ where γ is the nascent vortex strength and R∆θ is the ∂S R∆θ distance that blade moves at each time step. 29 Then the velocity WC becomes as: WC = 1 γ 2 R∆θ (3.23) Equating the equations (3.22) and (3.23) at r = rC gives the vortex core radius as: rC = R∆θ (3.24) π Thus the velocity induced by a discrete vortex at a point within its vortex core is: γ r Wvi = − 2 R∆θ rC for r ≤ rC (3.25) Vortex Sheet WC WC Shed Vortices dS WC WC Figure 3.11: Vortex shedding from a foil. From equations (3.22), (3.25), the velocity components induced by a vortex at a point p are: U vi = − π γ Y p − Yv 2 (R∆θ )2 For a point p in the vortex core π γ X p − Xv V =+ 2 (R∆θ )2 (3.26) i v U vi = − Y p − Yv γ 2 2 2π (X p − X v ) + (Y p − Yv ) For a point p out of the vortex core X p − Xv γ V =+ 2 2 2π (X p − X v ) + (Y p − Yv ) (3.27) i v Where (X p , Y p ) and ( X v , Yv ) are the coordinates of point p and vortex center respectively. 30 3.3.4: Angle of attack calculation for a blade of VAT The loads on a blade of vertical axis turbine are obtained in section 3.2. From equations 3.15 and 3.16, the tangential and normal forces on a blade are: 1 Ft = ρW 2 C C t 2 (3.28) π 1 Fn = ρW 2 C C n + ρC 2 ΩU ∞ (4 sin θ + TSR ) 4 2 (3.29) The tangential and normal directions are shown in Figure 3.12. The tangential force coefficient, C t , and the normal force coefficient, C n , in terms of lift and drag coefficients are as: C t = C l sin α − C d cos α , C n = C l cos α + C d sin α Note that in equations 3.28 and 3.29, the effects of wake which are modeled by discrete vortices, are included in the quasi-steady part of loads. W t α n RΩ ∑V i v ∑U i v R U∞ θ Figure 3.12: Velocity components for a blade. The local angle of attack at mid-chord of a blade is calculated as follows: Wn W α = sin −1 (3.30) The tangential and normal velocity components Wt , Wn , and the relative velocity W at each phase angle θ are: Wt = Rω + U ∞ + ∑ U vi cos θ + ∑ Vvi sin θ ( ( W = Wt 2 + Wn2 ) ) , Wn = U ∞ + ∑ U vi sin θ − ∑ Vvi cos θ ( 1 2 ) (3.31) (3.32) The velocities U iv and Vi v are the velocity components (in X and Y directions) induced by a discrete vortex or a vortex filament at the mid-chord of foil. 31 3.3.5: Effects of connection point between supporting arm and blade The aerodynamic loads obtained in section 3.2 are for the case that supporting arm is connected to the middle of blade. However, the connection point may be influential on flow dynamics over the blade and eventually on the loads applied to the blade. In order to take into account the connection point effect, here it is assumed that only the quasi-steady loads and circulation are influenced by the connection point. As shown in Figure 3.13, the mid-chord has a relative velocity Vκ relative to the connection point which is equal to: Vκ = (1. 2 − κ )CΩ Where κ is the position of connection point measured from the leading edge as a fraction of chord length. As this relative velocity is in normal direction to the blade, the tangential velocity does not change. But the normal velocity component at the mid-chord of blade changes to: Wn = U ∞ + ∑ U vi sin θ − ∑ Vvi cos θ + Vκ ( ) (3.33) Note that the induced velocity components U vi and Vvi are still calculated at the mid-chord. The loads on a blade connected to the supporting arm at a point other than the mid-chord are obtained by calculating the angle of attack from equation 3.30 using equation 3.33. In this case, the position of mid-chord is calculated from the following equations: X C / 2 = X κ + (1 / 2 − κ )C ⋅ cos θ X κ = − R sin θ where (3.34) YC / 2 = Yκ + (1 / 2 − κ )C ⋅ sin θ Yκ = R cos θ where C 2 Vκ κ Ω θ Figure 3.13: The blade connected to the arm at a point other than the mid-chord. 32 3.3.6: Calculation of turbine power coefficient As it is mentioned in section 3.2, it is preferable to work with force and moment at the quarter-chord of blade because the center of pressure for a stationary thin foil is usually close to its quarter-chord (see Fig 3.14). Therefore by placing the loads at the quarter-chord of blade, the torque produced by a blade of turbine at each phase angle θ is: Ti (θ ) = Ft ⋅ R + Fn ⋅ (κ − 1 4 )C + M C / 4 (3.35) Where κ is the connection point of arm and blade. The average torque of turbine with N blades is: N ∑T i T = (3.36) i =1 N where Ti is the average torque of i th blade over one revolution. The power coefficient for a turbine with N blades is defined as follows: CP = TΩ Pmax (3.37) where Pmax is the maximum attainable power from the free stream in the frontal area of turbine which is: Pmax = 0.5 ρU ∞3 A , A = 2 R ⋅ 1 C 4 M C/4 Ft κ Fn Ω Figure 3.14: The loads on a blade of VAT at the quarter-chord of blade. 33 3.3.7: Computational procedure for DVM A time-marching scheme is used for calculation of power coefficient in DVM. In this scheme the turbine rotates an angle ∆θ at each time-step ∆t , which is equal to ∆t = ∆θ (rad ) . The Ω convergence criterion is the difference of power coefficient at two successive rotations. The calculation continues until the convergence criterion is less than 0.001. Here, ∆θ in terms of degree is chosen to be 10 degrees. This value provides a good balance between the accuracy and the computational time. The flowchart of this scheme is shown below. Initial Data Calculate the flow velocity, the angle of attack and the circulation of blades Calculate the strength and position of nascent vortices Repeat last steps including nascent vortices Calculate induced velocities on all discrete vortices Calculate the tangential force, the normal force and the torque of blades Calculate the new positions of blades and the free vortices θ n +1 = θ n + ∆θ No End of Rotation? Yes End Yes ∆C p < 0.001 ? No Figure 3.15: The flowchart of computational procedure for DVM. 34 3.4: RANS simulation for the calculation of foil characteristics The lift and drag characteristics of the foils are required for the calculations in DVM. Due to the lack of experimental data for a wide range of angles of attack and Reynolds numbers, these data are obtained from RANS simulations. The commercial package FLUENT is used to obtain the required data. In order to get a reasonable result out of a RANS simulation, the domain of solution needs to be discretized into very small sections called cells or elements. The discretization of domain allows converting the Partial Differential Equations (PDE’s) to algebraic equations which can be solved numerically through an iterative process. The PDE’s that need to be discretized are continuity, momentum equation (in X and Y direction), energy equation (when using of a coupled solver), and turbulence equation(s). 3.4.1: Grid generation for the RANS simulations The GAMBIT software is used to generate the required mesh for the RANS simulation. An accurate prediction of lift and drag on the blades, especially for unsteady, separated flow requires a high quality fine grid around the foil. The size of grid required adjacent to the blades is dependent on the thickness of turbulent boundary layer. The boundary layer profile shown in Figure 3.16 can be divided into three different regions: the viscous sublayer, the buffer layer and the fully turbulent layer [61]. Figure 3.16: Profile of turbulent boundary layer [1]. The viscous sublayer is the innermost layer in which the viscous forces are dominant and the value of non-dimensional cell height, y+ is below 5 [61] where : 35 (3.38) In the buffer layer, the shear stress effects of viscosity and turbulence are of the same order of magnitude and the y+ value is between 5 and 60 [61]. The fully-turbulent layer, also called the log-law region, is where the effects of turbulence are dominant. In order to make sure that the boundaries far from the body have minimum effects on the solution, it is necessary to conduct the simulations in a large domain as shown in Figure 3.17. Thus, a very fine C-shape structured grid was used to resolve the flow around the foil and the wake structure as shown in Figure 3.18. Figure 3.17: Domain of solution in simulations. Figure 3.18: Fine grid around the foil in simulations. 3.4.2: Solver specifications for the RANS simulations The FLUENT solver is used for the solution of RANS equations in the domain described in the last section. A coupled, second order, implicit, and unsteady solver is used obtain a more accurate and a faster solution. Care was taken in the selection of time step as an improperly large time step affects the convergence of solution, while an overly large time step causes significant errors in the results. On the other hand, the choice of an extremely small time step can significantly increase the simulation time. 36 The Spalart-Allmaras turbulence model is used to model the turbulence in the flow. The main advantages of this turbulence model are the low computational cost, numerical stability and acceptable accuracy [61]. This model was first developed for aerospace applications but it is gaining popularity in the turbomachinery fields. However, as for many other turbulence models, it is incapable of perfectly modeling severe flow separations. The convergence criteria for the simulations were set to at least 10 −5 for the continuity, momentum and turbulence equations. 3.4.3: Results of the RANS simulations Since the blade profile of vertical axis turbine for experiments is NACA-63(4)-021 and regarding the sizing of turbine in experiments, the lift coefficient and drag coefficient of this foil are calculated at two Reynolds numbers: 200,000 and 500,000. The results of RANS simulations for this foil at angles of attack up to 40 degree are shown below. 1.2 1 1.0 0.8 0.8 0.6 0.6 CL CD 1.2 0.4 0.4 CL at Re=500,000 CL at Re=200,000 CD at Re=500,000 CD at Re=200,000 0.2 0.2 0 0.0 0 10 20 Angle of attack(degree) 30 Figure 3.19: The characteristics of NACA 63(4)-021. 37 40 Chapter Four Ducted Vertical Axis Turbine 4.1: Introduction Numerous numerical and experimental results [1,3,4] have shown that higher power can be extracted from a ducted turbine by concentrating the energy of a larger area into a smaller area using a venturi duct shape. However, in order for the ducted turbine to become economically viable, the duct structural costs must be more than offset by the advantages of increased power and annual energy production. A schematic of a ducted turbine is shown in Figure 4.1. In this chapter a new method is developed for performance calculation of ducted turbines which can be used for designing efficient and cost effective ducts. Ω Figure 4.1: Schematic of a ducted turbine. 38 4.2: Boundary element methods for internal-external flows Boundary element methods have been used for years to numerically predict flow fields around arbitrary three-dimensional bodies. In these methods, the geometry of the body is represented by a set of elements or panels. Singularities are distributed over these panels to create the flow field around the body. The singularities may be sources, vortices, or some combination of these [62]. If the singularities are distributed with constant strength over each panel, a low order boundary element method results. If the singularities are distributed with a linearly or quadratically varying strength over each element, a high order boundary element method results. High order boundary element methods tend to be more accurate than the low order boundary element methods because they model better a continuous singularity strength distribution over the body. However, high order methods require significantly longer computation times and care must be taken to ensure that all panels representing the body match up exactly to satisfy the requirement that the singularity strength must be continuous across panels. The long computation times required for high order methods can be a limiting factor for complicated geometries; hence, there is considerable interest in improving the accuracy of low order methods without sacrificing their shorter computation time. Typically, boundary element methods are used to model bodies such as an aircraft in an external uniform flow. But application of these methods to internal flow problems has been limited mainly to higher order methods. This is because most of the low order methods suffer severely from leakage when applied to internal flows. As discussed by Holt and Hunt [63], leakage manifests itself as a non-constant mass flux down the length of tunnel. Leakage occurs because the boundary conditions are applied only at the control point of each panel and the flow is free to leak everywhere else. The leakage problem is more pronounced [62] for geometries such as a wind tunnel that has a contraction area, a test section and a diffuser because the tunnel walls typically have to turn the flow and the contraction and diffuser impart large changes in velocity down the length of the wind tunnel. Following some improvements in the mathematical basis of boundary element methods, Ashby et al. [62] developed a low order boundary element method that shows a good mass continuity inside the duct. This method is used for the calculation of potential flow field around the duct in this research. The details of this method are explained in the following sections. 39 4.2.1: General solution for potential flow problem Consider a body with known boundaries S B submerged in a potential flow in a region R out of the body as shown below. n R y Body Coordinates l→∞ x n γn SB γk r r n p 2ε S∞ Figure 4.2: A body submerged in a potential flow. The governing equation for the potential flow in the body frame of reference, in terms of velocity potential φ is: ∇ 2φ = 0 (4.1) The boundary condition on the surface SB for the solution of equation 4.1 is: r ∇φ ⋅ n = V n (4.2) r Where Vn is the normal velocity to the surface, n is a vector normal to the surface pointing outwards in the region R , and ∇φ is measured in a fixed frame of reference ( x, y ) . For internal or external flows Vn is zero for all stationary solid boundaries, and is non-zero in the inlets and outlets. The disturbance created by the motion of body must decay at far away from the body, this means: r lim ∇φ − V f = 0 ( ) (4.3) l →∞ r where l represents a point away from the origin, and V f is the undisturbed velocity far away r from the body and can be equal to the free stream velocity V∞ . 40 By using the Divergence theorem [52], the equation 4.1 can be used to obtain: r ∫∫ ∇ φdV = ∫∫ ∇ ⋅ ∇φdV = ∫ ∇φ ⋅ ndS 2 R R (4.4) S where S is the boundary of domain R in the absence of any free vortices and is equal to S = S B + S ∞ . The boundaries S B and S ∞ are shown in Figure 4.2. Let the vector ∇φ be replaced by the vector φ1 ∇φ 2 − φ 2 ∇φ 1 , where φ1 , φ 2 are the two scalar functions of position satisfying Laplace’s equation 4.1. We can then write: ∇ 2φ = ∇ ⋅ ∇φ = φ1∇ 2φ 2 − φ 2 ∇ 2φ1 By replacing ∇ 2φ from above, the equation 4.4 becomes as follows: ∫∫ ∇ ⋅ ∇φdV =∫∫ (φ ∇ φ 2 R R 1 2 ( ) r − φ 2 ∇ 2φ1 )dV = ∫ φ1 ∇φ 2 − φ 2 ∇φ1 ⋅ ndS S (4.5) In particular choose the functions φ1 and φ 2 to be as follows: φ1 = ln r and φ 2 = φ where φ is the flow potential of interest in the region R and r is the distance from a point p (see Figure 4.2). When the source point p is outside of R, both φ1 and φ 2 satisfy Laplace equation 4.1, then the equation 4.5 becomes: r ∫ (ln r∇φ − φ∇ ln r ) ⋅ ndS = 0 (4.6) S However, the function φ1 is unbounded when point p is inside the region as r → 0 . Then the point p must be excluded from the region of integration as shown in Figure 4.2, where the point p is excluded by a small circle of radius ε . In the remaining of region R , the potential φ1 satisfy Laplace's equation as ∇ 2 (ln r ) = 0 . Similarly ∇ 2φ 2 =0, then the equation 4.6 becomes: r ( ln r∇φ − φ∇ ln r ) ⋅ ndS = 0 S + circle ε ∫ (4.7) 41 To evaluate the integral over circle with radius ε , we introduce a polar coordinate system at p r r ∂φ v 1 and since the vector n points inside the small circle, n = −er , n.∇φ = − and ∇ ln r = er . ∂r r Then the equation 4.7 becomes: r ∂φ φ − dS + ∫ (ln r∇φ − φ∇ ln r ) ⋅ n dS = 0 ln r circle ε S ∂r r −∫ Assuming that the potential φ and its derivatives (4.8) ∂φ are well-behaving functions as ε → 0 (i.e. ∂r they do not vary much in the small circle), the first term in the first integral vanishes, while the φ dS = 2πφ ( p ) . circle ε r second term yields ∫ Then from the equation 4.8, the potential φ at point p is as: , φ ( p) = − 1 2π r ∫ (ln r∇φ − φ∇ ln r ) ⋅ ndS (4.9) S Now consider potential function φi for a fictitious flow that occurs out of region R (inside the boundary for external flows and outside of boundary for internal flows). For this flow the point p which is in the region R is exterior to S B , then the equation 4.6 yields: − 1 2π ∫ (ln r∇φ SB i r − φ i ∇ ln r ) ⋅ ni dS = 0 (4.10) r r r Here ni points outward from S B , ni = −n . Since the boundary S is equal to S = S B + S ∞ as shown in Figure 4.2, then the potential at point p can be written as follows by adding the equations 4.9 and 4.10: φ ( p) = − 1 2π 1 − 2π r ∫ [ln r (∇φ − ∇φ ) − (φ − φ )∇ ln r ]⋅ ndS i i SB (4.11) r ∫ (ln r∇φ − φ∇ ln r ) ⋅ ndS S∞ 42 The contribution of integral around S ∞ in equation 4.11, when S ∞ is considered to be far from S B is defined as: φ∞ ( p ) = − 1 2π r ∫ (ln r∇φ − φ∇ ln r ) ⋅ nds (4.12) S∞ The potential φ ∞ usually depends on the selection of coordinate system. Then the equation 4.11 becomes: φ ( p) = − 1 2π r ∫ [ln r (∇φ − ∇φ ) − (φ − φ )∇ ln r ] ⋅ ndS + φ ( p ) i i ∞ (4.13) SB r It is more common to consider the positive direction of vector n into the region of interest. Then r r we replace vector n with vector − n . Also it is advantageous to reverse the direction of vector r. Note that by reversing the direction of r only the sign of ∇ ln r changes. Therefore the potential at point p can be written as: φ ( p) = 1 2π r ∫ [ln r (∇φ − ∇φ ) + (φ − φ )∇ ln r ] ⋅ ndS + φ ( p ) i i ∞ (4.14) SB The potential φ s = 1 ln r is the potential of a source with unit strength. Besides, the following 2π relation exists between the source and doublet potentials of unit strength [52]: φd = − ∂φ s r r 1 = −∇φ s ⋅ n = − ∇ ln r ⋅ n 2π ∂n The equation 4.14 can now be written in terms of source and doublet potentials of unit strength as follows: r φ ( p ) = ∫ φ s (∇φ − ∇φi ) ⋅ ndS − ∫ (φ − φi )φ d dS + φ ∞ ( p ) SB SB (4.15) The equation 4.15 represents the potential at point p, φ ( p ) , as a distribution of sources on the boundary S B with strength σ per unit length and doublets on the boundary with strength µ per unit length. 43 Using the equation 4.15, the source strength σ and the doublet strength µ are: r σ = (∇φ − ∇φi ) ⋅ n = ∂φ ∂φ i − ∂n ∂n (4.16) µ = (φ − φi ) Equation 4.16 shows that the source strength σ is equal to the jump of ∂φ on the boundary and ∂n the doublet strength µ is equal to the jump in φ on the boundary (see Figure 4.3). Figure 4.3: The potential and its gradient inside and outside of body. Therefore from equations 4.15 and 4.16, the potential at each point p at region R is: φ ( p) = 1 2π ∫ SB ∂ σ ln r + µ ∂n (ln r )dS + φ ∞ ( p ) (4.17) r The direction of µ is the same as direction of n (into the flow field of interest). The governing equation for a fictitious potential flow φi inside the body is: ∇ 2φ i = 0 (4.18) The same method described above gives the potential value for a point p, inside the body as: φi ( p ) = 1 2π ∫ SB ∂ σ ln r + µ ∂n (ln r )dS + φ ∞ ( p ) (4.19) 44 4.2.2: Boundary conditions The boundary condition given by equation 4.2 directly specifies a normal velocity component Vn on surface S B , this direct formulation is called the Neumann problem. It is also possible to specify φ on the boundaries, this indirect formulation is called the Dirichlet problem. A combination of the above boundary conditions is called a mixed boundary condition problem. 4.2.2.1: Neumann problem Since the solution given by equation 4.17 should satisfy the boundary condition of equation 4.2, we can write the boundary condition as: r 1 ∇φ ⋅ n = 2π r ∂ σ∇ ln r + µ∇ (ln r )dS + ∇φ∞ ( p ) ⋅ n = 0 SB ∂n ∫ (4.20) The equation 4.20 is the basis for many numerical solutions and should hold for every point on the surface S B . This equation reduces to a set of algebraic equations by writing it for a number of points (called collocation points) on the surface S B . Since the numerical formulation of equation 4.20 is not unique and many problems can be solved by using a preferred type of singularity or any linear combination of two singularities, in many applications additional considerations such as the Kutta condition, is required to obtain a unique solution. 4.2.2.2: Dirichlet problem For an impermeable boundary ∂φ = 0, as required by the boundary condition 4.2. One can show ∂n that the potential inside body φi remains constant if there is no internal singularity [52]. Then the Dirichlet boundary condition for an impermeable boundary can be written as: φi = const (4.21) Then from equations 4.17 and 4.21, the inner potential φi for a point p inside the boundary S B can be written as follows: φi = 1 2π ∂ σ ln r + µ (ln r )dS + φ∞ ( p ) = const SB ∂n ∫ 45 (4.22) The equation 4.22 is the basis for methods utilizing Dirichlet boundary condition. However, there are many differences between various numerical methods of solution, depending on the value selected for the inner potential φi . For example, by setting φi = 0 , the resulting singularity distribution from equation 4.22 will include φ∞ and the strength of singularities will be large. Other values for the inner potential (not necessarily constant) can be specified as well. In a special case when the inner potential is set to φi = φ ∞ , the equation 4.22 reduces to a simple form as given below: 1 2π ∫ SB ∂ σ ln r + µ ∂n (ln r )dS = 0 (4-23) The potential φ can be assumed to consist of two components, φi the potential of fictitious flow ~ and a perturbation potential φ as: ~ φ = φi + φ Then using the Neumann boundary condition, (4.24) ∂φ = 0 , we obtain: ∂n ~ ∂φ r ∂φ = − i = − n ⋅ ∇φ i ∂n ∂n (4.25) As the source strength σ is equal to the jump in normal derivative of velocity potential across the surface as given by equation 4.16, then: ~ ∂φ ∂φi ∂φ = − σ= ∂n ∂n ∂n (4.26) Comparing equations 4.25 and 4.26, the source strength σ based on the value of fictitious potential φi is: r σ = − n ⋅ ∇φ i (4.27) 46 4.2.3: Kutta condition The analysis up to this point described the method for solving the Laplace equation 4.1. This method consists of solving both equations 4.22 and 4.27 to find the strength of source and doublet at each point on the surface S B . However this solution is not unique. To obtain a unique solution, W.M. Kutta suggested considering the physics of flow by assuming a wake model, [52]. The simplest wake model assumes that the flow leaves the sharp trailing edge of a foil smoothly, with a finite velocity along the bisector line. This means that the pressure difference at the trailing edge is zero. In other words if the circulation around the foil is modeled by a vortex distribution, the Kutta condition requires that the vortex strength at the trailing edge should be zero [52]: γ T .E = 0 (4.28) In a two dimensional flow ∂µ ( x ) ∂S = −γ (S ) [52], where S is the line bisecting the trailing edge angle δ T . E . as shown in Figure 4.4. Then the Kutta condition can be rewritten as: ∂µ µ w − µ T . E . = = 0 → µ T .E. = µ w ∂S ∆S Where µ w is the strength of doublet in the wake which is assumed to be constant (see Fig. 4.4). And µ T .E . is the doublet strength at the trailing edge and is equal to µ T .E . = µU − µ L . This leads to the following relation for the doublets at the trailing edge: µU − µ L − µ W = 0 (4.29) Where µU and µ L are the corresponding upper and lower surface doublet strengths at the trailing edge as shown in Figure 4.4. Figure 4.4: Flow leaves the T.E. along the bisector line based on kutta condition. 47 4.2.4: Modeling of the vortices in the flow field Equation 4.17 gives the velocity potential resulting from a free stream at each point p under specified boundary conditions. In the case that there are some discrete vortices in the flow field as shown in Figure 4.2, the total velocity potential at each point p can be obtained using the principle of superposition as follows: Nv φt ( p ) = φ ( p ) + ∑ φ vl ( p ) (4.30) l =1 Here φt ( p ) is the total potential at point p. Also φ ( p ) and φ vl ( p ) are the potential at point p due to free stream and l th discrete vortex in the flow field respectively. The total potential φt ( p ) is still analytic in all of region R except at the position of vortices in the flow field; this means that the Laplace equation still holds: ∇ 2φ t ( p ) = 0 (4.31) By summation of the equations 4.17 and 4.30, the total potential at a point p in the region R can be expressed as: φt ( p ) = 1 2π Nv ∂ ( ) ( ) σ ln r µ ln r dS + φ p + φ vl + ∑ ∞ ∫SB ∂n l =1 (4.32) Similarly the total potential for a point p inside the body can be calculated by superposition as: (φt )i Nv = φ i + ∑ φ vl (4.33) l =1 Substituting the potential φi from equation 4.19 in 4.33, the total potential at a point p inside the body is: (φt )i ( p ) = 1 2π Nv ∂ ( ) ( ) σ ln r µ ln r dS + φ p + φ vl + ∑ ∞ ∫S B ∂n l =1 (4.34) In equations 4.32 and 4.34, the definitions of σ and µ remain the same as in equations 4.16. This implies that the amount of jump in the potential across the surface, µ , does not change when there are some discrete vortices in the flow as vortices induce equal amount of potential across the boundary surface, when the assumption of thin boundary surface is valid. 48 When some discrete vortices exist in the flow field, the boundary condition for normal velocity equation 4.2 changes as: r ∇φ t ⋅ n = 0 (4.35) From equations (4.16),(4.30) and (4.35), the source strength σ can be calculated as follows: r Nv r l =1 σ = (∇φ − ∇φi ) ⋅ n = ∇φt − ∑ ∇φvl − ∇φi ⋅ n r Nv r l =1 σ = (∇φ − ∇φi ) ⋅ n = − ∑ ∇φ vl − ∇φi ⋅ n (4.36) Using similar arguments leading to equation 4.21, for an enclosed boundary the gradient of potential in the normal direction is zero as required by the boundary condition in equation 4.35, the potential inside the enclosed boundary is constant. (φt )i = const (4.37) Equating the equations 4.34 and 4.37 gives the following equation that needs to be solved in order to find the potential in the domain: 1 2π Nv ∂ ( ) ( ) σ ln r µ ln r dS φ p φ vl = const + + + ∑ ∞ ∫SB ∂n l =1 (4.38) The value of the constant in equation 4.38 is crucial for the accuracy of results and for decreasing the numerical leakage. By choosing this constant as φ∞ , the strength of calculated singularities are much smaller and the numerical leakage decreases [62]. The source strength can then be calculated using the equation 4.36: Nv r l =1 σ = − ∑ ∇φ vl − U ∞ ⋅ n In summary, for the calculation of potential flow field when there are some discrete vortices in the flow field, we need first to calculate the source strength on S B from equation 4.36 and then calculate the doublet strength on S B using the equation 4.38. 49 By breaking the surface S B into a series of panel, the equation 4.38 can be written as follows for collocation point of each panel (see Figure 4.5). Np Np 1 1 ∂ Nv ∑ 2π ∫ σ ln rdS + ∑ 2π ∫ µ ∂n (ln r )dS + ∑ φ j j =1 j j =1 l v =0 (4.39) l =1 Figure 4.5: The body surface is represented by some panels. The doublet and source strength can be factored out of integrals in equation 4.38 by assuming a constant strength source and doublet distribution over each panel (making it a low order panel method). Then the equation 4.39 transforms to a system of algebraic equations for a series of panels as follows: Np ∑A σ ij j Np Nv j =1 l =1 + ∑ Bij µ j = ∑ C il j =1 i = 1,..., Np (4.40) Where the coefficients Aij , Bij and C il are: Aij = 1 ln rdS 2π ∫j Bij = 1 2π ∂ ∫ ∂n (ln r )dS (4.41) j C il = −φ vl = Y − Yl γv tan −1 i 2π X X − i l The coefficients Aij and Bij are integrals over the surface of panel j. They represent the influence of a constant distribution of sources and doublets with unit strength placed on panel j acting on panel i (Fig. 4.5). C il is the potential induced by lth vortex on the panel i. Note that ( X i , Yi ) and ( X l , Yl ) are in the global coordinate system. 50 4.2.5: Calculation of influence coefficients The potentials of a source and a doublet with unit strength are: φs = 1 ln r 2π (4.42) r 1 φ d = − ∇ ln r ⋅ n 2π Thus, the velocity potential at the collocation point of panel i (point p) due to a source on panel j using the local coordinates of panel j can be written as: φs = 1 ln (ζ − ζ 0 ) 2 + η 2 2π (4.43) Where (ζ 0 ,0) and (ζ ,η ) are the local coordinates of singularity and point p respectively as shown in Figure 4.6. η Panel i p r n r r χ Panel j r er ζ1 ζ 0 ζ2 ζ Figure 4.6: Panel i in the local coordinates of panel j . Using a polar coordinates system, gradient of ln r can be written as: ∂ ln r r 1 ∂ ln r r ∂ ln r r 1r er + er = e r eθ = r ∂θ r ∂r ∂r r Since the direction of vector r is towards the point p, then: ∇ ln r = r r er ⋅ n = cos χ The velocity potential of a doublet on panel j acting on a point p can be written as: φd = − 1 η 2π (ζ − ζ 0 )2 + η 2 (4.44) 51 4.2.5.1: Constant-strength source distribution By using equation 4.43, the velocity potential at a point p due to a distribution of sources with unit strength placed on panel j is (see Figure 4.6): Aij = 1 2π ζ2 ∫ζ ln (ζ − ζ 0 ) 2 + η 2 dζ 0 = 1 1 4π ζ2 ∫ζ ln((ζ − ζ 0 ) ) 2 + η 2 dζ 0 1 where ζ 1 and ζ 2 are the local coordinates of edges of panel j as shown in Figure 4.6. Calculation of this integral gives φ s in the local coordinates of panel j as follows: Aij = 1 (ζ − ζ 1 ) ln (ζ − ζ 1 ) 2 + η 2 − (ζ − ζ 2 ) ln (ζ − ζ 2 ) 2 + η 2 − 2(ζ 2 − ζ 1 ) 4π η η + 2η tan −1 − tan −1 ζ −ζ 2 ζ − ζ 1 [ ( ) ( ) (4.45) The velocity components induced by this source distribution at point p are: (ζ − ζ 1 ) + η 2 ∂ 1 φs = u= ln ∂ζ 4π (ζ − ζ 2 )2 + η 2 2 θ1 647 4θ 2 48 4 647 48 η η ∂ 1 −1 − tan −1 φs = v= tan ζ −ζ 2 ζ −ζ1 ∂η 2π (4.46) Note that these velocity components are at local coordinates of panel j. When the point p is on the panel ( r = 0 ), in this case the velocity components become: u (ζ ,0 ± ) = 2 1 (ζ − ζ 1 ) ln 4π (ζ − ζ 2 )2 The ζ velocity component diminishes at the center of panel and is infinite at panel edges. But for evaluating the η component of velocity when the point p is on the panel j, it is important to distinguish between the conditions when the panel is approached from its upper or from its lower side. When the point p is above the panel, θ1 = tan −1 η ζ −ζ1 = 0 and θ 2 = tan −1 These conditions are reversed on the lower side, then: v(ζ ,0 ± ) = ± 52 1 2 η ζ −ζ2 =π . 4.2.5.2: Constant-strength doublet distribution From equation 4.44, the velocity potential at a point p on panel i due to a doublet distribution with unit strength on panel j is: Bij = 1 2π ζ2 η ∫ζ (ζ − ζ ) 1 0 2 +η 2 dζ 0 Again (ζ 0 ,0) and (ζ ,η ) are the local coordinates of singularity and point p respectively, as shown in Figure 4.6. With calculation of this integral, the potential at point p due to this distribution of doublets becomes as follows: Bij = 1 2π −1 η η − tan −1 tan ζ −ζ 2 ζ − ζ 1 (4.47) The velocity components in ζ and η directions at point p induced by the doublet distribution are: u= v= ∂φ d 1 η η = − 2π (ζ − ζ 1 )2 + η 2 (ζ − ζ 2 )2 + η 2 ∂ζ (4.48) ∂φ d ζ −ζ1 ζ −ζ 2 1 − =− 2 2 2 2 2π (ζ − ζ 1 ) + η ∂η (ζ − ζ 2 ) + η When the point p is on the panel ( r = 0 ), the velocity component at ζ direction vanishes as η → 0: u (ζ ,0 ± ) = 0 The induced velocity component at η direction when the point p approaches the panel j as η → 0 is: 1 1 1 v(ζ ,0 ± ) = − − 2 ζ − ζ 1 ζ − ζ 2 It is clear from the above equation that the induced velocity component at η direction is singular at the panel edges. 53 4.3: Development of the BEM-DVM for ducted vertical axis turbines DVM model was developed for performance analysis of unducted vertical axis turbines in chapter three. In addition, BEM model was presented for the calculation of the potential flow field in internal and external flows in the first section of this chapter. Here, a new model referred to as “BEM-DVM” is developed for the performance analysis of a ducted vertical axis turbine based on BEM model and DVM model. The boundary layer of duct is not taken into account in this model. A schematic of modeling for a ducted turbine is shown in Figure 4.7. In the development of this model, the same assumptions as in DVM model are used (see section 3.3). The effects of ducting on turbine are taken into account using BEM model; however, it is assumed that ducting influences only the quasi-steady components of loads and circulation. Using this assumption simplifies adjusting the computational procedure of the unducted turbine for a ducted turbine. In each stage of computational procedure for BEM-DVM, the following steps are taken: • Calculation of sources strength, by including the effect of discrete vortices, vortex filaments as given in equation 4.36 • Calculation of doublets strength by solving the equation 4.39 • Calculation of the flow velocity, the angle of attack of blades by including the effect of sources and doublets • Calculation of induced velocity on each discrete vortex by including the effect of sources and doublets Figure 4.7: Schematic of modeling for a ducted vertical axis turbine. 54 In the computational procedure for ducted turbine, the difference of power coefficient based on the frontal area of turbine at two successive rotations is used as convergence criteria; the solution continues until the convergence criteria becomes less than 0. 001. At each time step ∆t , the turbine rotates by an angle ∆θ = 10 degree. The flowchart of calculation for the BEM-DVM is shown in Figure 4.8. Initial Data Calculate the strength of sources Calculate the strength of doublets Calculate the flow velocity, the angle of attack and the circulation of blades Calculate the strength and the position of nascent vortices Repeat last steps including nascent vortices Calculate the induced velocity on all free vortices Calculate the tangential force, the normal force and the torque of blades Calculate the new position of blades and free vortices θ new = θ old + ∆θ No End of Rotation? Yes End Yes ∆C p < 0.001 ? Figure 4.8: The flowchart of BEM-DVM. 55 No 4.4: A method for the calculation of tank wall effects on the turbine performance The experimental results from testing a turbine in a towing tank or wind tunnel include the wall effects. It is important to have a good estimation of wall effects for design purposes. One can calculate the tank wall effects using BEM-DVM model developed in the last section, however, the source strength needs to be calculated from the following relation: r Nv r l =1 σ = (∇φ − ∇φi ) ⋅ n = Vn + − ∑ ∇φ vl − ∇φi ⋅ n (4.49) Where Vn is the normal velocity on the boundaries which is zero for impermeable surfaces and non-zero for input and output of tank (Figure 4.9). φvl is the velocity potential of l th discrete vortex and φi is the velocity potential of fictitious flow. Figure 4.9: Schematic of a ducted turbine operating in a tank. In order to have kinematic similarity between the two-dimensional numerical model of a turbine operating in a tank and the three-dimensional experiments, it is required to conduct the numerical calculations in a tank with a width that give the same blockage ratio as that of experiments. The blockage ratio which is a characteristic of amount of interaction between the wall and the device is defined as follows: Λ = Ad , p ,a At .c.a here At .c.a is the cross section area of tank and Ad . p.a is the net projected area of device and its supports, Ad . p.a can be obtained by projecting the device on a plane perpendicular to the direction of free steam flow. In chapter 6, the blockage ratio used in the calculations is provided. 56 Chapter Five Duct Optimization 5.1: Introduction Numerical optimization methods have become more practical in engineering design due to increased efficiency of the optimization algorithms and the computational power of computers. Optimization methods are used to find a set of design parameters, which can in some way be defined as optimum. In a simple case, this might be the minimization or maximization of some system characteristics that depend on the design parameters x. In a more advanced formulation the objective function, f ( x ) , to be minimized or maximized, might be subject to constraints in the form of equality constraints, inequality constraints and some bounds on the parameters. A general description of an optimization problem is as: minimize x subject to g i (x ) = 0 g i (x ) ≤ 0 f (x ) (5.1) i = 1,..., me i = me + 1,..., m 57 Where the vector function g ( x ) returns a vector of length m , containing the values of equality and inequality constraints evaluated at x . In general there are two main classes of optimizers: • Search-based methods • Gradient based methods A brief explanation of these methods is provided below. A complete explanation of these methods can be found through references [64,65,66]. 5.1.1: Search-based methods Search-based methods, such as genetic algorithms, are able to handle a broad range of optimization problems, including global optimization. These methods are usually easy to implement; however, they are slow in terms of the number of flow field evaluations required. Furthermore, they are not efficient near an optimum solution. The details of this method are beyond the scope of this research. 5.1.2: Gradient-based methods Gradient-based methods are used far more in optimization routines than search-based methods. While they are efficient near an optimum solution, they still require many flow field evaluations and are not able to handle as broad range of optimization problems as search-based optimizers. The key issue in developing an efficient gradient-based optimizer is reducing the time required for computing the gradient. Usually the gradient is found using a finite-difference approach; however, the associated computational cost is high. Another method for the calculation of the gradient of a function is automatic differentiation which uses the chain rule for the calculation of the gradient of a function. This method gives more accurate calculation of the gradient and is more widely used then the finite difference approach. The optimization method used in this research, is a gradient-based method referred to as “Sequential Quadratic Programming”, SQP, which is explained briefly in appendix A. SQP methods represent the state of the art in optimization of nonlinear functions. Much more about these methods are available in references [64,65]. This optimization method uses a finite-difference approach for finding the gradient of objective function. 58 5.2: Duct optimization for ducted vertical axis hydro current turbines Here, a numerical optimization approach is developed in order to find the duct shape that maximizes the performance (i.e. Cp) of a turbine under a number of constraints. The details of this approach are explained in the following sections: 5.2.1: Duct representation Parametric representation of a model using a number of design parameters is essential in reducing the optimization effort while maintaining a high level of accuracy and flexibility in representing the model. Because of this parametric representation the Bezier curves and, in particular, the B-spline curves are widely used for optimization purposes. Essentially, these curves are piecewise polynomials which are explained in appendix B. A complete description of these curves can be found through references [67,68]. Since the turbine performance is optimized through modification of the duct shape, the design variables are selected to be the Y coordinate of control points which represent the duct shape. At each level of the optimization process, the duct is created by fitting B-splines to the control points. Two kinds of ducts are considered in the optimization study; symmetric ducts and asymmetric ducts. For creation of the duct at each stage of the optimization process, the bottom side of the duct is created, and then the upper side of the duct is obtained by mirroring the bottom side around the X axis. For symmetric ducts, one B-spline is used for defining the curved part of bottom side as shown in Figure 5.1, while the straight part is created using a straight line. Y WD X LD Figure 5.1: Schematic of a symmetric duct. 59 In the optimization process of asymmetric ducts, two B-splines are used for defining the duct’s bottom side as shown in Figure 5.2. Y WD X X LE LD Figure 5.2: Schematic of an asymmetric duct. After the creation of duct at each stage of optimization process, a number of panels are distributed over the boundaries of duct. The concentration of the panels is more in the leading edge and the trailing edge of the duct. 5.2.2: Objective function The power coefficient of the turbine with negative sign (–Cp) is used as the objective function; as the optimizers usually find the minimum of a function, while we want to maximize the power coefficient of the turbine. This power coefficient is still calculated based on the frontal area of the turbine. 5.2.3: Optimization constraints The constraints used in this optimization study are only geometrical; no other constraints such as environmental constraints or cost constraints are considered. These geometrical constraints depend on the width, WD, and the length, LD, of the duct, which are fixed at each optimization case. The constraints used for different optimization cases are provided in the next chapter. 60 5.2.4: Optimization methodology The commercial package MATLAB is used for this optimization study as it has the Spline and Optimization toolboxes and a user-friendly interface. Some subroutines were developed in MATLAB to make an automatic optimization process. In the optimization process, the difference of the objective function between two successive cycles is used as convergence criteria. The optimization continues until the convergence criteria is less than 0.001. The flow chart of the optimization process is shown in Figure 5.3. Initial Design Variables Calculate Splines (MATLAB Spline Toolbox) Put Panels on the Duct (180 panels on each side of duct) Calculate the Performance of Turbine (BEM-DVM Code) Yes ∆C p < 0.001 ? No Find New Set of Design Variables (MATLAB Optimization Toolbox) End Figure 5.3: Flowchart of the optimization process. 61 Chapter Six Results and Discussion 6.1: Introduction The numerical models developed for unducted and ducted vertical turbines were explained in Chapters 3 and 4. The methodology developed for duct shape optimization was also explained in Chapter 5. In this chapter, the results of numerical models are compared with the experimental and RANS simulation results, followed by the optimization results. 6.2: Experimental testing of a vertical axis turbine at UBC towing tank Researchers at the Naval Architecture and Maritime Laboratory of the University of British Columbia conducted three rounds of experimental studies on a vertical axis turbine in the towing tank at UBC; the author formed part of the team who conducted these experimental tests [3]. The effects of different turbine parameters as well as ducting effects on the turbine performance were investigated in these experiments. The results of these experiments have been used for the validation of numerical models. The characteristics of the turbine model tested in these experiments are given in Table 6.1. The characteristics of this turbine which is referred to as “UBC-Turbine” are similar to the model used in the experimental studies conducted at NRC Canada in the 1980’s [3]. 62 Table 6.1: The characteristics of the UBC-Turbine. PARAMETER DIMENSION / CHARACTERISTIC Diameter (across foil chord) of the turbine (D) 36 inch Number of blades (N) 3 Blade span (H) 27 inch Blade profile NACA 634-021 Blade chord length (C) 2.70 inch(ideal); 2.57 inch (manufactured) Outer shaft diameter 1.9 inch The dimensions of the towing tank at the University of British Columbia were 200 feet long, 12 feet wide and 8 feet deep. Figure 6.1 shows the towing tank and the test carriage spanning the tank width. The test carriage was used as a platform for mounting the turbine model and the required instrumentation for the experiments. Figure 6.1: The towing tank facility at UBC [3]. The experimental setup is composed of attaching the turbine blades to a central shaft that is supported by two ball bearings at the top where they are mounted to a force balance and the 63 other bearing located at the shaft bottom as shown in Figure 6.2. The force balance consists of two parallel plates which are rigidly connected together with a hole, where the shaft goes thorough as shown in Figure 6.3. The bottom bearing is constrained to a horizontallymounted bottom plate supported by two vertical rectangular beams forming a U frame. Each of these beams is bolted to the carriage and stiffened using three guy wires. The turbine assembly with the arms supporting the blades at the top and bottom locations as well as its position within the tank are shown in Figure 6.2. This figure also shows the supporting frame and the force balance for mounting the instrumentation. Shaft Arms U frame Blades Figure 6.2: The turbine assembly (all dimensions are in inches) [3]. Figure 6.3: Gearbox drive-train configuration [3]. 64 Figure 6.3 also shows the drive and control mechanism of the turbine shaft, done by an AC motor attached to a 20:1 gearbox [3]. The motor and the gearbox are mounted on the top plate of the force balance. A flexible coupling was used to couple the torque sensor and the gearbox. The encoder which provides the angular position of a reference point on the turbine shaft, was mounted directly around the main shaft above the bearing of the lower plate. The duct tested in these experiments is a venturi-type duct as shown in Figure 6.4. This duct is obtained from previous NRC trials [3]. Figure 6.4: Symmetric duct dimensions (all dimensions are in inches) [3]. Figure 6.5 illustrates a schematic view of the ducted turbine in the towing tank. The duct position within the tank, and the top and bottom plates are also shown in this figure. At the center of the top plate, a large Plexiglas window was installed to allow for installation and removal of the turbine while leaving the ducting in place, as well as to facilitate visualization of the flow and the turbine. As Figure 6.5 shows, the ducted turbine in the tank has a much larger blockage ratio than the unducted turbine the tank as shown in Figure 6.2. One can argue that this ratio is an indicator of the amount of interaction of the turbine assembly with the tank wall; the higher this ratio, the greater the interaction would be. 65 Figure 6.5: Schematic of ducted turbine in the tank (dimensions are in inches) [3]. 6.3: Comparison of the numerical results with the UBC experimental results Since the experimental results include the towing tank wall effects, the numerical results presented here include these effects using the method explained in section 4.3. The blockage ratios for unducted and ducted turbines during the experiments were 0.012 and 0.08 respectively. From these blockage ratios, the tank width of 12.5 meter is used in all numerical calculations for both unducted and ducted turbines. The tank length of 60 meter is also used to make sure all vortices are between the tank walls during the calculations, and the turbine is located 20 meters down the tank inlet as shown in Figure 6.6. 20 m Ω 60 m Figure 6.6: A schematic drawing of the ducted turbine operating in a tank. 66 To facilitate the comparison of numerical and experimental results, the experimental Cp values for the arms-only have been added to the numerically calculated Cp values. The dimensions of the supporting arms which are from profile NACA- 0012 are given below. Figure 6.7: The profile and dimensions of the arms (all dimensions are in inches) [3]. Figure 6.8 gives the experimental Cp values of the arms-only for velocities U∞=1.5 m/s and U∞=2 m/s. One can observe that the parasite drag of the arms substantially increases at high tip speed ratios. It is important to note that the vortex shedding of the shaft and the interactions between the arms, shaft and the blades are still ignored when we add the experimental Cp values for the arms to Cp values from the numerical calculations. Figure 6.8: Experimental Cp vs. TSR values for the supporting arms [3]. Figure 6.9 shows the comparison of the Cp results for unducted turbine at velocity U∞=1.5 m/s. The RANS simulation results in this figure are from the studies conducted by Nabavi[1]. One can observe that the BEM-DVM results are in good agreement with the experimental and RANS simulation results at TSR values above 2.5. The BEM-DVM also gives a good estimation of the TSR value at which the turbine shows its best performance. At TSR=2.75 67 where the turbine performs its best, Cp value of 0.29 is obtained from the experiment while Cp values of 0.31 and 0.32 are obtained from the BEM-DVM and the RANS simulations respectively. However, the BEM-DVM model under-predicts the turbine performance at TSR values less than 2.75. This is possibly due to lack of dynamic stall modeling which is dominant at this range of TSR values. 0.5 BEM-DVM FLUENT EXPERIMENT 0.4 Cp 0.3 0.2 0.1 0.0 2 2.25 2.5 2.75 3 3.25 3.5 TSR Figure 6.9: Results of Cp vs TSR for unducted turbine at U∞=1.5 m/s. The torque variations of the turbine over one revolution at TSR=3 and U∞=1.5 m/s are given in Figure 6.10. For the torque curves, a blade is considered to be at 0 degrees when it is headed directly into the flow, and is at 180 degrees when it is moving in the same direction as the flow (see Figure 3.5). As Figure 6.10 shows, the turbine has high torque fluctuations with three peaks over one revolution. The position of the maximum torque predicted by BEM-DVM model is in good phase agreement with the FLUENT prediction. However, the maximum experimental torque happens at slightly higher angles because of some issues related to phase measurement during the experiments. On the other hand, the torque values calculated by both numerical methods are below the experimental torque results. This is possibly due to lack of dynamic stall modeling in BEM-DVM model, and inability of RANS simulations to model the dynamic stall phenomenon accurately. 68 The torque of one blade over one revolution at TSR=3 is also shown in Figure 6.11. Higher torque is obtained in the upstream region of the turbine shaft than the downstream region. BEM-DVM under-predicts the RANS simulation results when the blade is in the upstream region of the turbine shaft, while it over-predicts the RANS simulation prediction when the blade is in the downstream region of the turbine shaft. Figure 6.10: The torque of unducted turbine at TSR=3 and U∞=1.5 m/s. Figure 6.11: The torque of one blade of the unducted turbine at TSR=3 and U∞=1.5 m/s. 69 Figure 6.12 shows the comparison of the Cp values for when the turbine is operating within the symmetric duct shown in Figure 6.4. The experimental and RANS simulation results were not available for ducted turbine at TSR values above 3. As this figure shows, BEMDVM gives a good prediction of the maximum performance point which is at TSR=3. However, it under-predicts the power coefficient of the turbine at TSR values less than 2.75. The maximum experimental Cp is increased to 48% for the ducted turbine, up from 29% percent for the unducted turbine. Figures 6.13 and 6.14 illustrate the torque of the ducted turbine and the torque of one blade at TSR=3 and U∞=1.5 m/s respectively. As the Figure 6.13 shows, the torque curve of the ducted turbine has less fluctuation than the torque curve of the unducted turbine. However, the torque curve obtained from BEM-DVM shows six peaks instead of three peaks for experimental and RANS simulation results, which lead to a higher prediction for the power coefficient. These extra peaks might come from the shaft vortex shedding which is not taken into account in BEM-DVM model. As shown in Figure 6.14, the RANS simulation predicts lower torque values for one blade in the downstream region of the turbine shaft. 1.0 BEM-DVM FLUENT EXPERIMENT 0.8 Cp 0.6 0.4 0.2 0.0 2 2.25 2.5 2.75 3 3.25 TSR Figure 6.12: Results of Cp vs TSR for ducted turbine at U∞=1.5 m/s. 70 3.5 Figure 6.13: The torque of the ducted turbine at TSR=3 and U∞=1.5 m/s. Figure 6.14: The torque of one blade of the ducted turbine at TSR=3 and U∞=1.5 m/s. 71 6.4: Experimental testing of a vertical axis turbine at IOT In October 2008, another experimental study was conducted by the researchers from the Naval Architecture and Maritimes laboratory of UBC in the ice tank of the institute of ocean technology (IOT); the author formed part of the team who conducted these experiments. As a part of this experimental study, an asymmetric duct designed by the author of this thesis was tested [4]. This duct, which is shown in Figure 6.15, was tested in two positions relative to the turbine; its initial position is shown with solid line, while the turbine is moved 0.3m forward in its secondary position (dashed line). The trailing edges of the duct tested in the experiments were truncated by 16 cm due to the manufacturing limitations. 1.65 m 1.83 m Figure 6.15: The asymmetric duct and its positions relative to the turbine. The ice tank in IOT is 90 meter long, 12 meter wide and 3 m deep. The carriage which was used as a platform for the turbine has a four-wheel synchronous motor-drive, adjustable frame, 80,000 kg mass, 745 kW of power, with a speed range of .0002 m/s to 4.0 m/s [69]. The Tank and the carriage are shown in Figure 6.16. The two main beams of the carriage which can be raised or lowered via a hydraulic motor are used as a platform for mounting the experimental setup as shown in Figure 6.17. 72 Figure 6.16: The tank and the carriage in IOT. For attaching the turbine and the duct to the two main beams of the carriage, the frame shown in Figure 6.18 was designed by one of engineers involved in these experiments. Figure 6.17: Main beams of carriage. Figure 6.18: Frame designed for experiments [4]. In these experiments, a new set of blades with the same profile given in Figure 6.7 and a new shaft with the same dimensions given in Table 6.1 were used. The torque and the phase angle were obtained using a 1000 N.m torque meter and an encoder. During the experiments, the 73 top end of the blades was 13.5 inches below the free surface level to reduce the free surface effects. Pictures of the turbine placed in two positions relative to the duct during the experiments are shown in Figures 6.19 and 6.20. Figure 6.19: The turbine in its initial position relative to the duct. Figure 6.20: The turbine in its secondary position relative to the duct. 74 6.5: Comparison of the numerical results with the IOT experimental results Figures 6.21 and 6.22 shows the comparison of Cp vs. TSR results for the turbine at its initial position at free stream velocities of U∞=1.5 m/s and U∞=2 m/s. The BEM-DVM results are in good agreement with the experimental results. BEM-DVM also predicts the maximum performance point at TSR=3, which matches well with the experimental prediction. 1.0 BEM-DVM Experiment 0.8 Cp 0.6 0.4 0.2 0.0 2.25 2.5 2.75 3 3.25 TSR Figure 6.21: Cp vs TSR comparison for the turbine in its initial position at U∞=1.5 m/s. 1 BEM-DVM Experiment 0.8 Cp 0.6 0.4 0.2 0 2.25 2.5 2.75 3 3.25 TSR Figure 6.22: Cp vs TSR comparison for the turbine in its initial position at U∞=2 m/s. 75 The comparison of the torque variations over one revolution are shown in Figures 6.23 and 6.24 at TSR=3 and free stream velocities of 1.5 m/s and 2m/s. Similar to the results for the symmetric duct, the torque curves from BEM-DVM have extra peaks which lead to a higher prediction for the power coefficient. As mentioned earlier, these extra peaks might come from the shaft vortex shedding which is not taken into account in BEM-DVM model. Figure 6.23: Torque vs. θ for the turbine in its initial position at TSR=3, U∞=1.5 m/s. Figure 6.24: Torque vs. θ for the turbine in its initial position at TSR=3 and U∞=2 m/s. 76 The comparison of the results for the turbine operating in the secondary position is given below. Figures 6.25 and 6.26 show the comparison of Cp vs. TSR values at free stream velocities of U∞=1.5 m/s and U∞=2 m/s. There is a good agreement between numerical and experimental results at TSR values above 2.75. However, larger discrepancy between the results is observed at TSR values less than 2.75. 1.0 BEM-DVM Experiment 0.8 Cp 0.6 0.4 0.2 0.0 2.25 2.5 2.75 3 3.25 TSR Figure 6.25: Cp vs TSR for the turbine in its secondary position at U∞=1.5 m/s. 1 BEM-DVM Experiment 0.8 Cp 0.6 0.4 0.2 0 2.25 2.5 2.75 3 3.25 TSR Figure 6.26: Cp vs TSR for the turbine in its secondary position at U∞=2 m/s. 77 The comparison of the torque curves for the turbine in its secondary position is shown in Figures 6.27 and 6.28, at TSR=3 and free stream velocities of U∞=1.5 m/s and U∞=2 m/s. Although in these cases only three peaks are obtained from the numerical calculations, the phase of the peaks is about 45 degrees different than the experimental torque peaks. Figure 6.27: Torque variations for turbine in its secondary position at TSR=3,U∞=1.5 m/s Figure 6.28: Torque variations for turbine in its secondary position at TSR=3,U∞=2 m/s 78 6.6: A discussion on the sources of discrepancies There are some discrepancies between the numerical results and experimental results presented in the last sections. The sources of these discrepancies are discussed briefly below: • The numerical results are obtained from two-dimensional models, while the experiments are three-dimensional. • The viscosity of the fluid is not taken into account in numerical models. • The lack of dynamic stall modeling in BEM-DVM model is possibly the most important source of discrepancy between the numerical and experimental results, especially for TSR values below 2.5 where the blade stall is very likely. • The interactions between the arms, the shaft and the blades have not been taken into account accurately, as the effects of the arms and shaft are considered by adding their experimental parasite drag to the numerical results. • The tank wall effects are brought into account through the blockage ratio concept. But as the three-dimensional blockage ratios are translated to two-dimensional blockage ratios in numerical calculations, there is some uncertainty about accuracy of this procedure. • The blades and ducts in numerical models have sharp trailing edges, while the trailing edges of the blades and ducts used in experiments are truncated due to manufacturing limitations. • The free surface of the tank, which might affected the experimental results, is not modeled in the numerical model. • There were some instrumentation issues in the experiments which may have affected the experimental results. The variation of the angular velocity at different phase angles was one of the main problems. This problem resulted from the delay in the response time of the motor control; there were also some issues with the encoder at the correct reading of the phase angle. The friction losses in the gearbox and the motor and bearings may have affected the results. An extensive discussion on the experimental errors is given by Rawlings [3]. • Some inevitable experimental issues such as the wake created by the struts of the frame during the experiments may also have affected the results. 79 6.7: A parametric study on the performance of a VAT Having a good understanding of the impact of different parameters on the performance of a vertical axis turbine is important in designing an efficient turbine. In this section, the results of a parametric study are presented. 6.7.1: Effects of TSR on the torque and the angle of attack (AOA) of a blade Figure 6.29 shows the AOA variations of a blade for the unducted turbine over a revolution at two different TSR values and U∞=1.5 m/s. The AOA of the blade at TSR=3.0 exceeds the static stall angle of attack (which is 15 degree for NACA 634-021, as shown in Figure 3.20) at some phase angles θ, while at TSR=3.5 the AOA of the blade remains below the static stall angle of attack at all phase angles θ. Also, the maximum AOA in the downstream region of the turbine shaft is always lower than the maximum AOA in the upstream region. The higher the TSR becomes, the smaller the maximum AOA in the downstream region of the turbine shaft becomes relative to the maximum AOA at the upstream region. In addition Figure 6.29 shows the torque variations of a blade of unducted turbine over one revolution. A lower maximum torque is produced in the downstream region of the turbine shaft relative to the upstream region at both TSR values 3 and 3.5. Also one can observe that the blade in the downstream region of the turbine shaft becomes less productive with increasing TSR. 45 60 AOA at TSR=3.0 AOA at TSR=3.5 Torque at TSR=3.0 Torque at TSR=3.5 40 AOA (deg) 20 15 0 -20 0 -40 -15 0 45 90 135 180 225 270 315 -60 360 Theta Figure 6.29: Variations of AOA and torque for a blade of turbine at U∞=1.5 m/s. 80 Torque (N.m) 30 6.7.2: Ducting effects on the torque fluctuations of a turbine As it is shown in Figure 6.30, the numerical torque curves for the unducted turbine fluctuate severely at all TSR values. However as it is shown in Figure 6.31, the numerical torque curves for the ducted turbine using the symmetric duct have lower fluctuations especially at higher TSR values. TSR=2.5 80 TSR=2.75 TSR=3.0 TSR=3.25 TSR=3.5 Torque(N.m) 60 40 20 0 0 45 90 135 180 225 270 315 360 Theta Figure 6.30: Torque vs. theta (θ) for the unducted turbine at U∞=1.5 m/s. TSR=2.5 120 TSR=2.75 TSR=3.0 100 TSR=3.25 TSR=3.5 Torque(N.m) 80 60 40 20 0 -20 0 45 90 135 180 225 270 315 Theta Figure 6.31: Torque vs. theta (θ) for the ducted turbine at U∞=1.5 m/s. 81 360 It is more convenient to define a torque fluctuation coefficient, CTF, which is defined as follows from values of the torque curve: CTF = Tmax − Tmin Tavg where Tmax , Tmin and Tavg are the maximum, minimum and average torque of the curve respectively. The CTF results at different TSR values for both unducted and ducted vertical axis turbine are shown in Table 6.2. The CTF values for ducted turbine decreases with increasing TSR value. They show a similar trend to experimental torque fluctuation coefficients [3]. Possibly, lower torque fluctuations for the ducted turbine are the result of flow acceleration within the duct which makes the downstream region of the turbine more productive. Table 6.2: CTF values for an unducted turbine and a ducted turbine at U∞=1.5 m/s. TSR CTF of Unducted Turbine CTF of Ducted Turbine 2.5 1.05 1.87 2.75 0.58 0.80 3 0.61 0.45 3.25 0.86 0.40 3.5 1.21 0.33 6.7.3: Tank wall effects on the power coefficient of a turbine A full-scale vertical axis turbine is expected to operate in the ocean, while normally the experiments of a turbine are conducted in a towing tank. It is important to quantify the wall effects on the turbine tested in a towing tank. Figure 6.32 shows the effects of different tank sizes over the performance of unducted turbine. One can observe that the performance of the unducted turbine drops by as much as 5% when the turbine is operating in a free stream condition, compared to the turbine operating in a tank 9.5 meters in width. 82 Also the Figure 6.33 shows the effect of different tank sizes over the performance of the ducted turbine using the symmetric duct shown in Figure 6.4. Similar to the trend for unducted turbine, the performance of ducted turbine drops in free stream condition. However, the drop in the performance is more pronounced in this case, as the ducted turbine in the tank has higher blockage ratio than the unducted turbine in the tank. 0.5 Free Stream Tank Width : 12.5 m Tank Width : 9.5 m 0.4 Cp 0.3 0.2 0.1 0 2.5 2.75 3 3.25 3.5 TSR Figure 6.32: Cp vs. TSR for the unducted turbine at different tank size at U∞=1.5 m/s. Free Stream Tank Width : 12.5 m Tank Width : 9.5 m 0.8 Cp 0.6 0.4 0.2 0 2.5 2.75 3 3.25 3.5 TSR Figure 6.33: Cp vs. TSR for ducted turbine at different tank size at U∞=1.5 m/s. 83 6.8: Results of duct optimization The theory and the methodology of duct optimization for a vertical axis hydro current turbine were presented in Chapter 5. The results of duct optimization are presented in this section. These results are obtained for TSR=3.0 and U∞=1.5 m/s. As mentioned in Chapter 5, –Cp calculated based on the frontal area of turbine is used as the objective function. The width and length of the duct are specified in each case (i.e. they are not design variables). The constraints used in this optimization study are only the geometrical constraints, and no other constraints such as environmental or cost constraints are considered. For calculation of the power coefficient for the ducted turbine, 180 panels are placed on each side of the duct which are concentrated more at the leading and trailing edges of the duct. 6.8.1: Symmetric duct optimization results A schematic of a symmetric duct is shown in Figure 6.34. For optimization of this duct, it is assumed that the trailing edge angle θ TE =30 degrees, and the length of the duct, LD=5HD, where HD=(WD/2-D/2). Eleven control points are used to represent the curved part of the duct as shown in Figure 6.34. The curved part of the duct is obtained by fitting a spline of degree 10 to these control points. HD LD=5HD Y WD X 3 4 5 6 7 8 9 2 10 1 θ TE Figure 6.34: Symmetric duct and its control points. 84 11 Since the duct is symmetric, only the following variables, bounds and constraints are used for the optimization of this duct: Table 6.3: The design variables and their bounds in symmetric duct optimization. Design Variables Bounds Y i (i = 3 : 6) − WD / 2 < Y i < −0.55D where Y i is the Y coordinate of the ith control point and D is the diameter of the turbine. The X coordinates of the control points are obtained from the following relation: X i (i = 1 : 11) = (i − 1) * dLD where dLD = 1 LD 10 Also, the constraints for the optimization are as follows: Y3 >Y2 Y4 >Y3 Y5 >Y4 Y6 >Y5 θ 2 < θ1 Where tan θ 2 = Y 3 −Y 2 0.1 ⋅ LD θ3 < θ2 θ4 < θ3 Where tan θ 4 = Y 5 −Y 4 0.1 ⋅ LD θ5 < θ4 Where tan θ 3 = Y 4 −Y3 0.1 ⋅ LD Where tan θ 5 = Y 6 −Y5 0.1 ⋅ LD The Y coordinate of the rest of control points are calculated as follows: Y 1 = Y 11 = −WD / 2 Y 2 = Y 10 = 0.1 ⋅ LD ⋅ tan(θ TE ) − WD / 2 Y3 =Y9 Y4 =Y8 Y5 =Y7 The initial profile for optimization is selected to be a parabolic curve which passes through points (-LD/2,-WD/2), (0,-0.55D) and its slope at leading edge of the duct is tan(30◦). 85 The results of duct optimization are shown in Figure 6.35. The Cp value of the turbine using the optimized symmetric ducts approaches 62% by increasing the duct width. As this figure shows, the predicted Cp value for the turbine using the optimized duct with width WD=1.8D is 59.6%, up from 56.9% for the turbine using the tested symmetric duct which has the same width (without adding the parasite drag of the arms). The comparison of the optimized duct and tested symmetric duct at WD=1.8D is given in Figure 6.38. The amount of improvement in power coefficient of the turbine using the optimized duct with width WD=1.8D is not high relative to the experimentally tested duct. This is because the tested duct which is similar tested in NRC trials, was highly optimized, and possibly that the gradient-based optimization methods only find local optimum. One needs an initial condition close to the global optima to find it which is difficult to guess what that initial condition should be. The optimized ducts at different duct widths are shown in Figures 6.36 to 6.39. One should keep in mind that the duct which gives the maximum boost in Cp value may not necessarily be the best duct. The appropriate duct for a turbine depends on the duct cost, drag on the duct and the site requirements. Therefore it is up to the designer to select the most suitable duct, considering all factors which play a role in the overall design. 0.8 Optimized Duct Tested Duct 0.60 0.6 Cp 0.55 0.61 0.62 0.57 0.45 0.4 0.2 1.2D 1.4D 1.6D 1.8D 2D 2.2D WD Figure 6.35: Results of symmetric duct optimization. 86 2.4D Figure 6.36: Optimized symmetric duct for WD=1.4D. Figure 6.37: Optimized symmetric duct for WD=1.6D. 87 Figure 6.38: Optimized symmetric duct for WD=1.8D. Figure 6.39: Optimized symmetric duct for WD=2.0D. 88 6.8.2: Asymmetric duct optimization results Similar to the symmetric duct optimization, it is assumed that the trailing edge angle θ TE =30◦ and the length of the duct, LD=5HD, where HD=(WD/2-D/2) as shown below: HD LD=5HD Y WD X X LE 2 3 4 5 6 7 1 L.E. θTE 8 T .E. Figure 6.40: Asymmetric duct and the control points. Two curves are used for the creation of an asymmetric duct, as shown in Figure 6.40. Each curve is created by fitting a spline of degree 7 to 8 control points. The design variables and their bounds for the asymmetric duct optimization are given below. Table 6.4: The design variables and their bounds at asymmetric duct optimization. Design Variables XLE: distance from leading edge to the shaft of the turbine Bounds − 3D < X LE < 3D − W D / 2 + 0.1H D < YU2 < −0.55 D YUi (i = 2 : 7) : The Y Coordinate of the − W D / 2 < YUi (i = 3,6) < −0.55 D upper control point i − W D / 2 − H D + 0.2 ⋅ LD ⋅ tan(θ T . E . ) < YU7 < −0.55 D YLi (i = 2 : 7) : The Y Coordinate of the − WD / 2 − H D < YL2 < −WD / 2 − 0.1H D lower control point i − WD / 2 − H D < YLi (i = 3,6) < −WD / 2 89 Due to the way the design variables and their bounds are selected, only the following constraint is used: YL7 < YU7 The X coordinates of the control points for this duct is obtained from the following relation: X U1 = X L1 = X LE X Ui (i = 3 : 4) = X Li (i = 3 : 4) = X Ui −1 + (i − 2) * dLD X Ui (i = 5 : 7) = X Li (i = 5 : 7) = X Ui −1 + 3 * dLD Where dLD = 1 LD 15 X T . E . = X L. E . + L D As shown in Figure 6.40, the control points are more concentrated around the leading edge of the duct. In order to obtain a rounded leading edge, the X coordinate of the first and second control points are selected to be the same for both the upper and the lower splines (see the control polygon property of the Spline in Appendix B). The Y coordinate of the trailing edge is calculated from the following relation: YTE = YU7 − 0.2 ⋅ LD ⋅ tan(θ T .E . ) The initial condition for optimization of the asymmetric duct is as follows: X i = [− LD / 5,−WD / 2 + 0.2 H D ,−WD / 2 + 0.4 H D ,−WD / 2 + 0.5H D ,−WD / 2 + 0.6 H D , − WD / 2 + 0.4 H D ,−WD / 2 + 0.2 H D ,−WD / 2 − 0.2 H D ,−WD / 2 − 0.3H D ,−WD / 2 − 0.4 H D , − WD / 2 − 0.2 H D ,−WD / 2 − 0.2 H D ,−WD / 2 − 0.3H D ] The results of the asymmetric duct optimization are given in Figure 6.41. The Cp value of the turbine using the optimized asymmetric ducts approaches 63% by increasing the duct width. Also, the predicted Cp value for the turbine using the optimized asymmetric duct with width WD=1.8D is 61%, up from 57% for the turbine using the tested asymmetric duct which has the same width (without adding the parasite drag of the arms). The shapes of the optimized asymmetric ducts at different widths are also shown in Figures 6.42 to 6.44. 90 0.8 Optimized Duct Tested Duct 0.61 0.6 0.62 0.63 0.57 0.57 Cp 0.50 0.4 0.2 1.2D 1.4D 1.6D 1.8D 2D 2.2D WD Figure 6.41: Results of asymmetric duct optimization. Figure 6.42: Optimized asymmetric duct at WD=1.4D. 91 2.4D Figure 6.43: Optimized asymmetric duct at WD=1.6D. Figure 6.44: Optimized asymmetric duct at WD=1.8D. 92 Chapter Seven Conclusions and Recommendations A new discrete vortex method, referred to as “DVM”, has been developed for performance calculation of unducted vertical axis turbines. The flow in this two-dimensional method is assumed to be potential. The unsteady wake of turbine is modeled by discrete vortices shed from the blades. The interactions between the blades and discrete vortices are taken into account by representing each blade with a vortex filament. The details of this model are explained in Chapter 3. A new numerical model, referred to as “BEM-DVM”, is developed for performance calculation of ducted vertical axis turbines. This model is developed based on DVM model and a boundary element method referred to as “BEM”. This boundary element method is used to calculate the two-dimensional potential flow field around the duct. In this research, BEM-DVM model has also been used to quantify the effects of towing tank walls on experimental results. Understanding the wall effects is important for the design of a full scale turbine operating in the ocean. The details of this model are explained in Chapter 4. 93 The numerical models for unducted and ducted vertical axis turbines were validated using a wide range of experimental and RANS simulation results for an unducted and a ducted turbine. The experimental results were obtained from the experiments conducted in the towing tank at UBC [3] and the towing tank at IOT [4]; the author formed part of the team who conducted these experimental tests. The RANS simulation results were taken from the simulations carried out by Nabavi [1]. A comparison between the numerical and experimental results is provided in Chapter 6. In addition, an optimization procedure was developed to find the duct shape that maximizes the power coefficient of a vertical axis turbine. The optimization method consisted of coupling the MATLAB optimization toolbox with the BEM-DVM model. The MATLAB spline toolbox was used to create the duct shape at each level of the optimization process. Both symmetric and asymmetric ducts were investigated in the optimization study, and new duct shapes were obtained with improved turbine power coefficients. The optimization procedure is explained in Chapter 5, and the optimization results are provided in Chapter 6. 7.1: Conclusions • It is shown that the loads and circulation of an unstalled blade of an unducted vertical axis turbine consist of a quasi-steady component and a pitching component. Refer to Chapter 3 for more details. • The comparison of the numerical and experimental results confirmed that the net output power of an unstalled unducted vertical axis turbine can be calculated using stationary foil characteristics through quasi-steady modeling of the vertical axis turbine. The comparison of results is provided in Chapter 6. • It is shown that lower torque is always produced in the downstream region of the turbine shaft relative to the upstream region. Also, the ratio of the maximum torque in the downstream region of the turbine shaft relative to the maximum torque in the upstream region decreases by increasing TSR value. Refer to section 6.7.1 for details. • The numerical results are in good agreement with the experimental data for TSR values above 2.5. The numerical results do not match well with the experimental results for TSR values below 2.5, mainly because the dynamic stall phenomenon is not taken into account. Refer to Chapter 6 for more details. 94 • The results confirmed that the effects of ducting on a vertical axis turbine can be taken into account by superposing the duct interference on the turbine performance by assuming a potential flow in the duct. The details of the model for ducted turbines are provided in Chapter 4, and the results are given in Chapter 6. • Ducting of a vertical axis turbine can significantly increase the power output of a turbine. The experimental results showed 70% increase in power coefficient of a turbine using a duct with width 1.8D. Refer to sections 6.3 and 6.5 for details. • Numerical and experimental results proved that ducting decreases the torque fluctuations of the vertical axis turbine significantly, especially at high tip speed ratios. The numerical results from the symmetric duct (see Figure 6.4) showed 26%, 53% and 73% reductions in torque fluctuation coefficient (CTF) at TSR values of 3, 3.25 and 3.5 respectively. Refer to section 6.7.2 for the table of results. • The results for the effects of towing tank walls showed that the power coefficient of a vertical axis turbine increases when it operates in a towing tank. This is especially true for a ducted turbine which has a higher blockage ratio than the unducted turbine in the towing tank. The numerical results showed that the amount of increase in power coefficient grows when the tank width decreases. This implies that a turbine has the lowest power coefficient when it operates in a place where there is no wall in its vicinity (like in the ocean). Refer to section 6.7.3 for the results at different tank widths. • Ducts developed from the optimization study improved the power coefficient of the vertical axis turbine relative to the experimentally tested ducts. The optimized symmetric duct with width 1.8D showed a 2.7% increase in power coefficient relative to the tested symmetric duct. The optimized asymmetric duct with width 1.8D showed a 4% increase in power coefficient relative to the tested asymmetric duct. Refer to section 6.8 for more details. • The optimization results showed that the power coefficient of optimized ducts approaches 62% for the symmetric ducts and 63% for the asymmetric ducts, as the duct width increases. This implies that there is an upper limit for the amount of increase in power coefficient due to ducting of the turbine. Refer to section 6.8 for more details. 95 7.2: Original contributions • A new numerical model was developed for the performance calculation of unducted vertical axis turbines. In this model, new formulations are developed for the loads on an unsteady blade of VAT. A new approach is also developed for calculation of the nascent vortex position relative to the blade. In addition, a new approach is developed to take into account the effects of connection point between the blade and the supporting arm. • A new numerical model was developed to predict the performance of ducted vertical axis turbines. This model is fast, compared to other means of numerical analysis, and it is able to quantify the effects of ducting on power coefficient of a turbine within the engineering range of accuracy. • Based on the method developed for ducted turbines, a new approach was developed to calculate the effects of towing tank walls on the experimental results. • Using a wide range of RANS simulation results and experimental results from testing in towing tanks at UBC and IOT, a good correlation is obtained between the experimental and numerical results. • An optimization procedure was developed to find the duct shape that maximizes the output power of a vertical axis turbine. This procedure is fast and has the ability to take into account different constraints on the length and width of the duct. 7.3: Recommendations for the future work The recommendations for future work are as follows: 7.3.1: Foil shape study Most of the studies conducted on foil selection for vertical axis hydro current turbines are limited to available NACA airfoils which are designed for performing in air. Designing foils specifically for vertical axis turbines operating in the water is a very important step towards the development of an efficient and reliable turbine. It is also important to note that the boundary condition on a moving foil is not exactly the same as for a stationary foil in a uniform flow. The geometry of a foil for vertical axis turbines needs to be adjusted to obtain the same amount of lift and drag as that of a stationary foil. 96 7.3.2: Inclusion of a dynamic stall model into the numerical models The numerical results were not in good agreement with experimental data at TSR values lower than 2.5. This problem may be solved by including a dynamic stall model to adjust the foil characteristics at angles higher than the static stall angle of attack. 7.3.3: Inclusion of a boundary layer method for duct calculations The boundary layer of the duct is not taken into consideration in the method developed in this research. However, including the boundary layer is necessary for the ducts that have large separation zones. 7.3.4: Inclusion of three-dimensional effects As a blade’s span is finite, there are always some losses from a blade in the form of vortices shed from blade’s tip. These losses might be significant depending on the aspect ratio of the blade and the way arms and the blade are connected. A better correlation between numerical and experimental results might be obtained by including these losses through extending the numerical models to three dimensions. 7.3.5: Using a high-order boundary element method The boundary element method used in this research showed a reasonable mass continuity for the flow in the duct. However, better mass flow continuity might be obtained by adopting a higher-order boundary element method. 97 References 1: Nabavi, Y., “Numerical Study of the Duct Shape Effect on the Performance of a Ducted Vertical Axis Tidal Turbine.” MASc thesis, University of British Columbia, Vancouver, B.C., Canada, 2007 2: Milne-Thompson, L., "Theoretical Hydrodynamics", MacMillan & Co., 1968 3: Rowlings, G.W., “Parametric characterization of an experimental vertical axis hydro turbine” MASc thesis, University of British Columbia, Vancouver, B.C., Canada, 2008 4: Nabavi, Y., Klaptocz, V., Rowlings, B., Alidadi, M., “The Experimental Study in IOT Tank”, Internal Report, Naval Architecture and Maritime Laboratory, Department of Mechanical Engineering, UBC, 2008 5: United Nations Environment Programme , “Global Trends in Sustainable Energy Investment”, 2007 6: Twidell, A. and Weir, T. “Renewable Energy Sources”, Taylor & Francis, 2nd Edition, 2006 7: Paraschivoiu, I., “Wind Turbine Design : With Emphasis on Darrieus Concept”, Polytechnic International Press, 2002 8: South, P., Rangi, R.S., “Preliminary Tests of a High Speed Vertical Axis Wind Mill Model,”, National Research Council of Canada, Report No, LTR-LA-74, 1971 9: Quite revolution Turbine, http://www.quietrevolution.co.uk 10: Betz, A., “Das Maximum der theoretisch möglichen Ausnützung des Windes durch Windmotoren”, Zeitscrift für das gesamte Turbinenwesen, Heft 26, Sept.26, 1920 11: Newman, B.G., “Actuator-disk theory for vertical-axis wind turbines”, J. Wind Eng. And Industrial Aerodynamics, 15/3, 1983, p.347-355 98 12: Kirke, Brian, “Developments in ducted water current turbines”, 2005 , http://www.cyberiad.net/library/pdf/bk_tidal_paper25apr06.pdf 13: Riegler, G., “Principles of energy extraction from a free stream by means of wind turbines”, Wind Engineering 7/2, 1983, p.115-126 14: Templin , R. J. “ Aerodynamic Performance Theory for the NRC Vertical-Axis Wind Turbine” , National Research Council of Canada Report LTR-LA-160, June (1974) 15: Wilson, R., Lissaman, P., “Applied Aerodynamics of Wind Powered Machines”, Technical Report NSF-RA-N-74-113, Oregon State University 16: Strickland , J. H. , “The Darrieus Turbine : A Performance Prediction Model Using Multiple Streamtubes”, Advanced Energy Projects Department, Sandia Laboratory , SAND 75-0431, Oct. 1975 17: Paraschivoiu, I., “Double-Multiple Streamtube Model for Darrieus Wind Turbines”, Second DOE/NASA Wind Turbines Dynamics Workshop, NASA CP-2186, Cleveland, Ohio, February 1981, pp 19-25 18: Islam, M., “Analysis of Fixed-Pitch Straight-Bladed VAWT with Asymmetric Airfoils”, PhD Thesis, University of Windsor, 2008 19: Larsen, H., “Summary of the Vortex Theory of the Cyclogiro”, Proceeding of the 2nd national conference on Wind Engineering Research, June 1975 20: Holme, O., “A Contribution to the Aerodynamic Theory of the Vertical Axis Wind Turbine”, Proceeding of the International Symposium on Wind Energy Systems, Cambridge, England, Paper C4, September 1976 21: Wilson , R.E., “ Vortex Sheet Analysis of the Giromill”, Journal of Fluids Engineering , Vol. 100 , September 1978 , pp. 340-342 99 22: Wilson , R.E., and Walker , S.N. , “ Fixed Wake Analysis of the Darrieus Rotor”, Journal of Fluids Engineering , Vol. 105 , December 1983 , pp. 389-393 23: Fanucci, J., Walters, R., “Innovative Wind Machines: The Theoretical Performance of a Vertical Axis Wind Turbine”, Proceedings of the VAWT Technology Workshop, Sandia Lab. Report SAND 76-5586, pp. III-61-93 24: Strickland, J. H., Webster, B. T.,Nguyen, T., “A Vortex Model of the Darrieus Turbine : An Analytical and Experimental Study” ASME Journal of Fluids Engineering , Vo1. 101, December 1979 25 : Li, Y., “Development of a Procedure for Predicting the Power Output from a Tidal Current Turbine Farm”, PhD thesis, University of British Columbia, Vancouver, B.C., Canada, 2008 26: Rajagopalan, R. G, “Viscous Flow Field Analysis of a Vertical Axis Wind Turbine”, Proceedings of the Intersociety Energy Conversion Engineering Conference, 1986, p 12421246 27: Ishimatsu, K., Kage, K., Okubayashi, T., “Numerical simulation for flow fields of Darrieus turbine”, Nippon Kikai Gakkai Ronbunshu, B Hen/Transactions of the Japan Society of Mechanical Engineers, Part B, v 61, n 587, Jul, 1995, p 2543-2548 28: Allet, A. (Ecole Polytechnique de Montreal), Brahimi, M.T., Paraschivoiu, I., “On the aerodynamic modeling of a VAWT”, Wind Engineering, v 21, n 6, 1997, p 351-365 29: Guerri, O., Sakout, A., Bouhadef, K., “Simulations of the fluid flow around a rotating vertical axis wind turbine”, Source: Wind Engineering, v 31, n 3, May, 2007, p 149-163 30: Jiang, Z. (College of Mechanical and Energy Engineering, Zhejiang University), Doi, Y., Zhang, S., “Numerical investigation on the flow and power of small-sized multi-bladed straight Darrieus wind turbine”, Journal of Zhejiang University: Science A, v 8, n 9, August, 2007, p 1414-1421 100 31: Wang, T., Coton, F.N., "Numerical simulation of wind tunnel wall effects on wind turbine flows", Wind Engineering , 2000; 3:135-148 32: Lawn, C.J. (Department of Engineering, University of London), “Optimization of the power output from ducted turbines”, Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, v 217, n 1, 2003, p 107-118 33: Werle, M.J. (FloDesign Inc.), Presz Jr., W.M., “Ducted wind/water turbines and propellers revisited”, Source: Journal of Propulsion and Power, v 24, n 5, September/October, 2008, p 1146-1150 34 : Gilbert, B.L. and Foreman, K.M., “Experiments with a diffuser-augmented model wind Turbine”, J. Energy Resources Technology, Trans ASME, vol.105, March 1983, p.46-53 35: Davis, Barry V., “Water Turbine Model Trials”, Nova Energy Limited for NRC Hydraulics Laboratory NEL-002, March 1980 36: Davis, B.V., D.H. Swan, and K.A. Jeffers, “Ultra Low Head Hydroelectric Power Generation Using Ducted Vertical Axis Water Turbines”, Nova Energy Limited for NRC Hydraulics Laboratory NEL-021, March 1981 37: Davis, B.V., D.H. Swan, and K.A. Jeffers, “Ultra Low Head Hydroelectric Power Generation Using Ducted Vertical Axis Water Turbines”, Nova Energy Limited for NRC Hydraulics Laboratory NEL-022, October 1983 38: Davis, B.V., J.R. Farrell, D.H. Swan, and K.A. Jeffers, “Research and Development of a 50kW to 100kW Vertical Axis Hydro Turbine for a Restricted Flow Installation”, Nova Energy Limited for NRC Hydraulics Laboratory NEL-038, March 1984 39: Davis, B.V., D.H. Swan, and K.A. Jeffers, “The Ducted Vertical Axis Hydro Turbine for Large Scale Tidal Energy Applications”, Nova Energy Limited for H. A. Simmons NEL-070, March 1984. 101 40: Davis, B.V. and D.H. Swan, “Commissioning and Testing of a 100kW Vertical Axis Hydraulic Turbine”, Nova Energy Limited for NRC Hydraulics Laboratory NEL-081, December 1985. 41: Gorlov, A. M., "Tidal Energy", Encyclopedia of Ocean Sciences, Academic Press, London, pp. 2955-2960, 2001 42: Shiono, M., Suzuki, K., Kiho, S., “An Experimental Study of the Characteristics of a Darrieus Turbine for Tidal Power Generation”, Electrical Engineering in Japan, Vol. 132, No. 3, 2000 43: “Cycloidal Tidal Power Generation – Phase 2.” QinetiQ Ltd. for the Department of Trade and Industry Contract T/06/00229/00/REP/2, 2004. 44: Setoguchi, T., Shiomi, N., Kaneko, K., “Development of two-way diffuser for fluid energy conversion system”, Renewable Energy 29 (2004) 1757–1771 45: Coiro, D.P., A. De Marco, F. Nicolosi, S. Melone, F. Montella. “Dynamic Behaviour of the Patented Kobold Tidal Current Turbine: Numerical and Experimental Aspects.” Acta Polytechnica Vol. 45 No. 3, pp. 77-84, 2005. 46: Guido, C., S. Francesco, L. Greco, A. Moroso, H. Eriksson. “An Experimental Investigation and a Theoretical and Computational Methodology to Study an Innovative Technology for Marine Current Exploitation: the Kobold Turbine.” Bolletino della Communita Scientifica in Australasia, December 2006 47: Ponta, F. and P. Jacovkis, “Marine-Current Power Generation by Diffuser-Augmented Floating Hydro-Turbines”, Renewable Energy, Elsevier, 2007 48: Noll, R., Ham, N. ,"Effects of dynamic stall on SWECS", Journal of Solar Energy Engineering, 104, pp. 96-101, 1982 49: Klimas, P., "Darrieus rotor aerodynamics", Journal of solar energy engineering , 104, pp. 102-105 , 1982 102 50: Paraschivoiu, I. Allet, A. , " Aerodynamic analysis of the Darrieus wind turbines including dynamic stall effects", Journal of Propulsion, 4(5), pp. 472-477, 1988 51: Jiang, D., Coton, F., Galbraith, M.R., "Fixed wake vortex model for vertical axis wind turbines including unsteady aerodynamics", Wind Engineering, 15(6), pp. 348-360, 1991 52: Katz, J., Plotkin, A., "Low Speed Aerodynamics", Cambridge University Press, Second Edition 53: Fung, Y., " An introduction to the theory of Aeroelastisity", John wiley & Sons, NY., 1955 54: Strickland, J., Webster, B., Nguyen, T., "A vortex model of the Darrieus turbine: an analytical and experimental study", Technical report SAND81-7017, Sandia National Laboratories, 1981 55: Wilson, R.E., Lissaman, P.B.S., James, M., McKie, W.R., "Aerodynamic loads on a Darrieus rotor blade", Journal of Fluids Engineering, Vol. 105, pp. 53-58, March 1983 56: McCroskey, W.J. “The Phenomenon of Dynamic Stall.”, March 1981, NASA Technical Reports Server, http://ntrs.nasa.gov 57: Gormont, R., "A mathematical model of Unsteady Aerodynamics and Radial Flow for Applications to Helicopters", Technical Report DAAJO2-71-C-0045, 1973 58: Johnson, W., "The Response and Airloading of Helicopter Rotor Blades due to Dynamic Stall", Technical Report TR 130-1, MIT Aeroelastic and Structure Research Lab, May 1970 59: Leishman, J., Beddoes, T., "A semi-empirical model for Dynamic Stall", Journal of the American Helicopter Society, 34(3), pp. 3-17, 1989 60: Ciffone, D.L., Orloff, K.L., "Far-field Wake-Vortex Characteristics of Wings", Journal of Aircraft, Vol. 12, No. 5, pp. 464-470, 1975 61: Fluent Inc. “User’s Guide,” Fluent 6.2.2, Ch 12, 2006 103 62: Ashby, D. L., Sandlin, D. R., "Application of a low order panel method to complex three dimensional internal flow problems," Ames Research Center, Contract NCC2-226, September 1986 63: Holt DR, Hunt B., "The use of panel methods for the evaluation of subsonic wall interference", AGARD-CP-335, 1982 64: Nocedal, J., Wright, S.J., “Numerical Optimization”, Springer Publication, 1999 65: Bertsekas, D.P., “Nonlinear Programming”, Athena Scientific Publication, Second Edition, 1999 66: MATLAB 6.2 documentation; optimization toolbox http://www.mathworks.com/access/helpdesk/help/toolbox/optim/ 67: Anderson, F. , "Bezier and B-spline Technology", www.cs.umu.se/education/examina/Rapporter/461.pdf 68: Farin, G., “Curves and Surfaces for CAGD”, Morgan Kaufmann Publication, 5th edition, 2001 69: National Research Council of Canada http://iot-ito.nrc-cnrc.gc.ca/facilities/it_e.html 104 Appendix A: Gradient-based optimization methods Several gradient-based algorithms are developed for different kinds of optimization problems, as an efficient and accurate solution of the optimization problem depends not only on the size of the problem in terms of the number of constraints and design variables but also on characteristics of the objective function and constraints. When both the objective function and the constraints are linear functions of the design variables, the problem is known as a Linear Programming (LP) problem. Quadratic Programming (QP) concerns the minimization or maximization of a quadratic objective function that is linearly constrained. For both the LP and QP problems, reliable solution procedures are readily available. More difficult to solve is the Nonlinear Programming (NP) problem in which the objective function and constraints can be nonlinear functions of the design variables. A solution of the NP problem generally requires an iterative procedure to establish a direction of search at each major iteration. This is usually achieved by the solution of an LP, a QP, or an unconstrained subproblem. The general idea for constrained optimization is to transform the problem into an easier subproblem that can then be solved and used as the basis of an iterative process. A characteristic of a large class of early methods is the translation of the constrained problem to a basic unconstrained problem by using a penalty function for constraints that are near or beyond the constraint boundary. In this way the constrained problem is solved using a sequence of parameterized unconstrained optimizations, which in the limit (of the sequence) converge to the constrained problem. These methods are now considered relatively inefficient and have been replaced by methods that have focused on the solution of the KuhnTucker (KT) equations. The KT equations are necessary conditions for optimality for a constrained optimization problem. If the problem is a so-called convex programming problem, that is, f ( x ) and g i ( x), i = 1,..., m are convex functions, then the KT equations are both necessary and sufficient for a global solution point. The Kuhn-Tucker equations is stated in equation (A-1). The solution of the KT equations forms the basis to many nonlinear programming algorithms. These algorithms attempt to compute the Lagrange multipliers directly. Constrained quasi-Newton methods guarantee superlinear convergence by accumulating second-order information regarding the KT equations using a quasi-Newton updating procedure. These methods are commonly referred 105 to as Sequential Quadratic Programming (SQP) methods, since a QP subproblem is solved at each major iteration. m ∇f x ∗ + ∑ λ∗i ⋅ ∇g i x ∗ = 0 ( ) λ g (x ) = 0 ( ) i =1 ∗ i ∗ (A-1) i = 1,..., me i λ∗i ≥ 0 i = me + 1,..., m Where the x ∗ is the local minimum and λ∗i is the i th Lagrange multipliers. A.1: Sequential quadratic programming (SQP) This method mimics Newton’s method for constrained optimization just as is done for unconstrained optimization. At every major iteration, an approximation is made of the Hessian of the Lagrangian function using a quasi-Newton updating method. This is then used to generate a QP subproblem whose solution is used to form a search direction for a line search procedure. The formulation of a QP subproblem is based on a Lagrangian function which is defined as follows: m L ( x, λ ) = f ( x ) + ∑ λ i g i ( x ) (A-2) i =1 • SQP procedure The SQP procedure consists of three main stages as follows: • 1- Updating the hessian matrix The hessian matrix H is the square matrix of second-order partial derivatives of a function. • 2- Quadratic programming solution A QP subproblem needs to be solved at every major iteration to find a search direction d • 3- Line search In this stage, the step length parameter β k is calculated to form a new iteration as follows: x k +1 = x k + β k d k 106 Appendix B: Introduction to spline B.1: Bezier curves The Bezier curves are made up from basic functions called Bernstein polynomials. In general form, these functions are defined as follows: n n −i Bi , n (u ) = u i (1 − u ) i (B-1) Where i is the index of the control point to be weighted and n is the degree of the curve. The parameter u ∈ (0,1) at equation (B-1) can be viewed as a point trying to interpolate the control points in discrete steps. As it moves along the trajectory, it is drawn towards control point i by an amount determined by Bi,n (u ) . In another words, the Bernstein polynomial represents the influence of the control point with index i . The definition of a generalized Bezier curve of degree n is defined as follows: n C (u ) = ∑ Pi Bi ,n (u ) (B-2) i =0 As the equation B-2 shows, there is no local control property. All control points affect the curve, and this can be a disadvantage when designing smooth shapes with the Bezier curve. A Bezier curve of degree three with four control points and an outlined control polygon is shown in Figure B-1. As in this figure is shown, the tangents to the Bezier curve at the end points are always along the line segments P0 P1 and Pn~ −1 Pn~ . Figure B-1: Bezier curve and control polygon. 107 Another fundamental property of the Bezier curve is that no matter how we manipulate the control points, we can always rely on the curve to lie within the control polygon composed from line segments P0 P1 , P1 P2 ,…, Pn−1 Pn called legs. This property is very beneficial at optimization process, because no point of the resulting curve gets beyond the control points which are used as design variables. By concatenating several Bezier curves we can obtain local control, this is done by sharing the first and last control points with the surrounding curves. Obviously, the first and last segments only share one of these. This will however introduce some discontinuities at the shared control point. A Bezier curve has C 1 continuity to its defined interval, but there is only C 0 continuity at shared vertices. Hence, we do not obtain a smooth transition between segments if we do not create these manually, which can be good or bad, depending on the application. For instance, a simple square can not be created with concatenated C 1 continuous curves unless we introduce a discontinuity at the corners of it, just as a circle requires smooth transitions and greater continuity. B.2: B-spline curves As it is discussed in the last section, the main drawbacks with Bezier curves are: • No real local control • Strict relation between Bezier curve degree and number of control points B-splines can solve these problems. In addition, they can be forced to interpolate any of its control points without repeating it, which is not possible with the Bezier curve. The B-Spline of degree p , for n + 1 control points P0 ,..., Pn is defines as follows: n C (t ) = ∑ Pi Bi , p (t ) (B-3) i =0 Where Bi , p is the basis function and p = m − n − 1 . Also t belongs to a non-decreasing sequence with t i ∈ [0,1] . 108 The basis functions are calculated using the following recursive algorithm which is referred as the Cox de Boor algorithm. Bi ,1 (t ) = 1.0 if t i ≤ t ≤ t i +1 , else 0.0 t i + p +1 − t t − ti Bi , p (t ) = Bi , p −1 (t ) + Bi +1, p −1 (t ) ti+ p − ti t i + p +1 − t i +1 (B-4) The main properties of the B-spline defined by equation B-3 are as follows: • It is p − k times differentiable at a point where k duplicate knot values occur. • The B-spline with no internal knots is a Bezier curve where the knots t p +1 ,..., t m − p −1 are called internal knots. 109
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Duct optimization for a ducted vertical axis hydro...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Duct optimization for a ducted vertical axis hydro current turbine Alidadi, Mahmoud 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Duct optimization for a ducted vertical axis hydro current turbine |
Creator |
Alidadi, Mahmoud |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | Turbines developed for the power production from the energy of tidal and river currents are commonly referred to as hydro current turbines. This research focuses on duct optimization for a ducted vertical axis hydro current turbine. As a first step, a discrete vortex method, referred to as “DVM”, is developed for the performance analysis of unducted vertical axis turbines. It is shown that the loads for an unstalled blade of a vertical axis turbine may be adequately approximated by a quasi-steady component and a pitching component. A new approach is developed for the calculation of nascent vortex position relative to the blade. Besides, a new approach is developed to take into account the effects of connection point between a blade and its supporting arm. Subsequently, a new numerical model, referred to as “BEM-DVM”, is developed for the performance calculation of ducted vertical axis turbines by combining a boundary element method referred to as “BEM” with the DVM model. Comparison of the results for unducted and ducted turbines showed that ducting significantly improves the power coefficient of a vertical axis turbine. In addition, the results showed that ducting decreases the torque fluctuations significantly, especially at high tip speed ratios. Based on the method developed for ducted turbines, a new approach is developed to calculate the effects of towing tank wall on experimental results. The results showed that the power coefficient of a turbine increases when it operates in a tank. This is especially true for a ducted turbine which has a higher blockage ratio than the unducted turbine in the tank. The results also showed that the amount of increase in power coefficient grows when the tank width decreases. An optimization study is conducted to find the duct shape that maximizes the output power of a vertical axis turbine. The optimized ducts results showed a higher power coefficient for the turbine relative to the results of experimentally tested ducts. It is also shown that there is an upper limit to the amount of increase in power coefficient due to ducting. |
Extent | 1462660 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-06-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067284 |
URI | http://hdl.handle.net/2429/9116 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_fall_alidadi_mahmoud.pdf [ 1.39MB ]
- Metadata
- JSON: 24-1.0067284.json
- JSON-LD: 24-1.0067284-ld.json
- RDF/XML (Pretty): 24-1.0067284-rdf.xml
- RDF/JSON: 24-1.0067284-rdf.json
- Turtle: 24-1.0067284-turtle.txt
- N-Triples: 24-1.0067284-rdf-ntriples.txt
- Original Record: 24-1.0067284-source.json
- Full Text
- 24-1.0067284-fulltext.txt
- Citation
- 24-1.0067284.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067284/manifest