HindawiPublishingCorporation JournalofRenewableEnergy Volume2013,ArticleID839319,17pages http://dx.doi.org/10.1155/2013/839319 Research Article Computational Study on the Aerodynamic Performance of s Wind Turbine Airfoil Fitted with Coand Jet H.Djojodihardjo,1,2M.F.AbdulHamid,2A.A.Jaafar,3S.Basri,3F.I.Romli,2 F.Mustapha,2A.S.MohdRafie,2andD.L.A.AbdulMajid2 1UniversitasAlAzharIndonesia,KebayoranBaru,Jakarta12110,Indonesia 2DepartmentofAerospaceEngineering,UniversitiPutraMalaysia,34400Serdang,Selangor,Malaysia 3FacultyofManufacturingEngineering,UniversityMalaysiaPahang,26300Kuantan,Pahang,Malaysia CorrespondenceshouldbeaddressedtoH.Djojodihardjo;[email protected] Received16January2013;Accepted2April2013 AcademicEditor:Wei-HsinChen Copyright©2013H.Djojodihardjoetal. This is an open access article distributed under the Creative Commons Attribution License,whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperly cited. Variousmethodsofflowcontrolforenhancedaerodynamicperformancehavebeendevelopedandappliedtoenhanceandcontrol thebehaviorofaerodynamiccomponents.TheuseofCoanda˘effectfortheenhancementofcirculationandlifthasgainedrenewed interest, in particular with the progress of CFD. The present work addresses the influence, effectiveness, and configuration of Coanda˘-jetfittedaerodynamicsurfaceforimprovingliftand𝐿/𝐷,specificallyforS809airfoil,withaviewonitsincorporation inthewindturbine.Asimpletwo-dimensionalCFDmodelingusing𝑘-𝜀turbulencemodelisutilizedtorevealthekeyelements thatcouldexhibitthedesiredperformanceforaseriesofS809airfoilconfigurations.Parametricstudyperformedindicatesthat theuseofCoanda˘-jetS809airfoilcanonlybeeffectiveincertainrangeoftrailingedgerounding-offradius,Coanda˘-jetthickness, andmomentumjetsize.ThelocationoftheCoanda˘-jetwasfoundtobeeffectivewhenitisplacedclosetothetrailingedge.The resultsarecomparedwithexperimentaldataforbenchmarking.Three-dimensionalconfigurationsaresynthesizedusingcertain acceptableassumptions.Atrade-offstudyontheS809Coanda˘configuredairfoilisneededtojudgetheoptimumconfigurationof Coanda˘-jetfittedWind-Turbinedesign. 1.Introduction to50×109MWh per year against the current total annual world electricity consumption of about 15 × 109MWh. The In line with the efforts to enhance the use of green energy global wind power capacity installed in the year 2004 was technology for energy extraction, conversion, and propul- 6614MW,anincreaseintotalinstalledgeneratingcapacityof sion, the fundamental principles and mechanisms that play nearly20%fromtheprecedingyear. keyrolesinthesetechnologieshavebeenthefocusinmany Wind energy conversion research and development currentresearcheffortsaswellasthepresentresearchwork, efforts have their origin on fluid physics, thermodynamics, since these have the eventual potential of nationaldevelop- andmaterialsciences.Learningfromvariousinnovationsin- mentsignificance.Sinceoneofthefirstattemptstogenerate troducedinaircrafttechnology,windenergytechnologycan electricitybyusingthewindintheUnitedStatesbyCharles takeadvantageofflowcirculationmodificationandenhance- Brushin1888[1],windenergyhasbeenutilizedforelectricity ment. generation using large Wind turbines in many wind farms Flowcontrolforenhancedaerodynamicperformancehas withpotentialwindenergyinlandaswellasoff-shore.Ithas takenadvantageofvarioustechniques,suchastheuseofjets been estimated that roughly 10 million MW of energy are (continuous,synthetic,pulsed,etc.),compliantsurface,vor- continuouslyavailableintheearth’swind[2].Thetechnical texcell,andthelike[3–16].Theseactivecontrolconceptscan potential of onshore wind energy is very large, 20 × 109 dramaticallyalterthebehaviorofaerodynamiccomponents 2 JournalofRenewableEnergy such as airfoils, wings, and bodies. Specifically, modified the rotor invariably stalls and experiences highly unsteady circulation and enhanced circulation have been utilized for aerodynamic forces and moments. It is noted that over the improvedaircraftandland-vehicleperformanceaswellasin pastseveralyears,therehasbeenanincreasedefforttoextend novelair-vehiclesliftandpropulsionsystemconcept,dueto theusableoperatingrangeofstall-regulatedWindturbines. their better energy conversion capabilities [17–23]. The po- Theconceptsofcirculationcontrolthroughtrailingedge(TE) tentialbenefitsofflowcontrolincludeimprovedperformance blowing,orCoanda˘-jet,thencameintoview. andmaneuverability,affordability,increasedrangeandpay- However, over the past decades, commonly used airfoil load,andenvironmentalcompliance.Activeflowcontrolmay familiesforhorizontalaxisWindturbinessuchastheNACA beappliedtodelayoradvancetransition,tosuppressoren- 44-,NACA23-,NACA63-,andNASALS(1)seriesairfoils hanceturbulence,ortopreventorpromoteseparation,that suffer noticeable performance degradation from roughness willresultindragreduction,liftenhancement,mixingaug- effects resulting from leading edge contamination [30]. mentation, heat transfer enhancement, and flow induced Annualenergylossesduetoleadingedgeroughnessaregreat- noise suppression [24], as appropriate. It could be men- estforstall-regulatedrotors.Thelossislargelyproportionalto tioned that more recently plasma-based flow control has the reduction in maximum lift coefficient (𝑐𝑙,max) along the alsobeenconsideredforwindturbineapplications[25–28]. blade.Highblade-roottwisthelpsreducethelossbykeeping Sucheffortsmayinvolvethedesignof“smart”windturbine theblade’sangleofattackdistributionawayfromstallwhich bladeswithintegratedsensor-actuator-controllermodulesto delays the onset of 𝑐𝑙,max. Roughness also degrades the improvetheperformanceofwindturbines,byenhancingen- airfoil’slift-curveslopeandincreasestheprofiledrag,which ergy capture, and reduction of aerodynamic loading and contributestofurtherlosses.Forstall-regulatedrotors,whose noise by way of virtual aerodynamic shaping. These are angle of attack distribution increases with wind speed, the beyondthescopeofthepresentwork. annual energy loss can be 20% to 30% where leading edge Inthisconjunction,andwiththeprogressofCFD,theuse contaminationfrominsectsandairbornepollutantsiscom- ofCoanda˘effecttoenhancelifthasattractedrenewedinterest mon. Variable-pitch toward stall would result in simi- [4–7, 10–14, 17, 29]. Tangential jets that take advantage of lar roughness losses as fixed-pitch stall regulation, while Coanda˘ effect to closely follow the contour of the body are variable-pitch toward feather would decrease the loss to consideredtobesimpleandparticularlyeffectiveinthatthey around 10% at the expense of the rotor being susceptible canentrainalargemassofsurroundingair.Thiscanleadto to power spikes in turbulent high winds. For variable-rpm increasedcirculationinthecaseofairfoils,ordragreduction rotors operating at constant tip-speed ratio and angle of (ordragincreaseifdesired)inthecaseofbluffbodiessuchas attackdistribution,thelossisminimalataround5%to10%. anaircraftfuselage. Tominimizetheenergylossesduetoroughnesseffectsand InthedesignofanovelWindturbineusingwiththelatest todevelopspecial-purposeairfoilsforhorizontalaxisWind aerospace engineering technology, including aerodynamics, Turbine(HAWT’s),theNationalRenewableEnergyLabora- new lightweight materials, and efficient energy extraction tory (NREL), formerly the Solar Energy Research Institute technology, highly efficient aerodynamic surface came as a (SERI),andtheAirfoilsInc.beganajointairfoildevelopment keyissue.Recenteffortsinthelastdecadehavealsoindicated effortin1984.Resultsofthisefforthaveproducedanewseries that Coanda˘ effect has been given new and considerable ofairfoilswhichareconsideredmoreappropriateforHAWT considerations for circulation control technique [17, 19, 20, applications.Thesebladesareprovidingsignificantincreases 29].Takingadvantageoftheserecentresults,thepresentwork inannualenergyproductionasaresultoflesssensitivityto looksintotheefforttooptimizeaerodynamicperformanceof roughness effects, better lift-to-drag ratios, and, in the case WindturbineandtocriticallyinvestigatetheuseofCoanda˘ ofstall-regulatedrotors,throughtheuseofmoresweptarea effect. foragivengeneratorsize.Becauseoftheeconomicbenefits The progress of high speed computers exemplified by providedbytheseairfoils,theyarerecommendedforretrofit the availability of new generations of notebooks has made bladesandmostnewWindturbinedesigns[30]. possible the use of first-principles-based computational ap- Annual energy improvements from the NREL airfoil proaches for the aerodynamic modeling of Wind turbine familiesareprojectedtobe23%to35%forstall-regulatedtur- blades,tonameanexample.AspointedoutbyXuandSankar bines,8%to20%forvariable-pitchturbines,and8%to10% [17], since these approaches are based on the laws of con- for variable-rpm turbines. Optimizing airfoil’s performance servationofmass,momentum,andenergy,theycancapture characteristics for the appropriate Reynolds number and much of the physics in great detail. These approaches are thickness provides additional performance enhancement in particularly helpful at high wind speeds, where appreciable therangeof3%to5%[31].Becauseoftheeconomicbenefits regionsofseparationarepresentandtheflowisunsteady. providedbytheseairfoils,theyarerecommendedforretrofit In the design of circulation control technique, one may bladesandmostnewWindturbinedesigns[30].Inorderto note the progress of stall-controlled and pitch-controlled obtain results that can be well assessed, especially with turbinerotortechnologyfordealingwithpowerproduction reference to extensive research carried out as described in at low and high wind conditions [30]. Of interest is the [17,19,20,29],theairfoilS809[31]isconsideredinthepresent deployment of Wind turbines at larger range of sites not work.Forthispurpose,inthisparticularworkwefocusona limitedtoarelativelysmallnumberofsiteswherefavorable particular problem, that is, the basic airflow idealization conditions exist. At the other extreme, precautions should overatwo-dimensionalairfoilinsubsonicflow.Twospecific be taken when the wind speeds become very high, when conditionsareconsidered.ThefirstisacleanS809airfoiland JournalofRenewableEnergy 3 the second is the GTRI Dual Radius CCW airfoil that has velocitycontributedbytherotationoftheWindturbineblade alreadybeeninvestigatedthoroughly,forbenchmarking.The (at a selected baseline radial position, say 0.7𝑅, with 𝑅 the Coanda˘effectstudiesarecarriedoutonS809airfoil. radiusofthewindturbinerotor)andthefree-streamvelocity oftheambientair. 2.ProblemFormulation 3.Objectives Circulationcontrolwing(CCW)technologyisknowntobe beneficialinincreasingtheboundcirculationandhencethe To arrive at desired design configurations, one is lead to sectional lift coefficient of airfoil. This technology has been comeupwithdesiredlogicalcauseandeffectlawsaswellas extensivelyinvestigatedbothexperimentallyandnumerically tofindwaystocarryoutoptimizationschemes.Itiswithsuch [3–22,29–32]overmanyyears.Circulationcontrolisimple- objectivesinmindthatusingnumericalanalyses,parametric mentedbytangentiallyblowingasmallhigh-velocityjetover studies can be carried out and may offer some clues on ahighlycurvedsurface,suchasaroundedTE. relevantparameterswhichmaybeutilizedinamultivariable Thiscausestheboundarylayerandthejetsheettoremain optimization (and to a larger scale, multidisciplinary opti- attached along the curved surface due to the Coanda˘ effect mization). The present work will investigate the influence, (i.e., a balance of the radial pressure gradient and centrif- effectiveness,andconfigurationofCoanda˘-jetfittedaerody- ugal forces) causing the jet to turn without separation. The namicsurface,inparticularS809airfoil,forimprovingitslift rearstagnationpointlocationmovestowardthelowerairfoil augmentation and L/D, with a view on its incorporation in surface,producingadditionalincreaseincirculationaround thedesignofwindturbine.Inaddition,parametricstudyis the entire airfoil. The outer irrotational flow is also turned carried out for obtaining optimum configuration of the substantially,leadingtohighvalueofliftcoefficientcompa- Coanda˘effectliftenhancedairfoil.Thecomprehensiveresults rabletothatachievablefromconventionalhighliftsystems, asaddressedintheobjectiveoftheworkwillbeassessedto as illustrated in Figure1(b) [4]. Tongchitpakdee et al. [19, establishkeyfindingsthatwillcontributetothephysicalbasis 20]observedthatearlyCirculationControlledWingdesigns of Coanda˘-jet for lift enhancement in aerodynamic lifting typically had a large-radius rounded TE to maximize the surfaces, and an analysis based on physical considerations lift benefit. These designs, however, also result in high drag will be carried out with reference to their application in penalty when the jet is turned off. One way to reduce this propeller(orhelicopterrotor)andwindturbinebladeconfig- dragistomakethelowersurfaceoftheTEaflatsurface,while urationfordesignpurposes.Thefeasibilityandpracticability keepingtheuppersurfacehighlycurved.Thiscurvatureon ofCoanda˘configuredairfoilforwindturbineapplicationsare theuppersurfaceproducesalargejetturningangle,leading assessed. tohighlift.Inthisconjunction,theabilityofcirculationcon- trol technology to produce large values of lift advantageous 4.ComputationalModelingandSalient forWindturbinedesignsinceanyincreaseinthemagnitude of the lift force (while keeping drag small and L/D high) FeaturesofCFDComputationalScheme will immediately contribute to a corresponding increase in inducedthrustandtorquewillbetakenintoconsideration. 4.1. Governing Equations: Reynolds-Averaged Navier Stokes. Numerical simulation carried out by Tongchitpakdee Two-dimensionalincompressibleReynolds-averagedNavier- et al. [19, 20] looked into two approaches of introducing Stokes(RANS)equationwillbeutilizedtoperformtheanaly- Coanda˘-jet,thatis,attheappropriatelychosenpointinthe sisandtoobtainnumericalsolutionsonacomputationalgrid vicinity of the trailing as well as leading edge. A leading surrounding a reference airfoil. The governing equation is edgeblowingjetwasfoundtobehelpfulinincreasingpower givenby(steadystateandignoringbodyforces) generationatleadingedgeseparationcases,whileaTEblow- 𝜕𝑢 ingmayproducetheoppositeeffect. 𝜌𝑢 𝑖 Withsuchbackground,thepresentworkisaimedatthe ⏟⏟⏟⏟⏟⏟𝑗⏟𝜕⏟⏟𝑥⏟⏟⏟𝑗⏟ search for favorable Coanda˘-jet lift enhanced configuration Changeinmeanmomentumoffluidelement forwindturbinedesigns,focusingontwo-dimensionalsub- sonicflow,byperformingmeticulousanalysisandnumerical ⏞⏞⏞⏞⏞V⏞⏞⏞i⏞s⏞c⏞⏞o⏞⏞u⏞s⏞s⏞⏞t⏞r⏞e⏞⏞s⏞s⏞⏞⏞⏞⏞⏞ simulation using commercially available CFD code. Results = 𝜕 [[ −𝑝𝛿 +𝜇(𝜕𝑢𝑖 + 𝜕𝑢𝑗)− 𝜌𝑢𝑢 ]], obtainedfromrecentresearcheswillbeutilizedfordirecting 𝜕𝑥 [ ⏟⏟⏟⏟⏟⏟⏟⏟𝑖𝑗⏟ 𝜕𝑥 𝜕𝑥 ⏟⏟⏟⏟⏟𝑖⏟⏟⏟𝑗⏟ ] 𝑗 𝑗 𝑖 the search as well as for validation. For example, best lift [Normalstress Reynoldsstress] enhancing configuration [19], that is, lower part of the TE 𝜕𝑢 surfacetobeflat,willbeincorporatedinthepresentnumer- 𝑖 =0. icalinvestigation.Casestudiesperformedhereusebaseline 𝜕𝑥𝑖 chord dimension of 1m airfoil chord length. Two values of (1) free-streamvelocitiesareinvestigated,thatis,5.77m/secand 14.6m/sec,correspondingtoReynoldsnumberof3.95×105 Computationalfluiddynamicscodeswillbeusedtoana- 6 and 10 , which should fall within the typical situations for lyze the flow around two-dimensional airfoil. Various CFD large wind turbines. The free-stream velocity in the present codes are available and are alternatively utilized. Careful two-dimensional study is taken to be equal to the resultant reviewandanalysisoftheapplicationofthefirstprinciplein 4 JournalofRenewableEnergy Inboard: circulation control wing/thrust deflector (a) Possible leading edge blowing Tangential blowing over rounded trailing edge surface Pressure/centrifugal force balance 𝐶 𝐿 Slot Δ𝐶 /𝐶 >80 𝐿 𝜇 Circulation control Boundary layer control Jet sheet (b) Figure1:(a)Circulationcontrolwing/uppersurfaceblowingSTOLaircraftconfiguration,from[4];(b)Coanda˘-jettangentialblowingover TEsurfaceforliftenhancement,from[5]. the numerical computationare carried out prior to the uti- finenessshouldbechosentoobtainefficientsolution,which lizationofcommerciallyavailableCFDcodesforthepresent isacombinationofcomputationalspeedandaccuracy.The study.Forvalidationofthecomputationalprocedure,resort secondistheturbulencemodeling,whichisbasedonphysical ismadetoappropriateinviscidcase,andonlineaerodynamic modeling of real flow. Associated with turbulence model- codesareutilized.Foratwo-dimensionalgeometrythemesh ing, the wall boundaries are the most common boundaries generator partitions the subdomains into triangular or encountered in fluid flow problems. Turbulent flows are quadrilateralmeshelements[33].Careshouldbeexercisedin significantlyaffectedbythepresenceofwalls.Thesuccessful thechoiceandgenerationofthecomputationalgrid,which prediction of wall bounded turbulent flows relies on the will be elaborated subsequently. The comprehensive results accuracy with which the flow in the near-wall region is as addressed in the objective of the work are assessed to represented. Experiments have shown that the near-wall establishkeyfindingsthatwillcontributetothephysicalbasis region can be subdivided into several layers, characterized of Coanda˘-jet for lift enhancement in aerodynamic lifting bythevalueof𝑦+,whichincludetheviscoussublayerwhere surfaces. An analysis based on physical considerations is thelaminarpropertyisdominant,thefullyturbulentregion carriedoutwithreferencetoitsapplicationinwindturbine whichconsistsoftheturbulentlogarithmiclayerforarange bladeconfigurationfordesignpurposes. of𝑦+values,andtheouterlayer,asschematicallyexhibitedin The basis of the computational approach in the present Figure2(b). The no-slip condition at stationary walls forces work is the two-dimensional incompressible Reynolds- themeanvelocitycomponentstoazeromagnitudeandcan averaged Navier-Stokes (RANS) equation (1) and its imple- alsosignificantlyaffecttheturbulencequantities.Ifthegrid mentation into the CFD equation. When pursuing further distribution is fine enough so as to satisfy the no-slip con- the numerical solution of dynamic fluid flow using CFD dition,near-walltreatmentisnotnecessary[34–42]. techniques, in order to arrive at accurate solution, at least two aspects should be given careful attention, with due 4.2.TurbulenceModelingConsiderations. Turbulencemodel- considerationofcomputationaluncertainties.Thefirstisthe inginCFDisveryessentialinthechoiceofgridfinenessand gridgeneration,andthesecondistheturbulencemodeling obtainingthecorrectsimulationofparticularflowfield.For of the particular flow. The first is based on computational thispurpose,onehastochooseaparticularmodeloutofa algorithmincludingthechoiceofappropriategrid;thegrid hostofavailableturbulencemodelsdevelopedtodate.Oneof JournalofRenewableEnergy 5 0.05 ismodeledbyusingtheKomolgorov-Prandtlexpressionfor turbulentviscosity 𝑘2 0.04 𝜇𝜏 =𝜌C𝜇 𝜀 , (2) whereC𝜇isamodelconstant.Inexpression(2)𝜇𝜏wasbased 0.03 ontheeddyviscosityassumptionintroducedforconvenience 𝑐 ht/ of further analysis by Boussinesq [41] to draw analogy g ei betweenthemomentumtransfercausedbyturbulenteddies H 0.02 andtheviscosityinlaminarflow. The dimensionless distance in the boundary layer, sub- layer scaled, representing the viscous sublayer length scale, 0.01 plays significant role in capturing relevant physical turbu- lencephenomenaneartheairfoilsurfacecommensuratewith the grids utilized in the numerical computation. In this regards,theflowfieldinthevicinityoftheairfoilsurfaceis 0 0.8 1 1.2 1.4 1.6 1.8 2 usuallycharacterizedbythelawofthewall,whichattempts 𝑈/𝑈∞ toidentifyintricaterelationshipsbetweenvariousturbulence SA 𝑘-𝜀 model scaling in various sublayers. For example, the wall functions approach (wall functions were applied at the first SST Experimental node from the wall) utilized by k-𝜀 turbulence model uses Figure 2: Typical 𝑦+ values in the turbulent boundary layer and empiricallawstomodelthenear-wallregion[35]tocircum- velocityprofilesinwallunitsforsubsonicflowoverflatplate,M= vent the inability of the k-𝜀 model to predict a logarithmic 0.2,Re=1×106,mediumgrid;nozzlepressureratio=1.4.(Source: velocityprofilenearawall.Thelawofthewallischaracterized Menter[23],SalimandCheah[45],andEconomonandMilholenII byadimensionlessdistancefromthewalldefinedas [46]). 𝑢 𝑦 𝑦+ = 𝜏 . (3) 𝜐 Subject to local Reynolds number considerations, the wall theseturbulencemodelsthatisconsideredtobeappropriate 𝑦+ is often used in CFD to choose the mesh fineness anduserfriendlyisthek-𝜀turbulencemodel,withoutdisre- requirements in the numerical computation of a particular gardingothermodelsthatmaybesuitableforthepurposeof flow. thepresentwork. Basedonpreliminaryattemptstochoosetheappropriate The k-epsilon turbulence model was first proposed by turbulence model which yields desirable results, in the HarlowandNakayama[34]andfurtherdevelopedbyJones presentwork,k-𝜀turbulencemodelischosen[34,35].Inthe and Launder [40]. Continuing extensive research has been framework of eddy viscosity based turbulence models, the carriedouttounderstandthenatureofturbulence,andthis flow of a turbulent incompressible fluid is governed by the theoryhasbeenfurtherelaboratedandmanyothertheories Reynolds-averaged Navier-Stokes (RANS) equation (1) for havealsobeenintroducedsuchaselaboratedin[43],which thevelocityuandpressurep. seem to be satisfactory for only certain classes of cases. Forpracticalimplementationpurposes,itisworthwhile TheturbulentviscositymodelsbasedonReynolds-averaged tointroduceanauxiliaryparameter𝛾 = 𝜀/𝑘tolinearizethe Navier-Stokes(RANS)equations(1)arecommonlyemployed equations of the k-𝜀 turbulence model using the following inCFDcodesduetotheirrelativeaffordability[44].However, equivalentrepresentation: sincethechoiceofturbulencemodelandassociatedphysical 𝜕𝑘 ] phenomena addressed are relevant in modeling and com- 𝜕𝑡 +∇⋅(𝑘𝑢− 𝜎𝑇∇𝑘)+𝛾𝑘=𝑃𝑘, (4a) putersimulationoftheflowsituationneartheairfoilsurface 𝑘 consideredhere,acloserlookatturbulencemodelsutilized 𝜕𝜀 ] bytheCFDcodechosenwillbemade.Althoughthepractical 𝜕𝑡 +∇⋅(𝜀𝑢− 𝜎𝑇∇𝜀)+𝐶2𝛾𝑘=𝛾𝐶1𝑃𝑘. (4b) implementationofturbulencemodel,especially,forthenear- 𝑘 walltreatment,hasbeensomesomewhatofamystery[35], The corresponding linearization parameter is given by 𝛾 = the numerical implementation of turbulence models has a C𝜇(𝑘/]𝑇).Hencethestandardk-𝜀turbulencemodelisbased decisive influence on the quality of simulation results. In particular,apositivity-preservingdiscretizationofthetrou- on the assumption that ]𝑇 = C𝜇(𝑘2/𝜀), where k is the blesome convective terms is an important prerequisite for turbulent kinetic energy and 𝜀 is the dissipation rate. 𝑃𝑘 = the robustness of the numerical algorithm [35]. The k-𝜀 (]𝑇/2)|∇𝑢+∇𝑢𝑇|2 and 𝜀 are responsible for the production modelintroducestwoadditionaltransportequationsandtwo anddissipationofturbulentkineticenergy,respectively.The dependentvariables:theturbulentkineticenergy,𝑘,andthe default values of the involved empirical constants are C𝜇 = dissipation rate of turbulence energy, 𝜀. Turbulent viscosity 0.09,𝐶1 =1.44,𝐶2 =1.92,𝜎𝑘 =1.0,and 𝜎𝜀 =1.3. 6 JournalofRenewableEnergy FollowingiterativeprocedureelaboratedbyKuzminetal. the wall (airfoil) surface, practically the meshing layer of [35],withprescribedhomogenousNeumannboundarycon- width y is removed from the computational domain, thus ditionandimpermeablewall,itremainstoprescribethewall reducing the total number of grids involved, and hence shear stress as well as the boundary conditions for 𝑘 and 𝜀 reducing the computational time. Kuzmin et al. [35] have on the wall. The equations of the k-𝜀 model are invalid in introducedawallfunctionthatisabletocontrolthe𝑦+value thevicinityofthewall,wheretheReynoldsnumberisrather suchthattheydistancewillnotfallintotheviscoussublayer, low and viscous effects are dominant. Therefore, analytical whichcouldjeopardizethewallfunctionassumptionmade. solutions of the boundary layer equations are commonly The smallest distance of 𝑦+ that can be defined, following employedtobridgethegapbetweentheno-slipboundaries GrotjansandMenter[47],correspondstothepointwherethe and the region of turbulent flow. The friction velocity 𝑢𝜏 is logarithmiclayermeetstheviscoussublayer. assumedtosatisfythenonlinearequation The distance y is automatically computed iteratively by solving𝑦+ =(1/𝑘)log𝑦++𝛽=(1/0.41)log𝑦++5.2,sothat |𝑢| 1 = log𝑦++𝛽 (5) 𝑦+ = (𝑢𝜏𝑦)/𝜐 = 𝜌𝑢𝜏(𝑦/𝜇), where 𝑢𝜏 = 𝜌C𝜇(𝑘2/𝜀) (which 𝑢 𝜅 𝜏 hasbeenfurtherderivedbyKuzminetal.[35]tobeequivalent valid in the logarithmic layer, where the local Reynolds to𝑢𝜏 = C1𝜇/4√𝑘)isthefrictionvelocityandisequalto11.06. number This corresponds to the distance from the wall where the 𝑢 𝑦 logarithmiclayermeetstheviscoussublayer.Hence,thevalue 𝑦+ = 𝜏 (6) of𝑦,asrepresentedbythechoiceofthecomputationalgrid ] size, can never become smaller than half of the height of ]isintherange11.06≤ 𝑦+ ≤300.Theempiricalconstant𝛽 the boundary mesh cell [48], as illustrated in Figure3(a). equals 5.2 for smooth walls. Note that the above derivation This means that 𝑦+ can become higher than 11.06 if the resultsintheboundaryvalueoftheturbulenteddyviscosity meshisrelativelycoarse.Thek-𝜀variablesarenotspecified ]𝑇tobeproportionalto].Indeed, a priory, but during the computation using CFD software, k-𝜀variablevaluesareimplicitlysettoobtainacceptable𝑦+ 𝑘2 ] =𝐶 =𝜅𝑢 𝑦=𝜅𝑦+], (7) valuesthatproviderelativelyrapidconvergenceandaccuracy 𝑇 𝜇 𝜀 𝜏 ∗ (asvalidatedwithappropriatedata). Therefore, care is exercised in the choice of grids in the where 𝑦∗+ = |𝑢|/𝑢𝜏. Based on such information, the search vicinityoftheairfoil,toobtainacertainacceptableerrortol- forappropriategridfines,inadditiontotheoutcomeofgrid erance(whichmayalsobeattributabletonumericalerrorand finenessratiostudy,isalsoguidedbyplausiblecomparisonof uncertainties). The plausibility of the numerical results will thevalidationofCFDcomputationalresultsofthetestairfoil be judged by comparisonto other established results in the withexperimentalvalues.HeretheCFDcodeutilizedusing literature (numerical or experimental) for specific cases. In 𝑦+ value of 11 and as suggested in [35] was found to yield thepresentstudy,itwasfoundthattheturbulenceintensity satisfactoryresultsasindicatedthere.Figure2belowshows of5%andlengthscaleof0.01myieldresultsthatagreewith theinfluenceof𝑦+ inrepresentingtheboundarylayernear benchmarkingdata.ItisnotedthatHowelletal.[49]alsouse thewall. thisvalueandturbulenceintensityof1%intheirnumerical Various studies [23, 35, 42, 47] have indicated that studyonVAWT. integrating of the k-𝜀 type models through the near-wall regionandapplyingtheno-slipconditionyieldunsatisfactory results.Therefore,takingintoaccounttheneedforeffective 5.ComputationalSet-UpandValidation choice of grid mesh compatible with the turbulence model, adoptedusingthecommercialCFDcode,aparametricstudy 5.1.ComputationalGrid. Inthepresentwork,thegridswere iscarriedoutontwoairfoilswhereeitherexperimentaldata constructed following free mesh parameter grid generation or computational results are available, to validate the com- usinganalgebraicgridgeneratorandvariedfromextremely putationalprocedure,turbulencemodelandgridsizechosen coarseuptoextremelyfinegrids.However,fortheS809airfoil and establish their plausibility for further application in case, near the surface of the airfoil, boundary layer based the present work. The idea is to place the first computa- meshingiscarriedoutthroughout.Gridsensitivitystudyon tionalnodesoutsidetheviscoussublayerandmakesuitable the lift force per unit span for clean airfoil at zero angle of assumptionsabouthowthenear-wallvelocityprofilebehaves, attackisplottedinFigure4. in order to obtain the wall shear stress. Hence, the wall Figure4showsthattheliftanddragarevirtuallyinsen- functionscanbeusedtoprovidenear-wallboundarycondi- sitive of the size of the grid generated. The grid sensitivity tionsforthemomentumandturbulencetransportequations, studyhasledtoacceptablegridsizesintherangeoffineup rather than conditions at the wall itself, so that the viscous toextrafinegridsize,sincewithinthisrangetheliftaswell sublayerdoesnothavetoberesolvedandtheneedforavery asthedragforcesdonotchangeappreciably.Theassociated finemeshiscircumvented.ThewallfunctionsinCOMSOL rangeoffreemeshparameters(rangingfromextrafinetofine [48] are chosen such that the computational domain is meshinggridsize)isdefinedinTable1. assumedtostartatadistanceyfromthewall(seeFigure3). Inordertocontrolthe𝑦+value(setat11.06,ontheairfoil By applying the wall function at the nodes of the first surface),themaximumelementsizewassetat0.005,whileat meshinglayerofthecomputationalgridatadistanceyfrom theTErounding-offsurfacethemaximumelementsizewas JournalofRenewableEnergy 7 𝑈/𝑈𝜏 Buffer layer Fully turbulent region Viscous or or sublayer blending region log-law region 1.2(1.2𝑦) 1.2𝑦 𝑦 ln(𝑈𝜏(𝑦)/(cid:13)) (a) (b) Figure3:Forwallfunctions,thecomputationaldomainstartsadistanceyfromthewall.(a)GridgenerationfollowedfromCOMSOL4.2 User’sGuide;(b)thestructureofturbulenceboundarylayer,adaptedfromMenter[23]. m) 11.4 m) 2.02 N/ 11.3 N/ 2 Lift force per unit span ( 1111111100001......2119876 0.0002 0.00005 Drag force per unit span ( 11111....1.99998.986428 0.0002 0.00005 Extra... Coarser Coarse Normal Fine Finer Extra... Extra... Coarser Coarse Normal Fine Finer Extra... (a) (b) Figure4:GridsensitivitystudyforcleanS809airfoilatzerodegreeangleofattack;(a)liftforceperunitspan;(b)dragforceperunitspan. Table1:Therangeoffreemeshparameters. Table2:Computationalgridmeshinglayerproperties. Parameter Range Parameter Range Maximumelementsize 0.156–0.42 Numberofmeshlayersintheboundarylayers 8 Minimumelementsize 0.018–0.12 Boundarylayermeshexpansionfactor 1.2 Maximumelementgrowthrate 1.08–1.13 Thicknessoffirstmeshlayer 0.00005–0.0002 Resolutionofcurvature 0.25–0.3 Thicknessadjustmentfactor 1 Resolutionofnarrowregions 1 setat0.001.AttheCoanda˘-jetoutlet,thesurfaceisdivided 0.0002m, commensurate with the curvature on the airfoil into a minimum of 10 grid meshes. To obtain acceptable surface.Havingchosenthethicknessoftheelementsofthe computationalresultsusingCFDcode,appropriateboundary first meshing layer adjacent to the airfoil surface, the fol- layer mesh properties should be carefully chosen and the lowingconsecutivelayerswereexpandedwithanexpansion thinboundarylayersalongtheno-slipboundariesshouldbe factorof1.2,foreightconsecutivelayers.Anotheroptioncan resolvedintomeshlayerscommensuratewiththegridfine- be taken to allow thickness adjustment of the first meshing ness.Aboundarylayermeshisaquadrilateralmesh(2Dflow) layer, by choosing the thickness adjustment factor [48]; for withdenseelementdistributioninthenormaldirectionalong thepresentstudy,thethicknessadjustmentfactorofonewas specificboundariestoresolvethethinboundarylayersalong chosen. theno-slipboundaries.Theboundarylayerpropertieschosen Thegridgeneratorissufficientlygeneralsothatonecan aredefinedinTable2. easilyvarythejetslotlocationandsize.Gridspacingandclus- Alongwiththerequirementofthe𝑦+valuetocapturethe teringcanhavesignificanteffectsonWindturbineloadand salientturbulentflowfieldsublayers,arangeofvaluesofthe performancepredictions[49].Theouterboundaryisplaced grid dimension could be chosen to represent the thickness farawayfromthebladesurface,atleastatsixchordlengths of the first meshing sublayers. In the present example the (6C) from the airfoil surface, to avoid significant influ- grid dimension was chosen to be between 0.00005m and ence from outer boundary into the interior domain and to 8 JournalofRenewableEnergy 6C 6C 6C 6C Boundary layer mesh Figure5:Computationalgridaroundatypicalairfoil,shownhereforS809Coanda˘configuredairfoil. allow the disturbances to leave the domain [20]. The com- bejustifiedintheresultspresentedinFigure5.Itshouldbe putationaldomainandthegridinthevicinityoftheairfoil noted,however,thatFigure5indicatesthatwithhigherangle surfaceandCoanda˘-jetslotareexhibitedinFigure5. ofattack,far-fieldboundariesshouldbeplacedfurtheraway fromtheairfoil. 5.2. Initial and Boundary Conditions. In general, the initial conditions are set to be equal to the properties of the free- 5.3.BaselineValidation. Priortoitsutilization,variousbase- stream flow condition [19]. The flow properties everywhere line cases have been tried out, with favorable results. As a insidetheflowfieldareassumedtobeuniformandsettothe caseinpoint,forbenchmarkingpurposesthecodehasbeen free-streamvalues.Thefollowinginitialvaluesareused: appliedtocalculatethelift-slopecharacteristicsofS809airfoil (cleanconfiguration,withoutCoanda˘)andGTRICCWDual 𝑢=𝑢 , V=V , 𝜌=𝜌 , 𝑝=𝑝 , ∞ ∞ ∞ ∞ Radius airfoil (with Coanda˘) and compared to the results (8) 𝑇=𝑇 . obtained by Somers [50], as shown in Figure6, and Englar ∞ et al. [5], as shown in Figure7, respectively. The results as Theouterboundaryisplacedfarawayfromthebladesur- shown in both Figures demonstrated the plausibility of the face,ataminimumof6chords(6C).Nonreflectingboundary presentnumericalcomputationalscheme. conditionsareappliedattheouterboundariesofthecompu- Computational results exhibited in Figure6 served to tational domain. The jet is set to be tangential to the blade validatethecomputationalprocedureassociatedwiththeuse surface at the Coanda˘-jet nozzle location. The jet velocity of CFD code in the numerical simulation and as a baseline profileisspecifiedtobeuniformatthejetexit.Ontheairfoil in furthering the computational study of Coanda˘-jet in the surface,exceptatthejetexit,no-slipboundaryconditionsare presentstudy.Figure6comparesthepresentcomputational applied.Togainbenefits,thejetvelocityisheredesignedto simulation using k-𝜀 turbulence model and Kuzmin’s [35] belargerthanthepotentialflowvelocityatthevicinityofthe wall function (the associated 𝑦+ value is 11.06) with the outeredgeoftheboundarylayer.Inaddition,thethicknessof experimental data of Somers [50]. The numerical results theCoanda˘-jetisdesignedtobelessthanthelocalboundary above simulated similar experimental conditions of Somers layerthickness. for clean S809 airfoil (with wall boundary conditions and Outer boundary conditions are chosen to minimize without Coanda˘-jet). It also shows that the lift and drag blockage effect, that is, significant discrepancies with free coefficientsatzeroangleofattackareincloseagreementwith flightsituations,whethertheoreticallyornumericallysimu- theexperimentaldata. latedorexperimental.Walleffectsblockagechangesstream- Similar numerical simulation was carried out for GTRI linecurvatureinteractionboundarylayersandtheseshould Dual Radius CCW airfoil with leading edge blowing and be minimized. The best measure is by carrying out com- comparedtotheexperimentaldata(inFigure7)fromEnglar putational experiments with reasonable outer boundaries et al. [5] with the same boundary conditions and jet con- of the computational grid and validating with both outer figuration using k-𝜀 turbulence model and Kuzmin’s wall potentialflowresultsandnumericalparametricstudies.This function[35].Ofsignificantinterestarethegridsize,thewall rationale has shown that, for the case considered, an outer function,and𝑦+valueutilizedfornumericalsimulation.The gridboundarydistancetothecomputationaldomaincenter- associated𝑦valueforthecellsadjacenttotheairfoilsurface lineofsixchordsproducesresultswithhighaccuracy,ascan wasdesignedtobenomorethan0.0002minordertosatisfy JournalofRenewableEnergy 9 0.05 1.2 0.045 1.1 1 0.04 0.9 nt 0.035 ent 0.8 cie 0.03 ci 0.7 ffi coeffi 00..65 g coe 00.0.0225 Lift 00..43 S809 Dra 0.015 0.01 0.2 0.1 0.005 0 0 −0.1−1 0 1 2 3 4 5 6 7 8 9 10 −1 0 1 2 3 4 5 6 7 8 9 10 Angle of attack Angle of attack Experiment, in closed conduit Experiment, in closed conduit COMSOL𝑘-𝜀, in closed conduit COMSOL𝑘-𝜀, in closed conduit COMSOL𝑘-𝜀, in infinite medium COMSOL𝑘-𝜀, in infinite medium (a) (b) Figure6:ComparisonforvalidationofS809airfoilofCFDcomputationalresultsusingCOMSOLk-𝜀turbulencemodelandexperimental valuesfromSomers[50]atRe=1E6;(a)liftcoefficient;(b)dragcoefficient. 5 appropriateCFDnumericalprocedures.Simplifiedtheoreti- calanalysisiscarriedoutpriortothemodelingandnumerical 4 nt computation,whilefurthertheoreticalanalysiswillbecarried e ffici 3 outinevaluatingthecomputationalresults. e co 2 ft Li 1 GTRI Dual Radius CCW 6.ResultsandDiscussions 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 6.1.DesignConsiderationsforCoanda˘ConfiguredAirfoil. For Momentum jet the purpose of assessing the influence and the effectiveness oftheCoanda˘enhancedL/Dandliftaugmentationonwind Experiment, in closed conduit turbine blade, a generic two-dimensional study is carried COMSOL𝑘-𝜀, in closed conduit out. The problem at issue is how the Coanda˘-jet can be COMSOL𝑘-𝜀, in infinite medium introduced at the TE of the airfoil, bearing in mind that Figure 7: Comparison of lift-curve slope for GTRI Dual Radius such design may recover any losses due to the possible CCWairfoilwithLEblowingonCFDcomputationusingCOMSOL inception of flow separation there. In addition, for effective k-𝜀turbulencemodelandexperimentalvaluesfromEnglaretal.[5] Coanda˘-jetperformance,acurvatureshouldbeintroduced. atRe=395,000. Furthermore,thethicknessoftheCoanda˘-jetasintroduced ontheairfoilsurfacecouldhaveaverycriticaleffectonthe intended lift enhancement function. For best effect, the lower surface near the TE should be flat, as suggested by the𝑦+valuerequirements,whichwerechosentobe11.06.The Tongchitpakdee et al. [19, 20]. The design of the Coanda˘- resultsexhibitedinFigure7forthiscaseisalsoencouraging configuredTEshouldalsoconsidertheoff-designconditions. duetothecloseagreementbetweenthepresentcomputation With all these considerations, a configuration suggested is andtheexperimentaldata. exhibitedinFigure8. The results exhibited in both examples above serve to indicate that the computational procedure and choice of turbulencemodelseemtobesatisfactoryforthepresentcom- 6.2.ComputationalResultsforS809Airfoil. Next,wewould putationalstudyandcouldlendsupporttofurtheruseofthe liketoinvestigatetheinfluenceofspecificallydesignedairfoil approachinthenumericalparametricstudy. geometryforWindturbineapplication,andforthispurpose Allcomputationsinthepresentstudywereperformedon atypicalS809,incleanandCoanda˘-jetequippedconfigura- alaptopcomputerwitha2.10GHzIntelCoreDuoprocessor, tions. S809 airfoil represents one of a new series of airfoils 4GBofRAM,and32-bitoperatingsystem.Typicalcompu- whicharespecificallydesignedforHAWTapplications[20]. tation time for the computation of the flow characteristics Forthepresentstudy,thenumericalsimulationwascarried around a two dimensional airfoil is in the order of 4 hours out at two free-stream velocities. These are 5.77m/s (corre- witharound300,000degreesoffreedombyusingstationary spondingtoRe=3.95×105)and14.6m/sec(correspondingto segregatedsolverink-𝜀turbulentmodelanalysis.Withsuch Re=1×106),whichrepresentlowandhighfree-streamcases, computationalenvironment,whichinconsequencelimitsthe respectively,whilethechord-lengthismaintainedatc=1m, computational grids, the present work will seek and adopt and density at 𝜌 = 1.225kg/m3. The baseline for assessing 10 JournalofRenewableEnergy beenappliedtotheproductionaircrafts.Manyoftheroad- TE jet thickness blocks have been associated with the engine bleed require- mentsandcruisepenaltiesassociatedwithblownbluntTEs. In addition, there is a tradeoff between the use of a larger Parallel Tangent line radius Coanda˘ configured airfoil for maximum lift and the useofasmallerradiusoneforminimumcruisedrag. Distance from LE TE radius In contrast to the needs of TE rounding-off radius, performancedegradationassociatedwithitalwaysstandsas an issue due to the drag penalty when the jet is in the off Figure8:TEconstructionoftheCoanda˘configuredairfoil(thejet mode. To overcome such drawback, the TE radii should be flowistangentialtotheroundedcircularsector). specificallyandcarefullydesigned.Forthatpurpose,simula- tionsatseveralTEradius(from10mmto50mm)havebeen performed,atafixedCoanda˘-jetmomentumcoefficient𝐶𝜇 (𝐶𝜇 = 0.003, considered just enough to fit with the Wind theadvantagesofCoanda˘-jetfromparametricstudyonS809 turbine application), and at a constant free-stream velocity airfoil is the computational result for the clean airfoil. The of 𝑉∞ = 14.6m/sec (Re = 1 × 106), to investigate the effect computationalresultforthisbaselinecasehasbeenvalidated of TE radius on the aerodynamic characteristics of Coanda˘ bycomparisonofthecomputationalvalueforthesameS809 configuredairfoils.ResultsexhibitedinFigure10showthata airfoiltotheexperimentalresultsbasedonwind-tunneltest. higher𝐿/𝐷canbeachievedwithasmallerTEradius(30mm) Itshouldbenotedthattheflowfieldboundariesaredifferent andthatthe𝐿/𝐷isdecreasingastheTEradiusisincreased fromthoseutilizedfortheparametricstudy. from 30mm to 50mm. The effect of TE radius on the lift Thetwo-dimensionalnumericalsimulationstudyforthe augmentationissignificant,asexhibitedbythedashedlinein S809 airfoil is carried out in logical and progressive steps. Figure10.However,whentheTEradiusisincreasedbeyond First,thenumericalcomputationisperformedontheclean certain value (in Figure10, ≫35mm), the TE rounding-off S809 airfoil, then on the Coanda˘ configured S809 airfoil seemstobeineffective,evendetrimental. withoutthejet(i.e.,afterappropriatemodificationduetoTE Moving the location of the Coanda˘-jet forward implies rounding-off and back-step geometry), and then finally on theincreaseoftheTEradius.Consequently,thecamberofthe Coanda˘configuredS809airfoilinitsoperationalconfigura- airfoilwillalsobechanged.Thisinturnwillproducechanges tion. intheangleofattackoftheairfoilcomparedtothebaseline To address three-dimensional wind turbine configura- case. Therefore, although the present study looks into the tions,particularlyfortheoptimumdesignofahorizontalaxis influenceofmodifyingTEradius,whileholdingotherairfoil Windturbine(HAWT),logicaladaptationshouldbemade, geometrical parameters constant, the zero angle of attack takingintoaccountthefactthatdifferentairfoilprofilesmay conditionisatbestpossibleonlyasfirstapproximation.The beemployedatvariousradialsections.Certainassumptions robustness of the computational procedure adopted in the havetobemadeinordertoprojecttwo-dimensionalsimu- presentcomputationalset-upmaybeabletotakeintoaccount lation results to the three-dimensional case, which may be all these changes. Such a case is indeed revealed in Figures necessarytoevaluatetheequivalentBetzlimit. 10(c)and10(d)wherebytheinfluenceofallthechangesofthe geometricalvariablesaccompanyingthepresentparametric 6.2.1.TheInfluenceofCoanda˘-JetontheFlowofftheTE. The study should have been incorporated in the computational flowfieldinthevicinityoftheTEforbothconfigurationsis scheme which is focused on the TE radius changes and showninFigures9(a)and9(b).Carefulinspectionofthese produceschangesinLandDcontributedalsobythechanges figuresmayleadtotheidentificationofthegeometryofthe intheangleofattack.Additionally,thesefiguresshowthatthe flowthatcouldcontributetoincreasedlift,insimilarfashion changes in lift is much larger than the changes in drag, as asthatcontributedbyflap,jetflap,orGurneyflap.Figure9(b) exhibitedinFigures10(c)and10(d). typifiestheflowfieldaroundCoanda˘configuredS809airfoil Variation of the Coanda˘-jet thickness from 0.5mm to without Coanda˘-jet (only with its back-step configuration), 3.0mm at a fixed 𝐶𝜇 (𝐶𝜇 = 0.005) and at a constant free- whichisusedheretogetinsighttotheactionoftheCoanda˘- stream velocity of 𝑉∞ = 14.6m/sec (Re = 1 × 106) is per- jet. formedtoinvestigatetheeffectofCoanda˘-jetthickness(also calledjetslotthickness)ontheaerodynamiccharacteristics 6.2.2. Parametric Study Results for Coanda˘ Configured S809 of Coanda˘ configured airfoils. From Figure11, it is found Airfoil. TheTEradiusplaysanimportantroleintheCoanda˘ thatahigherL/DcanbeachievedwithasmallerCoanda˘-jet configureddesignairfoil,sinceitmaypositivelyornegatively thickness (1.0mm) and that the L/D is decreased rapidly as influencethedownstreamflowbehavior.Tongchitpakdeeet theCoanda˘-jetthicknessisincreasedfrom1.5mmto3.0mm. al.[19,20]hadclaimedthatthelowersurfaceattheTEofthe A similar behavior is observed for the lift augmentation as applied Coanda˘-jet should be flat in order to minimize the exhibitedbythedashedlineinFigure11.However,generating drag when the jet is turned off. Also, as reported earlier by asmallerjetrequireshigherpressurethanalargeroneatthe AbramsonandRogers[51,52]inthelate1980s,inspiteofthe same momentum coefficient. Since higher lift with as low abilitytogeneratemoreliftthetechniquehasnotingeneral mass flow rate as possible is preferred, a thin jet is more
Description: