# Optimization of a winglet for improving the performance of an H-Darrieus turbine using CFD

László Daróczy, Gábor Janiga, Dominique Thevenin

► To cite this version:

László Daróczy, Gábor Janiga, Dominique Thevenin. Optimization of a winglet for improving the performance of an H-Darrieus turbine using CFD. 16th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Apr 2016, Honolulu, United States. hal-01879336

HAL Id: hal-01879336

https://hal.science/hal-01879336

Submitted on 23 Sep 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimization of a winglet for improving the performance of an H-Darrieus turbine using CFD

László Daróczy<fnref n="1" />*, Gábor Janiga<fnref n="1" />, Dominique Thévenin<fnref n="1" />

<figure><img src="image_2.png" /></figure>

ISROMAC 2016

International
Symposium on
Transport

Phenomena and
Dynamics of

Rotating Machinery

Hawaii, Honolulu

April 10-15, 2016

## Abstract

The importance of wind energy has progressed rapidly in the last years. Although Horizontal Axis Wind Turbines (HAWT) are most well-spread, there is an increasing interest in Vertical Axis Wind Turbines (VAWT), especially in the H-Darrieus concept, as these rotors are omni-directional and affordable. However, the physics of these rotors is more complex; they can only be analyzed using transient CFD simulations. Due to the finite aspect ratio of the rotors, a wingtip vortex is created, which generates losses. Optimizing the wingtip geometry could be advantageous for increasing the efficiency of the rotors: this can only be achieved with three-dimensional turbulent transient simulations.

For the optimization of winglets, the whole process (mesh generation, CFD computation, post-processing) has to be automated. This is achieved using the OPtimization Algorithm Library++ (OPAL++), a custom C++ code for the description of blended and canted winglets, coupled with a CD-Adapco StarCCM+ JAVA script for the automatization of the mesh generation and CFD computations. To check the viability of the present concept, two parameters have been varied in the simulations. As shown in what follows, an efficient automatic optimization of wind turbine wingtips can be implemented in this manner.

## Keywords

CFD — wind energy — optimization — winglet

## INTRODUCTION

As a result of the growing customer awareness and due to increasing fuel prices, the importance of wind energy has progressed rapidly in the last years. Energy companies are strongly encouraged to adopt sustainable energy generation processes. The gross production of wind energy in the EU27 increased from 80 PJ (2000) to 537 PJ (2010) [1], and an increase to 3900 PJ is prognosed for 2050 [2]. As a result, further increasing the efficiency of the rotors is a critical area, which can only be accomplished by computer-based optimization.

Although Horizontal Axis Wind Turbines (HAWT) are most well-spread, there is an increasing interest in Vertical Axis Wind Turbines (VAWT), as these rotors, especially the H-Darrieus concept are quite affordable due to a simple design. Being omni-directional, such rotors can be installed even in urban environments and may be still quite efficient at low wind speeds. However, these rotors may suffer from large wingtip losses at small aspect ratio, which justifies the following analysis and optimization of wingtips.

## 1. H-DARRIEUS ROTORS

The original patent of the Darrieus rotor from Georges Jean Marie Darrieus dates back to 1927, but its very complex aerodynamics are still not completely understood, as the operation of VAWTs is far more complex compared to Horizontal Axis Wind Turbines. The operation of wind turbines is usually characterized by the characteristic curve, or TSR-$C_p$ curve. The tip-speed-ratio (TSR = $\lambda = \frac{v_{tip}}{u} = \frac{R\omega}{u}$) is the ratio of the

highest velocity of the rotating blade compared to the wind speed ($u$), where $R$ is the maximal radius of the rotor and $\omega$ is the angular velocity of the rotor.

Throughout the rotations, the blades of the Darrieus rotors experience significantly varying angle-of-attacks, see Fig. 1. This leads to stall/dynamic stall effects and results in complex wake interactions. As a result, these rotors can only be analyzed using transient simulations. The efficiency of the energy

<figure><img src="image_3.png" /><figcaption>Figure 1. Angle-of-attack during rotor rotation</figcaption></figure>

generation of a rotor is described by the power coefficient ($C_p$), which is the ratio of the power generated by the wind

<footnote n="1">¹Laboratory of Fluid Dynamics and Technical Flows, Otto von Guericke University Magdeburg, Magdeburg, Germany</footnote>

<footnote>*Corresponding author: laszlo.daroczy@ovgu.de</footnote>
turbine compared to the available energy in the wind:

$$C_p = \frac{P_{mech}}{P_{total}} = \frac{\bar{T}\omega}{\frac{1}{2}\rho u^3 A} = \frac{\bar{T}\omega}{\frac{1}{2}\rho u^3 2RH}, \quad (1)$$

where A is the projected area of the rotor, ρ is the density, H is the height of the rotor and T is the (average) torque generated by the rotor. The moment or torque coefficient (CQ) is the indicator for the generated torque, and is computed as:

$$C_Q = \frac{\bar{T}}{\frac{1}{2}\rho u^2 RA} = \frac{C_p}{\lambda}. \quad (2)$$

## 1.1 Blade tips

Darrieus rotors are lift-driven rotors equipped with airfoils. Airfoils generate lift through the pressure difference on the two sides of the airfoil. All blades have a finite aspect ratio. As a result of the pressure difference, an airflow is induced from the higher pressure side to the suction side along the end of the blade and a wingtip vortex is generated, which results in losses. In wind turbines with small aspect ratio these losses can be large (even up to 25%;[3]), while at high aspect ratio corresponding losses become negligible. As shown in the study of Gosselin et al., the last 20-30% of the blade height is less effective [4]. Amet also came to a similar conclusion, with 22% loss of performance due to the tip vortices and blade/arm junctions [5], while Qin et al. reported 40% losses compared to 2D simulations [6].

This aerodynamic problem has been known for a long time, mostly in connection with airplanes, and several blade tip modifications were already proposed to overcome it, e.g., vortex diffuser, tip sails [7], wing-grids, raked tip [8] and many different winglet designs. Some of them are presented in Figure 3. Winglets are now widely used for airplane wings [9] and, to a less extent, for wind turbines.

Only few studies consider winglets on wind turbines in the scientific literature. Taborda et al. [10] analyzed a cross-flow vertical axis wind turbine based on 3D transient turbulent simulations using the Transitional SST model. The original rotor without winglet, with an asymmetric winglet and with a symmetric winglet was analyzed. They came to the conclusion, that the performance can be increased using winglets, the optimal solution being to use a symmetric winglet. Berlferhat et al. experimentally analyzed the polars for different winglets [11]. The aerodynamic efficiency of HAWTs with blade tip modifications was also discussed in the work of Gauanaa and Johansen [12] and in the work of Imamura et al. using a vortex lattice method [13].

Winglets on VAWTs have been even less considered. Amato et al. investigated Darrieus rotors equipped with different blade tips [14]. The original configuration, an aerodynamic bulkhead, elliptic termination and three different winglets were evaluated using numerical simulations, validated by the measurements of Castelli et al. [3]. This study indicated that using a correct wingtip configuration, the performance can be increased. However, a poorly chosen solution can further decrease it.

## 2. AIM OF THE PROJECT

The blade shape optimization of Darrieus wind turbines has not been often considered in the literature. The present project is even more challenging, since the evaluation of the winglet performance can only be achieved using three-dimensional, transient, turbulent, large-scale Computational Fluid Dynamics (CFD) simulations. As a result, only very few articles deal with the CFD simulation of winglets for VAWTs (e.g., [14]). The optimization of such winglets has never been discussed in the literature according to the best knowledge of the authors.

The first step of the current project is to derive a fully automatic procedure to optimize winglets for H-Darrieus rotors. To show the viability of this proof-of-concept, the general geometrical description of parametric winglets is first given, followed by a simplified model parameterized by only two parameters. Afterwards, a very simple Design-of-Experiment is created with 6 different configurations and CFD simulations are carried out for a single tip-speed-ratio.

Finally, results are analyzed and conclusions are drawn. The evaluation of the cases is executed in a completely automatic manner in an optimization framework, proving the viability of the developed approach.

## 3. PARAMETERIZATION

Winglets are essentially the composition of one or more wings, with (usually) blended center line and non-uniform cross-section. Thus, the geometry of the winglet can be defined based on several support sections and/or by varying all parameters along the span. This approach is described, e.g., in the work of Sobieczky [15].

### 3.1 General parameterization

In the current study variation of parameters along the span were defined. The different geometrical parameters of the winglet are defined along a non-dimensional parameter (ψ = 0...1). First, a guiding line (see Fig. 2(a)) of the winglet is defined by

$$\mathbf{x}(\psi) = (x(\psi), y(\psi), z(\psi)). \quad (3)$$

This is equivalent to the definition of the parameters along the span as the arc length can be simply computed as

$$s(\psi) = \int_{0}^{\psi} \sqrt{\left(\left(\frac{dx}{d\hat{\psi}}\right)^2 + \left(\frac{dy}{d\hat{\psi}}\right)^2 + \left(\frac{dz}{d\hat{\psi}}\right)^2\right)} d\hat{\psi}. \quad (4)$$

For each ψ a local coordinate system (x_L, y_L, z_L) is defined with (x(ψ), y(ψ), z(ψ)) as center and with α(ψ), β(ψ), γ(ψ) angle rotation around the x-, y- and z-axis, respectively. In this coordinate system, a profile is mounted on the (x_L, y_L) plane, defined by the equation C(x_L, y_L, ψ) = 0. Using this method, the contour of the winglet is completely defined.

This parameterization is implemented in C++ and the centerline positions and profiles are exported at the p(ψ) = (p_1(ψ), p_2(ψ), ..., p_n(ψ)) positions as support sections. Star-CCM+ can create the surface of the winglets as a NURBS based on these data.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Parameterization (a) general; (b) specific variant</figcaption></figure>

To illustrate the approach flexibility provided by the present model, a couple of winglets are presented in Fig. 3, as created by the current workflow. As one can see, the most common types of winglets can be described and fully parameterized.

## 3.2 Specific parameterization

Although the previously presented (and widely used) parameterization is very flexible, it is too general to be used for an optimization (as all parameters are arbitrary functions). As a result, a simplified 6 m was defined, see Fig. 2 (b). In this third the frontal projection of the centerline of the winglet is composed of a vertical segment with length $L_1$ and an inclined segment with $L_2$ length, the angle between the two segments is denoted by $\gamma$. If the two sections are connected by a circular arc with a radius $R_1$. Besides the 6 ntal projection, the side view of the winglet is defined so, that its height is equal to the length of the frontal projection: this way the same 6 m can be bended on the 6 ntal projection independently from the size of $\gamma$. The side view is defined by the trailing and leading edge, which are both created by two vertical segments connected by circular arcs. The lower camber has a length of $c_0$, the winglet ends with a camber of $c_{end}$.

In the current study $R_{1,2,3}$, $c_0$, $c_{end}$, $L_1$, $L_2$ were fixed and $\gamma$ and $\alpha_2$ were varied.

## 3.3 Analyzed configurations

Following the setup of the workflow, 8 different configurations were created and meshed. For ease of notation and meshing, the torus with $H_{half}=707+200$ mm height, 'Endplate' denotes the 'Base' configuration equipped with a NACA0030 endplate ($c = 150.15$ mm, $r_T = 4$ mm unding, $t = 2$ mm thickness).

Configurations A1, A2, B1, B2, C1 and C2 denote winglets with length $L_1 + L_2 = 707 + 200$ mm. For the transition $R_{1,2,3} = 100$ mm was applied. For the camber $c_0 = 85.8$

mm, $c_{end} = 85.8/3$ mm was retained, A, B and C denote respectively $\alpha_2 = 0^\circ, 15^\circ$ and $30^\circ$. Numbers 1 & 2 denote configurations with $\gamma = 60^\circ$ and $\gamma = 30^\circ$, respectively.

The configurations can be seen in Figure 4.

# 4. OPTIMIZATION WORKFLOW

For the optimization of winglets, the whole process (mesh generation, CFD computation, post-processing) has to be automated. This was achieved using the OPtimization Algorithm Library++ 2.4.5 (called simply OPAL++). OPAL++ is an object-oriented multi-objective optimization and parameterization framework developed at the University of Magdeburg "Otto von Guericke" [16]. The software has already been successfully applied to many different problems [16, 17] and is focused especially on CFD-based problems.

OPAL++ is able to operate in parallel on many PCs to speed up the optimization process. For a single configuration, the following steps are performed:

1. The geometry is created using a custom C++ code (app. 1200 lines), which is called by OPAL++ as an ODL (OPAL++ Dymic Libary, the surface geometry is exported in CSV format.

2. The geometry is read in CD-Adapco StarCCM+ and the surface is prepared using a NURBS surface. Afterwards, polyhedral mesh is generated. All operations are performed by a custom JAVA script (4000 lines).

3. The mesh quality is checked.

4. The simulation is started in parallel with CD-Adapco StarCCM+ (in the present study each simulation was performed with 24 cores using Intel(R) Xeon(R) E5-1650v3 3.5 GHz CPUs; each simulation required 1-2 weeks to finish). Upon completion, StarCCM+ exports the y+ values, forces and moments for the different parts of the rotor and for the whole rotor.

5. Post-processing is performed using a small bash script through OPAL++.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Example of different winglet designs: (a) Whitcomb winglet as a composition of two single centerline winglets; (b) Canted winglets; (c) Simple blended winglets; (d) Blended winglets with more complex leading edge</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Tested configurations from (a) isometric; (b) front; (c) top; (d) side view (gray='base'; black='endplate'; red='C2'; blue='B2'; purple='A2'; brown='C1', yellow='B1'; green='A1')</figcaption></figure>

6. Results are collected in a Microsoft Excel file for all the different winglets.

## 5. MESH GENERATION

After the geometry was imported into StarCCM+ v10, the mesh generation is executed in a completely automatic manner based on a custom JAVA script developed for the present case. Relevant parts of the geometry are retained and combined with the surroundings with the help of the 'Surface Wrapper'. Mesh for the stationary and rotating domain is created using polyhedral mesh, additional refinements are created with several volume controls, surface and edge controls to ensure an appropriate resolution. In the stationary domain, refinement is added for the wake, in the rotating domain for the blade, wingtip and trailing edge (see Fig. 5). Extreme care was given to resolve the flow around the wingtips. The most important settings are summarized in Table 1.

The inlet and outlet domains are extended to allow development of the flow and to ensure an appropriate distance from the boundaries with the help of 'Directed Mesh'.

## 6. CFD SETUP

### 6.1 Physical model

As a series of experimental validations based on two-dimensional CFD computations have revealed that k-ε Realizable and k-ω SST models could be appropriate for the simulation of the flow around an H-Darrieus rotor [18], only these two models were retained for the analysis. The flow around the

rotor was modeled as incompressible (this is a widely adopted simplification in the literature [3], as the relative speed of the blades stays well below $M = 0.3$ due to the relatively small rotation speed of small Darrieus rotors). Reference density $ρ = 1.225$ kg/m$^3$ was chosen according to DIN EN 61400 [19], and the dynamic viscosity was chosen correspondingly as $μ = 1.7894 · 10^{-5}$ Pa s.

The system of equations was solved in a Coupled manner, as the solver was revealed to be superior by earlier studies (according to the CD-Adapco StarCCM+ User Guide [20], the Coupled solver is always recommended when sufficient resources are available due to its increased robustness). Additionally, Cell Quality Remediation was activated. This model can efficiently improve stability locally if divergence is detected in a couple of deformed cells, which is unavoidable in an automatic 3D mesh generation process. To promote precision, second-order implicit temporal discretization is always applied in the current work together with second-order (central or upwind) derivatives.

### 6.2 Boundary Conditions

The sides of the computational domain are defined as velocity inlet, pressure outlet and symmetry boundaries (except for the validation in the wind tunnel, where no-slip conditions are applied except for the symmetry plane). Proper inlet conditions are clearly the most important ones: $u = 9$ m/s, $I = 0.1\%$ and $TVR = 10$ were retained for all configurations, in agreement with recommendations from the literature [21, 22] and as a result of preliminary studies [18].
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Cross-section of mesh at different locations</figcaption></figure>

For the computations with winglets, $\lambda = 2.6$ was used in all cases. For the experimental validation the angular velocity was varied while keeping the wind speed constant.

## 6.3 Multi-fidelity approach

In order to speed-up the CFD computation, a multi-level approach was adopted. In this multi-fidelity computational approach not only the temporal, but the mesh resolution was varied as well. The following three steps were defined for the computations:

1. Level 1 - Initialization of flow field: Coarse mesh resolution ($s_{act,i} = s_{calib,i} \cdot 1.5$), large time step ($\Delta\varphi = 5^\circ$), $N = 6$ revolutions, Solver Courant number=15.

2. Level 2 - Initialization of flow field near the blades: Fine mesh solution ($s_{act,i} = s_{calib,i}$), intermediate time step ($\Delta\varphi = 1^\circ$), $N = 1$ revolution, Solver Courant number=30.

3. Level 3 - Precise computation of torque: Fine mesh solution ($s_{act,i} = s_{calib,i}$), small time step ($\Delta\varphi = 0.5^\circ$), $N = 1$ revolution, Solver Courant number=30.

For all cases $N = 30$ inner iterations were performed, in order to ensure very low residuals (e.g., 6 r B1 all residuals stayed below $10^{-7}$ in absolute value). Although it would be advantageous to increase the number of revolutions, this is definitely not possible, as the detailed mesh and temporal

resolution already lead to huge requirements for the computational time. Additionally, 200 would be better suited for the solver Courant number, but this resulted in instability and divergence. With the current settings, an analysis of the coefficients revealed an appropriate convergence, and the results did not change significantly between Level 2 and Level 3, see Figure 6. Thus, additional revolutions were not considered necessary for the computations.

Mesh sizes were calibrated so that the number of cells stays in the acceptable range for the available computational system (below 12 million polyhedral cells).

## 7. VALIDATION OF THE MODEL

The experimental validation relies on measurements from the literature. In the work of Castelli et al. [3] a small H-Darrieus rotor was tested. The rotor has NACA0021 blades with 0.38 mm blunt edge, $c = 85.8$ mm chord length and was mounted at $0.5c$. The radius and height were 1030 mm and 1414 mm respectively. The geometry is summarized in Table 2. The measurement was conducted in the Milan-Bovisa low turbulence wind tunnel with 4 m x 3.8 m cross-section at 9 m/s wind velocity. In the original experiment the results were not corrected for wind tunnel blockage. Castelli et al. [3] found that the best agreement was delivered by the $k-\omega$ SST model, with $y^+ > 30$ and a resolution of 2 mm on the surface of the blades.
<table><caption>**Table 1.** Settings for mesh generation</caption><thead><tr><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>base size in stator</td><td>0.12</td></tr><tr><td>size on the rotating interface</td><td>14 mm</td></tr><tr><td>size in the wake of the rotor</td><td>55 mm</td></tr><tr><td>size on the edge of straight blade</td><td>1.25 mm</td></tr><tr><td>target maximal y+ value</td><td>100</td></tr><tr><td>base size in the rotating domain</td><td>14 mm</td></tr><tr><td>size around the blades</td><td>5.5 mm</td></tr><tr><td>size on the surface of the blades</td><td>2 mm</td></tr><tr><td>size near the trailing edge</td><td>0.5 mm</td></tr><tr><td>size around the wingtips</td><td>4 mm</td></tr><tr><td>size on the surface near the wingtips</td><td>0.5...0.75 mm</td></tr><tr><td>prism growth rate</td><td>1.3</td></tr><tr><td>volume growth rate</td><td>1.06</td></tr><tr><td>surface growth rate</td><td>1.075</td></tr><tr><td>volume blending</td><td>0.68</td></tr><tr><td>volume density</td><td>0.8</td></tr><tr><td>length of extruded inlet domain</td><td>13 m</td></tr><tr><td>length of extruded outlet domain</td><td>17 m</td></tr><tr><td>length of stator domain</td><td>4+6 m</td></tr><tr><td>width of the domain</td><td>12 m</td></tr><tr><td>height of domain</td><td>7.68 m</td></tr></tbody></table>

<table><caption>**Table 2.** Geometry for validation & optimization</caption><thead><tr><th>Parameter</th><th>Valid.</th><th>Optim.</th></tr></thead><tbody><tr><td>number of blades</td><td>3</td><td>3</td></tr><tr><td>airfoil</td><td>NACA0021</td><td>NACA0021</td></tr><tr><td>radius of the rotor</td><td>515 mm</td><td>515 mm</td></tr><tr><td>camber</td><td>85.8 mm</td><td>85.8 mm</td></tr><tr><td>mounting</td><td>0.5c</td><td>0.5c</td></tr><tr><td>(half)-length of blades</td><td>717 mm</td><td>707+200 mm</td></tr><tr><td>wind speed</td><td>9 m/s</td><td>9 m/s</td></tr><tr><td>TSR</td><td>2.05/2.35/2.6/3.1</td><td>2.6</td></tr></tbody></table>

To reduce the mesh size, only half of the rotor and wind tunnel were retained. The mesh was created with the same setup as the optimization process itself (except the wall of the wind tunnel, where prism layers were added and no-slip condition was set).

The recommendation of Castelli et al. was retained, i.e., $y^+ > 30$ and 2 mm cell size on the blades. The resulting mesh contained 7.05 million polyhedral cells. The simulations were performed with the $k-\epsilon$ Realizable and $k-\omega$ SST model for tip-speed-ratios for $\lambda = 2.05$, $\lambda = 2.35$, $\lambda = 2.6$, $\lambda = 3.1$ (near the optimal tip-speed-ratios; the curves between the operating points have been interpolated with cubic spline); the blunt edge was disregarded and modeled instead as being sharp. The characteristic curve of the complete rotor is compared in Fig. 7 to the measurement of Castelli et al.; this figure also contains the characteristic curve far away from the wingtip, near the symmetry plane.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Torque coefficient for C1 configuration for a single blade</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Experimental validation</figcaption></figure>

As one can see, the agreement in the location of the optimal tip-speed-ratios is quite good (exp.: $\lambda_{opt} = 2.7$; sim.: $\lambda = 2.5$), but the exact values are not reproduced. As struts were not modeled and gear losses were not computed, this is not really a surprise. Furthermore, the large deviation in the trend at small tip-speed-ratios is most likely due to the application of the Wall Function approach ($y^+ > 30$). In this domain the flow is very complex due to the presence of dynamic stall and a finer grid resolution would be required.

The two turbulence models deliver very similar results. This is again not a surprise: the $k-\omega$ SST model is a blended model, which falls back to $k-\epsilon$ in the free stream and to $k-\omega$ near the wall. Since $y^+ > 30$ was employed, the benefits offered by the $k-\omega$ model are replaced virtually with the Wall Functions. As a result, the $k-\epsilon$ Realizable model was retained for all further computations.
## 8. RESULTS AND DISCUSSION

The direct comparison of the performance of the different rotors is not possible, as the projected area (see Table 3) and height of the rotor changes from configuration to configuration. As a result, four different analysis were performed, as discussed in the followings.

<table><caption>**Table 3.** Projected area of the different configurations</caption><thead><tr><th></th><th>Proj. area [mm²]</th></tr></thead><tbody><tr><td>Base/Endplate</td><td>934210 mm²</td></tr><tr><td>[A/B/C]1</td><td>812814 mm²</td></tr><tr><td>[A/B/C]2</td><td>889168 mm²</td></tr></tbody></table>

### 8.1 Performance coefficient far away from the wingtip

In the first step the performance coefficient of the first 0.1 m blade section was analyzed, far away from the wingtip (near the symmetry plane). Although far away from the wingtip it is expected, that the different rotors will deliver a similar performance, this is not necessarily true. The wingtip vortex changes the flow structures near the end of the blades, and through 3D instabilities it can have an effect on the flow structure of the whole rotor, see Table 4. As one can see, the performance coefficients remain very similar to each other; only small differences can be seen. Nevertheless, more than 1% (relative) improvement is detected for the B2 and C2 configurations.

This can be either the effect of the decreased energy extraction of B2 and C2 configurations near the wingtip, which results in an increased energy flux, or could be a result of the different blockage, since the wingtips restrict the 3D expansion in the vertical direction, leading to an increased efficiency. For rotors with high aspect ratios (AR) this effect could provide an advantage.

<table><caption>**Table 4.** Performance coefficients far away from the wingtips</caption><thead><tr><th></th><th>C<sub>p</sub><sup>first</sup> [%]</th><th>Improv.</th></tr></thead><tbody><tr><td>Base</td><td>34.72%</td><td>—</td></tr><tr><td>Endplate</td><td>34.74%</td><td>0.0576%</td></tr><tr><td>A1</td><td>34.84%</td><td>0.3456%</td></tr><tr><td>A2</td><td>34.95%</td><td>0.662%</td></tr><tr><td>B1</td><td>35.02%</td><td>0.864%</td></tr><tr><td>B2</td><td>35.11%</td><td>1.12%</td></tr><tr><td>C1</td><td>34.73%</td><td>0.0288%</td></tr><tr><td>C2</td><td>35.20%</td><td>1.382%</td></tr></tbody></table>

### 8.2 Performance coefficient of the complete rotor

In the second step the total performance coefficient of the rotor was analyzed, non-dimensionalized by the projected area for

each case. The results can be seen in Table 5. As one can see, instead of the expected improvement, the winglets provide a decrease in the performance for the current case. However, in case of the endplate variant, a small improvement can be seen. The improvement is small, showing that the wingtip losses do not seem to be that significant in this case. Obviously, this statement does not necessarily hold for rotors with more significant wingtip effects.

<table><caption>**Table 5.** Performance coefficients for the whole rotor</caption><thead><tr><th></th><th>C<sub>p</sub><sup>rotor</sup> [%]</th><th>Improv.</th></tr></thead><tbody><tr><td>Base</td><td>32.47%</td><td>—</td></tr><tr><td>Endplate</td><td>32.68%</td><td>0.65%</td></tr><tr><td>A1</td><td>31.84%</td><td>-1.94%</td></tr><tr><td>A2</td><td>31.95%</td><td>-1.60%</td></tr><tr><td>B1</td><td>31.75%</td><td>-2.22%</td></tr><tr><td>B2</td><td>31.29%</td><td>-3.63%</td></tr><tr><td>C1</td><td>31.36%</td><td>-3.419%</td></tr><tr><td>C2</td><td>30.23%</td><td>-6.899%</td></tr></tbody></table>

### 8.3 Effect of the wingtip on the distribution of the energy generation

Afterwards, the energy generation of the different sections were analyzed and compared graphically. The results can be seen in Fig. 8. Section 0 denotes the first 0.1 m near the symmetry plane (far away from the wingtip), while Section 6 denotes the last section (0.6 m - to the end of the blade). All results are non-dimensionalized by the projected area of the section. As one can see, all configurations resulted in a decreased performance coefficient for the last section compared to the base geometry, except the endplate. This supports the importance of endplates. The decreased performance of the winglets is most probably the result of the decreased camber length and rotor radius. For an efficient improvement, smaller winglet configurations or rotors with more significant wingtip losses would be needed.

A2 proved to be the best configuration among the winglets for Section 6. The B2 and C2 configurations provided an increased performance coefficient far away from the wingtip, as discussed previously.

### 8.4 Size of wingtip vortex

Finally, the vortex structures were compared for a single angular position for all the different cases ($\varphi = 0°$, starting position), see Fig. 9. In the figure the isosurface for 75 $\frac{1}{s}$ vorticity magnitude is presented along with the vorticity magnitude in two cross-sections. As one can see, A1 and B1 configurations efficiently reduced the size of the vortices, while C1 and C2 provided similar sizes to the base and endplate configuration.

Additionally, a threshold was created for 75 $\frac{1}{s}$ vorticity magnitude and the volume of this threshold was computed to estimate the size of the vortex structures. The results are shown in Table 6. The values further confirm, that A1 and B1 provide the smallest vortex volumes, while C2 shows similar
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Performance coefficients of the different blade sections for the different configurations</figcaption></figure>

values to the base configuration. Interestingly, the rotor with the endplate did not show a decrease in the size of the vortex.

<table><caption>**Table 6.** Volume of fluid domain with vorticity magnitude larger than 75 $\frac{1}{s}$</caption><thead><tr><th></th><th>V<sub>vortex</sub><br>75 [m<sup>3</sup>]</th></tr></thead><tbody><tr><td>Base</td><td>4.437 &middot; 10<sup>-2</sup></td></tr><tr><td>Endplate</td><td>4.452 &middot; 10<sup>-2</sup></td></tr><tr><td>A1</td><td>2.706 &middot; 10<sup>-2</sup></td></tr><tr><td>A2</td><td>2.910 &middot; 10<sup>-2</sup></td></tr><tr><td>B1</td><td>2.903 &middot; 10<sup>-2</sup></td></tr><tr><td>B2</td><td>3.197 &middot; 10<sup>-2</sup></td></tr><tr><td>C1</td><td>3.298 &middot; 10<sup>-2</sup></td></tr><tr><td>C2</td><td>3.650 &middot; 10<sup>-2</sup></td></tr></tbody></table>

## 9. CONCLUSIONS

In the current manuscript a completely automatic method was developed and presented for the efficient optimization of winglets for VAWTs. Parameterization, geometry and mesh generation, and CFD simulation are efficiently realized in a fully automatic way.

Following a validation of the model using a rotor from the literature, a Design-Of-Experiment was executed with 8 different configurations to check the validity of the developed approach. Unfortunately, no configuration was found that leads to a performance improvement compared to the rotor with a simple endplate. However, two configurations show an increased performance far away from the wingtips, which could be advantageous for rotors with high aspect ratio. To understand this effect, further studies are needed. Some questions are still open, which will be the subject of the upcoming research. Among others (1) mesh independency is to be analyzed in detail; (2) Spalart-Allmaras, Transitional SST and

LES models are to be analyzed in 3D and (3) the validation is to be carried out for other rotor geometries as well.

After checking these issues, the optimization will be executed to identify if a better wingtip configuration can be identified.

## REFERENCES

[1] European Commission. ENERGY - Country Factsheets V.1.3. 2012.

[2] European Commission - Directorate-General for Energy. EU energy, transport, and greenhouse gas emissions trends to 2050. 2013.

[3] Marco Raciti Castelli, Guido Ardizzon, Lorenzo Battisti, Ernesto Benini, and Giorgio Pavesi. Modeling Strategy and Numerical Validation for a Darrieus Vertical Axis Micro-Wind Turbine (IMECE2010-39548). ASME 2010 International Mechanical Engineering Congress and Exposition, 7:409–418, 2010.

[4] Rémi Gosselin, Guy Dumas, and Matthieu Boudreau. Parametric study of H-Darrieus vertical-axis turbines using uRANS simulations (Paper CFDSC-2013 178). In 21st Annual Conference of the CFD Society of Canada, Sherbrooke, Canada, May 2013.

[5] Ervin Amet. Simulation numérique d'une hydrolienne à axe vertical de type Darrieus. PhD thesis, Polytechnical National Institute, France, Grenoble, 2009.

[6] Ning Qin, Robert Howell, Naveed Durrani, Kenichi Hamada, and Tomos Smith. Unsteady flow simulation and dynamic stall behaviour of vertical axis wind turbine blades. Wind Engineering, 35(4):511–528, 2011.

[7] Ilan Kroo. DRAG DUE TO LIFT: Concepts for Prediction and Reduction. Annual Review of Fluid Mechanics, 33(587-617), 2001.

[8] Joel E. Guerrero, Dario Maestro, and Alessandro Bottaro. Biomimetic spiroid winglets for lift and drag control. Comptes Rendus Mécanique, 340(1-2):67 – 80, 2012.

[9] Richard T. Whitcomb. A design approach and selected wind-tunnel results at high subsonic speeds for wing-tip mounted winglets (NASA TN D-8260). Technical report, NACA, 1976.

[10] M.A. Taborda, S. Laón, O.D. López, and Y.U. López. Design Improvement of Cross Flow Water Turbines using CFD. In International Conference on Renewable Energies and Power Quality (ICREPQ'15), 2015.

[11] Belferhat, S., Meftah, S.M.A., Yahiaoui, T., and Imine, B. Aerodynamic Optimization of a Winglet Design. EPJ Web of Conferences, 45:01010, 2013.

[12] Mac Gaunaa and Jeppe Johansen. Determination of the maximum aerodynamic efficiency of wind turbine rotors with winglets. Journal of Physics: Conference Series, 75(1):012006, 2007.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Vortex structures near the configurations (first line: base, endplate, A1, A2; second line: B1, B2, C1, C2; isosurface: 75 1/s)</figcaption></figure>

[13] Hiroshi Imamura, Yutaka Haseguwa, and Koji Kikuyama. Numerical analysis of the horizontal axis wind turbine with winglets. *JSME International Journal Series B*, 41(1):170–176, 1998.

[14] Federico Amato, Gabriele Bedon, Marco Raciti Castelli, and Ernesto Benini. Numerical Analysis of the Influence of Tip Devices on the Power Coefficient of a VAWT. *World Academy of Sci., Eng. and. Tech.*, 78:387–395, 2013.

[15] H. Sobieczky. *Notes on Numerical Fluid Mechanics*, volume 68, chapter Parametric Airfoils and Wings. Springer Vieweg, 1998.

[16] László Daróczy, Gábor Janiga, and Dominique Thévenin. Systematic analysis of the heat exchanger arrangement problem using multi-objective genetic optimization. *Energy*, 65:364–373, 2014.

[17] László Daróczy, M.H. Mohamed, Gábor Janiga, and Dominique Thévenin. Analysis of the effect of a slotted flap mechanism on the performance of an H-Darrieus turbine using CFD (GT2014-25250). In *ASME Turbo Expo Conference*, Düsseldorf, 2014.

[18] László Daróczy, Gábor Janiga, Klaus Petrasch, Michael Webner, and Dominique Thévenin. Comparative analysis of turbulence models for the aerodynamic simulation of H-Darrieus rotors. *Energy*, 2015. Accepted for publication.

[19] DIN Deutsches Institut für Normung e.V. und VDE Verband der Elektrotechnik Elektronik Informationstechnik e. V. Design requirements for small wind turbines, DIN EN 61400-2 (VDE 0127-2):2007-2. Technical report, 2007.

[20] CD-adapco. *USER Guide Star-CCM+ Version 9.06*. CD-adapco, 2014.

[21] Rosario Lanzafame, Stefano Mauro, and Michele Messina. 2D CFD modeling of H-Darrieus wind turbines using a transition turbulence model. *Energy Procedia*, 45:131–140, 2014.

[22] Robin Blair Langtry. *A Correlation-Based Transition Model using Local Variables for Unstructured Parallelized CFD codes*. PhD thesis, Universitaet Stuttgart, 2006.