UNIVERSIDADE ESTADUAL PAULISTA “JÚLIO DE MESQUITA FILHO” 

CÂMPUS DE ILHA SOLTEIRA 

 

 

 

 

 

 

GIUSEP ALEXANDER BACA BERNABE 

 

 

 

 

 

 

OPTIMIZATION OF A H-DARRIEUS TURBINE THROUGH PARAMETRIC 

GEOMETRY BASED ON JOUKOWSKY TRANSFORMATION FOR LOW-TSR 

OPERATIONAL RANGES 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Ilha Solteira 

2024



 

 

 
 

GIUSEP ALEXANDER BACA BERNABE 

 

 

 

 

 

 

 

OPTIMIZATION OF A H-DARRIEUS TURBINE THROUGH PARAMETRIC 

GEOMETRY BASED ON JOUKOWSKY TRANSFORMATION FOR LOW-TSR 

OPERATIONAL RANGES 

 
 
 
 

Dissertação de Mestrado apresentada ao 
Programa de Pós-Graduação de Engenharia 
Mecânica da Universidade Estadual Paulista 
– UNESP, como requisito para a obtenção 
do título de Mestre.  
 
Área de concentração: Ciências Térmicas 
 
Orientador: Prof. Dr. Leandro Oliveira 
Salviano 

 
 
 
 
 
 
 
 
 
 
 
 

 
 

Ilha Solteira 

2024



Baca BernabeOPTIMIZATION OF A H-DARRIEUS TURBINE THROUGH PARAMETRIC GEOMETRY BASED ON JOUKOWSKY TRANSFORMATION FOR LOW-TSR OPERATIONAL RANGESDissertação de Mestrado apresentada ao Programa de Pós-Graduação de Engenharia Mecânica da Universidade Estadual Paulista – UNESP, como requisito para a obtenção do título de Mestre.Ilha Solteira2024 71 Sim Dissertação (mestrado)Engenharia Mecânicaengenharia mecanica, ciências térmicasNão

Dissertação de Mestrado apresentada ao Programa de Pós-Graduação de Engenharia M

.

.

.

FICHA CATALOGRÁFICA

Desenvolvida pelo Serviço Técnico de Biblioteca e Documentação

Baca Bernabe, Giusep Alexander.
Optimization of a H-Darrieus turbine through parametric geometry based on 

Joukowsky transformation for low-TSR operational ranges / Giusep Alexander 
Baca Bernabe. -- Ilha Solteira: [s.n.], 2024

84 f. : il.

Dissertação (mestrado) - Universidade Estadual Paulista. Faculdade de 
Engenharia de Ilha Solteira. Área de conhecimento: Ciências Térmicas, 2024

Orientador: Leandro Oliveira Salviano
Inclui bibliografia

1. Tip speed ratio. 2. H-Darius turbine. 3. Metamodel of optimal prognosis.
4. Response surface. 5. Joukowsky transformation. 6. Gradient-based
optimization.

B116o

Elaborada por Raiane da Silva Santos - CRB-8/9999





 

 

 
 

ACKNOWLEDGMENTS 

To my family, friends, research colleagues and professors of the department of 

Mechanical Engineering. More specifically, I dedicate my thanks to: 

 My father Humberto and my mother Ruth, for the love, understanding 

and strength given in every step of my personal life and professional 

development, for their messages of encouragement, positive and 

constant support in all the needs that presented themselves. 

 My advisor Leandro Salviano, for the discussions, ideas to improve my 

work of research, clarification of concepts with which I was not familiar, 

disinterested friendship during mentorship on my masters and the help 

so that I could adapt in the city of Ilha Solteira to which I moved from 

my home in Lima, Peru since the beginning of my master's degree. 

 My colleagues, from the Computational Simulation Laboratory in 

Thermal Sciences, for their friendship and exchange of constant ideas 

to improve, not just academically, but also as a person. 

This work was carried out with the support of the Coordination of Improvement 

of Higher Education Personnel - Brazil (CAPES) - Financing Code 001. 

 



 

 

 
 

ABSTRACT 

 

The depletion of traditional fossil fuel resources due to increased energy usage has 

led to efforts to promote renewable energy sources such as wind energy. This study 

aims to identify the optimal profile for an operational range (TSR = 2.33, 2.64, and 

3.09) in terms of performance through the numerical integration of the Cp-TSR curve 

by combining Computational Fluid Dynamics (CFD) with the Metamodel of Optimal 

Prognosis (MOP) response surface approach. Three parameters (a/b ratio, m, and 

pitch angle) are analyzed using the Joukowsky transformation parametrization 

methodology for symmetrical airfoils. The optimized turbine, performed using 

Gradient-based approaches, achieved a higher average Cp at each TSR compared 

to the traditional NACA 0021. Furthermore, a comparative analysis of non-

dimensional velocity and pressure contour between the optimal and the original 

profile was conducted. The study demonstrated that the approach used in coupling 

CFD modeling and optimization tools is reliable, cost-effective, and appropriate for 

fluid flow phenomena for high-dimensional and computationally expensive CFD 

models at different TSR.  

 

Keywords: tip speed ratio; H-Darius turbine; metamodel of optimal prognosis; 

response surface; Joukowsky transformation; gradient-based optimization. 

 

 



 

 

 
 

RESUMO 

 

A depleção dos recursos tradicionais de combustíveis fósseis devido ao aumento do 

uso de energia tem levado a esforços para promover fontes de energia renovável, 

como a energia eólica. Este estudo visa identificar o perfil ótimo para uma faixa 

operacional (TSR = 2,33, 2,64 e 3,09) em termos de desempenho, por meio da 

integração numérica da curva Cp-TSR combinando a Dinâmica de Fluidos 

Computacional (CFD) com a abordagem de superfície de resposta do Metamodelo 

de Prognóstico Ótimo (MOP). Três parâmetros (relação a/b, m e ângulo de 

inclinação) são analisados utilizando a metodologia de parametrização da 

transformação de Joukowsky para aerofólios simétricos. A turbina otimizada, 

realizada usando abordagens baseadas em Gradiente, alcançou um Cp médio mais 

alto em cada TSR comparado ao tradicional NACA 0021. Além disso, foi conduzida 

uma análise comparativa do contorno de velocidade e pressão adimensional entre o 

perfil ótimo e o original. O estudo demonstrou que a abordagem utilizada no 

acoplamento da modelagem CFD e ferramentas de otimização é confiável, custo-

efetiva e apropriada para fenômenos de fluxo de fluidos para modelos CFD de alta 

dimensão e computacionalmente caros em diferentes TSR. 

 

Palavras-chave: razão de velocidade na ponta; turbina H-Darius; metamodelo de 

prognóstico ótimo; superfície de resposta; transformação de Joukowsky; otimização 

baseada em gradiente. 



 

 

 
 

LIST OF FIGURES 

Figure 1. Three bladed H-Type Darrieus turbine. ...................................................... 20 

Figure 2. Diagram depicting the aerodynamic coefficients of an H-Type Darrieus 

turbine. ...................................................................................................................... 21 

Figure 3. Fluid domain dimensions and boundary conditions .................................... 30 

Figure 4. Global mesh of the computational domain (outer and inner domain). ........ 32 

Figure 5. Inner domain, airfoil leading and trailing edge mesh. ................................. 32 

Figure 6. Average power coefficient 𝐶𝑝 for the first 90 revolutions ............................ 34 

Figure 7. Average power coefficient 𝐶𝑝 after 20 revolutions ...................................... 34 

Figure 8. Relative change of the average power coefficient 𝐶𝑝 with respect to the 𝐶𝑝 

at the 90th revolution .................................................................................................. 35 

Figure 9. History of power coefficient 𝐶𝑝 for various azimuthal increments ............... 36 

Figure 10. Comparison of experimental and numerical validation results for a small-

scale NACA 0021 Darrieus VAWT operating at TSR = 2.64 and 9 m/s free stream 

velocity. ..................................................................................................................... 39 

Figure 11. Comparison between the k-w model and the experimental data presented 

by (RACITI CASTELLI; ENGLARO; BENINI, 2011) .................................................. 39 

Figure 12. Flowchart representing the full plan for CFD modeling, Sensitivity Analysis, 

and Optimization ....................................................................................................... 40 

Figure 13. Joukowski transformation. ........................................................................ 42 

Figure 14. Contrast between different parameters configuration (𝑎, 𝑏 and 𝑚). .......... 44 

Figure 15. Schematic showing the pitch angle 𝛽 for a VAWT blade .......................... 45 

Figure 16. 3D response surfaces rendered for each individual TSR; where the x-axis 

(𝛽) represents the pitch angle of the profile, y-axis (𝑚) represents the relative 



 

 

 
 

thickness of the profile and the z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. ................................................................................................ 58 

Figure 17. 3D response surfaces rendered for each individual TSR; where the x-axis 

(𝛽) represents the pitch angle of the profile-axis (𝑎𝑏) represents the maximum 

thickness position of the profile and z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. ................................................................................................ 58 

Figure 18. 2D response surfaces rendered for each individual TSR; where the x-axis 

(𝛽) represents the pitch angle of the profile, y-axis (𝑚) represents the relative 

thickness of the profile and the z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. ................................................................................................ 59 

Figure 19. 2D response surfaces rendered for each individual TSR; where the x-axis 

(𝛽) represents the pitch angle of the profile-axis (𝑎𝑏) represents the maximum 

thickness position of the profile and z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. ................................................................................................ 59 

Figure 20. Parameter optimization process flowchart for the optimal profile at TSR = 

2.33. .......................................................................................................................... 60 

Figure 21. Parameter optimization process flowchart for the optimal profile at TSR = 

2.64 ........................................................................................................................... 60 

Figure 22. Parameter optimization process flowchart for the optimal profile at TSR = 

3.09. .......................................................................................................................... 61 

Figure 23. Optimal profile and NACA 0021 profile. .................................................... 63 

Figure 24. Instantaneous torque coefficient (𝐶𝑡) for the NACA 0021 and the 

optimized profile at TSR = 2.33. ................................................................................ 64 

Figure 25. Dimensionless velocity contour for the NACA 0021 at TSR = 2.33 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. ....................... 65 

Figure 26. Dimensionless velocity contour for the optimalprofile atTSR= 2.33 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡 ........................ 65 



 

 

 
 

Figure 27. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR=2.33) at the azimuthal position of minimum torque 

coefficient 𝐶𝑡.............................................................................................................. 66 

Figure 28. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR =2.33) at the azimuthal position of maximum torque 

coefficient 𝐶𝑡.............................................................................................................. 67 

Figure 29. Instantaneous torque coefficient (𝐶𝑡) for the NACA 0021 and the optimized 

profile at TSR = 2.64. ................................................................................................ 68 

Figure 30. Dimensionless velocity contour for the for the NACA 0021 at TSR= 2.64 at 

the azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. ................. 69 

Figure 31. Dimensionless velocity contour for the optimal profile at TSR= 2.64 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. ....................... 69 

Figure 32. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR = 2.64) at the azimuthal position of minimum torque 

coefficient 𝐶𝑡.............................................................................................................. 70 

Figure 33. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR=2.64) at the azimuthal position of maximum torque 

coefficient 𝐶𝑡.............................................................................................................. 71 

Figure 34. Instantaneous torque coefficient (𝐶𝑡) for the NACA 0021 and the optimized 

profile at TSR = 3.09. ................................................................................................ 72 

Figure 35. Dimensionless velocity contour for the NACA 0021 at TSR = 3.09 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. ....................... 73 

Figure 36. Dimensionless velocity contour for the optimal profile at TSR = 3.09 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. ....................... 73 

Figure 37. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR=3.09) at the azimuthal position of minimum torque 

coefficient 𝐶𝑡.............................................................................................................. 74 



 

 

 
 

Figure 38. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR =3.09) at the azimuthal position of maximum torque 

coefficient 𝐶𝑡.............................................................................................................. 75 

Figure 39. Comparison between the NACA 0021 and the optimized profile .............. 76 

 

 



 

 

 
 

LIST OF TABLES 

Table 1. Main geometrical features of the 3-bladed NACA 0021 H-Darrieus turbine. 23 

Table 2. Summary of simulation settings. .................................................................. 29 

Table 3. Meshing metrics for orthogonal quality and skewness. ............................... 31 

Table 4. Details of the sensitivity analysis meshing and discretization errors............ 37 

Table 5. Lower bound and upper bound for the input variables ................................. 51 

Table 6. Summary of Meta-models created for TSR = 2.33. ..................................... 55 

Table 7. Summary of Meta-models created for TSR = 2.64. ..................................... 56 

Table 8. Summary of Meta-models created for TSR = 3.09. ..................................... 56 

Table 9. Summary of the best Meta-models for each specific TSR. .......................... 57 

Table 10. Detailed information about results and validation of the optimization 

models. ...................................................................................................................... 62 

 

 

 



 

 

 
 

ACRONYMS 

AoA Angle of Attack 

BoI Body of Influence 

CAE Computer-aided Engineer 

CFD Computational Fluid Dynamics 

CoD Coefficient of Determination 

CoI Coefficient of Importance 

CoP Coefficient of Prognosis 

CST Class Shape Transformation 

DoE Design of experiments 

EA Evolution Algorithm 

EP Evolution Programming 

ES Evolution Strategies 

GA Genetic Algorithm 

GCI Grid Convergence Index 

HAWT Horizontal Axis Wind Turbine 

LES Large Eddy Simulation. 

LHS Latin Hypercube Sampling 

LRM Linear Regression Model 

MCS Monte Carlo Simulation 

MLS Moving Least Squares 

MoE Margin of Error 

MOP Metamodel of Optimal Prognosis 

NLP Non-linear Programming 

NLPQL Non-linear Programming by Quadratic Langrangian 

PSO Particle Swarm Optimization 



 

 

 
 

RANS Reynolds Average Navier–Stokes 

SST Shear Stress Transport 

TSR Tip Speed Ratio 

URANS Unsteady Reynolds Averaged Navier–Stokes 

VAWT Vertical Axis Wind Turbine 

 



 

 

 
 

NOMENCLATURE 

𝛥𝐴𝑖 Area of the 𝑖th mesh element [𝑚2] 

𝜖𝑖 Error term [-] 

�̂�𝑖(𝑥𝑖) Approximation values [-] 

𝜆 Tip Speed Ratio [-] 

𝜇 Viscosity [𝑃𝑎 ⋅ 𝑠] 

𝐶𝑝 Average power coefficient [-] 

𝜌 Air density [𝑘𝑔/𝑚3] 

𝜎 Turbine solidity [-] 

𝜏𝑖𝑗 Strain rate [-] 

𝛼 Argument of the Joukowsky transformation airfoil [rad] 

𝜚 Vector of unidentified vector coefficients [-] 

휁 Joukowsky transformation imaginary plane [-] 

𝐴 Cross sectional area [𝑚2] 

𝑎 Joukowsky transformation ellipse major axis [rad] 

𝑎(𝑥𝑖) Moving Least Square’s coefficients [-] 

𝑏 Joukowsky transformation ellipse minor axis [rad] 

𝐶 Blade chord [m] 

𝐶𝐽 Real positive number [-] 

𝐶𝑝 Instantaneous power coefficient [-] 

𝐶𝑡 Torque coefficient [-] 

𝐷 Rotor diameter [m] 

𝑓(𝑥𝑖) Kriging trend function [-] 

𝐹𝐷 Drag force coefficient [-] 

𝐹𝐿 Lift force coefficient [-] 



 

 

 
 

𝐹𝑁 Normal force coefficient [-] 

𝐹𝑇 Tangential force coefficient [-] 

𝐺𝑏 Generation term for turbulent kinetic energy due to buoyancy 
effects [m2/s3] 

𝐺𝑘 Turbulent kinetic energy production term [m2/s3] 

𝐺𝑤 Production term for specific dissipation rate [1/s] 

𝐺𝑤𝑏 Generation term for specific dissipation rate due to buoyancy 
effects [1/s] 

ℎ Representative grid size [-] 

𝑀 Instantaneous moment [𝑁 ⋅ 𝑚] 

𝑚 Abscissa of the circle center [m] 

𝑁 Total of number cells [-] 

𝑛 Ordinate of the circle center [m] 

𝑁𝑏 Number of blades [-] 

𝑃 Turbine power [W] 

𝑝(𝑥) Basis of the polynomial regression function [-] 

𝑟 Grid refinement factor [-] 

𝑅 Radius of the blade [m] 

𝑅2 Coefficient of determination [-] 

𝑅𝐽 Radius of the Joukowsky transformation circle [m] 

𝑅𝑌,𝑋
2  Coefficient of determination of the model [-] 

𝑆𝑆𝐸  Regression model unexplained variation [-] 

𝑆𝑆𝑅 Variation brought on by the regression model [-] 

𝑆𝑆𝑇 Overall variance of the output [-] 

𝑈∞ Free stream velocity [m/s] 

𝑤𝑏 Tangential velocity of the blades [rad/s] 

𝑋 Real plane ordinate [m] 



 

 

 
 

𝑋𝑙 Lower airfoil ordinate [m] 

𝑋𝑢 Upper airfoil ordinate [m] 

𝑥𝑖 Input parameter [-] 

𝑌 Real plane abscissa [m] 

𝑌𝑙 Lower airfoil abscissa [m] 

𝑌𝑢 Upper airfoil abscissa [m] 

𝑌𝑘 Turbulent kinetic energy dissipation term [m2/s3] 

𝑌𝑤 Dissipation term for specific dissipation rate [1/s] 

𝑦𝑖 Model output [-] 

�̂�𝑖(𝑥𝑖) Sum of the approximation value [-] 

𝑧 Joukowsky transformation real plane [-] 

𝑍 CFD simulation result [-] 

𝑍(𝑥𝑖) Kriging stochastic function [-] 

𝑍𝑖 Optimization model predicted value [-] 

𝑅𝑌,𝑋
2  Coefficient of determination of the complete model [-] 

𝑅𝑌,𝑋~𝑖
2  Coefficient of determination of the reduced model [-] 

 

 

 

 



 

 

 
 

CONTENTS 

1 INTRODUCTION ........................................................................................................ 17 

1.1 H-DARRIEUS TURBINE ............................................................................................. 20 

2 LITERATURE REVIEW .............................................................................................. 24 

3 METHODOLOGY ....................................................................................................... 27 

3.1 GOVERNING EQUATIONS ........................................................................................ 27 

3.2 TURBULENCE MODEL AND SOLVER ...................................................................... 28 

3.3 FLUID DOMAIN .......................................................................................................... 30 

3.4 COMPUTATIONAL MESH .......................................................................................... 31 

3.5 REVOLUTION ANALYSIS .......................................................................................... 33 

3.6 AZIMUTHAL INCREMENT ANALYSIS ....................................................................... 35 

3.7 GRID CONVERGENCE ANALYSIS ............................................................................ 36 

3.8 VALIDATION STUDY ................................................................................................. 38 

3.9 SENSITIVITY ANALYSIS ........................................................................................... 40 

4 RESULTS ................................................................................................................... 55 

4.1 RESPONSE SURFACES AND ASSESSMENT .......................................................... 55 

4.2 OPTIMIZATION PROCESS SPECIFICS AND RESULTS ........................................... 61 

4.3 COMPARISON OF THE OPTIMAL PROFILE AT TSR = 2.33. ................................... 63 

4.4 COMPARISON OF THE OPTIMAL PROFILE AT TSR = 2.64 .................................... 67 

4.5 COMPARISON OF THE OPTIMAL PROFILE AT TSR = 3.09. ................................... 71 

5 CONCLUSIONS ......................................................................................................... 77 

6 RECOMMENDATIONS .............................................................................................. 79 

 REFERENCE ............................................................................................................. 80 

 



17 

 

 
 

1.  INTRODUCTION 

 

The rapid ever-growing energy consumption in recent years has hastened the 

depletion of fossil fuel resources. As a result, significant efforts have been made in 

the last decade to develop renewable energy alternatives. Wind energy, in particular, 

has received a lot of attention due to its sustainability, global availability, and lack of 

pollution  (Bazmi; Zahedi, 2011; Liu et al., 2013). 

The Horizontal Axis Wind Turbine (HAWT) and Vertical Axis Wind Turbine 

(VAWT) are the two most common types of turbines used in current projects. In 

general, all large-scale wind farms use HAWTs, which are becoming more efficient 

and larger in order to maximize the wind energy extracted (Bianchini et al., 2017). On 

the other hand, because VAWTs can capture wind energy from any direction, no yaw 

system is required to turn the rotor toward the wind direction; this also makes them 

suitable for urban areas where wind speed and direction oscillate (Qin et al., 2011). 

The H-Darrieus turbine is the most well-known and suitable for urban applications of 

the various types of VAWTs (Ahmadi-baloutaki; Mojtaba, 2016; Brownstein; Wei; 

Dabiri, 2019; Dabiri, 2011; Hansen; Mahak; Tzanakis, 2021). Furthermore, because 

they have fewer moving parts, their maintenance costs are lower, and they are less 

noisy because their operational speed is relatively slow. 

The airfoil shape directly affects the aerodynamic and structural characteristics 

of the blade and determines the pressure distribution of the entire flow field, which is 

the basis and key of flow analysis. Therefore, the expression method of the airfoil 

shape is very important, which affects the efficiency of the subsequent calculation 

process and the comprehensive performance of the airfoil. Many methods of 

parametrization of the blade profile have been developed over the years like the 

spline curve method, Bezier curve method, Hicks-Henne method, PARSEC method 

and CST (Class Shape Transformation) method. The spline curve method (Zhang et 

al., 2021)  can locally control and smooth the curve shape, but if the control points 

are not selected properly, wavy airfoils will easily appear. The Bézier curve method 

has the advantages of intuitiveness and convexity, but the Bézier curve cannot be 



18 

 

 
 

modified locally, and changing any control point position will affect the whole fitted 

curve (Wen et al., 2019). The Hicks-Henne (Ulaganathan et al., 2009) method has 

good control over the airfoil shape, but the variable range should not be too large; 

otherwise, the airfoil will appear unsmooth, which is suitable for the small 

adjustments of the airfoil. PARSEC method (Fujii; Dulikravich, 1998) has good 

robustness and a wave shape is unlikely to appear, but the ability to describe the 

details is poor, so it is suitable for the rough airfoil design. Furthermore, the Class 

Shape Transformation (CST) method can be used to describe a larger design space, 

but this method is less robust. To this end, the present research project will use the 

Joukowsky transformation (Zhang et al., 2021), a method that derives the function 

equation for airfoils (symmetrical in the case of this study) based on the sectional 

expression function of the wind turbine. Thus, the Joukowsky transformation presents 

two geometric parameters to be studied: the ratio 
𝑎

𝑏
 that controls the position of the 

maximum thickness position and 𝑚 that controls the relative thickness of the airfoil. 

In addition to the effect of the airfoil shape, VAWT suffers from one complex 

aerodynamic problem to operate at low values of Tip Speed Ratio (TSR), TSR < 4, 

related to its self-start capabilities. It is crucial then to have a good understanding of 

the starting process in particular at low Reynolds number which is appropriate for the 

urban application (Mălăel; Moutet; Drăgan, 2015). This phenomenon is called 

dynamic stall and it is induced by the large variations in the Angle of Attack  (AoA) on 

the blades during each revolution of the turbine at low TSR which has a great impact 

on the VAWT load and power output (Ferreira et al., 2008). When the dynamic stall 

occurs, flow separation is created at the leading-edge of airfoils, this dynamic stall 

can introduce structural vibrations which results in a critical drop in efficiency 

(Buchner et al., 2018). 

A particular method to address those issues is to vary the pitch angle of the 

blades in order to maximize turbine torque and make easier the start-up (Zhao et al., 

2018). A reduction of the stall is the main mechanism whereby variable-pitch turbines 

produce significant improvement of the torque at start-up and low TSRs (Pawsey et 

al., 2002). The limited research available about low TSR VAWTs indicates that the 

aerodynamic characteristics of the Darrieus VAWT can be improved using some 

strategies such as blade pitching (Ferreira et al., 2014; Fiedler; Tullis, 2009; Klimas; 



19 

 

 
 

Worstell, 1981; Rezaeiha; Kalkman; Blocken, 2017a). According to the foregoing, 

prior research examining the impact of fixed pitch angle has discovered that it is very 

promising, especially given that it is relatively simple to use in practice and does not 

incur significant manufacturing, installation, or maintenance costs. Therefore, the 

optimization of the low TSR H-Darrieus Turbine will be based on three parameters: 

𝑎/𝑏 (position of the maximum thickness), 𝑚 (relative thickness), and pitch angle (𝛽) in 

order to enhance its performance in terms of power generation and subsequently 

increase its adoption as a standalone wind energy generation source.  

Maximizing aerodynamic performance is a common concern when optimizing 

wind energy conversion systems. However, low-fidelity, simplified methods are often 

used to assess Darrieus aerodynamics as 2D CFD analyses of airflow through 

VAWTs are resource-intensive (Chen; Chiu; Fuge, 2020). The Metamodel of Optimal 

Prognosis (MOP) is used to address this issue and lower computational expenses. 

MOP is one of the techniques for design exploration frequently employed in 

stochastic analysis and nonlinear optimization (Sacks et al., 1989; Simpson et al., 

2004). The implementation of MOP after the DoE study produces a response surface 

that can be assessed based on the determination (CoD) (Montgomery; Runger, 

2014a), importance (CoI) (Bucher, 2007), and prognostic coefficients (CoP) (Most; 

Will, 2009) of each prospective MOP. 

In this work, our primary objective is to optimize the performance of the H-

Darrieus Vertical Axis Wind Turbines (VAWTs) using three parameters (a/b, m, and 

β) under an operational range of (TSR = 2.33, 2.64, and 3.09). We aim to (1) 

investigate the significance of these parameters on the turbine's aerodynamic 

performance at different operational conditions (TSR), (2) apply a Metamodel of 

Optimal Prognosis (MOP) to reduce computational expenses, and (3) perform an 

optimization using gradient methods to find the ideal profile for the operational range 

in terms of performance through numerical integration of the Cp - TSR curve. Overall, 

the novelty of our study lies in its comprehensive approach to optimizing low TSR H-

Darrieus VAWTs considering 3 different TSR, simultaneously, using three 

parameters (a/b, m, and β). This approach has not been thoroughly explored in 

previous research. By investigating the significance of these parameters, applying a 

Metamodel of Optimal Prognosis (MOP) to reduce computational expenses, and 



20 

 

 
 

performing an optimization using gradient methods to find the ideal profile, we offer a 

more efficient method for optimizing turbine performance in various operational 

conditions. 

 

1.1.  H-DARRIEUS TURBINE 

 

The H-Darrieus turbine is one of the most well-known of the VAWT. Figure 1 

depicts the H-Darrieus turbine, which consists of two or more airfoil-shaped blades 

attached to a rotating vertical shaft. 

 

Figure 1. Three bladed H-Type Darrieus turbine. 

Source: Prepared by the author. 

 

The blades of the H-Darrieus turbine experience a flow velocity over a range 

of Angle of Attack (AoA) as these blades rotate. Figure 2 shows the convention used 

for the aerodynamic coefficients (Francisco; Trentin, 2021) 

The lift (𝐹𝐿) and drag (𝐹𝐷) forces produce a resultant force that has a normal 

component (𝐹𝑁) and a tangential component (𝐹𝑇), which provides the torque: 



21 

 

 
 

𝐹𝑇 = 𝐹𝐿 sin(𝛼) − 𝐹𝐷𝑐𝑜𝑠(𝛼) 
(1) 

𝐹𝑁 = 𝐹𝐿 cos(𝛼) + 𝐹𝐷𝑠𝑖𝑛(𝛼) 
(2) 

 

Figure 2. Diagram depicting the aerodynamic coefficients of an H-Type Darrieus 
turbine. 

Source: (Francisco; Trentin, 2021). 

 

The TSR is an important parameter for quantifying continuous changes in the 

blade AoA, which is calculated by multiplying the tangential velocity of the blade (𝑤𝑏) 

by the radius (𝑅) relative to the wind speed. The TSR can be expressed as follows 

(Carrigan et al., 2012): 

TSR =
𝜔𝑅

𝑈∞
 

(3) 

 



22 

 

 
 

The second dimensionless parameter is solidity (𝜎) which is defined as the 

ratio of projected rotor blade surface area (𝑁𝑏𝐻𝐶) to the swept area (𝐻𝐷) of the wind 

turbine and is given by the following expression (Shengmao Li; Yan Li, 2010): 

𝜎 =
𝑁𝑏𝐶

𝐷
 

(4) 

where 𝑁𝑏 is the number of blades, 𝐶 is the blade chord length in meters and 𝐷 is the 

diameter of the rotor in meters. 

 

The most frequently used parameters to evaluate the performance 

characteristics of VAWT are the instantaneous moment coefficient 𝐶𝑡 and the power 

coefficient 𝐶𝑃. The moment coefficient 𝐶𝑡 is defined as: 

𝐶𝑡 =
𝑀

1
2𝜌𝐴𝑅U∞

2
 (5) 

where 𝑀 is the instantaneous moment, 𝜌 is the air density, 𝐴 is swept area of the 

turbine (𝜋𝑅2), 𝑅 is the turbine radius, and 𝑈∞ the freestream velocity.  

 

The power coefficient 𝐶𝑃 is defined as: 

𝐶𝑃 =
𝑃

1
2𝜌𝐴U∞

3
 (6) 

where P is the power of the turbine, and it is defined as: 

𝑃 = 𝑀𝜔 
(7) 

𝐶𝑃 =
𝐶𝑡𝑅𝜔

𝑈∞
= 𝐶𝑚𝜆 

(8) 



23 

 

 
 

The equation (8) correlates the three most important parameters: TSR, torque 

coefficient (𝐶𝑡), and power coefficient (𝐶𝑝) in a simple equation. The rotor 

geometrical specifications are reported in Table 1. 

 

Table 1. Main geometrical features of the 3-bladed NACA 0021 H-Darrieus turbine. 

Characteristics Turbine 

Number of blades, 𝑛 [−] 3 

Diameter, 𝐷 [𝑚] 1.03 

Height, H [m] 1 

Swept area, 𝐴 [𝑚2] 1.03 

Solidity, 𝜎 [−] 0.5 

Airfoil NACA 0021 

Airfoil chord, 𝐶 [𝑚𝑚] 85.8 

Tip speed ratio, 𝜆 [−] 2.64 

Free stream velocity, [𝑚∕𝑠] 9 

Rotational speed, [𝑟𝑎𝑑∕𝑠] 46.14 

Max. thickness position, [−] 30% 

Profile thickness, [−] 21% 

Source: (Raciti Castelli; Englaro; Benini, 2011). 

 

 



24 

 

 
 

2.  LITERATURE REVIEW 

 

The application of computational simulation tools, such as CFD, has proven to 

be an excellent analysis option for H-Darrieus turbines (Ahmadi-Baloutaki; Mojtaba, 

2016; Brownstein; Wei; Dabiri, 2019; Dabiri, 2011; Hansen; Mahak; Tzanakis, 2021). 

Analysis for the Darrieus turbine with aerodynamic profile like NACA 4 digits 

family are very common, especially symmetrical profiles, and its experimental data 

it’s easy to find for numerical model validation (Ferreira et al., 2014; Fiedler; Tullis, 

2009; Klimas; Worstell, 1981; Rezaeiha; Kalkman; Blocken, 2017a). Extensive work 

has been put on the development on well-established recommendations for 

computational domain size, spatial and temporal discretization (Hansen; Mahak; 

Tzanakis, 2021; Mălăel; Moutet; Drăgan, 2015). 

It is possible to find works with the analysis of 2 and 3 blade turbines for 

different values of TSR and different NACA profiles. Nonetheless, these works don`t 

present details of the relationship between each individual input parameter over the 

average power coefficient 𝐶𝑝, as well as their possible interactions between input 

parameter, and limit themselves only to a comparative analysis between different 

profiles. (Raciti Castelli; Englaro; Benini, 2011) does a study for a 3 bladed H-

Darrieus turbine with a NACA 0021 profile (Further details of the turbine can be found 

on Table 1). Their work presents experimental data that is widely used by several 

authors for their CFD validation (Francisco; Trentin, 2021; Hashem; Mohamed, 

2018), and will be used as a based case for this thesis as well. 

It is common to find works in the literature that study a single configuration of 

symmetric H-Darrieus turbines like NACA 0012, 0015, 0018 and 0021 (Saqib 

Hameed; Shahid, 2012). (Roy et al., 2017) goes further and study 12 aerodynamic 

profiles, 4 of them being symmetrical and NACA 4 digits profiles, 4 NACA 5-digit 

profiles and 4 Selig profiles. (Hashem; Mohamed, 2018) perform one of the studies 

with the largest sample of 3bladed H-Darrieus turbines, taking into account the 

aerodynamic analysis with 24 different aerodynamic profiles.(Ghasemian; Ashrafi; 

Sedaghat, 2017) presents a study containing 20 possible airfoils configurations and 



25 

 

 
 

creates a discussion on the effects of thickness and asymmetry on the performance 

of H-Darrieus turbines. This discussion sets the way for deeper examinations of how 

geometry affects the performance of the. H-Darrieus turbines 

All previous works study the behavior of the symmetric NACA 4 digits, but the 

geometry of these NACA profiles are limited to one parameter: Relative thickness. 

The latest 2 digits of the NACA profile stablishes the profile thickness. Since only one 

parameter is necessary to build the airfoil, the maximum thickness position (with 

reference to the chord length) is always limited to be at 30% of the airfoil. Works like 

(Ma et al., 2018) introduces a Multi-Island Genetic Algorithm to optimize the airfoil. 

The results show that the maximum thickness of the optimized profile increases by 

12.8% compared to NACA0018 based profile and the position of the maximum 

thickness moves from 30 % to 27.5% of the airfoil chord that results in an increment 

20% in the average 𝐶𝑝 for a TSR no larger than 1. As it is shown in this study, the 

position of the maximum thickness of the airfoil is worth studying, nevertheless the 

parametrization method used in this study contemplates 12 design variables. In order 

to reduce the computational time and the design variable another parametrization 

method is required. 

Works like (Fujii; Dulikravich, 1998; Ulaganathan et al., 2009; Zhang et al., 

2021) present different airfoil shape parametrization methods, but since the objective 

is to have control over position of maximum thickness with the lower design variables 

possible the Joukowski transformation is quite pertaining. (Zhang et al., 2021) 

presents an airfoil parametrization method based on the Joukowski transformation 

that allows us to control the relative thickness and the position of maximum thickness 

of a symmetric profile with just the addition of 1 more variable. 

Another important variable that increases the average 𝐶𝑝  according to works 

like: (Fiedler; Tullis, 2009; Klimas; Worstell, 1981; Rezaeiha; Kalkman; Blocken, 

2017a) is the pitch angle 𝛽. A 5m H-Darrieus turbine's average 𝐶𝑝 was examined in 

an open-field experiment by (Klimas; Worstell, 1981) for pitch angles ranging from -7° 

to 3°, with -2° being the ideal value. Similar to this, in an open jet wind tunnel, 

(Fiedler; Tullis, 2009) examined the effects of the pitch angles: -7.8°, -3.9°, 0°, 3.9°, 

and 7.8°. The study found that a negative pitch angle of -3.9° provides the optimum 



26 

 

 
 

performance. In a different experiment, (Rezaeiha; Kalkman; Blocken, 2017a) 

investigated pitch angles ranging from -7° to +3°, the results show that a pitch angle 

of -2° and a TSR of 4 can result in a 6.6% increase in 𝐶𝑝. 

This work, in addition to studying aspects relative thickness that a NACA-4-

digit parameterization provides control over, also analyzes the position of maximum 

thickness of the symmetric airfoil trough the used of the Joukowski transformation. 

Another important parameter studied in this thesis is the pitch angle. Another novelty 

of this work concerns the sensitivity analysis of the three input parameters using the 

Meta-model of Optimal Prognosis (MOP) methodology, explained in latter sections. 

This methodology has not been used in studies of H-Darrieus turbines, and therefore 

needs further investigation.  

 

 

 



27 

 

 
 

3.  METHODOLOGY 

 

3.1.  GOVERNING EQUATIONS 

 

The continuity and momentum equations are the governing equations for 

resolving the flow field over the fluid domains (Wilson; Korakianitis, 2014). For 

incompressible and transient state flow, the continuity equation is expressed in its 

conservation form as follows: 

𝜕𝜌

𝜕𝑡
+

𝜕(𝜌𝑢𝑖)

𝜕𝑥𝑖
= 0 

(9) 

The Navier–Stokes’ equations, also referred to as the momentum equations, 

have the following form: 

𝜕(𝜌𝑢𝑖)

𝜕𝑡
+

𝜕(𝜌𝑢𝑖𝑢𝑗)

𝜕𝑡
= −

𝜕𝑝

𝜕𝑥𝑖
+

𝜕𝜏𝑖𝑗

𝜕𝑥𝑗
+ 𝑆𝑀 

(10) 

where the strain rate (𝜏𝑖𝑗) is related to the stress tensor by: 

𝜏𝑖𝑗 = 𝜇 [(
𝜕𝑢𝑖

𝜕𝑥𝑗
+

𝜕𝑢𝑗

𝜕𝑥𝑖
) −

2

3
𝛿𝑖𝑗

𝜕𝑢𝑗

𝜕𝑥𝑖

] 
(11) 

The Navier-Stokes’s equations are transformed to the Reynolds Average 

Navier Stokes (RANS) by decomposing the instantaneous velocity to a mean 

component. This decomposition gives us an extra term called Reynold stress. The 

Reynold stress is modeled using different turbulence models (Versteeg; 

Malalasekera, 2007). 

  

 



28 

 

 
 

3.2.  TURBULENCE MODEL AND SOLVER 

 

The k-ω model and the standard k-ε model are both eddy viscosity turbulence 

models commonly used to simulate turbulent flow in fluid dynamics. However, there 

are some differences between these two models. The k-ω model is specifically 

designed for flows with adverse pressure gradients, such as those caused by the 

suction effect of the blade in a wind turbine (Menter, 1992). In this type of flow, the 

pressure decreases as the velocity increases while the flow hits the edge of the 

blade, creating adverse pressure gradients. The standard k-ε model, on the other 

hand, tends to overpredict the level of turbulent shear stress in these types of flows 

(Versteeg; Malalasekera, 2007). The standard k-ε model is more sensitive to the 

free-stream value of turbulence properties, while the k-ω model tends to depend on 

free-stream turbulence properties. Additionally, the k-ε model requires a complex 

nonlinear damping function for near-wall modifications, while the k-ω model uses a 

simpler approach (Wilson; Korakianitis, 2014). The k-ω model is suggested to be 

more accurate in predicting the level of turbulent shear stress in adverse pressure 

gradient flows. 

Since both models have their strengths and weaknesses, a hybrid model was 

developed by (Menter, 1994). This model combines the most accurate regimes of the 

𝑘 − 휀 and the 𝑘 − 𝜔 model and resulted in the two-equation eddy viscosity shear 

stress transport (𝑆𝑆𝑇) model. In this hybrid model, the 𝑘 − 𝜔 model is used in the 

near-wall regions and transforms to a 𝑘 − 휀  in regions far from the walls and in fully 

turbulent regions. With this hybrid treatment for regions near and far from the walls, 

the 𝑆𝑆𝑇 𝑘 − 𝜔 model becomes an appropriate choice for adverse pressure gradients. 

The 𝑘 equation in the 𝑆𝑆𝑇 𝑘 − 𝜔 model, with the buoyancy term disregarded, is 

written as followed (Menter, 1994) : 

𝜕(𝜌𝑘)

𝜕𝑡
+

𝜕(𝜌𝑈𝑗𝑘)

𝜕𝑥𝑗
=

𝜕

𝜕𝑥𝑗

[(𝜇 +
𝜇𝑡

𝜎𝑘3
)

𝜕𝑘

𝜕𝑥𝑗

] + 𝑃𝑘 − 𝛽∗𝜌𝑘𝜔 
(12) 

where the production term 𝑃𝑘 is defined as: 



29 

 

 
 

𝑃𝑘 = (2𝜈𝑡𝑠𝑖𝑗𝑠𝑖𝑗 −
2

3
𝑘

𝜕𝑢𝑖

𝜕𝑥𝑗
) 

(13) 

And the 𝜔 can be written as: 

𝜕(𝜌𝜔)

𝜕𝑡
+

𝜕(𝜌𝑈𝑗𝜔)

𝜕𝑥𝑗
=

𝜕

𝜕𝑥𝑗
[(𝜇 +

𝜇𝑡

𝜎𝜔2
)
𝜕𝜔

𝜕𝑥𝑗
]+ (1 − 𝐹1)2𝜌

1

𝜎𝜔2𝜔 

𝜕𝑘

𝜕𝑥𝑗

𝜕𝜔

𝜕𝑥𝑗
+ 𝛼3

𝜔

𝑘
𝑃𝑘 − 𝛽3𝜌𝜔2 

(14) 

 

Finally, Table 2 summarizes the main simulation settings in light of what has 

been stated in this section. 

 

Table 2. Summary of simulation settings. 

Characteristics Turbine 

Fluid properties Incompressible 

Turbulence model k-𝜔 SST 

Solver Pressure based 

Gradient Least Square Cell Based 

Momentum Second Order Upwind 

Turbulent Kinetic Energy Second Order Upwind 

Specific Dissipation Rate Second Order Upwind 

Pressure-Velocity solution Coupled 

Maximum residuals 1e-04 

Source: Prepared by the author. 

 

 



30 

 

 
 

3.3.  FLUID DOMAIN  

 

The current study aims to investigate the performance of a low TSR H-

Darrieus wind turbine operating at a determined free stream velocity (as illustrated in 

Figure 3). To ensure accurate and reliable results, the mesh quality was thoroughly 

evaluated. The dimensions of the computational domain were determined based on 

well-established recommendations with respect to the inner diameter (𝐷) of the 

turbine (Hansen; Mahak; Tzanakis, 2021; Mălăel; Moutet; Drăgan, 2015). Figure 3 

also depicts the boundary conditions for the outer domain (stationary domain) and 

the inner domain (rotational domain). The adopted boundary conditions were 

constant, uniform velocity and 6.3% turbulence intensity at the inlet, constant gauge 

pressure at the outlet, symmetry at the bottom and upper face, and no–slip at blade 

walls. 

 

Figure 3. Fluid domain dimensions and boundary conditions  

Source: (HANSEN; MAHAK; TZANAKIS, 2021). 

 

 

 



31 

 

 
 

3.4.  COMPUTATIONAL MESH 

 

Meshing plays a significant role in CFD modeling and significantly affects 

computation time. Due to its easy mesh generation for twisted geometries, an 

unstructured mesh was selected for the 2-D H-Darrieus Turbine, as seen in Figure 4. 

The outer domain is composed of 2 Body of Influence (BoI), and a growth factor of 

1.2 was set for the refinement of these 3 surfaces as shown in Figure 4. The inner 

domain is the fluid area simulating the revolution of the H-Darrieus Turbine, and the 

mesh size on both sides of the interface (outer domain and inner domain) has the 

same size in order to achieve faster convergence in the numerical solution (Launder, 

1995; Raciti Castelli; Englaro; Benini, 2011). 

The mesh around the NACA 0021 (Figure 5.a) profile is set to be an inflation 

with 30 layers, a growth ratio of 1.2, and a first cell with a 𝑦+  ≈ 0.1 (Figure 5.c). 

Additionally, mesh quality was evaluated based on two criteria: orthogonal 

quality (closeness of the angle between adjacent element faces) and skewness 

(closeness to the ideal equilateral or equiangular face). Table 3 shows meshing 

metrics under the aforementioned criteria. Ansys software was used to carry out the 

numerical modeling for a total of 285 geometry samples. 

 

Table 3. Meshing metrics for orthogonal quality and skewness. 

Orthogonal quality Skewness 

Minimum 0.11345 Minimum 2.1423e-02 

Maximum 0.99457 Maximum 0.98231 

Average 0.95897 Average 1.1134e-2 

Standard deviation 5.345e-2 Standard deviation 2.035e-2 

Source: Prepared by the author. 



32 

 

 
 

Figure 4. Global mesh of the computational domain (outer and inner domain). 

Source: Prepared by the author. 

 

Figure 5. Inner domain, airfoil leading and trailing edge mesh. 

Source: Prepared by the author. 



33 

 

 
 

In mesh generation for computational fluid dynamics simulations, a tapered 

trailing edge refers to a type of mesh refinement that is applied to the trailing edge of 

an airfoil or blade. The trailing edge of a blade is typically sharp and thin, which can 

cause numerical instabilities and convergence issues in simulations. To avoid these 

problems, a tapered trailing edge mesh refinement is used to gradually increase the 

thickness of the mesh elements towards the trailing edge as shown in Figure 5d. 

 

3.5.  REVOLUTION ANALYSIS 

 

The simulation for the reference scenario (NACA 0021, TSR=2.64) is run for a 

total of 90 revolutions in order to determine how many revolutions are necessary to 

provide a statistically steady flow field. The change in average 𝐶𝑝 between two 

successive revolutions of the turbine falls below 0.28% after 20 revolutions, and the 

difference between the values at 20 revolutions and at 90 revolutions is 0.58%, 

according to an analysis of the time history of the average 𝐶𝑝 (Figure 6 and 8) 

calculated over 1 revolution of the turbine. After 20 revolutions (Figure 7), the 

simulations are judged to have reached a statically steady condition. While this 

number of revolutions may be higher than the value used in previous studies, such 

as (Trivellato; Raciti Castelli, 2014), it is consistent with the convergence criteria 

reported by (Balduzzi et al., 2016; Lam; Peng, 2016) findings. 



34 

 

 
 

Figure 6. Average power coefficient 𝐶𝑝
̅̅ ̅ for the first 90 revolutions 

Source: Prepared by the author. 

 

Figure 7. Average power coefficient 𝐶𝑝
̅̅ ̅ after 20 revolutions 

 

Source: Prepared by the author. 



35 

 

 
 

Figure 8. Relative change of the average power coefficient 𝐶𝑝
̅̅ ̅ with respect to 

the 𝐶𝑝
̅̅ ̅ at the 90th revolution 

 

Source: Prepared by the author. 

 

3.6.  AZIMUTHAL INCREMENT ANALYSIS 

 

The azimuth angle (Δ𝜃) represents the change in position of the wind turbine 

rotor as it rotates over time. The azimuth angle is directly related to the time step 

used in the computational fluid dynamics (CFD) simulation, with a lower azimuth 

angle corresponding to a smaller time step. This relationship between the azimuth 

angle and the time step is due to the fact that the turbine rotation is modeled using a 

moving reference frame in the CFD simulation. By varying the azimuth angle, the 

simulation captures the transient behavior of the turbine under different operating 

conditions, which is critical for accurately predicting its performance. Therefore, the 

plot in Figure 9 provides a way to analyze how the power output of the turbine 

changes over time and how this is influenced by changes in the azimuth angle. 

These results demonstrate that the moments on the blades cannot be largely 

predicted by simulations with Δ𝜃 values of 10° and 5°. Similarly, Δ𝜃 value 2.5° 



36 

 

 
 

frequently underestimates the value of average 𝐶𝑝, albeit by a considerably lesser 

margin. The difference in average 𝐶𝑝  for Δ𝜃 values of 1°, 0.5° and 0.25° is negligible. 

This comparison reveals that choosing Δ𝜃 = 1° at the specified TSR is a secure 

option. It is important to keep in mind that this value is restricted to flow conditions 

with weak flow separation at moderate TSRs (Rezaeiha et al., 2017a). 

 

Figure 9. History of power coefficient 𝐶𝑝
̅̅ ̅ for various azimuthal increments 

Source: Prepared by the author. 

 

3.7.  GRID CONVERGENCE ANALYSIS 

 

The NACA 0021 airfoil was taken as a reference case and examined for three 

different types of mesh: refined, medium, and course in order to gauge the 

independence of the results as a function of the degree of mesh refinement. This 

makes it possible to measure the discretization error using Celik’s suggested Grid 

Convergence Index (GCI) approach (Celik et al., 2008). 



37 

 

 
 

ℎ = [
1

𝑁
∑(∆𝐴𝑖)

𝑁

𝑖=1

]

1/2

 (15) 

 

where ℎ is the grid size, ∆𝐴𝑖 is the area of the 𝑖th mesh element and 𝑁 is the total 

number of cells used for computations.  

Then, simulations are run to determine the values of the objective variable 

average power coefficient (𝐶𝑝1, 𝐶𝑝2, 𝐶𝑝3), where the grid refinement factor (𝑟) must be 

greater than 1.3 based on empirical knowledge (Francisco; Trentin, 2021). Following 

these procedures, simulations were run to determine the values for the average 

power coefficient which are used to find 𝐺𝐶𝐼32 using Celik’s analytical procedure. 

Table 4 shows the meshing specifics obtained for the sensitive analysis and 

discretization error. The numerical error in the solution for the medium mesh 𝐺𝐶𝐼32 is, 

therefore, 0.1903%, as shown in Table 4. For this reason, using the medium mesh 

for CFD simulations with a small amount of uncertainty is acceptable. 

 

Table 4. Details of the sensitivity analysis meshing and discretization errors. 

𝑁1 395145 

𝑁2 224220 

𝑁3 119665 

𝑟21 1.328 

𝑟32 1.369 

𝐶𝑝1
̅̅ ̅̅̅ 3.0170 

𝐶𝑝2
̅̅ ̅̅̅ 3.0227 

𝐶𝑝3
̅̅ ̅̅̅ 3.1557 

𝐺𝐶𝐼32 0.1903% 

Source: Prepared by the author. 



38 

 

 
 

3.8.  VALIDATION STUDY 

 

The validation study focuses on a 3-bladed H-Darrieus wind turbine with a 

symmetric NACA 0021 airfoil. The turbine's geometrical specifications are listed in 

Table 1. The study investigates the turbine's performance over a range of tip speed 

ratios (TSR) between 1.44 to 3.29, with a free stream velocity of 9 m/s. As 

demonstrated in a previous study by (Belabes; Paraschivoiu, 2021) , flow turbulence 

intensity can have a significant impact on the performance of vertical axis wind 

turbines (VAWT). In this sense, for this study, a 6.3 % turbulent intensity at the inlet 

was used as it presents the best fit for the experimental data provided by (Raciti 

Castelli; Englaro; Benini, 2011). 

For different tip speed ratios, the validation study compares the calculated 

turbine average 𝐶𝑝 to the measured average 𝐶𝑝 values by (Raciti Castelli; Englaro; 

Benini, 2011). As shown in Figure 10 and 11, there is good agreement between the 

CFD results and the experimental data. The following observations can be made 

about the comparison: 

 The current results are much more consistent with the experimental data 

provided by (Raciti Castelli; Englaro; Benini, 2011) than previous studies 

like (Francisco; Trentin, 2021; Hashem; Mohamed, 2018; Raciti Castelli; 

Englaro; Benini, 2011). The lower observed discrepancies can be 

attributed to the fact that the current study adheres to the domain size 

guidelines provided by (Hansen; Mahak; Tzanakis, 2021) as well as the 

appropriate 𝐺𝐶𝐼32 analysis (section 3.7), the minimum requirements for 

the number of turbine revolutions and azimuthal increment shown in 

sections 3.5 and 3.6. As well as an appropriate inlet turbulent intensity. 

 Due to geometrical simplifications, an overestimation of calculated 

average  𝐶𝑝 (compared to measured values) is expected for low to 

moderate tip speed ratios (where dynamic stall is present (Ferreira et al., 

2008). However, there is a similar overestimation of calculated average 

 𝐶𝑝 for moderate to high TSRs. This can be attributed in part to the 



39 

 

 
 

absence of blade tip losses in 2D CFD simulations (Rezaeiha; Kalkman; 

Blocken, 2017a). 

 

Figure 10. Comparison of experimental and numerical validation results for a small-
scale NACA 0021 Darrieus VAWT operating at TSR = 2.64 and 9 m/s free stream 

velocity. 

Source: (Francisco; Trentin, 2021; Hashem; Mohamed, 2018; Raciti Castelli; Englaro; Benini, 2011). 

Figure 11. Comparison between the k-w model and the experimental data presented 

Source: (Raciti Castelli; Englaro; Benini, 2011) 



40 

 

 
 

3.9.  SENSITIVITY ANALYSIS 

 

The following methodology is used to optimize the VAWT. It is based on the 

input parameter identification, Design of experiments (DoE), the Metamodel of 

Optimal Prognosis (MOP), and gradient-based and nature-based optimization 

models. A flowchart outlining the methodology’s steps is shown in Figure 12. These 

steps are presented in the following subsections. 

 

Figure 12. Flowchart representing the full plan for CFD modeling, Sensitivity Analysis, 
and Optimization 

 

Source: Adapted from (Francisco; Trentin, 2021) 

 

3.9.1.  Geometry parameterization by means of the Joukowsky transformation 

 

The Joukowsky Transformation is a conformal mapping method that converts a 

circle into an airfoil shape in the complex number plane to convey various properties. 

A conformal transformation of this kind is one in which none of the analytic functions’ 



41 

 

 
 

derivatives are zero. 

The Joukowsky transformation has the following analytical equation: 

𝑧 = 휁 +
𝐶𝐽

2

휁
 (16) 

where 𝐶𝐽 is a real positive number represented by 𝐶𝐽 = 2(𝑎−𝑚). 𝑚 is the abscissa and 

𝑛 is the ordinate of the circle’s center in the second quadrant as shown in Figure 13.  

 

In order to create an airfoil with a specific relative thickness and relative camber 

on a plane 𝑧, the circle with radius 𝑅𝐽 is changed by the Joukowsky transformation on 

the plane 휁, as shown in Figure 13. 

Polar coordinates are used to indicate the contour points on the plane 휁: 

휁 = 𝑅𝐽𝑒
𝑖𝛼 (16) 

The airfoil line on the plane 𝑧 has the following parameter equation, according 

to (Zhang et al., 2021): 

𝑋 =

[
 
 
 

𝑅𝐽 +
(
𝐶𝐽

2 )
2

𝑅𝐽

]
 
 
 

𝑐𝑜𝑠(𝛼) 

(17) 

𝑌 =

[
 
 
 

𝑅𝐽 −
(
𝐶𝐽

2 )
2

𝑅𝐽

]
 
 
 

𝑠𝑖𝑛(𝛼) 

(18) 

 

The airfoil coordinate is (𝑋, 𝑌), the argument is 𝛼. The airfoil's relative thickness 

and camber are controlled by parameters 𝑚 and 𝑛, respectively (Figure 13). 

Changing 𝑛 widens the airfoil's relative camber, while changing 𝑚 increases its 



42 

 

 
 

relative thickness. However, the maximum thickness position is always at 25% of the 

chord and the maximum camber position is always at 50% of the chord, regardless of 

the parameter values. The positions of the relative thickness and camber are also 

constant. To address this, the circle is transformed into a spinning elliptical shape. 

The Joukowsky transformation creates an airfoil with variable thickness by 

shifting an ellipse along the real axis (𝑚). The position of the maximum thickness 

depends on the ratio of the major (𝑎) and minor (𝑏) axis of the ellipse. This allows for 

the creation of a symmetric airfoil with different relative thickness and maximum 

thickness positions by adjusting 𝑚. 

Figure 13. Joukowski transformation. 

 

Source: (Zhang et al., 2021) 

 

According to the established method in (ZHANG et al., 2021), the airfoil’s 

parametric equation for a symmetric profile is as follows: 

Ω = √(𝑎𝑐𝑜𝑠(𝜃))
2
+ (𝑏𝑠𝑖𝑛(𝜃))

2
 (19) 



43 

 

 
 

𝑋𝑢 =

[
 
 
 

Ω +
(
𝐶𝐽

2
)

2

Ω

]
 
 
 

cos(𝛼𝑢) , 𝛼𝑢 𝜖 [0, 𝜋] 

(20) 

𝑌𝑢 =

[
 
 
 

Ω −
(
𝐶𝐽

2 )
2

Ω

]
 
 
 

sin(𝛼𝑢) , 𝛼𝑢 𝜖 [0, 𝜋] 

(21) 

𝑋𝑙 =

[
 
 
 

Ω +
(
𝐶𝐽

2 )
2

Ω

]
 
 
 

cos(𝛼𝑙) , 𝛼𝑙 𝜖 [𝜋, 2𝜋] 

(22) 

𝑌𝑙 =

[
 
 
 

Ω −
(
𝐶𝐽

2 )
2

Ω

]
 
 
 

sin(𝛼𝑙) , 𝛼𝑙 𝜖 [𝜋, 2𝜋] 

(23) 

𝛼𝑢 = 𝜋 − 𝑎𝑟𝑐𝑐𝑜𝑠 (
𝑚 − 𝑎𝑐𝑜𝑠(𝜃)

Ω
) , 𝛼𝑢 𝜖 [0, 𝜋] 

(24) 

𝛼𝑙 = 𝜋 + 𝑎𝑟𝑐𝑐𝑜𝑠 (
𝑚 − 𝑎𝑐𝑜𝑠(𝜃)

Ω
) , 𝛼𝑙 𝜖 [𝜋, 2𝜋] 

(25) 

 

The airfoil's relative thickness and maximum thickness position depend on the 

ratio of the ellipse's major and minor axes, 𝑎/𝑏. Figure 14 shows that when 𝑏<𝑎, the 

relative thickness tends to shift towards the leading edge, while for 𝑏>𝑎, it shifts 

towards the trailing edge. The parameter 𝑚 determines the overall thickness of the 

profile, while 𝑟 controls the position of the maximum thickness. 

 



44 

 

 
 

Figure 14. Contrast between different parameters configuration (𝑎, 𝑏 and 𝑚). 

Source: Prepared by the author. 

 

3.9.2.  Pitch angle 

 

The pitch angle is an essential parameter in wind turbines, as it refers to the 

angle at which the rotor blades are twisted or set along their chord length as shown 

in Figure 15. It determines the orientation of the blades relative to the rotational plane 

of the turbine and has a significant impact on how the blades interact with the wind 

flow. The pitch angle is set during the turbine's design and manufacturing process 

and remains constant during operation. By adjusting the pitch angle, wind turbine 

designers can optimize the turbine's performance in different wind conditions. 

In the context of Vertical Axis Wind Turbines (VAWTs), the pitch angle has 

been recognized as a key factor in improving turbine performance. Researchers have 

explored the impact of pitch angle using various models and simulation techniques. 

The Double Multiple Stream tube Model (DMSM) was used by numerous writers to 

examine the impact of pitch angle on VAWTs (Paraschivoiu, 2002). However, 

(Ferreira et al., 2008) emphasized the flaws in this low-fidelity model’s ability to 

forecast the intricate aerodynamics and performance of VAWTs. Works like (Fiedler; 



45 

 

 
 

Tullis, 2009; Klimas; Worstell, 1981; Rezaeiha; Kalkman; Blocken, 2017a) on the 

other hand, demonstrated the high fidelity of CFD simulation models, 

In an open-field experiment, (Klimas; Worstell, 1981) investigated the impact 

of the fixed pitch using high-fidelity CFD simulations for angles spanning from -7° to 

3° on the 𝐶𝑝 of a 5𝑚 Darrieus turbine, with -2° being the ideal value. Similar to 

this,(Fiedler; Tullis, 2009) investigated the impact of the pitch angles -7.8°, -3.9°, 0°, 

3.9° and 7.8° on 𝐶𝑝 for a high-solidity H-Darrieus turbine in an open jet wind tunnel. 

According to the study, a negative pitch angle offers (-3.9°) the best performance. In 

a different experiment, (Rezaeiha; Kalkman; Blocken, 2017a) used Unsteady 

Reynolds Averaged Navier–Stokes (URANS) computations to study an H-Darrieus 

turbine with pitch angles varying from -7° to +3°,. Turbulence was modeled using the 

4-equation transition SST model. According to the findings, a pitch angle of -2° and a 

TSR of 4 can lead to an increase in 𝐶𝑝 of 6.6%. From the aforementioned, it can be 

inferred that prior research examining the impact of fixed pitch angle has found that it 

is very promising, especially given that it is relatively straightforward to use in practice 

and does not incur significant manufacturing, installation, or maintenance expenses. 

The pitch angle 𝛽 (Figure 15) is therefore the third parameter selected for the VAWT 

optimization. 

 

Figure 15. Schematic showing the pitch angle 𝛽 for a VAWT blade 

Source: Prepared by the author. 

 



46 

 

 
 

3.9.3.  Design of Experiments (DoE) 

 

The space of the design variables must be scanned by either deterministic 

sampling schemes or stochastic sampling schemes in order to conduct a global 

sensitivity analysis. The Monte Carlo Simulation (MCS) is a technique that is 

frequently used for stochastic sampling (Homem-de-Mello; Bayraksan, 2014). The 

design variables in this technique are presumptively distributed uniformly with 

specified lower and upper bounds. Small sample sizes frequently result in clusters 

and gaps, which raises the possibility of unexpected correlations between input 

variables. To solve these issues, the comprehensive Latin Hypercube Sampling 

(LHS), which accurately represents input distribution and input correlation even for a 

limited number of samples, is applied (Simpson et al., 2004). 

A minimum of 80 to 120 samples are used to provide adequate estimations of 

the importance of variance-based input variables. The initial numbers of samples are 

adjusted to the solver run duration and the number of input parameters (The number 

of failed designs should not exceed 20%). The sensitivity analysis has the advantage 

of reducing the number of design variables and identifying the parameters that are 

most effective (CoI) in improving the optimization aim. 

 

3.9.4.  Metamodel of optimal prognosis (MOP) 

 

One of the most often used methods for design exploration in stochastic 

analysis and nonlinear optimization is meta-modeling (SIMPSON et al., 2004). It is 

not always possible to simplify the physical models in real-world applications of virtual 

prototyping while still obtaining numerical models that are easily solvable. Every 

numerical simulation typically takes hours or even days. Although numerical methods 

and high-performance computing have advanced, in certain situations it is impossible 

to investigate several model configurations; therefore, effective surrogate models are 

needed (Sacks et al., 1989; Simpson et al., 2004). 



47 

 

 
 

Following the CoI study’s conclusion, the application of the MOP creates a 

response surface, where physical plausibility must be evaluated based on the 

determination (CoD), importance (CoI), and prognostic coefficients (CoP) of each 

prospective MOP. The aforementioned MOPs are based on moving least square, 

Kriging, and polynomial regression approximations (Most; Will, 2010). 

 

3.9.5.  Coefficient of determination 

 

A polynomial regression model’s accuracy as an approximation can be 

evaluated using the Coefficient of Determination (CoD) (Montgomery; Runger, 

2014b). This metric represents the proportion of variation that the approximation can 

account for: 

𝑅2 =
𝑆𝑆𝑅

𝑆𝑆𝑇
= 1 −

𝑆𝑆𝐸

𝑆𝑆𝑇
, 0 ≤ 𝑅2 ≤ 1 

(17) 

where 𝑆𝑆𝐸  quantifies the unexplained variation, 𝑆𝑆𝑅 represents the variation brought 

on by the regression, and 𝑆𝑆𝑇 equals the overall variance of the output 𝑌. 

 

3.9.6.  Coefficient of importance 

 

Using the CoD metric, (Bucher, 2007) created the Coefficient of Importance 

(CoI) to gauge the significance of the input variable. The CoI of one variable 𝑋𝑖, with 

regard to the response 𝑌, is defined as follows using a polynomial model that 

incorporates all variables under investigation. 

𝐶𝑜𝐼(𝑋𝑖 , 𝑌) = 𝑅𝑌,𝑋
2 − 𝑅𝑌,𝑋~𝑖

2  
(18) 



48 

 

 
 

where 𝑅𝑌,𝑋~𝑖
2  is the CoD of the reduced model, where all linear, quadratic, and 

interaction terms belonging to 𝑋𝑖 are eliminated from the polynomial basis, and 𝑅𝑌,𝑋
2  is 

the CoD of the complete model, which includes all terms of the variables in 𝑋. The 

same set of sampling points is utilized in both situations. A variable’s CoI is almost 

zero if it is of little significance because both the full and reduced polynomial 

regression models have comparable strengths. Since the CoD quantifies the 

explained variation of the polynomial approximation, the CoI is equivalent to the 

explained variation with regard to a single input variable. The overall effect sensitivity 

measure from equation (18) is thus estimated. 

 

3.9.7.  Coefficient of prognosis 

 

A model-independent measure to evaluate the meta-model quality was 

developed in (MOST; WILL, 2009). The Coefficient of Prognosis (CoP), which 

represents this metric, is described as follows: 

𝐶𝑜𝑃 = 1 −
𝑆𝑆𝐸

𝑃𝑟𝑒𝑑𝑖𝑐𝑡𝑖𝑜𝑛

𝑆𝑆𝑇
 

(19) 

where 𝑆𝑆𝐸  stands for the total of prediction errors squared. Based on cross-

validation, these inaccuracies are estimated. The set of support points is converted to 

𝑞 subsets during the cross-validation process. After that, the approximation model is 

created by eliminating point 𝑖 from the support points and utilizing the remaining point 

set to approximate the subset model output 𝑦�̂�. This indicates that only those points 

that are not used to develop the approximation model are used to estimate the model 

quality. This method works with regression models and even interpolation models 

because the prediction error is used instead of the fit. 

 

 



49 

 

 
 

3.9.8.  Polynomial regression 

 

Polynomial regression is a popular approximation method in which the model 

response is generally approximated by a polynomial basis function of linear or 

quadratic order with or without coupling terms. The sum of the approximation value 

�̂�𝑖(𝑥𝑖) and an error term 𝜖𝑖 can be used to represent the model output 𝑦𝑖 for a given 

set 𝑥𝑖 of the input parameters 𝑥 (Freeny; Box; Draper, 1988; Myers, 2018). 

𝑦(𝑥𝑖) = �̂�𝑖(𝑥𝑖) + 𝜖𝑖 = 𝑝𝑇(𝑥𝑖)𝜚 + 𝜖𝑖 
(20) 

where 𝑝(𝑥) is the basis of the polynomial and 𝜚 is a vector that contains the 

unidentified regression coefficients. Typically, these coefficients are computed using 

a series of sampled support points under the assumption that each point has 

independent errors with equal variance (Myers, 2016). 

 

3.9.9.  Moving Least Squares approximation 

 

The Moving Least Squares (MLS) method was developed by (Lancaster; 

Salkauskas, 1981) and is a continuation of polynomial regression. In a similar vein, 

the basis function can include any kind of function, though often only linear and 

quadratic terms are included. By finding the best local fit for the actual interpolation 

point, this basis function can be accurately represented. It is said that the 

approximation function is: 

𝑦(𝑥𝑖) = 𝑝𝑇(𝑥𝑖)𝑎(𝑥𝑖) 
(21) 

Using variable ("moving") coefficients 𝑎(𝑥𝑖) as opposed to the global 

polynomial regression coefficients. 

 



50 

 

 
 

3.9.10. Kriging approximation 

 

Kriging was originally designed to model spatial variations in geological 

models that had previously been measured (Matheron, 1976). These models are 

inherently random, with Gaussian variation around a constant mean value. 

Additionally, the use of Kriging models for approximating optimization and stochastic 

problems is gaining popularity. 

Because it can be used to approximate a wide variety of complex response 

functions, kriging is a reliable surrogate model for this kind of application. 

Additionally, because it interpolates the support data points and provides a 

confidence interval about a forecast result of the approximation model, it is also 

appropriate for deterministic models. The use of a set of parameters to define the 

model gives this approach its versatility, however, choosing the appropriate set of 

parameters for a Kriging model has a few limitations. Constrained iterative searching 

is required to find these parameters, which, under the assumption of a Gaussian 

random process, is a computationally expensive procedure (Martin; Simpson, 2012). 

A regression model 𝑓(𝑥𝑖)𝜚 and a stochastic process 𝑍(𝑥𝑖) make up the 

response surface that is produced by the Kriging approximation (Cressie, 1990; Li et 

al., 2016). 

𝑦(𝑥𝑖) = 𝑓(𝑥𝑖)𝜚 + 𝑍(𝑥𝑖) 
(22) 

The trend function is represented by 𝑓(𝑥𝑖) and the accompanying coefficient 

vector, 𝜚. 

 

3.9.11. Optimization statement 

 

Two sources provide the input parameters for the CFD simulations: Input pitch 

angle (𝛽) and geometric parameterization using the Joukowsky transformation. 



51 

 

 
 

Two parameters are identified for the geometric parametrization of the 

symmetrical profile: the shifted distance of the ellipse along the real axis (𝑚) and the 

ratio between the ellipse major and minor axes (
𝑎

𝑏
). The position of the profile’s 

maximum thickness is controlled by (
𝑎

𝑏
), while the airfoil’s relative thickness is 

controlled by 𝑚. The lower bound and higher bound for the inputs are shown in Table 

5. 

Table 5. Lower bound and upper bound for the input variables 

Input variable Lower bound Upper bound 

(
𝑎

𝑏
) 0.95 1.09 

𝑚 0.03 0.053 

𝛽 -7 4 

Source: Prepared by the author. 

The range for the position of maximum thickness (
𝑎

𝑏
) goes from 20% to 40% 

while the range for the relative thickness (𝑚) goes from 15% to 33%. Taking into 

consideration that the original NACA 0021 has a maximum thickness position of 30% 

and a relative thickness of 21%, the ranges for these two variables cover an 

adequate working space (Table 1). As for the pitch angle (𝛽), from previous studies 

of the behavior of pitch angle over H-Darrieus turbine (Fiedler; Tullis, 2009; Klimas; 

Worstell, 1981; Rezaeiha; Kalkman; Blocken, 2017a), a range of -7° to 4° is quite 

appropriate. 

This study finds the optimal profile for an operational range in terms of 

performance through the numerical integration of the 𝐶𝑝
̅̅ ̅ − 𝑇𝑆𝑅 curve. The chosen 

values for the 𝑇𝑆𝑅 values are 2.33, 2.64 and 3.09 since they provide the highest 𝐶𝑝
̅̅ ̅ in 

the original NACA 0021 profile and there is available experimental data to validate 

the results. Gradient-based methods were used for unrestricted optimization. The 

optimization statement is: 



52 

 

 
 

 Find 𝑋 = {
𝑎

𝑏
, 𝑚, 𝛽}

𝑇

that  maximizes the ∫ 𝐶𝑝
̅̅ ̅𝑇𝑆𝑅𝑓

𝑇𝑆𝑅0
. 

 

 

3.9.12. Gradient-based optimization methods 

 

To locate the following local optimum, gradient-based optimization techniques 

use local derivatives of the objective function. A convex function’s minimum can be 

found by looking for the location where the first derivative equals zero. A second-

order Taylor series can be used to approximate the objective function if the second 

derivative is known. This method is known as Non-linear Programming (NLP). 

Newton’s methods are first- and second order derivative-based optimization 

techniques. The derivatives must be estimated using numerical estimates if the 

computer-aided Engineer (CAE) solution is viewed as a black box (Kelley, 1999). 

The next subsections outline the gradient-based optimization techniques that 

will be used for the current investigation. 

Nonlinear Programming by Quadratic Langrangian (NLPQL). The NLPQL 

method uses a quadratic Langrangian to efficiently handle constraints. The approach 

searches for the subsequent local optimum from a given starting point and converges 

if the estimated gradients are below a predetermined tolerance. It is advised to utilize 

the best design of a global sensitivity analysis as a starting point since it is a local 

optimization method in order to locate the global optimum. In low dimensions with up 

to 20 design variables, the NLPQL is quite efficient. The computation of numerical 

derivatives gets increasingly expensive for larger dimensional situations, and 

alternative methods are more effective. 

The differentiation interval is essential when there are noisy model responses. 

If it’s taken too tiny, the solver noise severely distorts the calculated gradient, sending 

the NLPQL in the wrong direction. An enhanced differentiation interval could result in 



53 

 

 
 

a decent convergence of the optimization process if solver noise is not overriding the 

functional trend (Schittkowski, 1986). 

Downhill Simplex Method (Simplex). The Downhill Simplex optimization 

method is an iterative technique that monitors 𝑛 +  1 points in 𝑛 dimensions, where 𝑛 

is the number of parameters that need to be specified for the machine learning 

process. These locations are regarded as the simplex’s vertices (triangle in 2D, a 

tetrahedron in 3D) (Fortran et al., 1986; Salcedo; Cândido, 2001). 

Every iteration replaces a vertex with the target function’s worst value. 

Depending on the newly computed value of the target function, the simplex may 

subsequently move into a new position in a variety of ways. The Downhill Simplex 

technique has a set number of iterations, which reduces the amount of time it takes 

to run. In the present scenario, we can additionally force the evaluation of the target 

function to end if it takes longer than the worst time thus far. When the final iteration 

is complete, the Downhill Simplex methods come to an end. The Machine Learning 

algorithm’s parameter values are the variables’ ultimate values (Raska; Ulrych, 

2015). 

 

3.9.13. Margin of Error (MoE) 

 

Specific CFD simulations are done to confirm the optimal control points and 

the related results anticipated by each gradient-based model in order to evaluate the 

consistency of the results. The margin of error (MoE) is determined using the formula 

(37), and the model’s output is considered valid if the MoE’s absolute value is less 

than 1%. In contrast, the model output is disregarded if the MoE is higher than the 

aforementioned threshold. 

|
𝑍 − 𝑍𝑖

𝑍
| = 𝑀𝑜𝐸 

(23) 



54 

 

 
 

where 𝑍 is the CFD simulation result and 𝑍𝑖 is the optimization model predicted 

output. 

 

3.9.14. Optimal profile 

 

The average (𝐶𝑃) is used to choose an appropriate shape following 

computational validation of each optimized design point. As a result, the best-

performing optimization model is used to generate the optimal profile based on 

design points predicted by the model for each TSR. 

 



55 

 

 
 

4.  RESULTS 

 

4.1.  RESPONSE SURFACES AND ASSESSMENT 

 

After advanced LHS-based sampling of design points, CFD simulations are 

conducted for each point to collect essential data, including velocity and pressure 

distribution, for building suitable meta-models for each 𝜆. Ansys Fluent is used for 

CFD simulations, and OptiSlang is used for sensitivity analysis, both working in 

concert on the Ansys Workbench platform. 

Table 6 showcases the meta-models constructed for TSR = 2.33. These 

models include Linear Regression Models of first and second order, Moving Least 

Squares of first and second order, and Kriging. Based on the Coefficient of 

Determination (CoD) and Coefficient of Performance (CoP), the second-order Linear 

Regression Model (LRM) emerges as the best fit for the three input parameters (a/b, 

m, β) at TSR = 2.33. This conclusion aligns with the results obtained for TSR = 2.64 

and 3.09, as demonstrated in Tables 7 and 8. 

 

Table 6. Summary of Meta-models created for TSR = 2.33. 

Meta-models Approximations Coefficient of 
Determination 
(CoD) 

Coefficient of 
Prognosis (CoP) 

LRM of order 1  0.755571 0.770465 

LRM of order 2  0.999939 0.999962 

MLS of order 1 0.974276 0.944808 

MLS of order 2 0.985157 0.995001 

Kriging  0.979047 0.995323 

Source: Prepared by the author. 



56 

 

 
 

Table 7. Summary of Meta-models created for TSR = 2.64. 

Meta-models Approximations Coefficient of 
Determination 
(CoD) 

Coefficient of 
Prognosis (CoP) 

LRM of order 1  0.781969 0.792290 

LRM of order 2  0.999922 0.999082 

MLS of order 1 0.938770 0.926734 

MLS of order 2 0.951678 0.973970 

Kriging  0.950395 0.967435 

Source: Prepared by the author. 

 

Table 8. Summary of Meta-models created for TSR = 3.09. 

Meta-models Approximations Coefficient of 
Determination 
(CoD) 

Coefficient of 
Prognosis (CoP) 

LRM of order 1  0.766882 0.782424 

LRM of order 2  0.999289 0.999640 

MLS of order 1 0.979902 0.951639 

MLS of order 2 0.984791 0.986444 

Kriging  0.990001 0.997038 

Source: Prepared by the author. 

 

Table 9 provides a detailed overview of the three selected second-order Linear 

Regression Model (LRM) response surface meta-models. It includes the number of 

active and unsuccessful design points, CoD, CoP, the highest average Cp obtained 

in the sensitivity analysis, and the average Cp increment. These data points are 

contrasted with the results presented in the study by (Raciti Castelli; Englaro; Benini, 

2011). 



57 

 

 
 

Figure 16 illustrates 3D heat maps of response surfaces for each TSR, 

showcasing the range of two input parameters (β, m) and the corresponding 

simulated design points. Here, the x-axis denotes the profile pitch angle (β), the y-

axis signifies the relative thickness of the profile (m), and the z-axis represents the 

average Cp. In a similar vein, Figure 17 presents the outcomes for the range of two 

other input parameters (β, a/b), with the x-axis indicating the profile’s pitch angle (β) 

and the y-axis depicting the position of maximum thickness of the profile (a/b). 

Figures 18 and 19 provide the same information as Figures 16 and 17, 

respectively, but from a top-down, 2D viewpoint. Results demonstrate that the zone 

of maximum average Cp, referred to as the zone of interest and shown in red, is 

attained in one particular peak region, while none of the middle ground between 

these regions produces a discernible increase in average Cp. 

 

Table 9. Summary of the best Meta-models for each specific TSR. 

TSR Number 
of 
design 
points 

Number 
of failed 
points 

Coefficient of 
Determination 
(CoD) 

Coefficient 
of 
Prognosis 
(CoP) 

Highest 
Average 
𝐶𝑝 

Average 

𝐶𝑝 

increment 

2.33 97 3 0.999939 0.999962 0.3055 11.098% 

2.64 93 7 0.999922 0.999082 0.3512 14.364% 

3.09 95 5 0.999289 0.999640 0.3385 12.960% 

Source: Prepared by the author. 

 

 



58 

 

 
 

Figure 16. 3D response surfaces rendered for each individual TSR; where the x-axis 
(𝛽) represents the pitch angle of the profile, y-axis (𝑚) represents the relative 

thickness of the profile and the z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. 

Source: Prepared by the author. 

 

Figure 17. 3D response surfaces rendered for each individual TSR; where the x-axis 

(𝛽) represents the pitch angle of the profile-axis (
𝑎

𝑏
) represents the maximum 

thickness position of the profile and z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. 

Source: Prepared by the author. 

 



59 

 

 
 

Figure 18. 2D response surfaces rendered for each individual TSR; where the x-axis 
(𝛽) represents the pitch angle of the profile, y-axis (𝑚) represents the relative 

thickness of the profile and the z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. 

Source: Prepared by the author. 

 

Figure 19. 2D response surfaces rendered for each individual TSR; where the x-axis 

(𝛽) represents the pitch angle of the profile-axis (
𝑎

𝑏
) represents the maximum 

thickness position of the profile and z-axis represents the average 𝐶𝑝 at the 

aforementioned TSR. 

Source: Prepared by the author. 

 

The Coefficient of Importance (CoI) is illustrated in Figures 20-22. Each row 

represents the contribution values for each input parameters. As shown in Figure 20, 

the pitch angle β contributes most significantly (76% for TSR = 2.33), followed by the 

relative thickness m (18.9%) and the position of maximum thickness (a/b) (5.1%). As 

TSR increases, the individual contributions of pitch angle β and position of maximum 

thickness (a/b) also grow, as depicted in Figures 21 and 22. Conversely, the 

individual contribution of relative thickness m diminishes with increasing TSR, as 



60 

 

 
 

indicated in Figures 20-22. This pattern aligns with prior research such as the study 

by (Francisco; Trentin, 2021), which found that the pitch angle β contributes the most 

(over 58%) to the increase in average Cp in Vertical Axis Wind Turbines (VAWTs). 

 

Figure 20. Parameter optimization process flowchart for the optimal profile at TSR = 
2.33. 

Source: Prepared by the author. 

 
Figure 21. Parameter optimization process flowchart for the optimal profile at TSR = 

2.64 

Source: Prepared by the author. 



61 

 

 
 

Figure 22. Parameter optimization process flowchart for the optimal profile at TSR = 
3.09. 

 

Source: Prepared by the author. 

4.2.  OPTIMIZATION PROCESS SPECIFICS AND RESULTS 

We employed gradient-based optimization techniques, including NLPQL and 

the Downhill Simplex Method, to develop the Metamodel of Optimal Prognosis 

(MOP), which boasted a coefficient of prognosis over 99%. These techniques help 

identify the optimal input parameters (a/b, m, β) to maximize the numerical 

integration of the average 𝐶𝑝
̅̅ ̅- TSR curve. To ensure model accuracy, CFD 

simulations were performed to compare predicted values with the optimization 

models, requiring a margin of error below 1%. Table 10 presents the optimal design 

points for each optimization model, along with corresponding simulation results and 

model comparisons. 

 



62 

 

 
 

Table 10. Detailed information about results and validation of the optimization 

models. 

Opt. 
model 

Optima

l 
𝑎

𝑏
 

Optimal 
𝑚 

Optimal 
𝛽 

Predicted 

∫ 𝐶𝑝
̅̅ ̅𝜆𝑓

𝜆0
 

Simulated 

∫ 𝐶𝑝
̅̅ ̅𝜆𝑓

𝜆0
 

MoE Validation 

NLPQL 1.054 0.0407 -3.491 0.2545 0.2540 0.2150
% 

Passed 

SIMPLEX 1.053 0.0408 -3.491 0.2532 0.2523 0.3560
% 

Passed 

Source: Prepared by the author. 

 

As Table 10 shows, both models display high accuracy between the predicted 

and simulated results. However, the NLPQL model yields the smallest margin of 

error. The optimal value for the position of maximum thickness along the chord (a/b) 

is 1.054, which corresponds to a maximum thickness position at 31.21% of the chord. 

The relative thickness (m) is set at 0.0407, which equates to a profile's relative 

thickness of 22.3%, and the optimal pitch angle is -3.491°. 

For this study, the NACA 0021 profile serves as the baseline, with a maximum 

thickness position at 30% and a maximum thickness of 21% of the chord (Table 1). 

Figure 23 provides a comparison between the NACA 0021 profile and the optimized 

profile. In the chapters that follow, a comprehensive comparison between the optimal 

profile and the original NACA 0021 profile will be presented. 

 



63 

 

 
 

Figure 23. Optimal profile and NACA 0021 profile. 

 

Source: Prepared by the author. 

 

4.3.  COMPARISON OF THE OPTIMAL PROFILE AT TSR = 2.33. 

 

Figures 24, 29, and 34 compare the torque coefficients of the base NACA 

0021 profile and the optimal profiles. At TSR = 2.33, the optimal profile exhibits a 

smaller range of negative torque coefficients compared to the NACA 0021 profile 

(Figure 24). Both profiles display regions of high and low torque coefficients, with the 

extreme points occurring at different azimuthal positions (Figure 24) 



64 

 

 
 

Figure 24. Instantaneous torque coefficient (𝐶𝑡) for the NACA 0021 and the 

optimized profile at TSR = 2.33. 

Source: Prepared by the author. 

 

The optimized configuration demonstrates superior aerodynamic performance, 

achieving a marginally higher minimum torque coefficient, as illustrated in Figures 25 

and 26. These figures depict the dimensionless velocity at azimuthal positions 

corresponding to the highest and lowest torque coefficients, respectively. The NACA 

0021 profile exhibits drag zones at θ=240° (Figure 25a) and θ=196° (Figure 25b), 

resulting in flow separation, a turbulent wake region, and a decrease in both lift and 

the turbine’s power coefficient. In contrast, the optimized configuration shows no 

such drag zones and possesses a slightly higher minimum torque coefficient (Figure 

26). 

 

 

 

 



65 

 

 
 

Figure 25. Dimensionless velocity contour for the NACA 0021 at TSR = 2.33 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. 

Source: Prepared by the author. 

 

Figure 26. Dimensionless velocity contour for the optimalprofile atTSR= 2.33 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡 

Source: Prepared by the author. 

 

Figures 27 and 28 display the dimensionless pressure contours at azimuthal 

points of the lowest and highest torque coefficients for the NACA 0021 profile and its 

optimum configuration at TSR = 2.33. Figure 27a reveals a low-pressure zone at the 

trailing edges of blades 2 and 3, a consequence of flow separation and drag zones 

that decrease the lift and turbine’s power coefficient for the NACA 0021 profile. The 

optimized configuration (Figure 27b) does not exhibit these negative characteristics 



66 

 

 
 

to the same extent, resulting in a higher torque coefficient. Figure 28 shows that the 

highest torque coefficient of both profiles is strongly influenced by the azimuthal 

positions of blades 1 and 3, where a significant pressure difference exists between 

the two sides. This difference generates a substantial lift force, enhancing the 

turbine’s performance. 

 

Figure 27. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR=2.33) at the azimuthal position of minimum torque 

coefficient 𝐶𝑡 

 

Source: Prepared by the author. 

 



67 

 

 
 

Figure 28. Dimensionless pressure contours for the three blades of the NACA 0021 
and the optimal profile (TSR =2.33) at the azimuthal position of maximum torque 

coefficient 𝐶𝑡 

 

Source: Prepared by the author. 

 

4.4.  COMPARISON OF THE OPTIMAL PROFILE AT TSR = 2.64 

 

Figure 26 compares the instantaneous torque coefficient of the NACA 0021 

profile with the optimal profile at TSR = 2.64. The optimal profile surpasses the NACA 

0021 profile, as it does not exhibit negative torque coefficients at three different 

regions of its rotation. 



68 

 

 
 

Figure 29. Instantaneous torque coefficient (𝐶𝑡) for the NACA 0021 and the optimized 

profile at TSR = 2.64.  

Source: Prepared by the author. 

 

Figures 30 and 31 further investigate the differences between the base and 

optimal configurations by showing the velocity at the azimuthal positions of the 

highest and lowest torque coefficients. Figure 30a and Figure 30b indicate that the 

NACA 0021 profile exhibits distinct drag zones at θ = 242° and θ = 206°, leading to 

flow separation, a turbulent wake region, and negative torque coefficients. 

Conversely, the flow of the optimal configuration, as shown in Figure 31, is smoother, 

resulting in an overall enhancement of aerodynamic performance. 

 

 

 



69 

 

 
 

Figure 30. Dimensionless velocity contour for the for the NACA 0021 at TSR= 2.64 at 

the azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. 

Source: Prepared by the author. 

 

Figure 31. Dimensionless velocity contour for the optimal profile at TSR= 2.64 at the 
azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. 

Source: Prepared by the author. 

 

Figure 32a reveals a low-pressure zone at the trailing edges of blade 2, 

caused by the flow separation and drag zones displayed in Figure 30. This 

phenomenon decreases lift and subsequently reduces the torque coefficient of the 

turbine. In contrast, the optimal profile depicted in Figure 32b lacks these unfavorable 



70 

 

 
 

features, resulting in an increase in the torque coefficient. As shown in Figure 33, the 

azimuthal positions of blades 1 and 3 significantly influence the maximum torque 

coefficient of both the optimal and NACA 0021 profiles. The pressure difference 

between the two sides of blades 1 and 3 generates a substantial lift force, improving 

the turbine's performance. 

 

Figure 32. Dimensionless pressure contours for the three blades of the NACA 0021 

and the optimal profile (TSR = 2.64) at the azimuthal position of minimum torque 

coefficient 𝐶𝑡. 

 

Source: Prepared by the author. 

 



71 

 

 
 

Figure 33. Dimensionless pressure contours for the three blades of the NACA 0021 
and the optimal profile (TSR=2.64) at the azimuthal position of maximum torque 

coefficient 𝐶𝑡 

 

Source: Prepared by the author. 

 

4.5.  COMPARISON OF THE OPTIMAL PROFILE AT TSR = 3.09. 

 

Figure 34 presents the instantaneous torque coefficient of the optimal and 

NACA 0021 profiles at TSR = 3.09. Remarkably, neither profile displays negative 

torque coefficients. 

 



72 

 

 
 

Figure 34. Instantaneous torque coefficient (𝐶𝑡) for the NACA 0021 and the optimized 
profile at TSR = 3.09. 

 

Source: Prepared by the author. 

 

Figures 35 and 36 further delve into the differences between the base and 

optimal configurations at TSR = 3.09, displaying the dimensionless velocity at the 

azimuthal positions of the highest and lowest torque coefficients. Unlike Figures 25 

and 30, Figure 35 reveals that the NACA 0021 profile exhibits a smoother flow, 

resulting in an absence of negative torque coefficients. The same is true for the ideal 

configuration at TSR = 3.09. 

 

 

 

 

 

 



73 

 

 
 

Figure 35. Dimensionless velocity contour for the NACA 0021 at TSR = 3.09 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. 

 

Source: Prepared by the author. 

 

Figure 36. Dimensionless velocity contour for the optimal profile at TSR = 3.09 at the 

azimuthal position of (a) lowest and (b) highest torque coefficient 𝐶𝑡. 

 

Source: Prepared by the author. 

 

Figures 37 and 38 illustrate the dimensionless pressure contours at the 

azimuthal positions of the highest and lowest torque coefficients for the NACA 0021 



74 

 

 
 

profile and its optimal configuration at TSR = 3.09. At this TSR, the NACA 0021 

profile shows no low-pressure zone at the trailing edge. Similarly, the ideal profile has 

a consistent pressure distribution without low-pressure zones (Figure 37). Moreover, 

Figure 38 demonstrates that the pressure difference between the two sides of blades 

1, 2, and 3 significantly impacts the maximum torque coefficient for both the optimal 

configuration and the NACA 0021 profile. This difference generates a substantial lift 

force that enhances the turbine's performance. 

 

Figure 37. Dimensionless pressure contours for the three blades of the NACA 0021 
and the optimal profile (TSR=3.09) at the azimuthal position of minimum torque 

coefficient 𝐶𝑡. 

 

Source: Prepared by the author. 

 



75 

 

 
 

Figure 38. Dimensionless pressure contours for the three blades of the NACA 0021 
and the optimal profile (TSR =3.09) at the azimuthal position of maximum torque 

coefficient 𝐶𝑡. 

 

Source: Prepared by the author. 

 

Lastly, Figure 39 compares the NACA 0021 and the optimized profile. As 

shown in Figure 39, the optimized profile exhibits a higher average Cp in the range 

for which the profile was optimized (TSR =2.33-3.09). Moreover, for lower and higher 

TSR, the CFD simulation also indicates an increase in the average 𝐶𝑝. 

 



76 

 

 
 

Figure 39. Comparison between the NACA 0021 and the optimized profile  

 

Source: Prepared by the author. 

 

 

 



77 

 

 
 

5. CONCLUSIONS 

 

In this study, a low TSR Darrieus turbine was optimized using sensitivity 

analysis of three input parameters: a/b (position of maximum thickness), m (relative 

thickness), and pitch angle (𝛽). Gradient-based optimization methods were 

employed, and a turbulent, unsteady, 2-D incompressible flow regime was 

considered. The study focused find an optimal profile for an operational range in 

terms of performance through the numerical integration of the 𝐶𝑝
̅̅ ̅ − 𝑇𝑆𝑅 curve. 

Specific findings of the profile optimization are listed below: 

 The second-order linear regression was identified as the best-fit MOP 

for all TSR based on CoD and CoP (Table 9). 

 In the sensitivity analysis, Figure 20, 21, and 22 show the contribution 

of each of the input parameters to the increment in average Cp. Results 

indicate that pitch angle 𝛽 has the highest contribution (over 76%), 

followed by relative thickness m and maximum thickness position a/b. 

This finding is consistent with previous studies, such as (Francisco; 

Trentin, 2021), which also emphasized the significance of pitch angle in 

VAWT performance. 

 The optimal profile has a corresponding to a maximum thickness 

position at 31.21% of the chord, a relative thickness of 22.3% and a 

pitch angle of -3.491°. It should be noted that the original NACA 0021 

profile has a maximum thickness position of 30% of the chord and a 

maximum thickness position of 21% (Table 1). 

 The comparison of the optimal profile with the base NACA 0021 profile 

across different TSR values (2.33, 2.64, and 3.09) clearly demonstrates 

the superior performance of the optimized profile. This superiority is 

evident in the reduced negative torque coefficients, smoother flow 

patterns, and fewer low-pressure zones. The optimized profile also 



78 

 

 
 

shows a higher average power coefficient (Cp) across the optimized 

range, as well as at lower and higher TSR values. Therefore, the study 

concludes that profile optimization significantly improves the 

aerodynamic performance of the turbine, highlighting its importance for 

efficient turbine operation. 

 

In the recommendations section, information is provided about the potential 

directions for further study. 



79 

 

 
 

6.  RECOMMENDATIONS 

 

 It is crucial to extend the analysis over a full 3D geometry in order to 

include a more difficult but thorough investigation and avoid any 

overestimation due to geometrical simplification (2-D simulation of the 

VAWT) and the absence of blade tip loss considerations. 

 The goal of the current study was to increase the Darrieus turbine’s 

average at a specific wind speed of 9 m/s. Even though this wind speed 

gives us a clear grasp of the turbine’s aerodynamic behavior, it’s still 

vital to set a variety of inlet velocity speeds to forecast how the turbine 

will behave in a range of wind speeds related to dispersed generation. 

 Finally, it is crucial to expand the input variable (Inlet velocity, TSR) and 

extend research beyond the optimization of the average in order to fully 

comprehend the fluid-dynamic aspects behind the Darrieus turbine. 

Examples of this include incorporating higher TSRs (TSR> 4) and other 

related components for a more integrated wind energy system (power 

transmission and electric generator). 

 



80 

 

 
 

REFERENCE 

AHMADI-BALOUTAKI, M.; MOJTABA. Analysis and Improvement of Aerodynamic 
Performance of Straight Bladed Vertical Axis Wind Turbines. PhDT, p. 196, 2016.  
 
BALDUZZI, F. et al. Critical issues in the CFD simulation of Darrieus wind turbines. 
Renewable Energy, Oxford, v. 85, p. 419-435, 1 jan. 2016.  

 
BAZMI, A. A.; ZAHEDI, G. Sustainable energy systems: role of optimization modeling 
techniques in power generation and supply: a review. Renewable and Sustainable 
Energy Reviews, Oxford, v. 15, n. 8, p. 3480-3500, 1 out. 2011. 

  
BELABES, B.; PARASCHIVOIU, M. Numerical study of the effect of turbulence 
intensity on VAWT performance. Energy, Oxford, v. 233, p. 121139, 15 out. 2021.  
 
BIANCHINI, A. et al. Effectiveness of two-dimensional CFD simulations for Darrieus 
VAWTs: a combined numerical and experimental assessment. Energy Conversion 
and Management, Oxford, v. 136, p. 318-328, 15 mar. 2017.  
 
BROWNSTEIN, I. D.; WEI, N. J.; DABIRI, J. O. Aerodynamically interacting vertical-
axis wind turbines: performance enhancement and three-dimensional flow. Energies, 

Basel, v. 12, n. 14, p. 2724, 16 jul. 2019.  
 
BUCHER, C. Basic concepts for robustness evaluation using stochastic analysis. In: 
EFFICIENT METHODS FOR ROBUST DESIGN AND OPTIMIZATION, 482., 2007, 
London. Proceedings [...]. London: EUROMECH Colloquium, 2007. 
 
BUCHNER, A.-J. et al. Dynamic stall in vertical axis wind turbines: scaling and 
topological considerations. Journal of Fluid Mechanics, Cambridge, v. 841, p. 746-

766, 25 abr. 2018.  
 
CARRIGAN, T. J. et al. Aerodynamic Shape Optimization of a Vertical-Axis Wind 
Turbine Using Differential Evolution. ISRN Renewable Energy, New York, v. 2012, 

p. 1-16, 16 jan. 2012.  
 
CELIK, I. B. et al. Procedure for estimation and reporting of uncertainty due to 
discretization in CFD applications. Journal of Fluids Engineering, New York, v. 

130, n. 7, p. 0780011-0780014, 1 jul. 2008.  
 
CHEN, W.; CHIU, K.; FUGE, M. D. Airfoil Design Parameterization and Optimization 
using B\’ezier Generative Adversarial Networks. AIAA Journal, Reston, v. 58, n. 11, 

p. 4723-4735, 21 jun. 2020.  
 
CRESSIE, N. The origins of kriging. Mathematical Geology, New York, v. 22, n. 3, 
p. 239-252, abr. 1990.  
 
DABIRI, J. O. Potential order-of-magnitude enhancement of wind farm power density 
via counter-rotating vertical-axis wind turbine arrays. Journal of Renewable and 
Sustainable Energy, Maryland, v. 3, n. 4, p. 043104, 19 jul. 2011.  



81 

 

 
 

 
FERREIRA, C. S. et al. Comparison of aerodynamic models for Vertical Axis Wind 
Turbines. Journal of Physics: Conference Series, Bristol, v. 524, n. 1, p. 012125, 16 

jun. 2014.  
 
FERREIRA, C. S. et al. Visualization by PIV of dynamic stall on a vertical axis wind 
turbine. Experiments in Fluids, Secaucus, v. 46, n. 1, p. 97-108, 3 ago. 2008. 

  
FIEDLER, A. J.; TULLIS, S. Blade offset and pitch effects on a high solidity vertical 
axis wind turbine. Wind Engineering, London, v. 33, n. 3, p. 237-246, 1 maio 2009. 
DOI: http://dx.doi.org/10.1260/030952409789140955. 
 
FORTRAN, I. et al. Numerical recipes: the art of scientific computing. 2. ed. [S. l.]: 

Fortran Numerical Recipes, 1986.  
 
FRANCISCO, P. et al. Screening analysis and unconstrained optimization of a small-
scale vertical-axis wind turbine. Energy, Amsterdam, v. 240, p. 122782, 2022.  

 
FREENY, A.; BOX, G. E. P.; DRAPER, N. R. Empirical model building and response 
surfaces. Technometrics, Washington, v. 30, n. 2, p. 229, maio 1988.  
 
FUJII K., DULIKRAVICH G. S. (ed.). Parametric airfoils and wings. Notes on 
Numerical Fluid Mechanics, Heidelberg, v. 68, pp.71-88, 1998. 

  
GHASEMIAN, M.; ASHRAFI, Z. N.; SEDAGHAT, A. A review on computational fluid 
dynamic simulation techniques for Darrieus vertical axis wind turbines. Energy 
Conversion and Management, London, v. 149, p. 87-100, 1 out. 2017.  

 
HANSEN, J. T.; MAHAK, M.; TZANAKIS, I. Numerical modelling and optimization of 
vertical axis wind turbine pairs: a scale up approach. Renewable Energy, Oxford, v. 
171, p. 1371–1381, 1 jun. 2021.  
 
HASHEM, I.; MOHAMED, M. H. Aerodynamic performance enhancements of H-rotor 
Darrieus wind turbine. Energy, Amsterdam, v. 142, p. 531-545, 1 jan. 2018.  
 
HOMEM-DE-MELLO, T.; BAYRAKSAN, G. Monte Carlo sampling-based methods for 
stochastic optimization. Surveys in Operations Research and Management 
Science, Amsterdam, v. 19, n. 1, p. 56-85, 1 jan. 2014.  
 
KELLEY, C. T. Iterative methods for optimization. Philadelphia: SIAM, 1999.  
 
KLIMAS, P. C.; WORSTELL, M. H. Effects of blade preset pitch: offset on curved-
blade darrieus vertical axis wind turbine performance [SF2900a (8-81 ).] 
Albuquerque: Sandia National Laboratories, 1981.  
 
LAM, H. F.; PENG, H. Y. Study of wake characteristics of a vertical axis wind turbine 
by two- and three-dimensional computational fluid dynamics simulations. Renewable 
Energy, Amsterdam, v. 90, p. 386-398, 1 maio 2016.  
 



82 

 

 
 

LANCASTER, P.; SALKAUSKAS, K. Surfaces generated by moving least squares 
methods. Mathematics of computation, Providence, v. 37, n. 155, p. 141, jul. 1981.  
 
LAUNDER, B. E. Turbulence Modelling for CFD. By D. C. WILCOX. DC