Development of a BEM-CFD tool for Vertical Axis Wind Turbines based on the Actuator Disk Model

Available online at www.sciencedirect.com Available online at www.sciencedirect.com Available online at www.sciencedirect.com Available online at www.sciencedirect.com Availableonline onlineatatwww.sciencedirect.com www.sciencedirect.com Available

ScienceDirect ScienceDirect ScienceDirect ScienceDirect

Energy Procedia 106 (00 2016 ) 1 000–000 Energy Procedia (2018) Energy Procedia 106 (00 2016 ) 1010–1017 1000–000 Energy Procedia (2018) 000–000 Energy Procedia 148 Energy Procedia 00(2018) (2017)

www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia

73rd Conference of the Italian Engineering Association (ATI 2018), 12-14 The 1stst edition of the EnergyThermal EconomicsMachines Iberian Conference took place on 4-5 February 73rd Conference of the Italian Machines Engineering Association (ATI 2018), 12-14 The 1 edition of the EnergyThermal Economics Iberian Conference took place on 4-5 February September 2018, Pisa, Italy 2016 in Lisbon, Portugal. The event was organized by the Italy Portuguese Association for September 2018, Pisa, 2016 in Lisbon, Portugal. The event was organized by the Portuguese Association for

Energy Economics (APEEN) in close collaboration with its Spanish counterpart (AEEE). Development of aa(APEEN) BEM-CFD tool for Vertical Axis Wind Turbines Energy Economics in close collaboration with its Spanish counterpart (AEEE). Both are respectively, the Portuguese andtool the Spanish affiliates of IAEEAxis – International Development of BEM-CFD for Vertical Wind Both areThe respectively, the PortugueseSymposium and the Spanish IAEE – International 15th International onaffiliates DistrictofHeating and CoolingTurbines Association for Energy Economics. based on the Actuator Disk model Association for Energy Economics. based on the Actuator Disk model The Conference had 155 participants from 17 and 73 papers were presented, a,∗ b countries a b Assessing the155feasibility using the demand-outdoor The Conference had participants fromof 17 and 73 heat papers presented, B.ofRocchio , S. Deluca , countries M. V. Salvetti Zanforlin a , S. were b from a total 90 sent toa,∗be evaluated byb ,the Scientific Committee, after the initial B. Rocchio , S. Deluca M. V. Salvetti , S. Zanforlin from a total of 90and sent to beEngineering, evaluatedUniversity by the Scientific Committee, after, Pisa the (PI) initial Department of Civil Industrial of Pisa, Largo Lucio Lazzarino 56122, Itayly submission offunction 142 abstracts. for From a those 73 presented papers, 22 were selected to be temperature long-term district heat demand forecast Department ofofCivil and Industrial Engineering, of Pisa, LargoofLucio Lazzarino , Pisa (PI) 56122, 142Territory abstracts. From thoseUniversity 73 presented papers, 22 were to be Departmentsubmission of Energy, Systems, and Constructions Engineering, University Pisa, Largoselected Lucio Lazzarino, PisaItayly (PI) 56122, Italy a a b b Department

published in this issue of Energy Procedia.Engineering, University of Pisa, Largo Lucio Lazzarino, Pisa (PI) 56122, Italy of Energy, Systems, Territory and Constructions published in this issue of Energy Procedia. a a,b,c a b c I.The Andrić *, A. market Pina , experience P. Ferrãomay , J. be Fournier B. for Lacarrière , O. Le Correc Iberian energy valuable ., both other European The Iberian energy market experience may be valuable both for other European countries but also worldwide, namely for -Latin American countries linked to a IN+ Center for Innovation, Technology and Policy Research Instituto Superior Técnico, Av. closely Rovisco Pais 1, 1049-001 Lisbon, Portugal Abstract but balso worldwide, namely for Latin American countries closely linked to Abstract countries Portugal and Spain culture,&language and VeoliabyRecherche Innovation, 291business. Avenue Dreyfous Daniel, 78520 Limay, France Portugal andSystèmes Spain by culture, language and business. c work focuses The present on the numerical simulation of Vertical Wind 4Turbines means44300 of anNantes, “in-house” BEM-based Département Énergétiques et Environnement - IMTAxis Atlantique, rue Alfredby Kastler, France TheDefined present work focuses the numerical simulation of Vertical Axis Typical Wind Turbines means of an3D “in-house” BEM-based Both under security of supply perspective andFluent. under sustainability, the Iberian case User Function tothe beon used 39ith the CFD code ANSYS VAWT by unsteady and phenomena, such as Both under the security of supply perspective and under sustainability, the Iberian case User Defined Function to be used 39ith the CFD code ANSYS Fluent. Typical VAWT unsteady and 3D phenomena, such as dynamic stall,can flow andatip losses, are taken into an account by island originalinand literature-based sub-models. The presence of becurvature viewed as case study. It is still energy Europe which increases dynamic flow curvature andatip losses, aresuitable taken into account by island originalinand literature-based sub-models. The presence of be viewed as case study. It is still an energy Europe which increases the bladesstall, is can mimicked by replacing them with momentum sources. costs to consumers, reduces competitiveness and creates vulnerability in terms of the blades mimicked bythe replacing them withcompetitiveness suitable momentum sources. costs tosupply. consumers, reduces creates vulnerability terms of For the is present work, Actuator Cylinder Model hasconsumers been and employed. 3Dhigher analysis, of aininSANDIA rotor, are carried out in Abstract energy The energy bill of Iberian is 40% than France. A For the present work, the Actuator Cylinder Model has been employed. 3Dhigher analysis, of ainThe SANDIA rotor, are carried out in energy supply. The energy bill of Iberian consumers is 40% than France. A order to assess the accuracy of our model against numerical simulations and experimental data. current User Function real Iberian gas corridor isagainst a particularly urgent infrastructure needed – but not theUser Defined order to assess the accuracy of our model numerical simulations and experimental data. The current Defined Function is able to give a satisfactory agreement with the reference cases especially from a qualitative point of view, with a significant real Iberian is aaddressed particularly urgent infrastructure needed but notsolutions the District heating networks arecorridor commonly in the literature as one of the most –effective for decreasing the only one - togas correct this situation. is able to give a satisfactory agreement with referencetocases especially from a qualitative point of view, a significant computational time reduction to athethis factor of 10the compared the case with the bodies.which A typical wakewith behavior be one - to correct situation. greenhouse only gas emissions from building sector. These systems require highmoving investments are returned through can the heat

computational time reduction tothough a factor 10 compared to the case with bodies. A inherent typical wake can be noticed in ourtosimulations evenclimate its of recovery is strongly by the the moving lack ofheat turbulence thebehavior chosen sales. Due the changed conditions and buildinginfluenced renovation policies, demand in the to future could actuator decrease, Although Spain´s dependence on energy imports was still 70% in 2014, there was an noticed in our simulations even though its recovery is strongly influenced by the lack of turbulence inherent to the chosen actuator model. The torque and the power theenergy turbines are alsowas analyzed and compared against the reference cases, finding a prolonging the investment returncoefficient period. ofon Although Spain´s stillPortugal 70% in 2014, there an impressive reduction of 10%ofcompared toare2009. Also reduced its was energy model. The agreement. torque and the powerdependence coefficient the turbinesimports also analyzed and compared against the reference cases, finding a remarkable The main scope of this paper is to assess the feasibility ofto using the heat demand – outdoor temperature function for heat demand impressive reduction of 10% compared 2009. Also Portugal reduced its energy remarkable agreement. dependency from 90% in 2005 to 71% in 2014. In both countries, incentives to The model hasdistrict been successfully applied to predict the transient aerodynamic of anstudy. offshore 5MW troposkein turbineofsubforecast. The of Alvalade, located in Lisbon (Portugal), was used loads as a case The district is consisted 665 dependency from generation 90% in 2005 to crucial 71%transient into 2014. In both countries, incentives The to model has beenmotion successfully applied toThe predict the aerodynamic of in an offshore 5MWto turbine subrenewable energy were this performance. In 2015 real energy jected thethat pitching of its platform. operating conditions have beenloads chosen order tomedium, allow a troposkein qualitative comparison buildings vary in both construction period and typology. Three weather scenarios (low, high) and three district renewable energy generation were crucial to this performance. In in 2015 energy jected the pitching motion of its platform. The operating conditions have beenmotion chosen orderreal to allow a qualitative comparison generation from renewable sources represented respectively, 48,1% from electricity with a to floating 5MW horizontal axis turbine which performance under pitching available in literature. renovation scenarios were developed (shallow, intermediate, deep). To estimate theiserror, obtained heat demand values were generation from renewable sources represented respectively, 48,1% from electricity with a floating 5MW horizontal axis turbine which performance under pitching motion is available in literature. gross generation plus importing balance in Portugal and 38,4% inand Spain. compared with results from a dynamic heat demand model, previously developed validated by the authors. gross generation plus importing balance in Portugal and 38,4% in Spain. c  2018 The Authors. Published by Elsevier Ltd. The results showed that when only weatherLtd. change is considered, the margin of error could be acceptable for some applications © 2018 The Authors. Published by Elsevier c 2018  The Notwithstanding, Authors. Published by Elsevier Ltd.IEA thelicense This an access article the BY-NC-ND (https://creativecommons.org/licenses/by-nc-nd/4.0/) according to power generation remains the largest This iserror an open open accessdemand article under under the CC CCthan BY-NC-ND (https://creativecommons.org/licenses/by-nc-nd/4.0/) (theis in annual was lower 20% for license all weather scenariossector considered). However, after introducing renovation This is an and open access article under the CC BY-NC-ND license (https://creativecommons.org/licenses/by-nc-nd/4.0/) Notwithstanding, according to IEA the power generation sector remains the largest Selection peer-review under responsibility of the scientific committee of the the 73rd Conference of the the Italian Italian Thermalconsidered). Machines of 73rd Conference of Thermal Machines CO2 emitter Portugalupwith 15.8 MtCO2 in 2013 orweather 35.2% of total CO2 emissions. In scenarios, the error value in increased to 59.5% (depending on the and renovation scenarios combination Selection andAssociation peer-review under responsibility of the scientific committee of theof73rd Conference of the Italian Thermal Machines CO2 emitter in Portugal with 15.8 MtCO2 in 2013 or 35.2% total CO2 emissions. In Engineering (ATI 2018). The value of increased average within range CO2 of 3.8% up to 8%87 perMtCO2 decade,inthat corresponds to the theslope casecoefficient of (ATI Spain, the energyonsector is also thethe largest emitter with Engineering Association 2018). the case ofpower Spain, the hours energy sector is during also the largest CO2 emitter with 87onMtCO2 in decrease in the number of heating of 22-139h the heating season (depending the combination of weather and 2013. The generation sector accounted for 28.4% of total CO2 emissions, while Keywords: VAWT; Actuator Disk Model; Dynamic Stall; Flow curvature; Wake analysis. 2013. The power generation sector accounted for 28.4% of total CO2 emissions, while renovationVAWT; scenarios considered). On the other hand, function intercept increased for 7.8-12.7% per decade (depending on the Actuator Disk Model; Dynamic Stall; Flow curvature; Wake analysis. Keywords: other energy industries (including oil refining) accounted for a further 7.9%. other energy industries (including a further 7.9%. for the scenarios considered, and coupled scenarios). The values suggested couldoilberefining) used to accounted modify the for function parameters improve theTo accuracy of heat demand estimations. attain CO2 reduction goals, both Portugal and Spain have two alternatives – in fact,

To CO2 reduction goals,–both Portugal Spain have two alternatives – in fact, twoattain complementary policies which must and be implemented wisely but broadly: to

Corresponding author. Tel.: +39-050-2217286 ©∗∗2017 Thetwo Authors. Published by policies Elsevier Ltd. complementary – which implemented but broadly: to improve energy efficiency gains and must apply be a CO2 tax whichwisely revenues should be Corresponding author. Tel.: +39-050-2217286 E-mail address: [email protected] Peer-review under responsibility of the Scientific Committee of The 15th International Symposium on District Heating and improve energy efficiency gains and apply a CO2 tax which revenues should be allocated to energy related innovation. E-mail address: [email protected] Cooling. allocated to energy related innovation.

Keywords: Heat demand; Forecast; Climate change

1876-6102 2017The TheAuthors. Authors.Published PublishedbybyElsevier ElsevierLtd. Ltd. c©2018 1876-6102 1876-6102  © 2018 The TheAuthors. Authors. Published by Elsevier Ltd. c under 1876-6102 2018 Published byIEA Elsevier Ltd. Peer-review responsibility of theCC Scientific Committee of The2016 15th review. International Symposium This is an open access article under the BY-NC-ND license (https://creativecommons.org/licenses/by-nc-nd/4.0/) IEA (2016), Energy Policies of countries: Portugal OECD/IEA, 2016on District Heating and Cooling. This is an open access article under the CC BY-NC-ND license (https://creativecommons.org/licenses/by-nc-nd/4.0/) This is an open access article under the CC BY-NC-ND license (https://creativecommons.org/licenses/by-nc-nd/4.0/) IEA (2016), Energy Policies of IEA countries: Portugal 2016 review. OECD/IEA, 2016 Selection responsibility of the scientific committee of theof73rd of the Italian MachinesMachines Engineering Selection and andpeer-review peer-reviewunder under responsibility of the scientific committee the Conference 73rd Conference of the Thermal Italian Thermal Selection and peer-review under responsibility of the scientific committee of the 73rd Conference of the Italian Thermal Machines Engineering Association (ATI 2018). Engineering Association (ATI 2018). 1876-6102 © 2016 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license Association 2018).by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license 10.1016/j.egypro.2018.08.060 1876-6102 © 2016 (ATI Published http://creativecommons.org/licenses/by-nc-nd/4.0/). http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review by the scientific conference committee of EEIC | CIEE 2016 under responsibility of Guest Editors. Peer-review by the scientific conference committee of EEIC | CIEE 2016 under responsibility of Guest Editors. doi:10.1016/j.egypro.2016.12.099 doi:10.1016/j.egypro.2016.12.099

2

B. Rocchio et al. / Energy Procedia 148 (2018) 1010–1017 Author name / Energy Procedia 00 (2018) 000–000

1011

1. Introduction Vertical Axis Wind Turbines (VAWTs) are based upon the design patented in France in 1925 and in the US in 1931 by aeronautical French engineer G. J. M. Darrieus [1]. Historically, Horizontal Axis Wind Turbines (HAWTs) have been preferred for power generation mainly because of their higher efficiency and less fatigue-related issues. Nonetheless, in recent years VAWTs have proved to be more suitable than HAWTs for small-scale urban applications, thanks to lower noise and vibration levels, ability to work with skewed flows [2, 3, 4, 5, 6, 7], insensitivity to the wind direction and, therefore, absence of any yaw device. In addition, VAWTs are gaining growing interest for large-scale offshore floating applications because of their higher position stability that can help reduce the platform costs and more easily accessible drive train and generator [8]. Therefore, the investigation of the aerodynamic performances of VAWTs is crucial to completely understand this application. To the authors knowledge, up to now the only 3D CFD studies concerning the aerodynamics of VAWTs on floating platforms are those recently carried out by Lei and co-workers on the effects of periodic pitching [9] and surge [10] motions on the performance of a H-rotor with a diameter of 2m. Their results indicate both pitching and surge can widen the variation ranges of the aerodynamics forces and change the flow field around the rotor with respect to a fixed turbine. In particular, pitching can improve the power output, while surge is able to keep the floating wind turbines steadier. To the author knowledge, the only tool specifically proposed for floating VAWTs is the one recently implemented by Cheng and co-workers [11]. Given this rising importance, several mathematical models have been developed for VAWTs investigation based on different approaches and implying very different computational costs. Simply and fast models, i.e. Blade Element Momentum (BEM) theory, provide fast calculations when applied to the design of a conventional low-solidity and low-loaded turbine. To analyze fluid dynamic interactions between turbines and the evolution of turbine wakes, more sophisticated and computationally expensive model can be employed. Studies based on both Unsteady Reynolds Averaged Navier-Stokes Simulations (URANS) [12, 13] and Large Eddy Simulations (LES) [14] were shown to allow a deep insight in important physical mechanisms such as dynamic stall and tip losses, which govern the torque generation and wake evolution of single VAWTs. The numerical representation of the blade geometry would lead to simulation difficult to afford. Therefore, simplified models are need to avoid this obstacle without losing the solution accuracy, in terms of power production and wake evolution. To reduce the computational costs of these simulations, the Actuator Disk (AD) and Actuator Line (AL) models can be employed, in which the blade action on the fluid is modeled by introducing distributed volume forces. Hence, the solution of the boundary layer around the blades is avoided and simpler meshes with a significantly reduced number of elements can be employed. When using a rotating actuator disk/cylinder model, the aerodynamic forces are spread on the entire surface swept by the blades. Both RAD/AL URANS and LES models were successfully implemented for analyzing wind and tidal farms of HAWTs. However, Actuator models for VAWTs imply some complications with respect to HAWTs, i.e. the need for a more refined dynamic stall model, a flow curvature model, and the interaction of the blade with the wake shed by itself or other blades[15, 16, 17, 18, 19]. It is therefore proposed a function for ANSYS Fluent which employs different sub-models for the description of unsteady and 3D VAWT phenomena. A preliminary validation was carried out against a 3D URANS simulation of the SANDIA 17m Darrieus turbine. Finally, the present model is applied for the prediction of a floating wind turbine.

2. Physical modeling The flow around wind turbines is considered incompressible. A URANS approach to turbulence is used, together with the k − ω SST model to close the problem. The final goal of the present work is to develop a model, which is able to reproduce the effect of the rotating blades. The aerodynamic blade forces calculation requires the development of several sub-models, i.e. an “in-house” dynamic stall model and a virtual profile model for curved flow effects, based on the work by Migliore [20]. To ensure the real physical trend of the forces along the blade span, the tip loss coefficient is also implemented [21, 22].

B. Rocchio et al. / Energy Procedia 148 (2018) 1010–1017 Author name / Energy Procedia 00 (2018) 000–000

1012

3

2.1. Rotor model To reduce the computational costs of the simulations, The idea proposed by Sørensen [23] is to mimic the presence of the blades by replacing the effect that they have on the flow-field. Thus, the following boundary condition on the blades surface: Vrel = 0

(1)

is neglected and the body forces are considered in the URANS equations as source terms Fi . The blade is numerically discretized in finite elements and, for each, the forces are evaluated from the 2D airfoil coefficients. To avoid numerical instabilities, the source term has to be spread among the several grid-points. The spreading is done according to a Gaussian regularization kernel, whose standard deviation is the smoothing parameter. The Actuator Disk model (ADM) has been employed for the present work. Hence, the total force F is spread on the whole rotor area and on several cells in radial direction. The Gaussian kernel can be written as: 

η=e

 θ−θ¯ 2 2π/N

e



 r−R 2 r

(2)

where R is the turbine radius, θ¯ gives the instantaneous azimuthal position of the blade, r is the standard deviation of the Gaussian kernel along the radial direction, N is the number of the blades and P = (r, θ) is a generic grid-point. The spreading along the azimuthal direction depends on the number of turbine blades. Therefore, the forces are evaluated as follows: Fη fturb =   ηdA A

(3)

where A is the annular ring over which the force is spread. 2.2. Dynamic Stall model This phenomenon is characterized by a huge increment of the lift coefficient, also after the static stall angle, α st where α is the angle of attack. The reason behind this behavior is the growth of a concentrated vortex at the leading edge of the airfoil, the so-called Leading Edge Vortex LEV. During the upstroke (α(t) ˙ > 0), the boundary layer remains thin and attached when the angle of attack is lower than the static stall. By increasing the angle of attack the flow starts to separate from the trailing edge; the separation point moves upstream until it reaches a certain distance from the leading edge [24, 25, 26]. In this condition, a vortex forms at the leading edge. The suctions induced by the LEV cause an increase of C L also if the boundary layer is massively separated. By further increasing the angle of attack, the LEV starts to move downstream with an almost constant velocity [24, 27]. Due to the convection of the vortex away from the airfoil, the C L drops. Sometimes, under particular flow and motion condition, multiple vortexes ˙ < 0), may shed from the leading edge of the airfoil [28], leading to multiple peaks of C L . During the downstroke (α(t) the boundary layer remains separated until the angle of attack reaches a threshold value, which depends on the motion conditions, under which the boundary layer starts to reattach. Let us consider an airfoil under a pitching motion: α(t) = α0 + α¯ sin(ωt)

(4)

where f = ω/2π is the frequency of the motion. In this work the final lift coefficient is written as a sum of four terms: C L f in (t) = C Ld + C Lv (t) + C Lshedding +

παc ˙ 2U∞

(5)

where the first term is the dynamic lift coefficient, evaluated without the contribution of the leading-edge vortex, the second term is the first leading edge vortex, the third is the contribution of the vortex shedding and the last one is the unsteady effect [29]. All the terms are fully explained in [30]. The experimental data published by McCroskey in [28] for a NACA0012 airfoil are used to calibrate and test our model.

4

B. Rocchio et al. / Energy Procedia 148 (2018) 1010–1017 Author name / Energy Procedia 00 (2018) 000–000

1013

3. Model validation In the following section, the present mode is now validated by comparing its results with the solution obtained by Delafin [31] via ANSYS CFX and to the experimental data by Akins et Al. [32]. 3.1. Computational setup The 17m Darrieus turbine of the SANDIA National Laboratories was chosen because of the availability of data of instantaneous torque (Akins et Al. [32]), which allow a detailed comparison of the aerodynamic models employed. The shape of the blades approximates the Troposkein configuration through a central circular arc and straight extremities. The rotor has N = 2 blades with a NACA 0015 airfoil as a cross-section with c = 0.61m. The equatorial diameter is D = 17m. The tower of the turbine has a diameter of 0.51m and is taken into account in the simulations. Figure 1 shows the domain size: the z = 0 plane is at 13.5m from the ground and corresponds to the equatorial plane of the turbine. On this plane, the origin (0, 0, 0) is placed at the center of the turbine tower. The inlet is at x = −500m, the outlet at x = 500m, the lateral boundaries at y = ±500m, while lower and upper boundaries at z = −13.5m and z = 170.5m respectively.

z

184m

0

y x

1000m 500m

1000m 500m

Fig. 1. Computational domain for the 3D case.

As for boundary conditions, by following the computational set-up in [31], at the inlet a velocity profile has been imposed (U(z) = U∞ (z/zre f )a ), where U∞ is the reference velocity, zre f is the mid height of the turbine and a = 0.1 [31]. At the lateral and upper boundaries symmetry is imposed, while no-slip condition is considered for the lower boundary. The density of the air is imposed to be ρ = 1kg/m3 according to [31]. The computational grids are realized in ANSYS ICEM and they are made up of two separated parts, a cylindrical inner grid including the turbine and an exterior grid for the rest of the domain, joined together with an interface boundary condition. The overall grid has 1054880 cells: 536607 in the external zone and 518273 in the inner zone while the volume swept by the rotor blades is made up of nring = 92 cells along the azimuthal direction and nz = 92 cells for the vertical direction. Structured multi-block grids have been generated for both parts of the computational domain making an extensive use of the “O-grid” technique, where all the single blocks are structured (i.e. only made by hexahedral cells). The region colored in blue is where the momentum source term (fturb in Eq. 3) is placed. The k − ω SST model is employed to close the problem. The inflow turbulent features have been set according to [31]: turbulent intensity 10%, turbulence viscosity ratio µt /µ = 100. The SIMPLEC velocity-pressure coupling algorithm is used. The spatial discretization is set to Green-Gauss node-based for gradient. Second order schemes are used for pressure, momentum, turbulent kinetic energy (k) and specific dissipation rate (ω) formulations. Second order implicit scheme is also adopted for the temporal discretization. 3.2. Results The results for the optimal Tip Speed Ratio (T S R = 4.6) are now discussed. Figure 3a shows the color contours of the spread lift coefficient on the blades region. In Figure 3b the time-averaged stream-wise velocity is shown for

B. Rocchio et al. / Energy Procedia 148 (2018) 1010–1017

1014

Author name / Energy Procedia 00 (2018) 000–000

5

Fig. 2. Computational grid: (a) vertical cut-plane; (b) equatorial and vertical cut-planes. Cells belonging to the volume swept by the blades are colored in blue.

a vertical plane, i.e. x − z, which contains the turbine rotational axis. In this simulation the shaft of the turbine is taken into account and its presence modifies the aerodynamic loads of the blades. The recovery of the turbine wake is visible; in fact, the size of the wake decreases along the stream-wise direction, i.e. the x-direction, due to the turbulence produced by the rotation of the turbine. Indeed, the wake is turbulent, so the fluid particles inside of the wake may move outside of it and vice versa. The fluid particles outside of the wake have a higher velocity than the ones inside and the mixing induced by the turbulence cause recovery of the wake. Moreover, by increasing the turbulence intensity the wake recovery is found closer to the rotor.

-8.5

udm_8

3

14.5

26

37.5

-2.0

0.5

3

5.5

8.0

Fig. 3. (a) Lift color contourns [N] for the SANDIA-17m Troposkien turbine; (b) Time-averaged stream-wise velocity color contourns [m/s].

7000

8000

6000

7000

5000

6000 Torque [Nm]

Torque [Nm]

Figure 4a shows the one-blade torque variation with respect to the azimuthal position of the blade and Figure 4b shows the total turbine torque. The red line is the reference URANS carried out by Delafin [31], while the black triangles are the experimental data published in [32].

4000 3000 2000 1000

5000 4000 3000 2000 1000

0

0

-1000

-1000

0

50

100

150

200 θ [°]

250

300

350

0

20

40

60

80 100 θ [°]

120

140

160

180

Fig. 4. Torque behavior for T S R = 4.6: (a) One-blade Torque: UDF (—), CFX by Delagin et Al. (—); (b) Turbine Torque: UDF (—), CFX by Delagin et Al. (—), experimental data By Akins et Al. ().

6

B. Rocchio et al. / Energy Procedia 148 (2018) 1010–1017 Author name / Energy Procedia 00 (2018) 000–000

1015

This case appears satisfactorily described; the global trend of the torque curve is in agreement with the reference one. During the upwind the present model anticipates the maximum torque value with a difference of 500Nm respect the reference solution. In downwind the torque deficit due to the presence of the shaft is visible. When the blade passes through the shaft wake, the incoming velocity seen by each section of the blade decreases and the blade reaches a higher angle of attack. The blade works in off-design condition and the torque decreases. The trend of the curve described by our model is slightly different from the one obtained by Delafin and this is plausibly due to all the assumptions made in modeling the blades. In particular, the dynamic stall model affects a lot the general trend of the Torque curve, especially the position of the maximum Torque. Despite the differences, the mean integral quadratic error between the two curves is about 0.013. The trend described above can be recognized also in the total torque curve, see Figure 4b. Again, from a qualitative point of view the matching is satisfactory. Compared to the experimental data [32] our curve is almost shifted at high value of azimuthal position. The maximum value of the total torque agrees with the data given by Akins et al [32]. We believe that a finer tuning of the free parameters employed in our code might results in a better quantitative description of the torque produced. 0.4 0.35 0.3

CP

0.25 0.2 0.15 0.1 0.05 2

3

4

5 TSR

6

7

8

Fig. 5. Power coefficient of the turbine for different TSR. UDF results (-•-), CFX by Delafin et Al. (), experimental data By Akins et Al. ().

In [31] is also reported the (C P − T S R) curve. Thus, to complete our analysis other simulations are carried out for different TSRs. Figure 5 shows how the power production changes with respect of the turbine TSR. The greatest disagreements are at low TSR values where the effects of the dynamic stall are predominant. 4. Practical Application: floating platform In recent years detailed CFD analysis have been performed by researchers from different countries to predict the aerodynamic effects of the platform motion for the NREL 5MW baseline HAWT. Our goal is to compare the effect of different pitch amplitudes on the performances of the VAWTs and HAWTs of the same size (5MW). Inside the danish project Future Deep Sea Wind Turbine Technologies or Deep Wind, a new wind wind turbine was studied and tested. The geometry of this machine is reported in [33] with a NACA0018 as cross-section of the blades. For the sake simplicity the turbine shaft is no more modeled. Before analyzing the effect of the pitch motion on the power production, a preliminary analysis is carried out in order to verify the reliability of our model. For the purpose, different velocities are considered and the results are compared to the ones obtained in [33] with the DMST model. The first analyzed velocity is 9m/s corresponding to the rated wind velocity and to the maximum power coefficient. By increasing the incoming velocity, the turbine rotor is regulated by stall, so the rotor angular speed is fixed and the maximum power production is reached for 14m/s, when the blades are stalled. The last condition is 11m/s, because is the one used for the case with a pitch motion. The nominal angular speed for this turbine is set to be 5.26rpm. From Table 1, it can be seen that our model is in agreement with the results of DMST, especially if the case of the optimal TSR is considered. A pitch motion in the vertical plane (y − z plane) respect to the lower tip of the blades is then considered. The incoming velocity is 11m/s while the pitch amplitude and frequency are, respectively, 4◦ and 0.1Hz, which are the

B. Rocchio et al. / Energy Procedia 148 (2018) 1010–1017 Author name / Energy Procedia 00 (2018) 000–000

1016

7

Table 1. Studied conditions without pitch re f

re f

Condition

Pmodel mean [MW]

C Pmodel

Pmean [MW]

CP

difference(%)

9m/s 11m/s 14m/s

1.694 2.824 4.574

35.26 32.25 25.31

1.70 3.00 4.75

35.58 34.25 26.30

0.89 5.84 3.76

same conditions for the NREL turbine in [34]. It must pointed out that this frequency is close to the DeepWind one, 0.088Hz. To evaluate the effect of a different amplitude on turbine performances an amplitude of 8◦ is also simulated. 7x106 6x106 5x106

P[W]

4x106 3x106 2x106 1x106 0 -1x106

0

5

10

15

20

25 t[s]

30

35

40

45

50

Fig. 6. Comparison between the istantenous power production for different pitch amplitude: no pitch (—), A = 4deg.(—), A = 8deg(—).

The initial position of the blades is assumed to be 0◦ and 180◦ , while the rotor starts its sinusoidal motion from the vertical position. Figure 6 shows that the platform pitching leads to a slight increase and decrease of the instantaneous production. Since the pitching period is different from the revolution one, the unsteady trend of the power can be observed. Due to this difference, after two pitching period the turbine has done 1.75 revolutions. Thus, at the beginning of the third pitching period the blades are at 90◦ (heavily loaded) and 270◦ , while the shaft has its higher velocity value (i.e. vertical position). This is the condition that has to be analyzed to check the effect of the pitching motion on the turbine performances. The reference case (i.e. no pitch) gives a mean power output of 0.322. For the case of 4◦ the mean power production reaches 0.329, while by increasing again the pitch amplitude up tp 8◦ the mean power production is almost 0.324. The power ripple for the reference case is about 2.5 and it increase by introducing the oscillations. In particular, we obtain and increment of 2.96% for 4◦ and 4.63% for 8◦ . From [34], the mean power coefficient for the NREL turbine is around 0.48. Without the pitch motion the HAWT power output is constant and it start to oscillate by moving the turbine. Therefore, for and oscillation of 4◦ the power ripple passes from 0 to 1.88. 5. Conclusion A momentum source function which is able to represent the behavior of a real vertical axis wind turbine is presented. All the classical aerodynamic phenomena are taken into account (i.e. the dynamic stall, the tip loss correction, the virtual camber effect). Our source function is able to well predict the power production of VAWTs when they work in design condition, while differences are found at low speed ratio numbers due to the importance of the dynamic stall model. Also, the virtual camber model must be improved in order to obtain a better prediction for high TSR numbers. At this step of the development of this complete tool practical applications are possible with acceptable results. This source function implies a low computational cost and it can be used for practical applications. It has been successfully applied to an off-shore wind turbine under pitching conditions which appears to be less sensible to oscillation amplitude variations compared to a HAWT of the same size.

8

B. Rocchio et al. / Energy Procedia 148 (2018) 1010–1017 Author name / Energy Procedia 00 (2018) 000–000

1017

References [1] G.J.M. Darrieus, Turbine having its rotating shaft transverse to the flow of the current, US1835018 A, 1931. [2] Toja-Silva F, Colmenar-Santos A, Castro-Gil M. Urban wind energy exploitation systems: Behavior under multi directional flow conditions, Opportunities and challenges, Renewable and Sustainable Energy Reviews 2013;24:364-378 (August 2013). [3] Riegler H. HAWT versus VAWT: small VAWTs find a clear niche. Refocus 2003;4(4):44–46 (July 2003). [4] Mertens S, van Kuik C, van Bussel G. Performance of an H-Darrieus in the skewed flow on a roof. J. Sol. Energy Eng 2003; 125(4):433-440 433–440 (November 2003). [5] Ferreira CJ, van Bussel G, van Kuik G. An analytical method to predict the variation in performance of a H-Darrieus in skewed flow and its experimental validation. Proceedings of the European Wind Energy Conference & Exhibition, EWEC, Athens (Greece); 2006. [6] Orlandi A, Collu M, Zanforlin S, Shires A. 3D URANS analysis of a vertical axis wind turbine in skewed flows. J Wind Eng Ind Aerodyn 2015;147:77-84. (December 2015). [7] S. Zanforlin, S. Letizia. Improving the Performance of Wind Turbines in Urban Environment by Integrating the Action of a Diffuser with the Aerodynamics of the Rooftops, Energy Procedia 82 (2015) 774–781. [8] Borg M, Collu M. A comparison between the dynamics of horizontal and vertical axis offshore floating wind turbines. Phil. Trans. R. Soc. 2015; A 373: 20140076. (January 2015). [9] Lei H., Zhou D., Lu J., Chen C., Han Z., Bao Y. The impact of pitch motion of a platform on the aerodynamic performance of a floating vertical axis wind turbine. Energy 2017, 119 369-383. [10] Lei H., Zhou D., Bao Y., Chen C., Ma N., Han Z. Numerical simulations of the unsteady aerodynamics of a floating vertical axis wind turbine in surge motion. Energy 2017, 127 1-17. [11] Cheng Z., Madsen H.A., Gao Z., Moan T. A fully coupled method for numerical modeling and dynamic analysis of floating vertical axis wind turbines. Renewable Energy 2017, 107 604-619. [12] Castelli MR, Pavesi G, Battisti L, Benini E, Ardizzon G. Modeling Strategy and Numerical Validation for a Darrieus Vertical Axis Micro-Wind Turbine, Int. Mech. Eng. Congr. Expo. (2010). [13] F. Balduzzi, J. Drofelnik, A. Bianchini, G. Ferrara, L. Ferrari, M.S. Campobasso, Darrieus wind turbine blade unsteady aerodynamics: a three-dimensional Navier-Stokes CFD assessment, Energy. 128 (2017) 550–563. [14] C. J. Simao Ferreira, A. van Zuijlen, H. Bijl, G. van Bussel, G. Van Kuik. Simulating dynamic stall in a two-dimensional vertical-axis wind turbine: verification and validation with particle image velocimetry data. Wind Energy Journal. 13:1-17 (2010). [15] S. Shamsoddin, F. Port´e-Agel. A large-eddy simulation study of vertical axis wind turbine wakes in the atmospheric boundary layer. Energies 9.5 (2016) 366. [16] J. N. Sørensen and W. Z. Shen. Computation of Wind Turbine Wakes using Combined Navier-Stokes/Actuator-line Methodology, pages 156– 159. EWEA, 1999. [17] M. Torresi, B. Fortunato, and S. M. Camporeale. Efficient three-dimensional CFD analysis of a Darrieus rotor. Conference paper: 67 Congresso Nazionale ATI – Trieste, 11-14 Settembre 2012 [18] S. Shamsoddin, F. Port´e-Agel. Large eddy simulation of vertical axis wind turbine wakes. Energies 7.2 (2014): 890-912. [19] S. Letizia, S. Zanforlin, Hybrid CFD-source terms modelling of a diffuser-augmented vertical axis wind turbine, Energy Procedia 101 (2016) 1280-1287. [20] P.G. Migliore, W.P. Wolfe, J.B. Fanucci, Flow Curvature Effects on Darrieus Turbine Blade Aerodynamics, J. Energy. 4 (1980) 49–55. https://doi.org/10.2514/3.62459. [21] L. Prandtl. Vier Abhandlungen zur Hydrodynamik und Aerodynamik. Universit¨atsverlag G¨ottingen SUB G¨ottingen, G¨ottingen, 2010. [22] W.Z. Shen, R. Mikkelsen, J.N. Sørensen, C. Bak, Tip loss corrections for wind turbine computations, Wind Energy. 8 (2005) 457–475. [23] W. Z. Shen, J. H. Zhang, J. N. Sørensen. The actuator surface model: a new Navier–Stokes based model for rotor computations. Journal of Solar Energy Engineering 131.1 (2009): 011002. ¨ G¨ulc¸at. Modern subjects in “Fundamentals of Modern Unsteady Aerodynamics”, Springer, Singapore, pp. 259-340 (2015). [24] U. [25] J.Nedi´c, Jovan, J. Vassilicos. Vortex Shedding and Aerodynamic Performance of Airfoil with Multiscale Trailing-Edge Modifications. AIAA Journal, 53, pp. 3240-3250 (2015). [26] W. Sheng, R. Galbraith, F.N. Coton. A Modified Dynamic Stall Model for Low Mach Numbers. Journal of Solar Energy Engineering, 130, (2008). [27] J.G. Leishman, T.S. Beddoes. A semi-empirical model for dynamic stall. Journal of the American Helicopter Society, 34, pp.3-17 (1989). [28] W.J. McCroskey. The phenomenon of dynamic stall. NASA Technical Memorandum 81264, (1981). [29] J.W. Larsen , S.R.K. Nielsen, S. Krenk. Dynamic stall model for wind turbine airfoils. Journal of Fluids and Structures, 23, pp. 959–982 (2007). [30] B. Rocchio, C. Chicchiero, M. V Salvetti, S. Zanforlin, Towards a general dynamic stall model, in: AIDAA XXIV Int. Conf., 2017. [31] P.L. Delafin, T. Nishino, A. Kolios, L. Wang, Comparison of low-order aerodynamic models and RANS CFD for full scale 3D vertical axis wind turbines, Renew. Energy. 109 (2017) 564–575. [32] R.E. Akins, D.E. Berg, W.T. Cyrus, Measurements and Calculations of Aerodynamic Torques for a Vertical-Axis Wind Turbine, (1987). [33] L. Vita, T. F. Pedersen e H. A. Madsen, Offshore vertical Axis wind turbine with floating and rotating foundation, PhD thesis Technical University of DenmarkDanmarks Tekniske Universitet, Department of Electrical Engineering Institut for Elektroteknologi, (2011). [34] T. T. Tran and D. H. Kim, The aerodynamic interference effects of a floating offshore wind turbine experiencing platform pitching and yawing motions, Journal of Mechanical Science and Technology 29, pp. 549-561, (2015).