A Method for Preliminary Rotor Design-Part 1: Radially Independent Actuator Disk 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 model (RIAD) 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 3 losses: wake-rotation-loss, tip-loss and viscous-loss. The gradients for 5 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 at which it can be applied for rotor optimization, and especially load constraint power optimization as it is described in Loenbaek et al. (2020). The relationship between the RIAD input and the rotor chord and twist is established and it is validated against a BEM solver. Keyword: BEM, Rotor Optimization 10


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), Bottasso et al. (2010). 15 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 20 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 theory (BEM) (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 25 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 30 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 35 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.

40
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. As this paper presents an aerodynamic model to be used for optimization, this paper is Part 1 where the model is used for load constraint power optimization in Part 2 (Loenbaek et al., 2020).

45
In the following the Radially Independent Actuator Disc (RIAD) model is presented which start 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 is discussed. 50 2 https://doi.org/10.5194/wes-2020-98 Preprint. Discussion started: 6 October 2020 c Author(s) 2020. CC BY 4.0 License.

Relationship between global and local coefficients
Starting with the Fundamental theorem of calculus the following equation can be made 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. Introducing 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 ): Combining the equations the following equations can be found for C T and C P : Wherer is the Normalized-Radius (r = r/R). A sketch of the relationship between the local and global coefficients can be seen in figure 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 70 assumed to originate from a 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 made:

75
Where ∂D/∂r = 0 in equation 9 to exclude drag from induction. Using the common normalization of ∂L/∂r and ∂D/∂r together with equation 5 and 6 the following equations can be found: Where σ is the rotor solidity (σ = Bc/2πr) andṼ rel the normalized relative velocity. A sketch of the air-flow and forces can be seen in figure 2.
Combining equation 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, a p is the tangential induction factor.

90
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 BEM's. The closures are given as: Ning, 2014, p. 4 eq. 2) (17) a(1 − a) = λ 2r2 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 section 2.5. For now F ∈]0, 1] with F → 0 asr → 1. Combining equation 17, 18 with equation 16 an explicit equation for tan φ in terms of C LT can be found. This leads to an explicit equation 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 sketch showing the impact of different losses for a specific input 105 can be seen in figure 3. The 1D power is the classical 1D momentum theory result by Betz and Joukowsky (Okulov and van Kuik, 2012), but here applied for radially independent streamtubes 1 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, 110 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 figure 3 the wake-rotation-loss is seen to affect the root region of turbine. This is further investigated in figure 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 figure 4 b) the wake-rotation-loss is seen to be insignificant for λr > 5 which is approximately from the 115 mid-span and outwards for a modern utility scale turbine with λ ≈ 7 − 10. From the example in figure 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 captured in the tip-loss-factor (F ) which is further described in section 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 120 induced velocities coming from the vorticity released at the tip of the blade. The most important parameter for tip-loss, is the Figure 4. Significance of wake-rotation-loss. a) wake-rotation-factor vs. local-tip-speed-ratio (λr). Black dashed line is the limit at 1/2. b) difference between the wake-rotation-factor and the limit at 1/2 vs. local-tip-speed-ratio (λr). Notice that the y-axis is log-scale. Red dashed line is at 99% of the limit. Notice that the solid lines are for different values of CLT . Higher CLT leads to higher wake-rotation-loss.
tip-speed-ratio (λ), with the tip-loss getting smaller with increasing tip-speed-ratio. From figure 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 125 to a larger loss. Opposed 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 section 3.3. From figure 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 3, 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 BEM's, and although some different tip-correction has been proposed, the tip-loss model by Glauert is a common one to used and it is also the one used here. It is given as: Sørensen, 2016, p. 132 eq. 8.29) recursive problem due to the mutual dependence 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 140 F would be 1.
The Glauert tip-loss model, breaks the explicit relationship between C LT and C LP in equation 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 out of scope for 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 section 4.2.

Blade loading without drag in induction 150
Whether or not to include drag when computing the thrust loading (including drag in equation 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 it makes the mathematical derivation simpler as well as the resulting equation for local power (equation 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 equation 13, as the closure equation in equation 18 also needs an additional λr C d C l (a + a p ) to be 155 added on the right-hand-side, if it should be consistent with the BEM described in Ning (2014).
A consequence of excluding drag from equation 13 there arises a difference between the forces seen from the blade, and the forces seen by the air where equation 13 is the thrust force as seen by the air and where equation 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 equation 24 shows the relationship between the two local-thrust loading's: Where C LT 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 wrt. tip-speed-165 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 (equation 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 170 way to compute the gradients, without any loss of accuracy. It also has the benefit that little additional work needs to be done when equations 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 175 complex step (or perturbation) gives the following (taking equation 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 equation 25 and dividing by the step size the approximate gradient can be found as: 180 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 185 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 equation 19.
The problem of maximizing C P w.r.t. 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 bold-face 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 (equation 8) and the maximization can be made for C LP at each span location (r) independently, which can be stated as: Since equation 19 is a smooth function the optimization problem can be simplified as finding the root wrt. C LT for the following 200 equation: Where ∂C LP ∂C LT is computed with equation 26, and the outcome from the optimization is the C LT that maximizes C LP , which is 205 referred to as C LT,opt and C LP,opt respectively. The root for equation 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 noticed that to get the true gradient when including Glauert's tip-loss-factor (described in section 2.5), the complex step should be included when solving for the tip-loss-factor (equations 21, 22) otherwise the effect of the tip-loss on the gradient will not be included.
Applying the optimization with similar input as for figure 3 λ = 7, C l C d = 40 , the resulting C LT,opt and C LP,opt is shown 210 in figure 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 3 losses (wake-rotation-loss, tip-loss, viscous-loss) are highlighted as shaded region. C LP,opt is seen to have a local maximum atr = 0.31 from where viscous-loss and later tip-loss is 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 section 3.2 required two inputs λ, Cl Cd 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 Figure 5. Optimal local-thrust (CLT ) and optimal local-power (CLP ) vs. normalized radius (r), with λ = 7 and span-wise constant C l C d = 40 as input.
optimization problem can be stated as: Using the same assumption as for the loading optimization the optimization for C LT can be solved as it was described in section 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 from experience. 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 , equation 8 is used and the 11 https://doi.org/10.5194/wes-2020-98 Preprint. Discussion started: 6 October 2020 c Author(s) 2020. CC BY 4.0 License.
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 it 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 figure 6 the result of solving the optimization problem for optimal tip-speed-ratio with varying span-wise constant glide-235 ratio is shown. As expected, the optimal tip-speed-ratio is seen to increase as the glide-ratio increase due to the balance between Figure 6. Optimal tip-speed-ratio (λopt) vs. span-wise constant glide-ratio ( C l C d ). Both with and without tip-loss included in the optimization. For the case with tip-loss included, 5 points are highlighted, showing sketches of the local-thrust (CLT ) and the local-power (CLP . the viscous-losses (increases with increasing λ) and wake-rotation-losses as well as tip-losses (decreases with increasing λ).
It is interesting to notice the significant impact of including tip-loss on λ opt , with the inclusion of tip-loss leading to a larger λ opt . 4 points along the λ opt curve in figure 6 is highlighted, showing sketches of the local-thrust and local-power with the difference to the Betz-Joukowsky limit shown by shaded red region. The associated C P,opt is shown in figure 7. The 3 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 values 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%.
The novelty is the ease at which these results can be obtained and the generality at 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 to find 250 a solution without the use of an optimizer. Common to them all is the exclusion of drag for the induced velocity as it 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 equation 19 more complicated, but the optimization methods presented in section 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 255 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 it is shown in Part 2 of this paper (Loenbaek et al., 2020).
4 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 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).

270
An equation for twist can be found in much the same way, by using equation 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 equations 36 and 275 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.
Running the optimization for λ opt as described in section 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 = 50m is used.

285
The resulting planform design can be seen in figure 8. Using the planform design as input for CCBlade as well as the other inputs, the resulting local-thrust and local-power was found from CCBlade. A comparison between RIAD and CCBlade can be seen in figure 9. a) shows the values for C LT and C LP for both the solvers. b) 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 is seen to agree to machine precision, but with an error that is growing towards the tip and reaching a 290 difference of 10 −4 (which is still 3 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 section 2.6. Figure 9. a) local-thrust (CLT ) and local-power (CLP ) for RIAD (line) and CCBlade (dots). b) the difference between the two methods for local-thrust (∆CLT ) and local-power (∆CLP ). The difference is seen to increase towards the tip, but the agreement is still within 3 significant digits and is therefore thought to be insignificant. The difference is likely a small implementation difference within tip-loss modeling.
A rotor performance model called Radially-Independent-Actuator-Disc model (RIAD) 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 equation 19, from which different physical effects 300 can easily be interpreted, such as wake-rotation-loss, tip-loss and viscous-loss.
A method to computing gradients for RIAD was presented, through the use of the complex-step-method, which allows to compute the gradient to machine precision with a minimum of additional work required.
The gradients were used for classical power-coefficient (C P ) maximization, which was first applied for loading optimization (C LT along the span) for a given tip-speed-ratio and glide-ratio. The optimization was then extended for combined optimization 305 of 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 span-wise constant glide-ratio it was shown that viscous-loss is the most significant loss, regardless of the value of glide-ratio. The optimization results by themself have been done before, but the novel development is the ease at which the optimal result can be achieved and the generality at which the method can be applied. But the real strength of using RIAD for optimization is for load-constraint rotor optimization as it is described in Part 2 (Loenbaek et al., 2020).

310
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 3 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.

Rotor Local variables
Variables that are scalars at a given radius location (r). Bold-face variables indicates it is a function or vector changing with radius.

Symbol
Description Unit x Bold face local variables symbolizes a function or vector changing with the rotor radius (r) - 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 Axial induction factor a p Tangential induction factor -C l Lift coefficient - Airfoil glide ratio -