WIND ENERGY
Wind Energ.
2017; 20:909–926
Published online 5 December 2016 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/we.2070
RESEARCH ARTICLE
Aerodynamic shape optimization of wind turbine
blades using a Reynolds-averaged Navier–Stokes
model and an adjoint method
Tristan Dhert
1
, Turaj Ashuri
2,†
and Joaquim R. R. A. Martins
1
1
University of Michigan, Ann Arbor, MI 48109 , USA
2
University of Texas at Dallas, Richardson, 75080 Texas, USA
ABSTRACT
Computational fluid dynamics (CFD) is increasingly used to analyze wind turbines, and the next logical step is to develop
CFD-based optimization to enable further gains in performance and reduce model uncertainties. We present an aerodynamic
shape optimization framework consisting of a Reynolds-averaged Navier Stokes solver coupled to a numerical optimiza-
tion algorithm, a geometry modeler, and a mesh perturbation algorithm. To efficiently handle the large number of design
variables, we use a gradient-based optimization technique together with an adjoint method for computing the gradients of
the torque coefficient with respect to the design variables. To demonstrate the effectiveness of the proposed approach, we
maximize the torque of the NREL VI wind turbine blade with respect to pitch, twist, and airfoil shape design variables
while constraining the blade thickness. We present a series of optimization cases with increasing number of variables, both
for a single wind speed and for multiple wind speeds. For the optimization at a single wind speed performed with respect to
all the design variables (1 pitch, 11 twist, and 240 airfoil shape variables), the torque coefficient increased by 22.4% rela-
tive to the NREL VI design. For the multiple-speed optimization, the torque increased by an average of 22.1%. Depending
on the CFD mesh size and number of design variables, the optimization time ranges from 2 to 24h when using 256 cores,
which means that wind turbine designers can use this process routinely. Copyright © 2016 John Wiley & Sons, Ltd.
KEYWORDS
wind turbine design; aerodynamic shape optimization; computational fluid dynamics; NREL VI; numerical optimization; adjoint
methods
Correspondence
Joaquim R. R. A . Martins, University of Michigan, Ann Arbor, MI 48109, USA
E-mail: jrram@umich.edu
†
Present address: Arkansas Tech University, Russellville, AR 72081, USA
Received 22 November 2015; Revised 4 July 2016; Accepted 3 October 2016
1. INTRODUCTION
Wind energy has experienced a considerable reduction in the levelized cost of energy over the last few years. However,
further cost reduction is needed for wind energy to become more competitive with traditional energy resources such
as coal and natural gas. Besides decreasing the turbine capital cost (TCC) and balance of station, the annual energy
production (AEP) has to increase to further reduce costs. This reduction can be achieved by upscaling wind turbines,
developing new optimization methods to further optimize wind turbine designs, and introducing new wind turbine concepts
and components.
The complexity of this challenge increases with the size of wind turbines,
1
and placing wind turbines offshore
further exacerbates this problem.
2
Given the multidisciplinary nature of wind turbines, their design hinges on opti-
mizing the trade-off between aerodynamic performance, structural efficiency, and operational and manufacturing cost.
3
Therefore, multidisciplinary design optimization
4
could become an important tool in the design of the next generation of
high-performance wind turbines by maximizing the levelized cost of energy while accounting for all relevant disciplines
and interactions.
5
Copyright © 2016 John Wiley & Sons, Ltd.
909
Aerodynamic shape optimization of wind turbine blades T. Dhert, T. Ashuri and J. R. R. A. Martins
Previous studies have tackled several areas in wind turbine design, ranging from the optimization of rotor design
6–10
to
support structure design.
11–13
These studies used various levels of fidelity in the models, ranging from low fidelity (blade
element momentum-based fidelity)
14–20
to higher fidelity (free and prescribed vortex-based fidelity).
21–27
Currently, most wind turbine blade designs are based on low fidelity models because of their ease of implementation,
lower computational cost, and fast convergence to feasible solutions. However, these models do not accurately represent
compressibility, viscosity and the three-dimensional (3D) effects present in wind turbines, which result in suboptimal
designs. These important effects can only be captured accurately by using 3D computational fluid dynamics (CFD) that
include viscous terms.
One approach adopted to incorporate high-fidelity simulations into wind turbine design is to use two-dimensional
Reynolds-averaged Navier–Stokes (RANS) CFD solvers to optimize wind turbine airfoils, neglecting 3D effects.
28, 29
Another effort considered 3D CFD to optimize turbine blade winglets with respect to two design variables.
30
However,
CFD-based aerodynamic shape optimization for wind turbine blades has not yet been done because of the compounding
challenges of high cost simulations and large numbers of shape design variables. In recent work, Economon et al.
31
used
a RANS solver and a continuous adjoint approach to optimize the NREL VI wind turbine blade with respect to 50 airfoil
shape variables. Their optimization resulted in a 4% increase in torque, but only three design optimization iterations were
completed. Large-eddy simulation (LES) provides even higher fidelity than RANS for massively separated flow conditions,
but the computational cost is prohibitive for inclusion in an optimization cycle.
These challenges have been addressed in the aerospace engineering community, where CFD-based design optimiza-
tion has been used for decades, with applications to wing design
32–35
and helicopter blade design.
36–40
Given that wind
turbine blades are in essence rotating wings, and share some of the physics in helicopter blades, much of the design
optimization techniques developed by the aerospace community can be adapted to the design optimization of wind
turbine blades.
In the present work, we address the challenges of CFD-based wind turbine optimization by developing a robust and
efficient aerodynamic shape optimization framework. We use an arbitrary Lagrangian–Eulerian RANS model for the CFD
analysis and a discrete adjoint method that efficiently computes the gradients of the performance metrics with respect to the
design variables. In contrast to lower fidelity models, the RANS CFD used herein captures the complex flow phenomena
present in wind turbines, resulting in a more reliable performance analysis and ultimately a more refined design.
Unsteady effects such as tower shadow and dynamic stall are neglected in this work, since performing design optimiza-
tion based on unsteady CFD is currently too costly. While we perform some verification and validation of our CFD model
against the NREL VI experimental results, the focus of this work is on demonstrating the first CFD-based wind turbine
optimization that considers a comprehensive set of airfoil shape design variables. The validation of CFD models, including
the validity of turbulence models, is still the subject of ongoing research in the wind energy community.
To demonstrate the proposed approach, we present a series of torque coefficient maximization results starting with the
NREL VI as the baseline design. These optimizations range from a single wind speed case with pitch as a design variable
to a case that optimizes 251 design variables for multiple wind speeds within the rated region.
The remainder of this paper is organized as follows: First, we present the aerodynamic shape optimization methodology
for rotating flow problems. Next, we apply the methodology to the NREL VI wind turbine blade as a benchmark case.
41
This requires a problem statement fitted for this particular case. Finally, we present the results of the NREL VI wind turbine
blade optimization, followed by the conclusions.
2. METHODOLOGY FOR AERODYNAMIC SHAPE OPTIMIZATION
We now describe the various aerodynamic shape optimization framework components used in this work. In addition to the
CFD solver, the following components are required: geometry parametrization, mesh deformation, gradient computation,
and numerical optimization algorithm.
2.1. Geometric parametrization
To optimize a given shape, we need to parametrize it by using a set of shape parameters that the optimizer ultimately
controls. This parametrization must result in smooth shape changes to ensure that our gradient-based optimizer is effective.
In addition, the geometry parametrization should be differentiable, and the gradients of the surface coordinates with respect
to the shape parameters should be computed accurately and efficiently.
In this work, we use the free-form deformation (FFD) volume approach to parametrize the blade geometry by embedding
the design inside a hexahedron.
42
This approach is known to have an efficient and compact set of geometry design variables,
which facilitates the manipulation of complex geometries. Once the local coordinates of the vertices of the geometry are
expressed in the FFD volume by doing a Newton search, FFD control points on the surface of the volume are used to
deform the baseline design. When the lattice structure within the volume deforms by moving the surface control points, the
Wind Energ.
2017; 20:909–926 © 2016 John Wiley & Sons, Ltd.
910
DOI: 10.1002/we
T. Dhert, T. Ashuri and J. R. R. A. Martins Aerodynamic shape optimization of wind turbine blades
Figure 1. Parametrization of FFD volume shape for the NREL VI rotor. [Colour figure can be viewed at wileyonlinelibrary.com]
FFD method deforms the embedded object as well. We use an FFD volume composed of B-spline tensor patches in three
dimensions; this is called a trivariate B-spline volume. B-spline volumes are chosen because derivatives of any point inside
the volume are easily computed.
43
Two FFD volumes are constructed for the two-bladed NREL VI rotor: one for each blade. The control-point parameters
of the two FFD volumes are linked to ensure that the two blades have exactly the same shape. Figure 1 shows the FFD
volumes for the original NREL VI rotor. The control points are fixed in the cylindrical piece near the root because this
section is not shaped for aerodynamic performance.
Although the basic shape variables correspond to the movement of the control points at the surface of the FFD volume,
we also implement composite shape variables corresponding global shape changes, such as blade pitch and the spanwise
variation of twist, which rotate the blade sections about its 25% chord line. The twist can change in the 11 outboard
sections (the root section is set to match the blade pitch). The spanwise interpolation of the twist distribution is done by
using a cubic spline within the FFD. The airfoil is shaped at the 12 outboard sections (i.e., the whole blade excluding
the cylindrical piece next to the hub). Each airfoil section is shaped by moving 20 FFD control points (10 on the lower
surface and 10 on the upper surface), resulting in a total of 240 airfoil shape variables. Our previous work on wing design
shows that this number of shape design variables provides a good trade-off between computational effort and quality of
the optimum.
35
The geometry parametrization can also handle chord and span variation, but we do not exploit them in this work because
this would require introducing other disciplines, such as structures, to make the appropriate optimization trade-offs.
44
2.2. Mesh deformation
Once a new blade surface shape is determined from the set of shape variables obtained by using FFD, surface deformations
are applied to the CFD surface mesh, and then propagated to the entire CFD volume mesh.
To propagate the surface-shape deformations into the volume mesh, we use the hybrid approach developed by Ken-
way et al.,
43
which combines algebraic and linear elasticity methods. The main idea of the hybrid approach is to
apply a linear elasticity method on a coarser approximation of the mesh for large, low-frequency perturbations, while
using an algebraic approach to attenuate small, high-frequency deformations. The hybrid approach consists of the
following steps:
1. Select a subset of nodes for each mesh edge. These nodes form a coarser approximation of the complete mesh.
2. Deform this coarser approximation of the initial mesh by using the linear elasticity method.
3. Regenerate each block with linear or cubic-Hermite spine interpolation to obtain the full mesh approximation.
4. Find the remainder of the deformed mesh by using the algebraic technique.
As an example, we compare the original and deformed mesh of the NREL VI S809 airfoil in Figure 2.
Wind Energ.
2017; 20:909–926 © 2016 John Wiley & Sons, Ltd.
DOI: 10.1002/we
911
Aerodynamic shape optimization of wind turbine blades T. Dhert, T. Ashuri and J. R. R. A. Martins
Figure 2. Mesh deformation of the NREL VI blade at
z
D 3.3 m.
2.3. Computational fluid dynamics solver
In this research, we use the RANS-based CFD solver ADflow, which is an extension of the SUmb solver
45
that adds
adjoint and optimization capability.
46
ADflow uses parallel computing to solve the Euler, laminar Navier–Stokes and
RANS equations in either steady, unsteady or time spectral modes.
47
Depending on the size of the computational domain,
a cluster of 16 to 256 processors is used for the cases presented in this paper. ADflow is a finite-volume, cell-centered
multi-block CFD solver that provides options for a variety of turbulence models with one, two and four equations, and
a variety of adaptive wall functions. The Roe scheme
48
with no limiter and a four-stage Runge–Kutta scheme
49
with a
segregated Spalart–Allmaras turbulence model
50
are used for the design optimizations presented herein. It is a common
practice to use the k turbulence model for wind energy applications; however, because we are simulating the rotor in
the attached flow conditions, the Spalart–Allmaras and k turbulence models yield similar results.
When considering an isolated rotor rotating at a constant angular velocity ! D
!
x
, !
y
, !
z
T
, where the rotor plane
is perpendicular to a uniform steady wind in a fixed reference frame, the flow field is intrinsically unsteady. However, to
simplify the CFD analysis and make design optimization tractable, we approximate the flow as a steady-state solution by
solving the problem in a reference frame co-rotating with the rotor. Unsteady flow can still exist in the presence of unsteady
phenomena, such as flow separation. By using this approach, the computational domain does not rotate, but the flow field
incorporates a rotational term. Assuming that both reference frames have the same origin, the velocity u
!
of the co-rotating
reference frame can be related to the velocity u in a fixed reference frame as follows:
u D u
!
C ! x,(1)
where x D Œx x
0
, y y
0
, z z
0
T
is the position vector pointing from the origin of the reference frame .x
0
, y
0
, z
0
/ to a
point .x, y, z/ in the flow domain. The time derivative of a time-dependent vector function f .t/ in the co-rotating reference
frame is
df .t/
dt
D
df .t/
dt
!
C ! f .t/.(2)
The conservation of momentum of an arbitrary time-independent control volume with permeable boundary @ is
governed by the Navier–Stokes equations. By using equations (1) and (2) in the momentum equations for a fixed reference
frame, conservation of momentum can be expressed in the integral conservation form with co-rotating motion (assuming
no external forces):
@
@t
Z
ud C
Z
@
u
Œu .! x/
T
n
dS D
Z
@
.pn
n/ dS
Z
.! u/d,(3)
where is the local density, p is the static pressure,
is the viscous stress tensor (assuming a Newtonian fluid), u D
Œu, v, w
T
is the local velocity of the fluid, and n D
n
x
, n
y
, n
z
T
is the outward pointing unit normal vector. When using the
co-rotating reference frame, the momentum equations gain an extra body force that represents the combined effect of the
centrifugal and Coriolis forces. Since the forces are perpendicular to the direction of motion, no work is produced and thus
no additional term appears in the energy equation.
Wind Energ.
2017; 20:909–926 © 2016 John Wiley & Sons, Ltd.
912
DOI: 10.1002/we
T. Dhert, T. Ashuri and J. R. R. A. Martins Aerodynamic shape optimization of wind turbine blades
2.4. Adjoint formulation
As previously mentioned, RANS-based aerodynamic shape optimization with respect to a large number of design variables
requires a gradient-based optimization algorithm and an efficient approach for computing these gradients.
One of the most popular methods for computing derivatives is the finite-difference method. However, the computational
cost of this method does not scale well with the number of design variables, and it is subject to large numerical errors
because of subtractive cancellation. The complex-step method is more accurate but does not scale any better.
51
When using
the adjoint method, however, the cost of computing derivatives in a gradient is practically independent of the number of
design variables. For a review on derivative computation methods, see Martins and Hwang,
52
or for a more specialized
review focused on CFD, see Peter and Dwight.
53
A function of interest (I) for an aerodynamic shape optimization (which could be an objective or constraint function) is
dependent on the design variables describing the wind turbine blade geometry (x), and the flow state variables (w). Such a
function depends both directly on the design variables and indirectly through the solution of the governing equations. The
direct dependence assumes that the state variables remain constant, while the indirect dependence is due to the change in
the state variables required to satisfy the governing equations for a given set of design variables. Hence, the total derivative
of the function of interest can be written as
dI
dx
D
@I
@x
C
@I
@w
dw
dx
.(4)
The total derivative of the flow solution with respect to the design variables (dw=dx) can be found by observing that the
total derivative of the CFD residuals R.x, w.x// D 0 with respect to x remains zero for a feasible solution, and thus,
dR
dx
D
@R
@x
C
@R
@w
dw
dx
D 0.(5)
We can use the preceding equation to write a linear system whose solution is the total derivative of the flow solution with
respect to the design variables, i.e.,
@R
@w
dw
dx
D
@R
@x
.(6)
Substituting the solution of this linear system into equation (4), we obtain
dI
dx
D
@I
@x
@I
@w
@R
@w
1
@R
@x
.(7)
Noting that the matrix @R=@w can also be factorized with @I=@w, we can first solve the system of linear equations,
@R
@w
T
D
@I
@w
T
.(8)
where are called the adjoint variables. In this equation, the design variables x do not appear explicitly, and the adjoint
vector does not depend on the design variables. This constitutes the major advantage of the adjoint method: the cost of
computing the gradient dI=dx is independent of the size of the vector of design variables x. After solving for the adjoint
variables, we substitute them into the total derivative (4) that we ultimately want to compute, yielding
dI
dx
D
@I
@x
C
T
@R
@x
.(9)
The equations derived above represent the discrete adjoint approach for computing derivatives, where we are linearizing
the discrete system of equations R D 0, and then solving the linear equation. This approach contrasts with the continuous
approach, where the partial differential equations are linearized first and then discretized.
54
In the adjoint equation (8) and total derivative equation (9), the partial derivatives must be computed. We derive the
code to compute these partial derivatives by using the algorithmic differentiation tool Tapenade,
55
which automatically
parses source code and creates additional code that computes derivatives. Algorithmic differentiation relies on a systematic
application of the chain rule to each line in the source code. In this paper, we use the ADflow adjoint implementation of
Lyu et al.
46
Wind Energ.
2017; 20:909–926 © 2016 John Wiley & Sons, Ltd.
DOI: 10.1002/we
913