The Pennsylvania State University, Study Guides, Projects, Research of Designs and Groups

Download The Pennsylvania State University and more Study Guides, Projects, Research Designs and Groups in PDF only on Docsity! THE PENNSYLVANIA STATE UNIVERSITY SCHREYER HONORS COLLEGE DEPARTMENT OF AEROSPACE ENGINEERING BLADE ELEMENT MOMENTUM THEORY APPLIED TO HORIZONTAL AXIS WIND TURBINES THOMAS R. PURCELL Spring 2011 A thesis submitted in partial fulfillment of the requirements for a baccalaureate degree in Aerospace Engineering with honors in Aerospace Engineering Reviewed and approved* by the following: Dr. Dennis K. McLaughlin Professor of Aerospace Engineering Thesis Supervisor and Honors Advisor Dr. Susan W. Stewart Research Associate Faculty Reader Dr. George A. Lesieutre Department Head of Aerospace Engineering * Signatures are on file in the Schreyer Honors College. i ABSTRACT Research efforts in the field of small wind turbines at Penn State has led to the need of an open source Blade Element Momentum Theory (BEMT) code for preliminary performance analysis. PSUWTA, a MATLAB code, has been developed in hopes of fulfilling this need. When checked against experimental wind turbine data and WT_Perf, the National Wind Technology Center‟s non-open source BEMT code, the results of PSUWTA came to within engineering accuracy of both of these data sources. PSUWTA was then used to analyze Southwest Windpower‟s Whisper 500 turbine, which Penn State currently uses, as well as to analyze the two and three-bladed versions of the Carolus turbine, an in-house design currently under construction. The Whisper 500 was found to have sub-par performance due to poor design while the Carolus turbine‟s performance was found to be a substantial improvement. The source code for PSUWTA, as well as operating instructions, is also included. iv LIST OF FIGURES Figure 2.1: Blade Element Diagram (Hansen 2008) .................................................... 4 Figure 2.2: Aerodynamic Forces Acting on a Blade Element (Moriarty 1993) .......... 5 Figure 3.1: The Wortmann FX 60-126 Airfoil ............................................................ 9 Figure 3.2: Cl vs. Angle of Attack for the Wortmann FX 60-126 ............................... 10 Figure 3.3: Wortmann FX 60-126 Drag Polar ............................................................. 10 Figure 3.4: Expanded Cl and Cd vs. Angle of Attack for the Wortmann FX 60-126 .. 11 Figure 3.5: Cd vs. Angle of Attack for the Wortmann FX 60-126 Airfoil ................... 12 Figure 3.6: The S822 Airfoil ........................................................................................ 13 Figure 3.7: Cl vs. Angle of Attack for the S822 .......................................................... 13 Figure 3.8: S822 Drag Polar ........................................................................................ 14 Figure 3.9: Expanded Cl and Cd vs. Angle of Attack for the S822 ............................. 15 Figure 3.10: Cd vs. Angle of Attack for the S822 ....................................................... 15 Figure 3.11: Whisper Blade Local Pitch Angle vs. Radial Location ........................... 16 Figure 3.12: Carolus Local Pitch Angle vs. Radial Location ...................................... 17 Figure 3.13: Whisper Blade Chord vs. Radial Location .............................................. 18 Figure 3.14: Carolus Turbine Chord vs. Radial Location ............................................ 18 Figure 4.1: Overall Power Output at 250 RPM Comparison ....................................... 23 Figure 4.2: Overall Power Output at 350 RPM Comparison ....................................... 24 Figure 4.3: Whisper Turbine Experimental Power Production ................................... 25 Figure 4.4: Whisper Turbine Experimental RPM and Wind Speed Correlation ......... 26 Figure 4.5: Overall Power Coefficient Comparison .................................................... 27 v Figure 4.6: Experimental Wind Speed and Power Coefficient vs. Tip Speed Ratio ... 28 Figure 4.7: Comparison of Radial Power Distribution Comparison............................ 29 Figure 4.8: Radial Distributions of the Axial Induction Factor Comparison .............. 30 Figure 4.9: Reynolds Number Distribution Condition 1 Comparison ......................... 31 Figure 4.10: Reynolds Number Distribution Condition 2 Comparison ....................... 32 Figure 5.1: Power Output vs. RPM and Wind Speed .................................................. 34 Figure 5.2: Power Coefficient vs. RPM and Wind Speed ........................................... 35 Figure 5.3: Experimentally Measured RPM vs. Wind Speed ...................................... 36 Figure 5.4: Power Coefficient vs. Tip Speed Ratio ..................................................... 37 Figure 5.5: Experimentally Measured Power Coefficient vs. Tip Speed Ratio ........... 38 Figure 5.6: Radial Power Distribution of the Whisper Turbine, Wind Speed = 8 m/s ......................................................................................................................... 39 Figure 5.7: Power Coefficient vs. Radial Location, Wind Speed = 8 m/s ................... 40 Figure 5.8: Axial Induction Factor vs. Radial Location, Wind Speed = 8 m/s ............ 41 Figure 5.9: Angle of Attack vs. Radial Location, Wind Speed = 8 m/s ...................... 42 Figure 5.10: Thrust vs. Wind Speed and RPM ............................................................ 43 Figure 5.11: Thrust Force vs. Radial Location, Wind Speed = 8 m/s.......................... 44 Figure 5.12: Overall Torque vs. Wind Speed .............................................................. 45 Figure 6.1: Power Output vs. RPM and Wind Speed .................................................. 46 Figure 6.2: Power Coefficient vs. Wind Speed and RPM ........................................... 47 Figure 6.3: Power Coefficient vs. Tip Speed Ratio ..................................................... 48 Figure 6.4: Power Output vs. Radial Location, Wind Speed = 8 m/s .......................... 49 Figure 6.5: Power Coefficient vs. Radial Location, Wind Speed = 8 m/s ................... 50 Figure 6.6: Axial Induction Factor vs. Radial Location, Wind Speed = 8 m/s ............ 51 Figure 6.7: Angle of Attack vs. Radial Location, Wind Speed = 8 m/s ...................... 52 vi Figure 6.8: Overall Thrust vs. Wind Speed and RPM ................................................. 54 Figure 6.9: Thrust Force vs. Radial Location, Wind Speed = 8 m/s ........................... 55 Figure 6.10: Overall Torque vs. Wind Speed and RPM .............................................. 56 Figure 7.1: Power vs. Wind Speed and RPM .............................................................. 57 Figure 7.2: Power Coefficient vs. Wind Speed and RPM ........................................... 58 Figure 7.3: Power Coefficient vs. Tip Speed Ratio ..................................................... 59 Figure B.1: Airfoil Data Storage Function .................................................................. 65 Figure B.2: Airfoil Data Function Call ........................................................................ 66 Figure B.3: Airfoil Data Selection Code...................................................................... 66 Figure B.4: WT_Perf File Name Entry ........................................................................ 68 ix the blade element δr = Width of the blade element λ = Tip Speed Ratio µ = Air Viscosity ρ = Air Density σ = Solidity of the annulus containing the blade element ϕ = Local Inflow Angle Ω = Rotational rate x ACKNOWLEDGEMENTS To begin, I first want to thank Dr. McLaughlin for stepping up as my advisor when I transferred to the College of Engineering from DUS and began looking for an honors thesis topic. If it weren‟t for him, I would not have known about this project. His guidance in this project made it possible. I also want to thank Brian Wallace, who also really helped to guide my research and work with me on troubleshooting PSUWTA with me. If it weren‟t for him, I don‟t if I ever would have been able to get his code working properly. He also provided all of the experimental data on the Whisper Turbine from his own thesis, and I believe the information really helps add depth to my research here. I would also like to thank now 2 nd LT Jacob Marsh, USAF, for creating the operating code for PSUWTA that has made organizing all of the data much easier and allowed PSUWTA to expand its capabilities exponentially. I also want to thank my parents for always supporting me in my school work and everything else that I‟m involved with. Knowing that I have that support from home really makes life much easier here. I also want to thank my two brothers, Kevin and Nate, my cousin Bubba, and my girlfriend Alexis for always being there to go hang out, have fun with, and get my mind off of all my work for a while. You guys kept me sane. The same goes for the rest of my extended/crazy friends and family. You‟re all awesome. Thank you. Chapter 1 Introduction Looking towards the future, green energy sources, such as wind, solar, geothermal, and nuclear, are all being looked at to one day help alleviate our dependence on fossil fuels. In President Obama‟s most recent “State of the Union” speech, he called for America to generate 80% of its electricity from clean sources such as these by the year 2035 (Obama 2011). Each of these energy resources will need to be harnessed in order to meet this goal. Wind energy has already shown its worth in countries such as Denmark, where, as of 2009, wind power provides 19% of the country‟s electricity (Walsh 2009). The majority of this power comes from wind farms, large tracts of land used for holding hundreds of massive wind turbines, among other uses. This is what most people tend to think of when they think of wind energy. However, there is another avenue through which the wind‟s energy can be harvested. This lays in the usage of small scale wind turbines in isolated or distributed use for specific purposes. Rather than using a massive wind farm to power a whole town or city, small scale wind turbines can be used, for example, to supplement a home‟s electrical needs and help lower electricity bills. Using wind turbines in this sense seems to have been overshadowed since the easiest way to increase the power generated by wind turbines is to simply increase their overall size. However, for small wind turbines, where increasing the size is not an option, their performance must be maximized in other ways. From an aerodynamics perspective, this can be done by carefully designing the wind turbine blades. This is done by borrowing much theory from helicopters and propellers. In those cases, the goal of the designer is generally to use the minimum amount of power to achieve the maximum amount of thrust. However, for wind turbines, it nearly the opposite, with the goal being to achieve the maximum possible power generation with whatever wind conditions are present. 4 In this figure, r represents the distance from the hub to the blade; δr represents the width of the blade element; Ω, also labeled as ω on occasion, represents the rotational rate of the rotor; U, sometimes labeled as U∞, represents the free-stream velocity as unaffected by the turbine; a represents the axial induction factor; and a’ represents the tangential induction factor. From here each blade element can be examined more closely. Figure 2.2 shows a much more detailed description of an individual blade element and the aerodynamic forces acting on it: Figure 2.1: Blade Element Diagram (Hansen 2008) 5 Along with the previously defined variables, β represents the local pitch angle of the blade with respect to the rotor plane; α represents the angle of attack; W represents the total local relative wind velocity; ϕ represents the inflow angle, which is the angle between the relative wind velocity and the rotor plane; L represents the lift force acting on the blade element; and D represents the drag force acting on the blade element. With the turbine properly divided into blade elements defined in the terms of BEMT, it is now possible to derive the iterative process used to calculate the overall forces acting on the turbine, and therefore the power production of the wind turbine. The process starts with a guess for the axial and tangential inductions factors, typically one and zero respectively, and an assumed rotational rate and wind speed under which the turbine is operating. It is also assumed that the blade geometry, defined as the chord distribution, local pitch angle distribution, and airfoil used, is known. From here equations 2.1, 2.2 and 2.3 can be used to calculate the local velocity of the wind and the angle of attack on the blade element: √ (2.1) (2.2) (2.3) Figure 2.2: Aerodynamic Forces Acting on a Blade Element (Moriarty 1993) W W 6 With the angle of attack known, the coefficients of lift and drag, Cl and Cd, can be found by looking up the airfoil data. If these coefficients vary widely with Reynolds number, the local Reynolds number of the blade element, Re, can be calculated using the local chord length of the element, c, and the local relative velocity, W. Equation 2.4 defines the local Reynolds number: (2.4) Here ρ is the air density and µ is the viscosity of the air. Coefficients of lift and drag can then be looked up based on both Reynolds number and angle of attack if necessary. From here, the lift and drag coefficients can be transformed into force coefficients normal and tangential to the plane of the rotor, Cn and Ct respectively. They are defined here in equations 2.5 and 2.6: (2.5) (2.6) Next, the solidity, σ, of the blade element must be defined before the new induction factors can be calculated. Equation 2.7 shows the definition of this factor: (2.7) In the above equation, N is the number of blades on the turbine. With these factors calculated, new values for the axial and tangential induction factors can be found using equations 2.8 and 2.9: 9 Chapter 3 Wind Turbine Blade Geometries The geometry of the blades of a wind turbine is obviously crucial to the performance of the turbine as a whole. The three primary factors that determine the geometry of a turbine blade are the types of airfoils used, chord length vs. radial position, and local pitch angle vs. radial position. Section 3.1: Airfoil Selections Both the two bladed and three bladed Carolus turbines as well as the Whisper turbine use a single airfoil along the whole length of the blade. The Whisper turbine is a two bladed turbine that uses a Wortmann FX 60-126 Airfoil, shown in figure 3.1 below. This airfoil has a Cl,max of about 1.7 and stalls around a 13 o angle of attack, depending on the Reynolds number of the flow. It also has a maximum thickness of 12.6% and camber of 3.6% (Wortmann 2011). The Wortmann FX 60-126 airfoil has been around since the late 1960s and was originally designed to have fairly forgiving stall characteristics for sailplanes (Wortmann 2011). Figure 3.2 shows the lift coefficients of the Wortmann FX 60-126 airfoil for various Reynolds numbers: Figure 3.1: The Wortmann FX 60-126 Airfoil -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y/ c x/c Wortmann FX60-126 Airfoil 10 As can be seen in Figure 3.2, the Wortmann FX60-126 airfoil stalls around 13 o regardless of Reynolds number, and from there the lift coefficient trails off rather slowly with increasing angle of attack, hence the forgiving stall characteristics. It also has a maximum lift coefficient of about 1.6-1.7 depending on Reynolds number. Figure 3.3 shows the drag polar for the Wortmann airfoil: Figure 3.2: Cl vs. Angle of Attack for the Wortmann FX 60-126 Figure 3.3: Wortmann FX 60-126 Drag Polar -1.5 -1 -0.5 0 0.5 1 1.5 2 -15 -10 -5 0 5 10 15 20 25 C l Angle of Attack (Degrees) Cl vs. Angle of Attack Re = 300,000 Re = 500,000 Re = 700,000 Re = 900,000 Re = 1,100,000 Re = 1,300,000 Re = 1,500,000 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 -1.5 -1 -0.5 0 0.5 1 1.5 2 C d Cl Cd vs. Cl Re = 300,000 Re = 500,000 Re = 700,000 Re = 900,000 Re = 1,100,000 Re = 1,300,000 Re = 1,500,000 11 The FX60-126 has a rather wide drag bucket ranging from a lift coefficient of almost negative one all the way to close to 1.7 depending on Reynolds number. Increasing the Reynolds number at which this airfoil operates at significantly increases the size of the drag bucket in the range of lift coefficients greater than about 1.3. Figure 3.4 below shows Cl and Cd data as a function of angle of attack and expanded to the full range of angles of attack from -180 o to 180 o necessary for wind turbine analysis; however, the expansion method used goes beyond the scope of this thesis. With this range of angles of attack, it is difficult to discern any major differences in curve shapes between Reynolds numbers. However, figure 3.5, shows the coefficient of drag as a function of angle of attack for a limited range of angles to better illustrate the differences in the curves among Reynolds numbers: Figure 3.4: Expanded Cl and Cd vs. Angle of Attack for the Wortmann FX 60-126 -1.5 -1 -0.5 0 0.5 1 1.5 2 -200 -150 -100 -50 0 50 100 150 200 D im e n si o n le ss C o e ff ic ie n t V al u e Angle of Attack (Degrees) Cl and Cd vs. Angle of Attack Cl, Re = 300,000 Cd, Re = 300,000 Cl, Re = 500,000 Cd, Re = 500,000 Cl, Re = 700,000 Cd, Re = 700,000 Cl, Re = 900,000 Cd, Re = 900,000 14 As can be seen here, the S822 airfoil tends to stall around 15 o regardless of Reynolds number and the performance of the airfoil improves greatly with Reynolds number. As designed, the stall characteristics are also rather forgiving given the shallow drop in Cl past the stall angle (Somers 1993). Figure 3.8 shows the S822‟s drag polar: The S822 also has a wide drag bucket ranging from a lift coefficient of about negative one to 1.3 depending on Reynolds number. The drag bucket also drops significantly with Reynolds number. Once again, the airfoil data needs to be expanded out all the way from -180 o to 180 o . Figure 3.9 shows this data here: Figure 3.8: S822 Drag Polar 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -1.5 -1 -0.5 0 0.5 1 1.5 C d Cl Cl vs. Cd Re = 100,000 Re = 300,000 Re = 500,000 Re = 700,000 15 Once again, it is difficult to see the majority of the differences in lift and drag coefficients based on Reynolds number for this range of angles of attack. Figure 3.10 shows the drag coefficient as a function of angle of attack since the lift coefficient has already been presented in Figure 3.7: Figure 3.9: Expanded Cl and Cd vs. Angle of Attack for the S822 Figure 3.10: Cd vs. Angle of Attack for the S822 -1.5 -1 -0.5 0 0.5 1 1.5 -200 -150 -100 -50 0 50 100 150 200 D im e n si o n le ss C o e ff ic ie n t V al u e Angle of Attack (Degrees) Cl and Cd vs. Angle of Attack Cl, Re = 100,000 Cd, Re = 100,000 Cl, Re = 300,000 Cd, Re = 300,000 Cl, Re = 500,000 Cd, Re = 500,000 Cl, Re = 700,000 Cd, Re = 700,000 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 -20 -15 -10 -5 0 5 10 15 20 25 C d Angle of Attack (Degrees) Cd vs. Angle of Attack Re = 100,000 Re = 300,000 Re = 500,000 Re = 700,000 16 As was the case with the Wortmann FX60-126 airfoil, drag decreases with increased Reynolds number. However, drag does not tend to increase drastically until an angle of attack of about 15 o as compared to 10 o for the Wortmann FX60-126 airfoil. Section 3.2: Chord and Pitch Distributions Aside from the airfoil design, which is constant along the whole blade for each turbine, the pitch angle vs. span-wise location is critical in determining how the blade performs, because at various RPMs and wind speeds, the turbine blade must be kept from stalling or producing too much drag. Figure 3.11 below shows the pitch distribution for the Whisper turbine blades. The Whisper blade uses a simple linear pitch distribution along the span of its blades. Only performance analysis of radial distributions of various factors will indicate the results of this. It would appear that this was done to make manufacturing simple and not necessarily to optimize performance; however, this claim remains unsubstantiated until a thorough analysis can be conducted. Meanwhile, Figure 3.12 below shows the span-wise pitch distribution for the Carolus turbine blade: Figure 3.11: Whisper Blade Local Pitch Angle vs. Radial Location y = -7.2222x + 16.25 -2 0 2 4 6 8 10 12 14 16 18 0 0.5 1 1.5 2 2.5 Lo ca l P it ch A n gl e ( D e gr e e s) Radial Location (m) 19 to find out how much truth this statement bears. The same can be said for the Carolus turbine though since that design is currently under construction and has yet to be tested experimentally. Chapter 4 PSUWTA In order to continue Penn State‟s in house wind turbine research, it is necessary to have a proper open source BEMT code. This was the goal behind developing PSUWTA. As stated previously, BEMT is the most basic form of wind turbine analysis and yields fairly accurate results for its simplicity. Currently, the only code available for use is WT_Perf, a BEMT code developed and used by the National Wind Technology Center (NWTC), the wind energy division of the National Renewable Energy Laboratory (NREL). The problem here is that WT_Perf is not an open source code. In other words, the user does not have access to the source code nor knows how the code runs, what assumptions and equations it uses, etc. Although some of this can be gleaned from the input file necessary to run WT_Perf, it is necessary to have an open source code that can be edited and updated so that the user knows exactly how the code works and what assumptions it uses. Section 4.1: Code Specifics The purpose of this section is to document what additions and corrections to Blade Element Momentum Theory are being used in our code, PSUWTA. The code uses both an axial and a tangential induction factor to calculate the induction ratio, Prandtl hub and tip loss factors, and the Glauert correction for highly loaded elements. The code was originally written by Leo Albanese, a graduate student at Penn State, for large wind turbines and has since been modified to run for smaller wind turbines. PSUWTA begins by following the standard BEMT method presented in Chapter 2. The first of the corrections used in the code come during the calculation of the 21 induction ratios and are the Prandtl hub and tip loss factors. PSUWTA uses the following formulas for calculating the tip loss factor (Moriarty 1993): (4.1) ( ) (4.2) ( √ ) (4.3) In the above equations, is the number of blades on the turbine, is the non- dimensional radial location of the element and is the inflow angle. Similarly, the following equations are used in the calculation of the hub loss factor (Moriarty 1993): (4.4) (4.5) ( √ ) (4.6) Here, is the non-dimensional radial location of the hub, sometimes called the root cutout value. In order to calculate the total loss factor, the hub and tip loss factors are simply multiplied together (Moriarty 1993): (4.7) Once the local hub and tip loss factor is calculated, the code calculates the thrust coefficient for the blade element using equation 4.8 (Hansen 2008): (4.8) 24 For an RPM of 350, PSUWTA stays reasonably accurate relative to WT_Perf up to a wind speed of about 12 m/s. Of course, all of these predictions mean very little if they can‟t be verified against real world data. Research has been conducted at Penn State on the Whisper blades and experimental data has been collected; however, since the turbine is free spinning in a turbulent flow field where wind speed is constantly changing, it is difficult to compare for a specific RPM and wind speed combination. Nonetheless, comparisons can still be drawn. Figure 4.3 shows the raw power produced by the wind turbine at various wind speeds: Figure 4.2: Overall Power Output at 350 RPM Comparison 2 4 6 8 10 12 14 -2000 0 2000 4000 6000 8000 10000 Windspeed (m/s) P o w e r (W ) Power vs. Windspeed Whisper Blades 8 Reynolds Numbers - 15 Segments RPM = 350 PSUWTA WT_Perf 25 The connection between the data contained in Figure 4.3 and the comparisons of predicted power output in Figures 4.1 and 4.2 lies in the correlation between wind speed and RPM. Figure 4.4 shows this correlation as measured experimentally: Figure 4.3: Whisper Turbine Experimental Power Production 0 5 10 15 20 25 0 500 1000 1500 2000 2500 3000 3500 Raw Data Power Curve Wind Speed (m/s) A C P o w e r (W ) Raw Measured Power Manufacturer 26 With this in place, it should now be possible to draw comparison. For example, typically the Whisper turbine runs at 250 RPM for a wind speed of 4 m/s. As predicted by WT_Perf in Figure 4.1, the turbine should produce about 140 Watts; while according to PSUWTA, the turbine should produce about 175 Watts. Both of these predicted power outputs at 250 RPM fall well within the range of experimentally measured power outputs for the turbine at a wind speed of 4 m/s as shown in Figure 4.3. To check another case, the Whisper turbine tends to run at an RPM of 350 for a wind speed of about 7 m/s. According to WT_Perf and PSUWTA in Figure 4.2, the turbine should produce 1.35 kW and 1.26 kW respectively. Both of these predicted power outputs fall well within the range of experimentally measured power outputs and land almost right on the Manufacturer‟s line on Figure 4.3. Beyond the total power output, the power coefficient makes for an excellent tool when measuring the aerodynamic efficiency of a wind turbine. According to the Betz limit, the theoretical maximum power coefficient for a wind turbine should be about Figure 4.4: Whisper Turbine Experimental RPM and Wind Speed Correlation 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 220 240 260 280 300 320 340 360 380 400 Ensemble Averaged RPM Performance Wind Speed (m/s) R o ta ti o n a l R a te , R P M 29 The radial distribution of power predicted by PSUWTA is nearly exact compared to that predicted by WT_Perf. The differences towards the tip of the blade are negligible when integrated to find the total power produced by the wind turbine under this combination of RPM and wind speed. BEMT does not directly calculate Cp; rather, it iterates to find a set of core variables for each blade element, such as the inflow angle and induction ratios. It then uses these core variables to calculate things such as the thrust, power, and torque produced by each blade element. Of all the core variables, one of the most important is the axial induction factor. The reason for this is that it determines how much the wind slows while passing through the turbine and in essence, dictates the maximum amount of power the turbine can theoretically extract from the wind. Once again, according to momentum theory as applied to an actuator disk and the Betz limit, the theoretical maximum power coefficient of 0.593 occurs at an axial induction ratio of 1/3. In theory Figure 4.7: Comparison of Radial Power Distribution Comparison 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 Non-Dim. Radial Position P o w e r (W ) Power vs. Non-Dim. Radial Position Whisper Blades 8 Reynolds Numbers - 15 Segments Windspeed = 8 m/s RPM = 250 PSUWTA WT_Perf 30 then, the radial distribution of induction factors should integrate to produce an overall value of 1/3 for an ideal wind turbine. BEMT has the capability to predict the radial distribution of axial induction factors. Figure 4.8 shows these distributions as predicted by both PSUWTA and WT_Perf: As can be seen in Figure 4.8 above, both PSUWTA and WT_Perf predict almost identical distributions of the axial induction ratio. There is some discrepancy towards the hub of the blade, but these differences even out when integrated over the length of the blade. Finally, there is reason to justify using airfoil data for eight different Reynolds numbers. Figure 4.9 shows the radial Reynolds Number distribution for an RPM of 250 and a wind speed of 8 m/s: Figure 4.8: Radial Distributions of the Axial Induction Factor Comparison 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Non-Dim. Radial Position In d u c ti o n R a ti o ( a ) Induction Ratio vs. Non-Dim. Radial Position Whisper Blades 8 Reynolds Numbers - 15 Segments Windspeed = 8 m/s RPM = 250 PSUWTA WT_Perf 31 Both PSUWTA and WT_Perf predict Reynolds Numbers in the range of 100,000- 300,000 along the span of the blade. At such low Reynolds Numbers, airfoil performance can change greatly. Now compare that range to that of a turbine running at 11 m/s and 450 RPM: Figure 4.9: Reynolds Number Distribution Condition 1 Comparison 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 x 10 5 Non-Dim. Radial Position L o c a l R e y n o ld s N u m b e r Loc. Reynolds # vs. Non-Dim. Radial Position Whisper Blades 8 Reynolds Numbers - 15 Segments Windspeed = 8 m/s RPM = 250 PSUWTA WT_Perf 34 As can be seen in the above figure, the rotation rate that produces the most power varies with wind speed. For wind speeds less than 4 m/s, a rotation rate of 150 RPM produces the most power. From 4 m/s to 7 m/s, 250 RPM produces the most power; from 7 m/s to 8.5 m/s, 350 RPM produces the most power; from 8.5 m/s to 11 m/s, 450 RPM produces the most power, and so on. As will be shown in Figure 5.2, this is due to the fact that these rotational rates produce the highest power coefficients. Figure 5.2 shows the comparison of power coefficient to wind speed for various RPMs: Figure 5.1: Power Output vs. RPM and Wind Speed 2 4 6 8 10 12 14 -2000 0 2000 4000 6000 8000 10000 12000 14000 Windspeed (m/s) P o w e r (W ) Power vs. Windspeed PSUWTA.m Airfoil: Whisper Blades 8 Reynolds Numbers - 15 Segments 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 35 It is also predicted here that the maximum power coefficient for the Whisper blades should be around 0.4. However, in reality the turbine is free spinning and therefore might not be spinning at the optimal rotational rate for any given wind speed. Figure 5.3 here displays the experimentally measured correlation between RPM and wind speed for the whisper turbine and should help provide some insight as to where the turbine is actually operating. Figure 5.2: Power Coefficient vs. RPM and Wind Speed 36 For a wind speed of 4 m/s, the Whisper turbine is actually operating right around 240 RPM; however, PSUWTA shows that the optimal rotational rate for this wind speed would probably be about 200 RPM. For a wind speed of 7 m/s, the Whisper turbine is running at just under 350 RPM; however, based on the results of PSUWTA, it can be inferred that a rotational rate of 300 RPM would be ideal for maximum power production. Based upon this analysis, it would seem that the Whisper turbine is not running at the ideal settings for maximum power production. The data presented in Figure 5.2 can be further non-dimensionalized into power coefficient versus tip speed ratio. Figure 5.4 shows the fully non-dimensionalized data: Figure 5.3: Experimentally Measured RPM vs. Wind Speed 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 220 240 260 280 300 320 340 360 380 400 Ensemble Averaged RPM Performance Wind Speed (m/s) R o ta ti o n a l R a te , R P M 39 produces. Figure 5.6 shows the radial power distribution of the whisper turbine for a wind speed of 8 m/s: For every case of rotational rate, the outboard elements tend to produce more power than the inboard elements. This should come as no surprise since these elements are moving at a much faster velocity relative to the wind. In order to analyze each blade element, it is helpful to look at the power coefficient of each element to determine which elements are the most efficient and which are not. Figure 5.7 shows the radial distribution of power coefficients: Figure 5.6: Radial Power Distribution of the Whisper Turbine, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -50 0 50 100 150 200 250 Non-Dim. Radial Position P o w e r (W ) Power vs. Non-Dim. Radial Position PSUWTA.m Airfoil: Whisper Blades 8 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 40 As can be seen here, the power coefficients remain fairly constant for much of the span of the turbine. Near the hub and the tip, these drop to zero due to hub and tip losses modeled as vortexes being shed. The rotational rates of 250 RPM to 450 RPM tend to produce the best overall distribution of power coefficients as they remain relatively high and constant for the whole span of the blade. When compared to Figure 5.2, it is easy to understand why these rotational rates generate the highest overall power coefficients for a wind speed of 8 m/s. For a rotational rate of 150 RPM, the turbine is simply not spinning fast enough to extract enough energy from the wind, hence the lower power coefficients. For rotational rates greater than 450 RPM, the turbine is spinning too fast and the forces generated by the wind on the turbine will not produce as much power as desired. According to momentum theory, the axial induction factor is directly correlated to how much energy is extracted from the flow. Therefore, it is important to look at the radial distribution of this performance factor since it encompasses all of the aerodynamic Figure 5.7: Power Coefficient vs. Radial Location, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 Non-Dim. Radial Position P o w e r C o e ff ic ie n t (C p ) Power Coefficient vs. Non-Dim. Radial Position PSUWTA.m Airfoil: Whisper Blades 8 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 41 forces experienced by the each element and links this to the power produced by each element. Figure 5.8 shows the radial distribution of induction factors for the whisper turbine. Ideally, the axial induction ratio should be about 1/3 (0.333) regardless of radial location since this corresponds to the Betz limit and should produce the maximum theoretical power coefficient of 0.593. This is not the case however for the Whisper turbine as can be seen above in Figure 5.8. Having too high or too low of an induction factor will result in a less than ideal power coefficient for any element. Equation 5.1, shown below, has been derived from equations 2.3, 2.5, 2.7, and 2.8. This has been done to show the direct relationship of each aerodynamic factor to the axial induction factor: Figure 5.8: Axial Induction Factor vs. Radial Location, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Non-Dim. Radial Position In d u c ti o n R a ti o ( a ) Induction Ratio vs. Non-Dim. Radial Position PSUWTA.m Airfoil: Whisper Blades 8 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 44 Whisper turbine tends to operate in wind speeds ranging from 4 m/s to 8 m/s at Penn State, a range of about 100 N to 400 N (20-90lbs.) can be expected. More importantly, the radial distribution of thrust force needs to be analyzed to ensure the root bending moment of the turbine is not too great. Figure 5.11 here shows the radial distribution of the thrust force for a wind speed of 8 m/s: Overall, the majority of the force is concentrated towards the tip of the blade, but even so, assuming that all of the thrust force were concentrated at the tip, this would result in an overall bending moment of 400N * 2.25m = 900 N-m (~660 ft-lbs) under normal operating conditions. One last structural issue needs to be analyzed and that is the overall torque that the shaft must handle in order to allow the turbine to operate. Figure 5.12 shows the overall torque produced by the turbine for various wind speeds and rotational rates. Figure 5.11: Thrust Force vs. Radial Location, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 300 350 400 Non-Dim. Radial Position T h ru s t p e r L e n g th ( N /m ) Thrust per Length vs. Non-Dim. Radial Position PSUWTA.m Airfoil: Whisper Blades 8 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 45 Overall, the torque experienced by the shaft would probably not exceed about 75 N-m under usual operating conditions (4 – 8 m/s, 150 – 350 RPM). Of course, none of this structural analysis includes radial forces due to centrifugal forces caused by the rotation of the turbine, twisting or bending moments experienced by the blade, nor does it include dynamic analysis with respect to resonance frequencies and changing wind speeds and rotational rates. It is simply a baseline structural analysis that can easily be derived from BEMT. Overall, the whisper turbine appears to have much room for improvement. Currently, as can be inferred by the data generated by PSUWTA, the turbine is not operating at the ideal rotational rate for any given wind speed. Furthermore, the blade design used by the whisper turbine leaves plenty of room for improvement with respect to power producing capabilities. Figure 5.12: Overall Torque vs. Wind Speed Chapter 6 Carolus Turbine Analysis – Two-Bladed Using the blade geometries laid out in Chapter 3 for the Carolus turbine blades, an analysis of these blades in the two-bladed configuration was done using PSUWTA. The same methodology applied to the Whisper turbine will be applied here. Section 6.1: Power Output and Performance Once again, the most logical place to start is the overall power production. Figure 6.1 shows the power output for various rotational rates: Figure 6.1: Power Output vs. RPM and Wind Speed 2 4 6 8 10 12 14 -2000 0 2000 4000 6000 8000 10000 12000 14000 16000 Windspeed (m/s) P o w e r (W ) Power vs. Windspeed PSUWTA.m Carolus Blades 4 Reynolds Numbers - 15 Segments 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 49 shallow. At optimal tip speed ratios, the Carolus turbine is predicted to outperform the Whisper turbine; however, the Carolus blade is much more susceptible to worse performance at off design conditions than the Whisper turbine. In general, the Carolus turbine‟s power coefficient drops much more rapidly with increasing tip speed ratio than does the power coefficient for the Whisper turbine. How this will affect overall power output in reality is impossible to tell from this analysis without experimental data. It is possible that in a turbulent flow field where the wind speed is constantly changing, the gains made by optimizing the Carolus turbine will be offset by the amount of time spent in off-design conditions. Aside from overall performance, it is important to look at how the Carolus turbine blades perform along the length of the blade in order to identify which parts contribute the most and least to the overall performance and to see where possible improvements could be made. Figure 6.4 shows the radial power distribution produced by the Blades when run through PSUWTA: Figure 6.4: Power Output vs. Radial Location, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 Non-Dim. Radial Position P o w e r (W ) Power vs. Non-Dim. Radial Position PSUWTA.m Carolus Blades 4 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 50 Once again, the majority of the power produced by the Carolus blades comes from near the blade tip and drops off near the hub. As for where the power production becomes negative in some of the higher RPM cases, this due to the fact that the turbine actually requires energy to keep these elements moving at this rate rather than be driven by the wind. When compared to the Whisper turbine, the power production near the tip of the blade is actually slightly less. However, the Carolus Blade has the advantage in the middle and near the hub of the blade. In these areas, the power produced by each element in the Carolus turbine is greater than that of the Whisper Turbine, and in the end, this small advantage distributed over the majority of the blade causes the Carolus turbine to outperform the Whisper turbine in power production. Looking at the radial distribution of power coefficients is once again useful for determining how efficient each blade element is. Figure 6.5 shows this distribution: Figure 6.5: Power Coefficient vs. Radial Location, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 Non-Dim. Radial Position P o w e r C o e ff ic ie n t (C p ) Power Coefficient vs. Non-Dim. Radial Position PSUWTA.m Carolus Blades 4 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 51 For a wind speed of 8 m/s, the optimal rotational rate falls right between 250 and 350 RPM, hence the curves for these rotational rates having the highest values of power coefficient across the span of the blade. Outside of these rotational rates, the power coefficient values tend to fall off rather quickly, especially near the tip of the blade. Compared to the whisper turbine, around the optimal rotational rate, the Carolus blade has consistently higher Cp values. However, as conditions move away from optimal, the Whisper blades begin to outperform the Carolus blades. Beyond the power coefficient, the axial induction factor is a good indicator of how to attempt to modify the blade if need be since it shows whether the wind passing through the turbine slows down too much or too little. Once again, the optimal value for the axial induction factor is one third. Figure 6.6 shows the radial distribution of this factor: Figure 6.6: Axial Induction Factor vs. Radial Location, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Non-Dim. Radial Position In d u c ti o n R a ti o ( a ) Induction Ratio vs. Non-Dim. Radial Position PSUWTA.m Carolus Blades 4 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 54 As would be expected, with increasing wind speed and rotational rate comes a higher thrust force. For every combination of wind speed and rotational rate shown here, the Carolus turbine experiences a higher thrust force than the Whisper turbine under the same condition. The thrust curves for both follow the same basic trend and shape overall. Aside from just affecting what the turbine is mounted on, the thrust force can possibly bend the blades out of the rotor plane. In order to withstand this, it is necessary to understand how the thrust force is distributed along each blade. Figure 6.9 shows the thrust force per unit length as a function of radial location on the Carolus blades: Figure 6.8: Overall Thrust vs. Wind Speed and RPM 2 4 6 8 10 12 14 0 200 400 600 800 1000 1200 1400 1600 1800 Windspeed (m/s) T h ru s t (N ) Thrust vs. Windspeed PSUWTA.m Carolus Blades 4 Reynolds Numbers - 15 Segments 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 55 Following the same trend as the overall thrust force, the incremental thrust force on the Carolus blade is higher for every combination of wind speed and rotational rate. And once again, the majority of the thrust force is concentrated around the blade tip, leading to higher root bending moments that must be accounted for when it comes time to decide how to construct the blades. The last piece of structural information that can be easily obtained is the overall torque produced by the turbine that must be handled by the shaft connected to the electrical generator. The torques shown here are solely those resulting purely from aerodynamic forces and does not include things such as electrical forces, gravity, or friction. Figure 6.10 shows the overall torque experienced by the turbine: Figure 6.9: Thrust Force vs. Radial Location, Wind Speed = 8 m/s 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 300 350 400 450 Non-Dim. Radial Position T h ru s t p e r L e n g th ( N /m ) Thrust per Length vs. Non-Dim. Radial Position PSUWTA.m Carolus Blades 4 Reynolds Numbers - 15 Segments Windspeed = 8 m/s 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 56 Since power is simply the product of torque multiplied by rotational rate, it comes as no surprise that the torque experienced by the Carolus blade is higher than that experienced by the Whisper turbine. This applies only at the optimal combinations of wind speed and rotational rate that produce a maximum power coefficient. Outside of these optimum combinations though, the Whisper turbine may experience a higher torque since it sometimes will be producing more power in many cases. No matter what though, these torques must accounted for when designing the structural parts of the turbine. Figure 6.10: Overall Torque vs. Wind Speed and RPM 2 4 6 8 10 12 14 0 50 100 150 200 250 300 Windspeed (m/s) T o rq u e ( N -m ) Torque vs. Windspeed PSUWTA.m Carolus Blades 4 Reynolds Numbers - 15 Segments 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 59 Here, the optimal tip speed ratio appears to be about seven, while for the two bladed version of the turbine, it was closer to nine. However, the biggest change is clearly what happens away from the optimal tip speed ratio. Operating the three-bladed Carolus turbine anywhere outside of tip speed ratios ranging from about five to ten results in a massive loss in performance and power output. This is a much smaller range than that for the two-bladed Carolus turbine or even the Whisper turbine. As far as the radial distribution of performance parameters such as power output and power coefficient go, they followed nearly the exact same trends as with the two- bladed Carolus turbine with nothing really worthy of note to be shown. The same also happened with the radial distribution of axial induction ratios as well as angles of attack. Figure 7.3: Power Coefficient vs. Tip Speed Ratio 0 2 4 6 8 10 12 14 16 18 20 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Tip Speed Ratio C p Cp vs. Tip Speed Ratio PSUWTA.m Airfoil: Whisper Blades 8 Reynolds Numbers - 15 Segments 150 RPM 250 RPM 350 RPM 450 RPM 550 RPM 650 RPM 60 Section 7.2: Structural Loads For structural purposes, the overall thrust experienced by the turbine increased with the addition of a third blade. Even though the blades have been pitched by 6 o to keep the solidity the same, the fact there is now a third blade exposed to the wind increases the thrust experienced by the turbine regardless. As far as the radial thrust distribution goes, the majority of the thrust force is still concentrated at the tip. However, since the overall thrust force is now distributed over three blades instead of two, the root bending moment of each blade is less than that of the two-bladed Carolus turbine. Finally, the overall torque increases drastically for the three bladed-Carolus turbine when compared to the two-bladed version. This comes as no surprise since with a three bladed setup, higher power outputs can be obtained from lower rotational rates. Since power is the product of rotational rate multiplied by torque, an increase in power, and a decrease in rotational rate, results in a much greater torque that the turbine will experience. 61 Chapter 8 Conclusion PSUWTA, a Blade Element Momentum Theory code written in MATLAB, was developed in order to be used as a starting point for analyzing small wind turbines at Penn State. When compared to the output of WT_Perf, NWTC‟s BEMT code, PSUWTA continually comes within engineering accuracy of WT_Perf‟s results, all of which was checked against experimental data collected from the Whisper turbine used in both programs. This code can be used for preliminary analysis of new wind turbine designs as well as to verify the output of other more complicated wind turbine analysis codes. PSUWTA was then used to analyze the Performance of Southwest Windpower‟s Whisper 500 wind turbine as well as the Carolus wind turbine, an in-house design by Dr. Thomas Carolus, a visiting professor from the Institute for Fluid Dynamics and Thermodynamics at the University of Siegen in Siegen, Germany. The Carolus turbine comes in two versions, a two bladed or a three bladed turbine, both of which were analyzed. Both the Whisper and Carolus turbines are the same size, with their differences lying in airfoil selection, chord distribution, and twist distribution. The results for the Whisper Turbine indicate that its design, which uses linear distributions for both the chord length and local pitch angle, was made with simplicity in mind and not for aerodynamic performance. Even at the optimal combinations of wind speed and rotational rate where the turbine is producing the maximum amount of power at a maximum power coefficient, the inboard 20% of the blade is stalled. Even worse, it‟s overall power coefficient tends to only reach a maximum of about 0.39, which falls well short of the theoretical maximum of 0.593 set by the Betz limit. At lower wind speeds and rotational rates, this maximum falls as low as 0.32 in some cases. One redeeming quality though is that outside of optimal tip speed ratios, performance does degrade too quickly. 64 Appendix B PSUWTA Code Section B.1: Code Setup Instructions PSUWTA, as it stands at the time of writing, is not as cleaned up as it should be. In order to run it, the user must have a good foundation with MATLAB since many things are hard coded. The following items are needed in order to fully define the turbine in the program: 1) The number of blades on the turbine 2) The radius of the turbine in meters 3) The radius of the hub in meters (also called the root cutout value) 4) The number of blade elements the user wishes to use 5) Airfoil Data covering the full spectrum of angles of attack from -180 o to 180 o 6) The local pitch angle in radians as a function of non-dimensional radial location (r/R) 7) The chord length in meters as a function of non-dimensional radial location (r/R) 8) WT_Perf Performance Data (Optional – Only if the user desires to use it for comparison) Items 1-4 are simply be entered into lines 14-17 of PSUWTAOp as well as in lines 23 and 24 of PSUWTA. From here, these instructions will describe how to input the rest of this information into the code. Section B.1.1: Airfoil Data Entering the airfoil data is the most difficult part of running PSUWTA. First of all, the user must select the airfoil used by the wind turbine. Currently, the code only runs 65 using one airfoil for the entire span of the blade. The user can obtain the airfoil data for various Reynolds numbers from any source such as experimental data or XFoil. From here, the data must be expanded to cover the whole spectrum of angles of attack from -180 o to 180 o . This can be done using a Microsoft Excel spreadsheet called AirfoilPrep_v2p2.xlsm. It can be found online here: http://wind.nrel.gov/designcodes/preprocessors/airfoilprep/ Once this data has been obtained for however many Reynolds numbers the user would like to use, it must be imported into the program. The easiest way to do this, at least in the opinion of the author, is to create a function to store all of the data. An example function used for the S822 airfoil is provided below: In this function, the variable storage is used to hold all of the airfoil data. It is a three dimensional array, defined as storage(x,y,z). Here, x selects which angle of attack to extract data from. y can have values of 1, 2, or 3 where 1 corresponds to the value of the angle of attack itself; 2 corresponds to the lift coefficient; and 3 corresponds to the drag coefficient. Finally, z selects which Reynolds number to select data from. The user can set this up to have as many data sets for as many Reynolds numbers as desired. Since this function requires no input and simply returns all of the airfoil data, it Figure B.1: Airfoil Data Storage Function function ren = S822data storage(:,:,1) = [α1 Cl Cd α2 Cl Cd α3 Cl Cd]; storage(:,:,2) = [α1 Cl Cd α2 Cl Cd α3 Cl Cd]; storage(:,:,3) = [α1 Cl Cd α2 Cl Cd α3 Cl Cd]; ren = storage; end 66 essentially acts as a global variable and can be simply called. The function is called in PSUWTAOp at line 19. The user simply inserts the name of the function as shown below at that line: Now, that the airfoil data has been updated, it is now necessary for the user to update the Reynolds number data selection code. Open up the function PSUWTA and go to line 162. This is where the experience with MATLAB becomes necessary. An „if‟ structure must be constructed for the function to select which airfoil data set to use based on Reynolds number and should look like this when complete: The variable names used in this „if‟ structure should be very straightforward: Re represents Reynolds number, ren(:,:,z) is the variable containing all of the airfoil data with the z index indicating Reynolds number for the data set and airfoil is the variable that stores the airfoil data set for the Reynolds number that the current element is being analyzed at. The selection structure uses airfoil data for a specific Reynolds number for Figure B.2: Airfoil Data Function Call Figure B.3: Airfoil Data Selection Code ren = Function_Name; if Re <= 200000 airfoil = ren(:,:,1); elseif Re <= 400000 airfoil = ren(:,:,2); elseif Re <= 600000 airfoil = ren(:,:,3); else airfoil = ren(:,:,4); end 69 2) Do you wish to plot the radial power, torque, and thrust? [Y/N] : 3) Do you wish to plot the overall power? [Y/N] : 4) Do you wish to plot the overall thrust? [Y/N] : 5) Do you wish to plot the overall torque? [Y/N] : Question one will plot six figures: angle of attack, lift coefficient, drag coefficient, local velocity, local Reynolds number, axial induction ratio, and tangential induction ratio. Each of these variables will be on its own plot against non-dimensional radial location. Question two will also plot six figures: power per unit length, thrust per unit length, torque per unit length, power coefficient, thrust coefficient, and torque coefficient. Once again, each will appear on its own plot against non-dimensional radial location. Questions three through five each produce three figures a piece. Question number three plots power vs. wind speed, power coefficient vs. wind speed, and power coefficient vs. tip speed ratio. Each of these is on its own plot. Questions four and five produce the same plots except for the overall thrust or torque respectively. Finally, the program asks the user if he or she would like to select a new RPM and wind speed and plot again. This plotting loop can be run as many times as the user desires without the program having to recalculate anything. 70 iagram Block Di Section B.2 eJep a4} S}O}d (OZ pe|q B43 JO Auptjos jed0} ayy se anjea Wdy }xoU BY} JJ BT-g Sdais suny “anjen Way salepdn (6T SUSIDIJJ9OD IBY} SE jam se samod pue ‘anbio} ‘ysn4y} ‘Besp “I4I] [e}0} pue Jequatua.9U! a}ejND]e9 0} pT days Woy senjeA sasp (ST 1]2M se Aujawoad apejq ay} jo JO. paseq Juawaja ay} Jo ajSue youd Jeo] pue puoyp sazejnoje5 (ET fi “LT - § sdays ul pasn Wdy 404 E}ep [Ie S8s0}s pue ‘Sayep!]Osuod ‘sana!ay (ST “paeds pum yxeu 342 404 9T - OT Sdays suny “paads pum sayepdy (ZT “quawala apelq 3x9U 4 Joy ST - €T sdays suny uawaya apelq sayepdn (9T. Adessavau J! UO!}991J09 Juane|D, ay} se [Jam se sasso] diy pue qn} |ypuerg sapnjou| daquuiny s,pjouday - AyDOIeA |e90]_— - ypene jo ajsuy — - GGaws :aye|n9]e9 Ajsnoaueynuuls 0} doo} aanesoy! ue suny (pT Vy 85N 0} JUaWA|a ape|g S}a185 (ZT pazAjeue si Juawela Spelq yoiym Alen 0} doo| sulgag (6 asn 0} paads pulm syajas (OT sisAjeue 343 Ul pasn ads puim Aen 0} doo] suigag (6 poo SuluuNs UO 0} passed aq 0} sao}eUW T pue pT sdays wou sanjea jje sa1035 (ST = uonouny ay} Jo uns Siu 40} eyep 21035 0} papaau saouu}ew souaz pue $3739 (g (Aue J!) youd BAIIaI}OD —- oney paads diy - Wd - Ko paads pum - :Buyuyeyuod x1JjewW eyep sayeas9 (Z ones MOU! SuIje]NI]e9 JOj aoUeJa]OL —- aie Jo Aysoosiy anjea no yndyooy—- ajSuy Bu0D-alg = - snipey apejg — - Aysuag ar isuajetuesed yndul Sulmoj|o By} sias (9 ANIEA Wd 5192195 (S “9J0}S 0} papaau Sar1J}eW 34} ||e SOlezZ pue sa}easD (E ejep ay} |e JO/d pue ‘adeuew pazAjeue aq 0} anjea Wd 84] sayepdn yey} doo] e seis (b “papinoad j! a1ay eyep Jad 1M ul peas osje ued “susquunu s,pjouAay snouen ze jloJse 343 405 Beg pue J! JO SyUa!DYJa09 a4} Sa10ys pue Ul speay (Z | azAjeue 0} spaads pulm joaguey = - BUIGIN} ay} UO sapelq Jouaquiny = - asn 0} sjuawwala apejqjoseaquinn = - azAjeue OJ SanjeAWidy —- ssupyatuesed ndu ay} asooyD (T (uonouny sishjeue AJoayy wnqUaWOW JUaLA!a apelg) VLMINSd = PUNoISy2eg UB—I9 || (apoo juaWaseueW eyep Bul uns quos4) dOV.LMNSd = punosdy9eg anig 71 Section B.3: PSUWTA Source Code In the following subsections, the actual MATLAB source code for PSUWTA has been pasted in. In these copies, PSUWTA is set up to run an analysis on the two-bladed version of the Carolus turbine. Section B.3.1: BEMT Function Code “PSUWTA_v7_B.m” function [Clar dLa Cdar dDa Qty1 Qty2 dQ Cta Cpa P dP T dT Uinfo yar thetaa phia alphaa locvel locRe dCt dCp indrat tindrat dCq pn pt Q Cqa] = PSUWTA_v7_B(Nb,rpm,seg,ren,minwind,maxwind) % Program: PSUWTA_v7_B.m % Author: Leo Albanese, See performance.m for original version of code % Note: performance.m is heavily based off a code called WT_Perf produced % by NREL % Modified by Thomas Purcell 12/8/09 % Modified by Jacob Marsh and David Faust Spring 2010 % Modified by Jacob Marsh Fall 2010 % Modified by Thomas Purcell Spring 2011 % % About the Program: % % This function is designed to analyze the performance of a wind turbine % using Blade Element Momentum Theory. % % Please note that if the turbine uses collective pitch that you must enter % a function that defines collective pitch in degrees as a function of the % variable 'wind' (the wind speed in m/s). % Input R=2.25; % Assigns blade radius in meters Ro = 0.137; % Assigns root cut out value in m a_tol=1e-6; % Assigns the tolerance for inflow ratio convergence mu = 1.814826290763033E-5; %Viscosity of air at 60 deg F rho=1.225; % Assigns air density in kg/m^3 % Calculations 74 ca(q,i) = c; thetaa(q,i)=theta*180/pi; % Stores theta value into a global array ca(i)=c; % Stores chord value into a global array sig=Nb*c/(2*pi*y*R); % Calculates solidity at this radial position %Calculates induction factor a=0.0; % Initial value of axial induction factor del_a_O=0.0; % Stores previous a a_step=0.25; % Iteration step size of a del_a_test=1; % The delta a that is tested against aprime=0.0; % Initial Tangential Induction Factor ItSLSC=0; % Stores Iterations Since Last Sign Change of a iter=1; while (iter<5000 && abs(del_a_test)>=a_tol) % Stops performing loop when iterations Vi_tan=X*V*y*(1+aprime); % get too high or the del_a is below the Vi_norm=V*(1-a); % specified tolerance Vi=sqrt(Vi_tan^2+Vi_norm^2); %Local relative velocity Re = rho*Vi*c/mu; %Calculates Local Reynolds Number phi=atan(Vi_norm/Vi_tan); % Calculates inflow angle phi alpha=(phi-theta)*180/pi; % Calculates angle of attack alpha % Calculates Prandtl Losses at hub and tip fhub=(Nb/2)*((y-yo)/(yo*sin(phi))); Fhexp=exp(-fhub); Fhub=(2/pi)*atan(sqrt(1-Fhexp^2)/Fhexp); ftip=(Nb/2)*((1-y)/(y*sin(phi))); Ftexp=exp(-ftip); Ftip=(2/pi)*atan(sqrt(1-Ftexp^2)/Ftexp); F=Fhub*Ftip; k=1; % Counter variable for interpolating Cl and Cd Cl=0; % Initializing Cl and Cd Cd=0; 75 if ((alpha < -180) || (alpha > 180)) alpha=-10; end if Re <= 200000 airfoil = ren(:,:,1); elseif Re <= 400000 airfoil = ren(:,:,2); elseif Re <= 600000 airfoil = ren(:,:,3); else airfoil = ren(:,:,4); end % Loop to interpolate Cl and Cd while (Cd==0) if ((airfoil(k,1)<=alpha) && (airfoil(k+1,1)>=alpha)) Cl=airfoil(k,2)+(alpha- airfoil(k,1))*(airfoil(k+1,2)-airfoil(k,2))/(airfoil(k+1,1)- airfoil(k,1)); Cd=airfoil(k,3)+(alpha- airfoil(k,1))*(airfoil(k+1,3)-airfoil(k,3))/(airfoil(k+1,1)- airfoil(k,1)); else k=k+1; end end Cn = Cl*cos(phi)+Cd*sin(phi); % Calculates Cn to be used below Ct = Cl*sin(phi)-Cd*cos(phi); CT = sig*((1-a)^2)*Cn/(sin(phi))^2; CT = min(max(CT,-2),2); if CT > 0.96*F, if CT*(50-36*F)+12*F*(3*F-4) < 0, a_new = 0; else a_new = (18*F-20-3*sqrt(CT*(50-36*F)+12*F*(3*F- 4)))/(36*F-50); end else a_new = 1/(4*F*(sin(phi))^2/(sig*Cn)+1); end aprime = 1/((4*F*sin(phi)*cos(phi)/(sig*Ct))-1); del_a =a_new-a; 76 if ((del_a_O ~= 0.0) && (del_a/del_a_O < 0.0)) % Changes step size a_step=0.5*a_step; % of iterations ItSLSC=0; % accordingly elseif (ItSLSC == 10) a_step=2.0*a_step; ItSLSC=0; else ItSLSC=ItSLSC+1; end del_a_test=a_step*del_a; a=a+del_a_test; del_a_O=del_a_test; iter=iter+1; end indrat(q,i)=a; % Stores value of a into a global array tindrat(q,i)=aprime; % Stores value of a' into a global array phia(q,i)=phi*180/pi; % Stores phi value into a global array alphaa(q,i)=alpha; % Stores angle of attack in a global array locvel(q,i)=Vi; % Stores local velocity in a global array locRe(q,i) = Re; % Stores Local Reynolds number into a global array Clar(q,i)=Cl; % Stores Cl in global array Cdar(q,i)=Cd; % Stores Cd in global array dL = 1/2*rho*(Vi^2)*c*Cl; % Calculates incremental lift dLa(q,i)=dL; % Stores incremental lift dD = 1/2*rho*(Vi^2)*c*Cd; % Calculates Incremental Drag dDa(q,i)=dD; % Stores incremental drag Qty1(q,i)=dL*sin(phi); % Incremental tangential lift Qty2(q,i)=dD*cos(phi); % Incremental tangential drag pt(q,i)=Qty1(q,i)-Qty2(q,i); % Incremental tangential force pn(q,i)=dL*cos(phi)+dD*sin(phi); % Incremental normal force A = (pt(q,i)-pt(q,i-1))/(R*yit); B = (pt(q,i-1)*R*y-pt(q,i)*R*(y-yit))/(R*yit); dM = (1/3)*A*((R*y)^3-(R*(y-yit))^3)+(1/2)*B*((R*y)^2-(R*(y- yit))^2); dQact = Nb*dM; 79 PSUWTA_Ct = zeros (numwin,maxi); %PSUWTA Total Thrust Coefficient PSUWTA_dCt = zeros (numwin,seg,maxi); %PSUWTA Incremental Thrust Coefficient PSUWTA_T = zeros (numwin,maxi); %PSUWTA Total Thrust in Newtons PSUWTA_Cq = zeros (numwin,maxi); %PSUWTA Total Torque Coefficient PSUWTA_dCq = zeros (numwin,seg,maxi); %PSUWTA Incremental Torque Coefficient PSUWTA_Q = zeros (numwin,maxi); %PSUWTA Total Torque PSUWTA_Cp = zeros (numwin,maxi); %PSUWTA Total Power Coefficient PSUWTA_dCp = zeros (numwin,seg,maxi); %PSUWTA Incremental Power Coefficient PSUWTA_P = zeros (numwin,maxi); %PSUWTA Total Power in watts PSUWTA_Uinfo = zeros (numwin,4,maxi); %page 1 = wind speed (m/s) %page 2 = tip speed ratio %page 3 = collective pitch (deg) %page 4 = RPM PSUWTA_TSR = zeros (numwin,maxi); %PSUWTA Tip Speed Ratio Matrix, %row= wind speed, column = RPM PSUWTA_yar = zeros (maxi,seg); %Global array for plotting dimensionless radius PSUWTA_Theta = zeros (numwin,seg,maxi); %PSUWTA theta PSUWTA_Phia = zeros (numwin,seg,maxi); %PSUWTA phi PSUWTA_AoA = zeros (numwin,seg,maxi); %PSUWTA angle of attack PSUWTA_LocVel = zeros (numwin,seg,maxi); %PSUWTA Local Velocity PSUWTA_LocRe = zeros (numwin,seg,maxi); %PSUWTA Reynolds Number PSUWTA_dT = zeros (numwin,seg,maxi); %PSUWTA Radial Thrust PSUWTA_dP = zeros (numwin,seg,maxi); %PSUWTA Radial Power PSUWTA_indrat = zeros (numwin,seg,maxi); %PSUWTA Induction Ratio (a) PSUWTA_tindrat = zeros (numwin,seg,maxi); %PSUWTA Tangential Induction Ratio (a') PSUWTA_pn = zeros (numwin,seg,maxi); %PSUWTA Normal Force per length at each radial station PSUWTA_pt = zeros (numwin,seg,maxi); %PSUWTA Tangential Force per length at each radial station WTPerf_totP = zeros (67,19); %WT_Perf initialization WTPerf_totCp = zeros (67,19); %WT_Perf initialization WTPerf_totT = zeros (67,19); %WT_Perf initialization WTPerf_totCt = zeros (67,19); %WT_Perf initialization WTPerf_totQ = zeros (67,19); %WT_Perf initialization WTPerf_totCq = zeros (67,19); %WT_Perf initialization WTPerf_TSR = zeros (67,19); %WT_Perf initialization i = 1; % each i count is another RPM value run 80 while i <=maxi [PSUWTA_Cl(:,:,i)... PSUWTA_dL(:,:,i)... PSUWTA_Cd(:,:,i)... PSUWTA_dD(:,:,i)... PSUWTA_Qty1(:,:,i)... PSUWTA_Qty2(:,:,i)... PSUWTA_dQ(:,:,i)... PSUWTA_Ct(:,i)... PSUWTA_Cp(:,i)... PSUWTA_P(:,i)... PSUWTA_dP(:,:,i)... PSUWTA_T(:,i)... PSUWTA_dT(:,:,i)... PSUWTA_Uinfo(:,:,i)... PSUWTA_yar(i,:)... PSUWTA_Theta(:,:,i)... PSUWTA_Phia(:,:,i)... PSUWTA_AoA(:,:,i)... PSUWTA_LocVel(:,:,i)... PSUWTA_LocRe(:,:,i)... PSUWTA_dCt(:,:,i)... PSUWTA_dCp(:,:,i)... PSUWTA_indrat(:,:,i)... PSUWTA_tindrat(:,:,i)... PSUWTA_dCq(:,:,i)... PSUWTA_pn(:,:,i)... PSUWTA_pt(:,:,i)... PSUWTA_Q(:,i)... PSUWTA_Cq(:,i)... ]... = PSUWTA_v7_B_Carolus... ... (Nb,RPMval(i),seg,ren,minwindspeed,maxwindspeed); PSUWTA_TSR(:,i) = PSUWTA_Uinfo(:,2,i); i = i + 1; end %{ %Read in WT_Perf data filename = 'Data/Test06_SWWT_ManuBlades_15.bid'; [WTPerf_windspeed WTPerf_RPM WTPerf_LocVel WTPerf_Re WTPerf_Loss WTPerf_AxialInd WTPerf_TangInd WTPerf_AoA WTPerf_AlfaD WTPerf_Cl WTPerf_Cd WTPerf_Ct WTPerf_Cq WTPerf_Cp WTPerf_ThrustLen WTPerf_TorqueLen WTPerf_Power] = WTPerfLoad(filename); i = 1; while i <= 67; y = 1; while y <= 19; 81 WTPerf_totP(i,y) = sum(WTPerf_Power(:,i,y)); WTPerf_totCp(i,y) = WTPerf_totP(i,y)/(0.5*rho*(WTPerf_windspeed(1,i)^3)*pi*R^2); WTPerf_totT(i,y) = sum(WTPerf_ThrustLen(:,i,y)*(segwidth*Nb)); WTPerf_totCt(i,y) = WTPerf_totT(i,y)/(0.5*rho*(WTPerf_windspeed(1,i)^2)*pi*R^2); WTPerf_totQ(i,y) = sum(WTPerf_TorqueLen(:,i,y)*(segwidth*Nb)); WTPerf_totCq(i,y) = WTPerf_totQ(i,y)/(0.5*rho*(WTPerf_windspeed(1,i)^2)*pi*R^3); w = WTPerf_RPM(y)*2*pi/60; % Calculates angular velocity of blade in rad/s vtip=w*R; WTPerf_TSR(i,y) = vtip/WTPerf_windspeed(i); y = y + 1; end i = i + 1; end WTPerf_Re = WTPerf_Re*1.0e+06 ; %Makes WT_Perf and PSUWTA have the same units %} Plotagain = 'y'; while Plotagain =='y' || Plotagain =='Y'; close all clc plotted = 0; %Prompt for Windspeed and RPM for graphs wspd = input('Input desired Windspeed for graphs (>=2,<=15): '); wspdval = (2:15); while wspd == all(wspdval) wspd = input('Windspeed entered not valid, please enter one of the following (>=2,<=15 INTEGERS ONLY): '); end rpmselect = input('Input desired RPM for graphs (150,250,350,450,550,650): '); while rpmselect == all(RPMval) rpmselect = input('RPM entered not valid, please enter one of the following (150,250,350,450,550,650): '); end %Calculations for what value to call when looking for desired RPM or %Windspeed value in PSUWTA_PSU and WT_Perf 84 title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/locRevsNDR.jpg') saveas(gcf,'Output_v3D/locRevsNDR.fig') %Displays a graph of Radial Station Induction Ratio (a) figure(6) plot(yar,PSUWTA_indrat(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_AxialInd(:, WTPerf_wind,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Induction Ratio (a)','FontWeight','bold') LABEL(1) ={'Induction Ratio vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/avsNDR.jpg') saveas(gcf,'Output_v3D/avsNDR.fig') %Displays a graph of Radial Station Tangential Induction Ratio (a') figure(7) plot(yar,PSUWTA_tindrat(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_TangInd(:, WTPerf_wind,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Tangential Induction Ratio (a prime)','FontWeight','bold') LABEL(1) ={'Tangential Induction Ratio vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/aprimevsNDR.jpg') saveas(gcf,'Output_v3D/aprimevsNDR.fig') end ques1b = input('Do you wish to plot the radial power, torque, and thrust? [Y/N] : ', 's'); if isempty(ques1b) ques1b = 'N'; end if ques1b =='N' || ques1b =='n' %Does nothing and skipps plots 8-15 else %Executes the Plot functions plotted = 1; %Displays a graph of Coefficient of Thrust (Ct) figure(8) plot(yar,PSUWTA_dCt(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_Ct(:,WTPerf_wi nd,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Thrust Coefficient (Ct)','FontWeight','bold') 85 LABEL(1) ={'Thrust Coefficient vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/CtvsNDR.jpg') saveas(gcf,'Output_v3D/CtvsNDR.fig') %Displays a graph of Torque Coefficient (Cq) figure(9) plot(yar,PSUWTA_dCq(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_Cq(:,WTPerf_wi nd,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Torque Coefficient (Ct)','FontWeight','bold') LABEL(1) ={'Torque Coefficient vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/CqvsNDR.jpg') saveas(gcf,'Output_v3D/CqvsNDR.fig') %Displays a graph of Power Coefficient (Cp) figure(10) plot(yar,PSUWTA_dCp(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_Cp(:,WTPerf_wi nd,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Power Coefficient (Cp)','FontWeight','bold') LABEL(1) ={'Power Coefficient vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/CpvsNDR.jpg') saveas(gcf,'Output_v3D/CpvsNDR.fig') %Displays a graph of Power (W) figure(11) plot(yar,PSUWTA_dP(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_Power(:,WTPerf_ wind,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Power (W)','FontWeight','bold') LABEL(1) ={'Power vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/PowervsNDR.jpg') saveas(gcf,'Output_v3D/PowervsNDR.fig') %Displays a graph of Thrust per Length (T) figure(12) plot(yar,PSUWTA_dT(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_ThrustLen(:,WTP erf_wind,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Thrust per Length (N/m)','FontWeight','bold') LABEL(1) ={'Thrust per Length vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') 86 legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/ThrustvsNDR.jpg') saveas(gcf,'Output_v3D/ThrustvsNDR.fig') %Displays a graph of Torque per Length (Q) figure(13) plot(yar,PSUWTA_dQ(PSUWTA_wind,:,PSUWTA_rpm),yar,WTPerf_TorqueLen(:,WTP erf_wind,WTPerf_rpm)','LineWidth', 2) xlabel('Non-Dim. Radial Position','FontWeight','bold') ylabel('Torque per Length (N-m/m)','FontWeight','bold') LABEL(1) ={'Torque per Length vs. Non-Dim. Radial Position'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Best'); grid on saveas(gcf,'Output_v3D/TorquevsNDR.jpg') saveas(gcf,'Output_v3D/TorquevsNDR.fig') end ques1c = input('Do you wish to plot the overall power? [Y/N] : ', 's'); if isempty(ques1c) ques1c = 'N'; end if ques1c =='N' || ques1c =='n' %Does nothing and skipps plots 8-15 else %Executes the Plot functions plotted = 1; LABEL(4) = cellstr(['RPM = ',num2str(rpmselect)]); %Displays a graph of Power vs. Wind figure(14) plot(wspdval,PSUWTA_P(:,PSUWTA_rpm),WTPerf_windspeed,WTPerf_totP(:,WTPe rf_rpm),'LineWidth', 2) xlabel('Windspeed (m/s)','FontWeight','bold') ylabel('Power (W)','FontWeight','bold') LABEL(1) ={'Power vs. Windspeed'}; title (LABEL,'FontWeight','bold') legend(legendlab,'Location','Southeast'); %,'Location','Best' grid on xlim([2 15]) saveas(gcf,'Output_v3D/PowervsWind.jpg') saveas(gcf,'Output_v3D/PowervsWind.fig') %Displays a graph of Cp vs. Wind figure(15) plot(wspdval,PSUWTA_Cp(:,PSUWTA_rpm),WTPerf_windspeed,WTPerf_totCp(:,WT Perf_rpm),'LineWidth', 2)