- Articles & preprints
- Peer review
- Editorial board
- Manuscript tracking
Correspondence: Kenneth Loenbaek (firstname.lastname@example.org)
We present an analytical model for assessing the aerodynamic performance of a wind turbine rotor through a different parametrization of the classical blade element momentum (BEM) model. The model is named the Radially Independent Actuator Disc (RIAD) model, and it establishes an analytical relationship between the local thrust loading and the local power, known as the local-thrust coefficient and the local-power coefficient respectively. The model has a direct physical interpretation, showing the contribution for each of the three losses: wake rotation loss, tip loss and viscous loss. The gradient for RIAD is found through the use of the complex step method, and power optimization is used to show how easily the method can be used for rotor optimization. The main benefit of RIAD is the ease with which it can be applied for rotor optimization and especially load constraint power optimization as described in Loenbaek et al. (2021). The relationship between the RIAD input and the rotor chord and twist is established, and it is validated against a BEM solver.
Wind turbine rotors are with their increasing size subject to continuous optimization, with the overall objective of reducing the cost of energy. Such optimizations are very complex because both the aerodynamic and the structural performance need to be included in the optimization setup. Combining both aerodynamics and structural performance has shown very promising trends indicating that a further cost reduction is possible; see, e.g., Perez-Moreno et al. (2016), Zahle et al. (2015) and Bottasso et al. (2012). However, these optimization studies build on existing aerodynamic and aeroelastic tools, include numerous design variables and constraints, and can be very complex. Thus, it is challenging to make a very general optimization study to map the design space. That is the reason for investigating alternative methods and models to ease the exploration of optimal rotor designs.
The development of aerodynamic models for wind turbines is closely linked to that of propellers and helicopters. The first theoretical model for predicting the aerodynamic performance of a rotor was the so-called 1D momentum theory developed by Betz (Betz, 1926) and Joukowsky, which resulted in the famous maximum power extraction limit of 59.3 % known as the Betz–Joukowsky limit or often just as the Betz limit (Okulov and van Kuik, 2012). The model assumed constant loading along the rotor radius in the flow direction.
Later Glauert developed the blade element momentum (BEM) theory (Durand and Glauert, 1935; Sørensen, 2016), which is an extension where momentum theory is used for radially independent stream tubes. It also included correction models for tip loss and highly loaded rotors. The model has since then been extended with multiple correction models to account for yaw misalignment, shear profiles, turbulent inflow, etc.
In this paper, we present an aerodynamic rotor performance model which we refer to as the Radially Independent Actuator Disc (RIAD) model. It establishes a direct analytical relationship between the local thrust loading and the local power, which is a useful simplification for rotor optimization. The model is equivalent to BEM but reduces the rotor design space to only two independent variables at each radial station, i.e., the local-thrust coefficient (CLT) and the glide ratio () as well as the global tip speed ratio (λ). The equivalent input for a BEM at each radial station is the design lift coefficient (Cl), the design drag coefficient (Cd) and the rotor solidity (σ) with the same global parameter. Most BEM formulations do not compute the local power directly, which is often an important optimization objective. At the same time, rotor optimization constraints are often formulated in terms of loading. Both the objective and constraints are outputs from the BEM, and their equations are usually fairly convoluted. Using RIAD, the local thrust loading (CLT) is the independent design parameter, and the local power is computed explicitly through a single equation. It makes it easy to recast the optimization problem, which generally requires a robust optimization algorithm, into a straightforward root-finding problem, which makes the optimization faster and more robust.
The paper starts by presenting the derivation of the RIAD model, which then leads into computing the gradients for RIAD. These gradients are used for power optimization, leading to a simple optimization method. In the end, the relationship between RIAD inputs and blade chord and twist is then established, and RIAD is validated against a BEM solver. This is Part 1 of a two-part paper. Part 1 describes an aerodynamic model for a wind turbine rotor and the use of the model for power optimization. Part 2 is described in Loenbaek et al. (2021), where the model is applied for load-constrained power optimization.
In the following the Radially Independent Actuator Disc (RIAD) model is presented, which starts by establishing the relationship between global and local parameters for a wind turbine rotor as well as introducing normalization. The relationship between the local forces is then established, leading to an implicit equation for the local power. A set of approximate closure equations is then used to establish an explicit equation. The physical interpretation of the different factors and terms is then presented, and at the end some details regarding the tip loss factor and the exclusion of drag from the induced velocity are discussed.
Starting with the fundamental theorem of calculus the following equation can be created for the global values for thrust (T) and power (P):
where is the thrust loading density and is the power density.
The classical non-dimensional relations for thrust (CT) and power (CP) as well as the local equivalent which is introduced as the local-thrust coefficient (CLT) and local-power coefficient (CLP) are introduced as follows:
Combining the equations, the following equations can be found for CT and CP:
where is the normalized radius (). A diagram of the relationship between the local and global coefficients can be seen in Fig. 1.
To establish a relationship between the local coefficients CLT and CLP, the forces are assumed to be aerodynamic forces where the lift force is orthogonal to the local flow and where there is loss from viscous drag in the local flow direction. The force is assumed to originate from rotating wind turbine blades with rotational speed ω. The tangential force density is then given as . Introducing the local lift density and drag density as well as the local flow angle ϕ, the following equations can be created:
where in Eq. (9) to exclude drag from induction. Using the common normalization of and ,
where σ is the rotor solidity () and the normalized relative velocity. A diagram of the airflow and forces can be seen in Fig. 2.
This is a general equation for the relationship between the CLT and CLP, but it is implicit since tan ϕ depends on CLT.
where a is the axial induction factor and ap is the tangential induction factor.
Equation (16) introduces two new variables (a, ap). In order to make an explicit equation for tan ϕ as a function of CLT two equations relating CLT to a and ap respectively are needed. These are the closure equations.
There does not exist a general set of model closure equations, but different approximate closures have been proposed. The most widely used set is referred to as the Glauert closure, which is an implicit assumption made for most BEMs. The closures are given as
where F is the tip loss factor, which is further described in Sect. 2.5. For now with F→0 as . Combining Eqs. (17) and (18) with Eq. (16) an explicit equation for tan ϕ in terms of CLT can be found. This leads to an explicit equation for CLP:
Equation (19) is the main result of this section. It states an explicit relationship between the local power and the local loading.
Equation (19) has a straightforward physical interpretation which is also highlighted through under-bracing of factors and terms, with an additional loss coming from the tip loss factor (F). A diagram showing the impact of different losses for a specific input can be seen in Fig. 3.
The wake rotation loss is the power loss that originates from the conservation of angular momentum. When extracting power from the rotational motion, the force that is rotating the blades leads to an opposing and equal magnitude force on the fluid, resulting in the fluid rotating in the wake of the turbine. Since there needs to be conservation of power, the potential power that can be extracted is lowered, leading to wake rotation loss. From Fig. 3 the wake rotation loss is seen to affect the root region of turbine. This is further investigated in Fig. 4, where the wake rotation loss factor is plotted as a function of the local tip speed ratio () for different values of the local loading (CLT).
The effect of changing the local loading is seen to have a limited effect. From Fig. 4b the wake rotation loss is seen to be insignificant for , which is approximately from the mid-span and outwards for a modern utility scale turbine with λ≈7–10. From the example in Fig. 3 the wake rotation loss is seen to have the smallest impact on the global power with .
The tip loss is the power loss associated with the rotor having a finite number of blades and not acting as an actuator disc with an infinite number of blades. This effect is captured in the tip loss factor (F), which is further described in Sect. 2.5. The tip loss model that is used in this paper captures the impact on the induced velocities in the rotor plane, with the additionally induced velocities coming from the vorticity released at the tip of the blade. The most important parameter for tip loss is the tip speed ratio (λ), with the tip loss getting smaller with increasing tip speed ratio. From Fig. 3 the tip loss is seen to affect the power at the tip as one might expect (hub or root loss was not included but easily could have been).
The viscous loss is simply the loss associated with the viscous drag from the airfoil profile. The viscous loss is found to be linear in inverse glide ratio (), loading (CLT) and local tip speed ratio (), with a larger value for each of them leading to a larger loss. In contrast to both wake rotation loss and tip loss, a larger tip speed ratio is found to result in larger losses. As a consequence, there will exist an optimal tip speed ratio since either extreme (λ→0, λ→∞) will lead to 0 or negative power. This is further discussed in Sect. 3.3. From Fig. 3 the loss is seen to increase towards the tip since the local tip speed ratio () is increasing. Viscous loss is seen to be the most significant loss of the three, although it should be noted that a glide ratio of is a fairly low value for a realistic modern rotor design and is here chosen to make the loss easily visible for the figure.
The tip loss factor is commonly implemented for BEMs, and although some different tip correction has been proposed, the tip loss model by Glauert is a common one to use, and it is also the one used here. It is given as
where B is the number of blades (which for simplicity is set to three throughout the paper). The Glauert tip loss model leads to a recursive problem due to the mutual dependence between CLT, sin ϕ and F. It is not possible to find an explicit equation for F in terms of CLT, but it is possible to find a iterative scheme that can solve for F. The iterative scheme is given as
where an acceptable tolerance for F is reached with at most 30 iterations () where a good initial guess for F would be 1.
The Glauert tip loss model breaks the explicit relationship between CLT and CLP in Eq. (19), since F needs to be solved through iterations. An explicit relationship could be obtained by using the Prandtl tip loss model instead, which is given as
which does not depend on CLT. Using Prandtl's tip loss model, the effect on the tip loss is found to be larger compared to Glauert's tip loss model, but investigating it further here is outside the scope of this paper.
Although the Prandtl tip loss model is much simpler and easier to implement, the Glauert model is used throughout this paper as it is the one used for the BEM validation later in Sect. 4.2.
Whether or not to include drag when computing the thrust loading (including drag in Eq. 13) is still a standing question for the derivation of BEM, and excluding drag in the derivation here should not be seen as an argument in favor of this approach but merely as a consequence that makes the mathematical derivation simpler as well as the resulting equation for local power (Eq. 19). If drag should be included for the computation of the induced velocities, it should be noted that it is not just a matter of including drag in Eq. (13), as the closure equation in Eq. (18) also needs an additional to be added on the right-hand side if it should be consistent with the BEM described in Ning (2014).
As a consequence of excluding drag from Eq. (13), there arises a difference between the forces seen from the blade and the forces seen by the air, where Eq. (13) is the thrust force as seen by the air and where Eq. (14) is the tangential force seen by the blade. Often in optimization it is the thrust load at the blade that is of interest, and Eq. (24) shows the relationship between the two local thrust loadings:
where CLT is the local thrust as seen by the air.
In this section, a method for computing the gradients for RIAD is presented. The gradients are then used for power optimization. First, it is applied for loading optimization for maximum power, and it is then further extended for optimization with respect to tip speed ratio and loading. In the end, a discussion of how optimization with RIAD fits within the current state of the art is given.
The local-power equation (Eq. 19) is an analytical expression (with some complications for the tip loss factor), and in principle it is straightforward to compute the gradient for any of the input variables. But it is tedious and error-prone, and the tip loss factor makes it fairly complicated. The complex step method (Martins et al., 2000) is therefore found to be an easier way to compute the gradients, without any loss of accuracy. It also has the benefit that little additional work needs to be done when Eq. (19) is implemented to compute the gradients. The conceptual idea behind the complex step method is fairly simple, and for the sake of making the reader familiar with it, it is summarized here. For a proper description, the reader is referred to Martins et al. (2000).
The complex step method is based on the observation that the Taylor series expansion of an analytical function with a complex step (or perturbation) gives the following (taking Eq. (19) as an example with a step in CLT):
where i is the complex unit and h the step size. Taking the imaginary component of Eq. (25) and dividing by the step size, the approximate gradient can be found as
Equation (26) is seen to have some similarity to computing the gradient through finite difference. But the key difference is that finite difference requires the difference between two function evaluations, which leads to a rounding error. For finite difference there is an optimal step size (h) where the combination of the truncation error (𝒪(h2)) and the rounding error is as small as possible. This is not the case for the complex step method, where the rounding error is eliminated by not computing a difference between two function evaluations and the step size can be arbitrarily small. Using a step size of the truncation error () is found to be smaller than machine precision (10−16), and the gradient is therefore accurate to machine precision. This method applies to any analytical expression, but special care should be taken with functions that might lead to undesirable effects for complex numbers like the absolute value function or similar. This is not of concern for Eq. (19).
The problem of maximizing CP with respect to CLT is a classic problem that will be used here to demonstrate how easily RIAD can solve this problem. It is thought to show the strength of using RIAD for optimization as opposed to using a regular BEM since the solutions is fairly easy to find.
The CP maximizing problem can be stated as
where the boldface signifies that it is a function/vector changing with span (). Since the model assumes radial independence the maximization can be moved within the integration for CP (Eq. 8), and the maximization can be used for CLP at each span location () independently, which can be stated as
Since Eq. (19) is a smooth function, the optimization problem can be simplified as finding the root with respect to CLT for the following equation:
where is computed with Eq. (26), and the outcome from the optimization is the CLT that maximizes CLP, which is referred to as CLT,opt and CLP,opt respectively. The root for Eq. (29) can be found with common root solving algorithms like bisection or Brent's method, eliminating the need for an optimizer to solve the problem. It should be noted that to get the true gradient when including Glauert's tip loss factor (described in Sect. 2.5), the complex step should be included when solving for the tip loss factor (Eqs. 21 and 22); otherwise the effect of the tip loss on the gradient will not be included.
Applying the optimization with similar input as for Fig. 3 , the resulting CLT,opt and CLP,opt are shown in Fig. 5. CLT,opt is seen to be mostly affected at the root and tip of the rotor, compared to the Betz–Joukowsky optimum. At the root CLT,opt is tending to a value of and at the tip going towards 0. The mid-span is seen to have a slightly decreasing slope. For CLP,opt the three losses (wake rotation loss, tip loss, viscous loss) are highlighted as the shaded region. CLP,opt is seen to have a local maximum at from where viscous loss and later tip loss are seen to grow for increasing . For decreasing the wake rotation loss is seen to increase.
The loading optimization in Sect. 3.2 required two inputs to maximize CP, but in this section the optimization will be extended to also include the tip speed ratio as an optimization parameter, leaving only the glide ratio as an input. The optimization problem can be stated as
Using the same assumption as for the loading optimization the optimization for CLT can be solved as described in Sect. 3.2 assuming that CLT,opt is for a fixed λ. A nested optimization can therefore be stated as
The solution to the above optimization problem can be found by solving
where the bounding region is found by observing that λopt has an approximate proportional behavior of , and the limits are simply determined to contain the optimal solution. The outcome from the optimization is the tip speed ratio that maximizes the power coefficient. They are referred to as λopt and CP,opt. To compute the gradient of CP, Eq. (8) is used, and the complex step is applied as follows:
where for practical implementation the problem is discretized along the span, and the integration can be performed using the trapezoidal rule. As was the case for loading optimization, the problem can be solved by the use of root-solving algorithms like bisection or Brent's method.
In Fig. 6 the result of solving the optimization problem for optimal tip speed ratio with varying spanwise constant glide ratio is shown. As expected, the optimal tip speed ratio is seen to increase as the glide ratio increases due to the balance between the viscous losses (increases with increasing λ) and wake rotation losses as well as tip losses (decreases with increasing λ). It is interesting to note the significant impact of including tip loss on λopt, with the inclusion of tip loss leading to a larger λopt. Four points along the λopt curve in Fig. 6 are highlighted, showing diagrams of the local thrust and local power with the difference to the Betz–Joukowsky limit shown by the shaded red region.
The associated CP,opt is shown in Fig. 7. The three power loss contributions are shown as shaded regions, and the viscous loss is seen to be the most significant regardless of the glide ratio. Tip loss is seen to be the second most significant loss, at least for , which anyway would be a very low value for a modern utility scale wind turbine. The slope of CP,opt is seen to become flat for large values of the glide ratio, and the improvement in CP,opt from to is %, and that is with a glide ratio improvement of %.
The results of CP maximization presented in Figs. 5–7 are not in themselves novel results. Similar results have been shown by Wilson et al. (1976, Sect. 3.1–2), Manwell et al. (2010, Sect. 3.9), Sørensen (2016, chap. 5) and Jamieson (2018, Sect. 1.9). The novelty is the ease with which these results can be obtained and the extent to which this method can be applied. In all of the mentioned works, the optimization method relies on excluding mechanisms like rotational effects, drag, or tip loss or an assumption of constant axial induction to find a solution without the use of an optimizer. For Wilson et al. (1976), Manwell et al. (2010) and Sørensen (2016) drag is excluded for the induced velocity as was also done within this paper. However, the reason to exclude drag from the induced velocity within this paper is for the derivation of RIAD to be simpler. Including drag in the induced velocity will make Eq. (19) more complicated, but the optimization methods presented in Set. 3.2 and 3.3 would still be applicable. A large part of why RIAD is easy to use and implement for optimization is the use of the complex step method, which is arguably not an invention of RIAD, and it could as well be applied for a regular BEM, with the same result, although the optimization would be more convoluted. RIAD is established on a better BEM parametrization dedicated to optimization when solving load-constrained power optimizations as is shown in Part 2 of this paper (Loenbaek et al., 2021).
Section 2.3 presented the connection between inputs, such as local thrust, tip speed ratio and glide ratio, for rotor power performance, and Sect. 3.2 and 3.3 presented the power optimization. But the presented methods only contain information about the loading and power at the actuator disc, where this section establishes the connection between the RIAD inputs and the blade planform, such as blade chord and twist. The blade planform is then used as an input for a BEM solver to validate that RIAD and BEM are equivalent formulations.
To compute the chord it is seen that additional inputs are required, such as lift coefficient (Cl), the number of blades (B) and rotor radius (R).
An equation for twist can be found in much the same way by using Eq. (16) for tan ϕ and applying the closure equations. Combining it with , the following equation for the twist can be found:
To show that RIAD is an equivalent formulation of the BEM equations, a planform design is created through Eqs. (36) and (37) and evaluated with a BEM solver. The BEM solver used for the validation is CCBlade (Ning, 2014).
Running CCBlade requires an airfoil polar (Cl, Cd vs. α), and to keep it as simple as possible a single airfoil polar is used all along the blade span. The airfoil is taken to be FFA-W3-301 (Bjorck, 1990) with the aerodynamic data from Bak et al. (2013) (). The design point for the polar was for simplicity taken as the angle of attack with maximum glide ratio, resulting in the following:
Running the optimization for λopt as described in Sect. 3.3 gave λopt=8.4, which is the tip speed ratio used for the design. To give the design some real dimensions, a rotor radius of R=50 m is used.
The resulting planform design can be seen in Fig. 8. Using the planform design as input for CCBlade as well as the other inputs, the resulting local thrust and local power were found from CCBlade. A comparison between RIAD and CCBlade can be seen in Fig. 9. Figure 9a shows the values for CLT and CLP for both the solvers. Figure 9b shows the difference between RIAD and CCBlade for both CLT and CLP using a log scale on the y axis, both as a function of the normalized radius. In the root region the two methods are seen to agree to machine precision but with an error that is growing towards the tip and reaching a difference of 10−4 (which is still three significant digits). The growing error is found to disappear if tip loss is excluded (agrees to machine precision), and the difference is also seen to disappear if the drag is included for the induction as well as including the tip loss. The difference is therefore attributed to some small implementation difference regarding the tip loss, but as the difference is anyway insignificant, the difference is not investigated any further. It should be noted that for the comparison with CCBlade, it is CLT,blade that is used, where CLT,blade was discussed in Sect. 2.6.
A rotor performance model called the Radially Independent Actuator Disc (RIAD) model was presented. It is a different parametrization of the blade element momentum (BEM) equations which is found to be better for wind turbine optimization. The model relates the local rotor power output (local-power coefficient – CLP) to the local rotor loading input (local-thrust coefficient – CLT) at a given radial station (). The model is a simple equation, shown in Eq. (19), from which different physical effects can easily be interpreted, such as wake rotation loss, tip loss and viscous loss.
A method to compute gradients for RIAD was presented, through the use of the complex step method, which allows us to compute the gradient to machine precision with minimum additional work required.
The gradients were used for classical power coefficient (CP) maximization, which was first applied for loading optimization (CLT along the span) for a given tip speed ratio and glide ratio. The optimization was then extended for the combined optimization of the tip speed ratio and loading, leading to a nested optimization for CP which only requires the glide ratio along the span as an input. Using a spanwise constant glide ratio it was shown that viscous loss is the most significant loss, regardless of the value of the glide ratio. The optimization results by themselves have been done before, but the novel development is the ease with which the optimal result can be achieved and the extent to which the method can be applied. But the real strength of using RIAD for optimization is for load constraint rotor optimization as described in Part 2 (Loenbaek et al., 2021).
The relationship between local thrust along the span and the blade chord and twist was presented, and they were used to create the input for validation with a BEM solver (Ning, 2014). The difference between the two methods was found to agree to three significant digits with a likely small implementation difference for the tip loss modeling. In this way, it was shown that RIAD and BEM are equivalent, with the difference being the parametrization.
Code is not publicly available and can not be shared.
The validation data are from https://www.hawc2.dk/Download/HAWC2-Model/DTU-10-MW-Reference-Wind-Turbine (last access: August 2020) (HAWC2, 2020) (version 9.1).
KL came up with the concept and main idea, as well as performed the analysis. All authors have interpreted the results and made suggestions for improvements. KL prepared the paper with revisions from all co-authors.
The authors declare that they have no conflict of interest.
We would like to thank Innovation Fund Denmark for funding the industrial PhD project which this article is a part of.
We would like to thank all the former employees at Suzlon Blade Sciences Center for being a great source of motivation with their interest in the results.
We would like to thank Mads Holst Aagaard Madsen from DTU Risø for the inspiration to use the complex step method.
This research has been supported by Innovation Fund Denmark (grant no. 7038-00053B).
This paper was edited by Alessandro Bianchini and reviewed by Peter Jamieson and two anonymous referees.
Bak, C., Zahle, F., Bitsche, R., Yde, A., Henriksen, L. C., Nata, A., and Hansen, M. H.: Description of the DTU 10 MW Reference Wind Turbine, DTU Wind Energy Report-I-0092, DTU, Denmark, Rosklilde, 1–138, https://doi.org/10.1017/CBO9781107415324.004, 2013. a
Betz, A.: Wind-Energie, und ihre aus- nutzung durch Windmuehlen (Wind energy and their utilization by windmills), Vandenhoeck, Gottingen, 1926. a
Bjorck, A.: Coordinates and Calculations for the FFA-Wl-xxx, FFA-W2-xxx, FFA-W2-xxx and FFA-W3-xxx Series of Airfoils for Horizontal Axis Wind Turbines, Tech. rep., The Aeronautical Research Institute of Sweden, Bromma, Sweden, 1990. a
HAWC2: DTU 10-MW Reference Wind Turbine, available at: https://www.hawc2.dk/Download/HAWC2-Model/DTU-10-MW-Reference-Wind-Turbine, last access: August 2020. a
Loenbaek, K., Bak, C., and McWilliam, M.: A method for preliminary rotor design – Part 2: Wind turbine Optimization with Radial Independence, Wind Energ. Sci., 6, 917–933, https://doi.org/10.5194/wes-6-917-2021, 2021. a, b, c, d
Manwell, J. F., McGowan, J. G., and Rogers, A. L.: Aerodynamics of Wind Turbines, in: Wind Energy Explained, 2, John Wiley & Sons, Ltd, Chichester, UK, 91–155, https://doi.org/10.1002/9781119994367.ch3, 2010. a, b
Martins, J. R., Kroo, I. M., and Alonso, J. J.: An automated method for sensitivity analysis using complex variables, in: 38th Aerospace Sciences Meeting and Exhibit, January 2000, Reno, USA, https://doi.org/10.2514/6.2000-689, 2000. a, b
Okulov, V. L. and van Kuik, G. A.: The Betz-Joukowsky limit: on the contribution to rotor aerodynamics by the British, German and Russian scientific schools, Wind Energy, 15, 335–344, https://doi.org/10.1002/we.464, 2012. a, b
Perez-Moreno, S. S., Zaaijer, M. B., Bottasso, C. L., Dykes, K., Merz, K. O., Réthoré, P. E., and Zahle, F.: Roadmap to the multidisciplinary design analysis and optimisation of wind energy systems, J. Phys.: Conf. Ser., 753, 062011, https://doi.org/10.1088/1742-6596/753/6/062011, 2016. a
Zahle, F., Tibaldi, C., Verelst, D. R., Bak, C., Bitsche, R., and Blasques, J. P. A. A.: Aero-Elastic Optimization of a 10 MW Wind Turbine, in: Proceedings – 33rd Wind Energy Symposium, January 2015, Kissimmee, Florida, 201–223, 2015. a