A method for preliminary rotor design – Part 1: Radially Independent Actuator Disc model

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.


Introduction
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 (C LT ) and the glide ratio (C l /C d ) as well as the global tip speed ratio (λ). The equivalent input for a BEM at each radial station is the design lift coefficient (C l ), the design drag coefficient (C d ) 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 (C LT ) 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 loadconstrained power optimization.

The RIAD model
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.

Relationship between global and local coefficients
Starting with the fundamental theorem of calculus the following equation can be created for the global values for thrust (T ) and power (P ): where ∂T ∂r is the thrust loading density and ∂P ∂r is the power density.
The classical non-dimensional relations for thrust (C T ) and power (C P ) as well as the local equivalent which is introduced as the local-thrust coefficient (C LT ) and local-power coefficient (C LP ) are introduced as follows: Combining the equations, the following equations can be found for C T and C P : wherer is the normalized radius (r = r/R). A diagram of the relationship between the local and global coefficients can be seen in Fig. 1.

Relationship between local coefficients
To establish a relationship between the local coefficients C LT and C LP , 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 1 ωr ∂P ∂r . Introducing the local lift density ∂L ∂r and drag density ∂D ∂r as well as the local flow angle φ, the following equations can be created: 1 ωr where ∂D/∂r = 0 in Eq. (9) to exclude drag from induction. Using the common normalization of ∂L/∂r and ∂D/∂r, together with Eqs. (5) and (6), the following equations can be found: where σ is the rotor solidity (σ = Bc/2π r) andṼ rel the normalized relative velocity. A diagram of the airflow and forces can be seen in Fig. 2. Combining Eqs. (13) and (14) an equation for C LP can be found as This is a general equation for the relationship between the C LT and C LP , but it is implicit since tan φ depends on C LT .

Explicit equation for C LP (closure relation between C LT and tan φ)
The local flow angle φ is given from the induced velocities (Sørensen, 2016, p. 101, Eq. 7.3) as where a is the axial induction factor and a p is the tangential induction factor. Equation (16) introduces two new variables (a, a p ). In order to make an explicit equation for tan φ as a function of C LT two equations relating C LT to a and a p 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 (Ning, 2014, p. 4 Eq. 2 a(1 − a) = λ 2r 2 a p 1 + a p (Sørensen, 2016, p. 50 Eq. 4.36), where F is the tip loss factor, which is further described in Sect. 2.5. For now F ∈]0, 1] with F → 0 asr → 1. Combining Eqs. (17) and (18) with Eq. (16) an explicit equation for tan φ in terms of C LT can be found. This leads to an explicit equation for C LP :  Equation (19) is the main result of this section. It states an explicit relationship between the local power and the local loading.

Physical interpretation and input sensitivity
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 1D power is the classical 1D momentum theory result by Betz and Joukowsky (Okulov and van Kuik, 2012), but here it is applied for radially independent stream tubes, noting the following relation: 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 (λr) for different values of the local loading (C LT ). 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 λr > 5, 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 C P = −0.01.
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 cap-tured 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 (C d /C l ), loading (C LT ) and local tip speed ratio (λr), 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 (λr) 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 C l C d = 40 is a fairly low value for a realistic modern rotor design and is here chosen to make the loss easily visible for the figure.

Tip loss factor
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 Sørensen, 2016, p. 132 Eq. 8.29), 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  (b) Difference between the wake rotation factor and the limit at 1/2 vs. local tip speed ratio (λr). Note that the y axis is log scale. The red dashed line is at 99 % of the limit. Note that the solid lines are for different values of C LT . Higher C LT leads to higher wake rotation loss.
between C LT , sin φ and F . It is not possible to find an explicit equation for F in terms of C LT , 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 (|F i+1 − F i | < 10 −9 ) where a good initial guess for F would be 1.
The Glauert tip loss model breaks the explicit relationship between C LT and C LP 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 C LT . 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.

Blade loading without drag in induction
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 λr C d C l (a + a p ) 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 C LT is the local thrust as seen by the air.

Gradients for RIAD and power optimization
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.

Gradients with complex step
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 C LT ): 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 (O(h 2 )) 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 h = 10 −9 the truncation error (O(h 2 ) ≈ 10 −18 ) 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).

Loading optimization for max power
The problem of maximizing C P with respect to C LT 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 C P maximizing problem can be stated as where the boldface signifies that it is a function/vector changing with span (r). Since the model assumes radial independence the maximization can be moved within the integration for C P (Eq. 8), and the maximization can be used for C LP at each span location (r) 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 C LT for the following equation: where ∂C LP ∂C LT is computed with Eq. (26), and the outcome from the optimization is the C LT that maximizes C LP , which is referred to as C LT,opt and C LP,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 λ = 7, C l C d = 40 , the resulting C LT,opt and C LP,opt are shown in Fig. 5. C LT,opt is seen to be mostly affected at the root and tip of the rotor, compared to the Betz-Joukowsky optimum. At the root C LT,opt is tending to a value of 3/4 and at the tip going towards 0. The mid-span is seen to have a slightly decreasing slope. For C LP,opt the three losses (wake rotation loss, tip loss, viscous loss) are highlighted as the shaded region. C LP,opt is seen to have a local maximum at r = 0.31 from where viscous loss and later tip loss are seen to grow for increasingr. For decreasingr the wake rotation loss is seen to increase.

Tip speed ratio optimization for maximum power
The loading optimization in Sect. 3.2 required two inputs λ, C l C d to maximize C P , 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 C LT can be solved as described in Sect. 3.2 assuming that C LT,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 √ C l /C d , 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 C P,opt . To compute the gradient of C P , 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 rootsolving 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 C P,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 C l C d > 25, which anyway would be a very low value for a modern utility scale wind turbine. The slope of C P,opt is seen to become flat for large values of the glide ratio, and the improvement in C P,opt from C l C d = 100 to C l C d = 150 is C P,opt = 3.5 %, and that is with a glide ratio improvement of C l C d = 50 %.

Compared to other work
The results of C P 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).

From RIAD to rotor blade planform
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.

Equations for chord and twist
An equation for chord can be found from Eq. (13) while applying the closure equations (Eqs. 17 and 18) for cos φ, resulting in the following: To compute the chord it is seen that additional inputs are required, such as lift coefficient (C l ), 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 φ = α + θ twist + θ pitch , the following equation for the twist can be found:

Validation with BEM
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 (C l , C d 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) (Re = 10 × 10 6 ). 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 C LT and C LP for both the solvers. Figure 9b shows the difference between RIAD and CCBlade for both C LT and C LP 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 C LT,blade that is used, where C LT,blade was discussed in Sect. 2.6.

Conclusion
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 -C LP ) to the local rotor loading input (local-thrust coefficient -C LT ) at a given radial station (r). 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 (C P ) maximization, which was first applied for loading opti-   The difference between the two methods for local thrust ( C LT ) and local power ( C LP ). The difference is seen to increase towards the tip, but the agreement is still within three significant digits and is therefore thought to be insignificant. The difference is likely a small implementation difference within tip loss modeling. mization (C LT 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 C P 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. Drag loading density N m −1 r Normalized rotor radius variable (r = r R ) -C LT Local-thrust coefficient (normalized ∂T /∂r -taken as the loading seen by the air) -C LT,blade Local-thrust coefficient as seen by the blade (including drag) -C LT,opt Local-thrust coefficient that maximizes C LP for a given λ and C l C d -C LP Local-power coefficient (normalized ∂P /∂r) -C LP,opt Optimal local-power coefficient for given λ and C l C d a

A1 Rotor local variables
Axial induction factor a p Tangential induction factor -C l Lift coefficient - Airfoil glide ratio -  Rotor power coefficient with λ opt and C LT,opt λ Rotor tip speed ratio λ = ωR V λ opt Rotor tip speed ratio that maximized C P for a given C l C d -B Number of blades -