Non-Empirical BEM Corrections Relating to Angular and Axial Momentum Conservation

The Blade-Element Momentum (BEM) model for Horizontal-Axis Wind Turbines (HAWTs), although extremely useful, is known to be approximate due to model formulation insufficiencies, for which add-ons and corrections have been formulated over the past many decades. Scrutiny of the axial and azimuthal momentum conservation properties reveals momentum simplifications and absence of momentum sources not included in momentum theory underlying the standard BEM. One aspect relates to azimuthal momentum conservation, the wake swirl. This correction can be expressed analytically. Another aspect relates to axial momentum conservation, the wake expansion. This correction is not analytically quantifiable. The latter correction term is therefore quantified from postprocessing a large number of axisymmetric Actuator Disk (AD) Navier-Stokes computations with systematic variation of disk loading and tip-speed ratio. The new momentum correction terms are then included in the BEM model, and results benchmarked against references. The corrected BEM is derived by re-visiting the governing equations. For a disk represented by a constant-circulation set of blades, the corrected BEM contains no approximation to the underlying conservation laws. The study contributes by bridging the gap between BEM and the axisymmetric AD method for all disk load levels and tip speed ratios relevant for a wind turbine. The wake swirl correction leads to higher power efficiency at lower tip-speed ratios. The wake expansion correction causes a redistribution of the potential for power extraction, which increases on the inner part of the rotor and decreases on the outer part of the rotor. The overall rotor-averaged value of Betz limit is unaffected by the new corrections, but exceeding Betz locally on the innerand mid-section of the rotor is shown to be possible. The two corrections significantly improve the axi-symmetric static BEM modelling accuracy for the radial distributions as well as for the rotor-integrated quantities, by reducing errors, approximately one order of magnitude. The relevance of these corrections for modern multi-MW rotors is quantified and discussed.


Introduction
The disturbance of a flow field from a spinning rotor, e.g., a propeller or a turbine, is caused by the rotor adding energy to the passing fluid (propeller) or retrieving energy from it (turbine).Both situations are of significant engineering interest.This flow disturbance is termed "induction".The most used induction model for a turbine is the Blade-Element Momentum (BEM) model.The axial momentum theory contained herein was initially formulated by Rankine [1] in 1865, and the concept of an actuator disk introduced by Froude [2] in 1889.Joukowski [3] 1918, Betz [4] 1920, and Glauert [5] 1926 further improved what came to be the classic BEM formulation by adding wake swirl to the model.The great virtue of the BEM model is its ability to compute the induction flow field in a simple yet accurate way that would otherwise require use of complex numerical tools, e.g., Computational Fluid Dynamics (CFD).Although the present study has an analytical and generic approach, it relates specifically to wind turbines.Upon a first encounter with the BEM model, it is tempting to underestimate its modelling accuracy, since the list of assumptions is long:

•
Annular element independence is assumed.

•
Uniform steady inflow is assumed.

•
Turbulence is assumed zero.

•
Axial momentum conservation is simplified.
Despite these shortcomings, BEM has proven to be the aerodynamic induction model of choice, and an indispensable workhorse for all stages of wind turbine design: Initial blade layout and static analysis, aero-elastic simulations, stability analysis, loads and performance for product certification, etc.The applicability, accuracy, and computational efficiency of the BEM can be attributed to two key issues:

•
The annular element independency yields just two non-linear equations to be solved iteratively for each of the annuli.The BEM computational efficiency, both in memory and cpu, is therefore many orders of magnitude higher than CFD-based induction models (full CFD, actuator disk, actuator line, vortex lattice methods, or other) which are all globally coupled and much more complex.Besides, the annular element independency assumption can be shown to be exact in the special case where the rotor/blade is completely in-plane, the blade is uniformly lightly loaded, and tip-losses are negligible [6].

•
Add-on correction models have been developed that expands the BEM applicability to practically all wind turbine operating, idling, and stand-still conditions.These contributions are: Yaw-correction formulation devised by Glauert [5] in 1926, wake stall correction also by Glauert [7] in 1935, Prandtl's tip-loss model [7,8], "frozen turbulence" approach by Taylor [9] in 1938 (later to be used with transient BEM for wind turbines), dynamic inflow model by Øye [10] in 1990.Corrections accounting for transient/stall effects of the blade-element part of the BEM include the dynamic stall model by Beddoes-Leichmann [11] in 1989 and the 3D-correction for a rotating blade in stall by Snel [12] in 1993.Note that other variations of most of these corrections have emerged later.The complete list is long and not presented here.There exists no simple BEM model correction that accounts for rotor/blade out-of-plane effects, but the degree of out-of-plane deflection on wind turbine performance is generally assumed to be limited.The impact from wind turbine blade coning on aerodynamic induction has been investigated, see e.g., [13].
Two BEM assumptions for which corrections have not been commonly established will be the focus of the present work.They relate to angular and axial momentum conservation.The swirl formulation contained in the BEM is simplified due to non-enforcement of angular momentum conservation with the purpose of keeping the annular elements mathematically uncoupled.The swirl of the wake scales inversely with the tip speed ratio.Modern wind turbines operating at high tip speed ratios (λ > 6) therefore have little swirl, so the simplified swirl contained in the classic BEM is sufficiently accurate for many practical purposes.Madsen [14] devised an improved BEM formulation that included the far-wake pressure drop term Glauert purposefully neglected in the original model, although still without imposing angular momentum conservation.This caused a one-way coupling between the annular elements, requiring the numerical solution to start with the innermost annulus and then proceed "outwards" element by element toward the tip.
The simplified axial momentum conservation contained in the BEM is the result of an annular element control volume analysis where the axial momentum is balanced by the far upstream inlet flux, the far wake outlet flux, and the actuator disk axial forcing on the annular element.Inviscid flow is assumed, so no viscosity-driven transfer of momentum is exchanged between adjacent annuli.In the trivial case of an unloaded rotor, there is no pressure drop across the disk and hence no induced pressure field surrounding the disk.However, in the normal case of a loaded rotor, there is a pressure drop across the disk and an induced spatially varying pressure field surrounding the disk.Furthermore, the disk loading will cause an annular element expansion, taking place in the upstream and downstream vicinity of the disk.The inside and outside borders of an expanding annulus will have non-zero axial surface projections where the disk-induced pressure field will exert different pressure forces, thereby adding axial momentum to the annular element.
Glauert commented in [7] that the simplified axial momentum conservation "may imply the neglect of the mutual interference between the various annular elements", but also that the modeling error was "believed to be extremely small in general".Goorjian proved in [15] that the simplified axial momentum conservation, i.e., Equation (1) with zero ∆F s and zero ∆p w , is incorrect when pressure gradients induced by the rotor disk are taken into consideration.The inclusion of ∆F s , named ∆X in their work, was later devised for infinite tip-speed-ratio and constant spanwise loading by Sørensen and Mikkelsen [16] and for a range of tip-speed-ratios, Sørensen [17].Inclusion of the missing momentum conservation terms has also been analyzed and discussed by van Kuik et al. [18], Branlard [19], and van Kuik [20].Xiros [21] reasoned that the mutual momentum impact between adjacent annuli would cause increased potential for power extraction for a highly loaded rotor and even exceed Betz limit.Madsen [14] quantified the assumed influence of annular element axial momentum interexchange by fitting the classical BEM to actuator disk results.It was shown that the obtainable local power efficiency C p indeed can exceed Betz limit on the inner part of the rotor, but that the local power efficiency drops towards the rotor tip, such that the overall obtainable power efficiency C p_rotor seems to be unaffected.
The scope of the two momentum conservation corrections must be clarified.Figure 1 shows, as mentioned, that the BEM model has many constituting parts.The "1D momentum theory" is the starting point, which with the addition of angular momentum conservation, BET, and annular element discretization becomes the "axi-symmetric static BEM model with aligned uniform inflow".This model (shade highlighted on Figure 1) is the one formulated in Section 2 and to which the proposed momentum conservation corrections are applied.It later evolves into the "add-on improved dynamic BEM" upon incorporation of the BEM add-ons as shown on the figure.The "add-on improved dynamic BEM" is the BEM used in aero-elastic simulators throughout the wind turbine industry.It contains the ideal correct model from each of the constituting parts plus the modelling errors from each of these constituting parts.The dominating errors will be case dependent: BET airfoil polar errors, tip-loss errors, yaw correction errors, etc.In many cases will the errors from omitting the two proposed momentum conservation corrections be non-dominating.However, in some cases they might be dominating, especially with the continuous improvement over time of the add-on parts and subsequent reduction of associated modelling errors.
Benchmarking of the axi-symmetric static BEM model is done by comparing against swirled Actuator Disk CFD with zero viscosity and without a turbulence model, i.e., the inviscid incompressible Navier-Stokes governing equations.It is important to realize that with the strict enforcement of momentum conservation from the two proposed corrections, the axi-symmetric static BEM method should compare very well with AD CFD method, and discrepancies should be attributed mainly to discretization errors inherently present in both methods.Benchmarking of the axi-symmetric static BEM model against full 3D CFD, experimental scaletests, or full-scale wind turbine measurements is not fully appropriate in this case, since the noncontainment of all the BEM add-ons (see Figure 1) in the axi-symmetric static BEM model would cause comparison discrepancies from the lack of tip-loss, root-loss, etc. to contaminate the soughtafter discrepancy.While the add-on improved dynamic BEM model should be compared against experimental/physical test results, the analytically formulated axi-symmetric static BEM is void of the more empirically founded add-ons, and therefore better suited for comparison against an equivalent numerical model like the AD CFD, see e.g.[13,14,16,17].
The present paper revisits the set of governing equations beneath the axi-symmetric static BEM model, identifies the missing terms, quantifies the missing terms either analytically or numerically with no fitting involved, and benchmarks the corrected axi-symmetric static BEM against highresolution axi-symmetric actuator disk reference results.Quantification of the error involved from not imposing strict momentum conservation is a main finding of this work.The spanwise error -and correction term -distributions are quantified for all relevant load levels and tip-speed ratios.
The Blade Element Theory (BET) part of the model in Equations ( 1)-( 9) is described for each annular element by a lift force  , a drag force  , and the ratio   ⁄ between the two forces, Equation (4).The relation between  ,  and the axial and azimuthal forces  ,  are described by the schematic of forces on Figure 2 right.Infinite   ⁄ is assumed throughout the investigation for the purpose of comparison with analytical references and other benchmark results.Blade layout distributions for chord, twist, and airfoil type, as well as airfoil specific polar data has been omitted in the BEM formulation, since focus is on fundamental momentum conservation for each annulus.Instead, constant blade circulation, equivalent to constant spanwise rotor loading  has been applied for all BEM computations covering a thrust coefficient range from low to high, i.e. approximately 1.The similarly constant AD CFD load levels are specified in Table 1.
The outline of the paper is as follows: the governing equations of the axi-symmetric static BEM are presented in Section 2. The angular momentum correction which includes Glauerts neglected term is presented in Section 3. Formulation of the axial momentum correction, i.e. the axial momentum added to the annular elements from the disk-induced pressure field is presented in Benchmarking of the axi-symmetric static BEM model against full 3D CFD, experimental scale-tests, or full-scale wind turbine measurements is not fully appropriate in this case, since the non-containment of all the BEM add-ons (see Figure 1) in the axi-symmetric static BEM model would cause comparison discrepancies from the lack of tip-loss, root-loss, etc. to contaminate the sought-after discrepancy.While the add-on improved dynamic BEM model should be compared against experimental/physical test results, the analytically formulated axi-symmetric static BEM is void of the more empirically founded add-ons, and therefore better suited for comparison against an equivalent numerical model like the AD CFD, see e.g., [13,14,16,17].
The present paper revisits the set of governing equations beneath the axi-symmetric static BEM model, identifies the missing terms, quantifies the missing terms either analytically or numerically with no fitting involved, and benchmarks the corrected axi-symmetric static BEM against high-resolution axi-symmetric actuator disk reference results.Quantification of the error involved from not imposing strict momentum conservation is a main finding of this work.The spanwise error -and correction term -distributions are quantified for all relevant load levels and tip-speed ratios.
The Blade Element Theory (BET) part of the model in Equations ( 1)-( 9) is described for each annular element by a lift force F L , a drag force F D , and the ratio L/D between the two forces, Equation (4).The relation between F L , F D and the axial and azimuthal forces F x , F θ are described by the schematic of forces on Figure 2 right.Infinite L/D is assumed throughout the investigation for the purpose of comparison with analytical references and other benchmark results.Blade layout distributions for chord, twist, and airfoil type, as well as airfoil specific polar data has been omitted in the BEM formulation, since focus is on fundamental momentum conservation for each annulus.Instead, constant blade circulation, equivalent to constant spanwise rotor loading F L has been applied for all BEM computations covering a thrust coefficient range from low to high, i.e., approximately 1.The similarly constant AD CFD load levels are specified in Table 1.
Section 4 and quantified numerically in Section 5. Comparisons between results from classical axisymmetric static BEM, corrected axi-symmetric static BEM, and actuator disk CFD is presented in Section 6.The two corrections' performance impact on a modern multi-MW rotor is discussed in Section 7. Summary and perspectives in Section 8 and conclusions in Section 9 finalize the paper.The appendix describes the actuator disk model used for mapping of the annular element pressure momentum term ∆ and benchmarking.

Axi-Symmetric Static BEM Formulation
We assume a rotor disk perpendicular to a passing flow, see Figure 2 (left).The common assumptions apply to the undisturbed flow: steady, uniform, inviscid and incompressible.
Axial momentum conservation over an annular element control volume: Two terms, ∆ and  ∆ are neglected in the classical (Glauert) BEM formulation.∆ is the far wake sustained pressure drop due to the fully developed swirl.It will be quantified in Section 3. The other new term, ∆ , accounts for the area-integrated axially projected pressure forcing exerted on an annular element by the two neighboring annuli, the inner and the outer.It will be presented and quantified in Sections 4 and 5.
Axial forcing exerted by the disk onto the annular element is: Azimuthal (angular) forcing exerted by the disk onto the annular element is: The expression in (3) for  follows from the conservation of angular momentum across the actuator disk for an annulus control volume in which  is the external angular forcing.Note the factor 2 on the RHS of (3).It is simply a consequence of defining  as the average azimuthal flow velocity upon passage of the rotor disk.The exchange of angular momentum between the disk and the fluid passing the disk is due to an assumed constant angular forcing, so the angular acceleration during disk passage will also be constant.When  is the average azimuthal flow velocity over the disk, then the final azimuthal flow velocity, when the flow has just passed the disk, is 2 .
Geometry relation between axial and azimuthal blade element forcing, see Figure 2 (right): The outline of the paper is as follows: the governing equations of the axi-symmetric static BEM are presented in Section 2. The angular momentum correction which includes Glauerts neglected term is presented in Section 3. Formulation of the axial momentum correction, i.e., the axial momentum added to the annular elements from the disk-induced pressure field is presented in Section 4 and quantified numerically in Section 5. Comparisons between results from classical axi-symmetric static BEM, corrected axi-symmetric static BEM, and actuator disk CFD is presented in Section 6.The two corrections' performance impact on a modern multi-MW rotor is discussed in Section 7. Summary and perspectives in Section 8 and conclusions in Section 9 finalize the paper.The appendix describes the actuator disk model used for mapping of the annular element pressure momentum term ∆F s and benchmarking.

Axi-Symmetric Static BEM Formulation
We assume a rotor disk perpendicular to a passing flow, see Figure 2 (left).The common assumptions apply to the undisturbed flow: steady, uniform, inviscid and incompressible.
Axial momentum conservation over an annular element control volume: Two terms, ∆F s and A w ∆p w are neglected in the classical (Glauert) BEM formulation.∆p w is the far wake sustained pressure drop due to the fully developed swirl.It will be quantified in Section 3. The other new term, ∆F s , accounts for the area-integrated axially projected pressure forcing exerted on an annular element by the two neighboring annuli, the inner and the outer.It will be presented and quantified in Sections 4 and 5.
Axial forcing exerted by the disk onto the annular element is: Azimuthal (angular) forcing exerted by the disk onto the annular element is: The expression in (3) for F θ follows from the conservation of angular momentum across the actuator disk for an annulus control volume in which F θ is the external angular forcing.Note the factor 2 on the RHS of (3).It is simply a consequence of defining U θd as the average azimuthal flow Energies 2019, 12, 320 6 of 28 velocity upon passage of the rotor disk.The exchange of angular momentum between the disk and the fluid passing the disk is due to an assumed constant angular forcing, so the angular acceleration during disk passage will also be constant.When U θd is the average azimuthal flow velocity over the disk, then the final azimuthal flow velocity, when the flow has just passed the disk, is 2U θd .
Geometry relation between axial and azimuthal blade element forcing, see Figure 2 (right): Equation ( 4) says that the ratio between axial and azimuthal forcing equals the ratio between the azimuthal and the axial wind speed in the rotating rotor disk reference.This relation only holds when the forcing from the rotor on the wind is perpendicular to the wind, i.e., when the drag force is zero.Non-zero drag is of course possible but leads to a more complicated version of Equation ( 4), see [22] for details, where the drag contribution is contained in the lift-to-drag ratio term L/D: Also note that any of Equation ( 4a) and (4b) conveniently constitutes a simple blade element formulation, where the blade/disk forcing is easily expressed, either as a uniform or a non-uniform distribution over the blade span.The zero-drag version of Equation ( 4) is used throughout the present investigation.The axial and azimuthal annulus forcing from the disk on the passing flow are projections of the summed forces over N B blade elements, Angular momentum conservation downstream of the actuator disk for an annular element control volume: In the axi-symmetrical context no azimuthal wake forcing can exist, so the LHS and RHS of Equation ( 5) are straightforward expressions of the wake's azimuthal momentum influx and outflux respectively.Note that Equation ( 5) is not contained in the classical (Glauert) BEM where the wake expansion is not quantified.Instead, classical BEM uses the simplification 2U θd = U θw , thereby implicitly neglecting wake expansion in the swirl by assuming that r w = r for any annulus.The RHS presence of the far-wake annulus radius r w introduces a one-way coupling between the annuli, where an annulus far-wake radius quite naturally depends on the accumulated far-wake radii of the (inner) annuli it surrounds.The coupling is weak for high tip-speed ratios but will become significant with decreasing tip-speed ratio and increasing swirl.
Mass conservation in annular element control volume upstream of actuator disk: Mass conservation in annular element control volume downstream of actuator disk: Mechanical energy conservation along an annulus streamline upstream of the actuator disk (Bernoulli): U rd in ( 8) is the radial velocity component of the flow passing the disk.Mechanical energy conservation along an annulus streamline downstream of actuator disk (Bernoulli): Energies 2019, 12, 320 The radial velocity components in Equations ( 8) and ( 9) must be included for correctness.The rotor disk is confined to a plane perpendicular to the free stream direction and hence cannot exert any forcing in radial direction.Therefore, the radial velocity component is weakly coupled to the forcing terms driving the dynamics, and mathematically cancels out early during the reduction of Equations ( 1)- (9).Equations ( 1)-( 9) constitute the governing set of equations for the BEM model, and contain no approximation from the underlying conservation laws.This means that, provided well-posedness of the equations, errors should become arbitrarily small with increasing numerical resolution.The reduction of variables and solution algorithm for the BEM model is not included for sake of brevity.Except for the new annular element forcing term ∆F s , a similar set of BEM equations was presented in [22] along with a proposed solution algorithm.

Angular Momentum Correction
The fully developed far wake is purely parallel.The azimuthal (angular) swirl velocity component of the steady inviscid incompressible Navier-Stokes equations in polar coordinates reduces to: The far wake term for the pressure drop due to swirl, ∆p w , can then be expressed as: This expression for the wake swirl-induced pressure drop is also given by Madsen [14] although he evaluates the swirl pressure drop in the disk plane and avoids coupling of the far wake angular velocity U θw to the system of equations through downstream angular momentum conservation (5).The neglected wake pressure term ∆p w was originally mentioned, although only implicitly, by Glauert [7] p. 211, stating that the (classical) BEM axial momentum equation is not correct, as it "neglects the fact that reduced pressure occurs in the wake owing to the rotational motion".

Ideal Swirl Correction
The far wake swirl velocity term, U θw , in Equation ( 11) is computed from the BEM governing set of equations.No swirl diffusion is assumed, although in principle a wake center vortex singularity will form if the disk loading is maintained all the way to the disk center.We term the swirl correction which uses Equation (11) for quantification of ∆p w the "ideal swirl correction".

Solid Body Swirl Correction
However, by choosing the "ideal swirl correction" approach, it is implicitly assumed that the impact of viscosity is negligible, also in the far wake state.Viscous shear stresses will be present in the real wake, primarily (1) in the border zone between the wake and the surrounding free flow, and (2) towards the center of the wake where azimuthal swirl velocities increase steeply, due to trailing root vorticity.The latter viscous shear stress contribution will cause a viscous driven exchange of angular momentum between the annular elements in the wake, tending to level out element differences in angular velocity.For angular momentum decay of trailing vorticity, see e.g., [23].
In a slightly extended approach we assume that angular momentum exchange with the surrounding free flow is negligible, and that internal swirl diffusion will cause the developed swirl in the far wake to be of solid-body-rotation type, such that: Energies 2019, 12, 320 8 of 28 ω solid is constant across all annuli in the far wake state and can be expressed from the condition for angular momentum conservation in the wake.Index i refers to the annular element indexation: Equations ( 12) and ( 13) replace Equation ( 5) in the set of governing BEM equations and Equation (11) still apply.As will be shown in Section 5, this solid-body far wake rotation assumption can lead to a more realistic wake development, (1) by eliminating the ideal far wake vortex core singularity, and (2) by further reducing the discrepancy between BEM and Actuator Disk CFD results for highly swirled rotor flows.
When introducing the impact of viscosity acting primarily on the swirl core in the wake, then from a broader perspective's view, it is clear, that viscosity gradually levels out velocity gradients everywhere as the wake propagates downstream, also radial gradients.So, to clarify, we assume a wake condition, e.g., a few disk diameters downstream where the large swirl gradients have been modified by diffusion, and where the wake expansion has developed (though theoretically, it is an infinite asymptote), but where wake decay and mixing up with the surrounding flow is yet for the most part to happen further downstream.

Axial Momentum Correction
The added correction term ∆F s in the axial momentum conservation Equation (1) accounts for the axially projected area-integrated pressure difference between the inner and outer surface of the annular element control volume.This correction term will vanish in the absence of wake expansion and associated pressure gradients but will be finite for the general case of a loaded HAWT rotor with wake expansion.It is convenient to express ∆F s as a non-dimensional local thrust coefficient times the stagnation pressure times the annulus cross-sectional area at the disk: Wake expansion is largely governed by the rotor thrust coefficient, C t−rotor .Variations in local thrust coefficient across the annuli are less influential.Therefore, and for the sake of simplicity, it is preferred to express the local annular element force coefficient ∆C ts for any element as a function of the rotor thrust coefficient, the relative radius at the disk, and the tip-speed ratio: Rearranging Equation ( 14) and expressing ∆F s in integral form: s_in and s_out are illustrated with dotted blue and full black line respectively on Figure 3 (left).The streamline and pressure contour CFD visualization on Figure 3 (right) seeks to illustrate what the wake expansion correction term is: Axial projection of area-integrated annular element pressure forcing.Using the axi-symmetry property of an annulus leads to: Energies 2019, 12, 320 9 of 28 n•n x will vanish on the undisturbed parallel part of the annulus far upstream and on the fully developed parallel wake part of the annulus far downstream, and likewise for the axial pressure gradients.The indefinite integral limits in (17) can therefore be replaced by finite axial position values upstream and downstream where the flow is parallel and void of axial pressure gradients within reasonable tolerance, e.g., x = ±10R.16) will be non-zero.

Quantification of Axial Momentum Correction Term
Mapping of the annular element force coefficient ∆ , as a function of the rotor thrust coefficient, the relative radius of the tube element at the disk, and the tip-speed ratio, will be computed by brute numerical force: A 2D axi-symmetric Navier-Stokes Actuator Disk model with swirl.Each AD computation consists of four steps:


Step 1: Computation of the AD converged solution on a high-resolution mesh.


Step 2: Post-processing of the annuli streamlines upstream and downstream of the disk, see e.g. Figure 3 right.


Step 3: Re-computation of the AD converged solution on a high-resolution mesh containing the annuli streamlines as open internal boundaries.


Step 4: Post-processing of the annuli force coefficients ∆ by integration around each annulus as expressed in Equation (17).
Nineteen AD computations have been performed.Constant circulation i.e. blade loading, along the blades (azimuthally averaged as an axi-symmetric disk) is applied for all cases.The tip-speed and disk loading combinations are shown in Table 1.Note the difference between the blade loading coefficient,  , and the rotor loading coefficient,  .The former is the blade loading.In the asymptotic limit of  ∞ the blade loading becomes axially aligned and will equal the rotor loading.Otherwise, when  is finite, the blade loading will have an azimuthal component as described by Equation ( 4).The rotor loading,  , is the axial projection of the overall blade loading, so a constant blade loading will lead to decreasing rotor loading when the tip-speed ratio is decreased, because a larger part of the blade loading will point in the azimuthal direction.
The 19 distributions of ∆ versus relative radius are shown on Figure 4 for tip-speeds ratios and rotor loadings.The ∆ magnitude increases with the disk loading, as one would expect.the  16) will be non-zero.
The local force coefficient is a function of spanwise position r R , rotor loading C t−rotor , and tip-speed ratio λ as specified in Equation (15).A large number of actuator disk CFD simulations have been performed and post-processed, such that a complete mapping of the annular element axial force coefficient ∆C ts has become available.Description of the actuator disk method and the post-postprocessing is presented in Appendix A.

Quantification of Axial Momentum Correction Term
Mapping of the annular element force coefficient ∆C ts , as a function of the rotor thrust coefficient, the relative radius of the tube element at the disk, and the tip-speed ratio, will be computed by brute numerical force: A 2D axi-symmetric Navier-Stokes Actuator Disk model with swirl.Each AD computation consists of four steps:

•
Step 1: Computation of the AD converged solution on a high-resolution mesh.

•
Step 2: Post-processing of the annuli streamlines upstream and downstream of the disk, see e.g., Figure 3 right.

•
Step 3: Re-computation of the AD converged solution on a high-resolution mesh containing the annuli streamlines as open internal boundaries.

•
Step 4: Post-processing of the annuli force coefficients ∆C ts by integration around each annulus as expressed in Equation (17).
Nineteen AD computations have been performed.Constant circulation i.e., blade loading, along the blades (azimuthally averaged as an axi-symmetric disk) is applied for all cases.The tip-speed and disk loading combinations are shown in Table 1.Note the difference between the blade loading Energies 2019, 12, 320 10 of 28 coefficient, C t−blade , and the rotor loading coefficient, C t−rotor .The former is the blade loading.In the asymptotic limit of λ = ∞ the blade loading becomes axially aligned and will equal the rotor loading.Otherwise, when λ is finite, the blade loading will have an azimuthal component as described by Equation ( 4).The rotor loading, C t−rotor , is the axial projection of the overall blade loading, so a constant blade loading will lead to decreasing rotor loading when the tip-speed ratio is decreased, because a larger part of the blade loading will point in the azimuthal direction.The 19 distributions of ∆C ts versus relative radius r R are shown on Figure 4 for tip-speeds ratios and rotor loadings.The ∆C ts magnitude increases with the disk loading, as one would expect.the annular element force coefficient peaks at around 5% of C t−rotor on the inner part of the rotor for infinite tip-speed ratios.The sustained low pressure of a swirling wake center will cause ∆C ts to decrease.This explains qualitatively the trend towards a decreased-even negative-rotor center ∆C ts values with decreased tip-speed ratio λ.
The ∆C ts trend towards the disk perimeter is the reverse: negative values.The zero-crossing location is around 90% relative radius, and the negative ∆C ts peak at the tip is around −20% of C t−rotor .This rather sharp drop of ∆C ts can be regarded as a disk perimeter downwash effect characterized by a large radial flow component.However, it should not be regarded as a tip-loss effect.Tip-loss in the Prandtl sense is the physical consequence of a rotor having a finite number of blade tips shedding discrete tip vortices.Therefore, an axi-symmetric disk has zero tip-loss.The AD simulations are set up to equal the BEM flow, i.e., no turbulence model, and no viscosity.Therefore, the absence of eddy diffusivity (turbulence model) will cause the AD simulations to diverge upon entering a reversed-flow wake state at high disk loading conditions.The specific AD simulation for λ = 2 with the highest loading is locally at the tip very close to wake flow reversal.This probably explains the different shape of the radial ∆C ts distribution for this computation.1) and tip-speed ratio λ.

Results
Spanwise distributions and rotor/disk integrated quantities will be presented for the baseline AD method, the classical BEM method, and 4 different variations of the corrected BEM.These denotations are used: • AD: The Actuator Disk baseline method against which the BEM methods are benchmarked, see Appendix A for details.
• BEM-E-S12: A hybrid between BEM-E-S1 and BEM-E-S2 where: Wake expansion correction term ∆F s given by Equations ( 14)-( 16).Far wake swirl correction term ∆p w and angular momentum conservation is a weighted average of S1 (60%) and S2 (40%) All results displayed use the color codes for disk loading given in Table 1.Focus is on the normalized axial flow velocity at disk passage, where a is the axial induction at the disk.A u d discrepancy corresponds roughly to an angle-of-attack discrepancy for an actual blade through the simple relations: This leads to: Equation (19) states that the relative error on the axial flow velocity at disk passage equals the relative error of the corresponding angle-of-attack.
The entire set of AD computed radial distributions of the normalized axial flow velocity at disk passage u d = U d /U ∞ is displayed in Figure 5.These velocity distributions constitute benchmark data for subsequent error quantification of the BEM results in Figures 6-10.Equation (19) states that the relative error on the axial flow velocity at disk passage equals the relative error of the corresponding angle-of-attack.
The entire set of AD computed radial distributions of the normalized axial flow velocity at disk passage   =    ∞ ⁄ is displayed in Figure 5.These velocity distributions constitute benchmark data for subsequent error quantification of the BEM results in Figures 6-10.  1) and tip-speed ratio .
The classical BEM is benchmarked against the AD results in Figure 6.The relative error in axial flow velocity at disk passage versus relative radius r/R is displayed for the 19 configurations listed in Table 1.The relative error is defined as diff(  ) = ( _ −  _ )  _ ⁄ .The AD method yields up to 5% higher   than the classical BEM at infinite tip-speed ratio at the inner and mid part of the rotor and down to −30% lower   towards the tip.The 0-error crossing location is at / ≈ 0.83 for all disk loadings.When lowering  to normal levels of 5 and 10 the classical BEM underestimation  1) and tip-speed ratio λ.
The classical BEM is benchmarked against the AD results in Figure 6.The relative error in axial flow velocity at disk passage versus relative radius r/R is displayed for the 19 configurations listed in Table 1.The relative error is defined as diff(u d ) = (u d_AD − u d_BEM )/u d_AD .The AD method yields up to 5% higher u d than the classical BEM at infinite tip-speed ratio at the inner and mid part of the rotor and down to −30% lower u d towards the tip.The 0-error crossing location is at r/R ≈ 0.83 for all disk loadings.When lowering λ to normal levels of 5 and 10 the classical BEM underestimation of u d increases further on the mid part to 6-7% and up to 20-30% on the inner part of the rotor.The 0-error crossing location moves outwards to around 90% r/R.At very low tip-speed ratio of 2 the tendencies continue: Up to +10% classical BEM underestimation of u d on the mid part and a much wider inner part of the rotor where discrepancies increase up to +30%.The 0-error crossing location moves further outboard to around 95%.As a rough indication, Equation ( 19) tells us that 5-10% error on u d corresponds to 5-10% error on the BEM-predicted angle-of-attack (AoA).This would for a typical blade operating at 7-8 degrees AoA correspond to around 0.5 degrees systematic error on AoA.The AoA errors on the inboard part of the rotor is several times higher, but the swept area is limited here, so the impact on power is also very limited, at least for normal λ values in the range 5-10.The AoA errors on the tip section also peaks at roughly 2-3 degrees, being overestimated by the classical BEM.The inclusion of tip-loss effects on a real 3-bladed rotor causes unloading of the blade towards the tip, so the significant AoA error at the tip is therefore less influential on the power for a real wind turbine.
Energies 2018, 11, x FOR PEER REVIEW 13 of 28 towards the tip, so the significant AoA error at the tip is therefore less influential on the power for a real wind turbine.The BEM-E benchmark results are displayed in Figure 7.The error levels have diminished significantly.BEM-E contains the axial momentum correction, and angular momentum conservation, but not the far-wake swirl correction pressure term.The best agreement is therefore seen in the swirlfree configurations for  = ∞, where diff(  ) is within a fraction of a percent on the inboard and mid part of the rotor.Towards the tip some wiggles are seen, and diff(  ) is bounded within +/− a few percent.Significant error reduction is also observed the finite  cases, where at least part of the discrepancy must be attributed to the lack of the far-wake swirl correction pressure term in the BEM-E method.
The BEM-E-S1 benchmark results are displayed in Figure 8.This BEM method contains both the axial momentum correction and the angular momentum correction in its original inviscid form (S1). Error levels are unchanged in the swirl-free case of  = ∞, but have diminished further for the finite  cases, especially for  = 10.For the cases  = 5 and  = 2 the improved agreement with the AD reference is confined to the mid and outer part of the rotor, but the inboard (center) part of the rotor The BEM-E benchmark results are displayed in Figure 7.The error levels have diminished significantly.BEM-E contains the axial momentum correction, and angular momentum conservation, but not the far-wake swirl correction pressure term.The best agreement is therefore seen in the swirl-free configurations for λ = ∞, where diff(u d ) is within a fraction of a percent on the inboard and mid part of the rotor.Towards the tip some wiggles are seen, and diff(u d ) is bounded within +/− a few percent.Significant error reduction is also observed the finite λ cases, where at least part of the discrepancy must be attributed to the lack of the far-wake swirl correction pressure term in the BEM-E method.The center-wake diffusion issue is not only an AD CFD model numerical effect.In real life any strong vortex will decay over time.Vortex decay is the result of laminar and turbulent (eddy) The BEM-E-S1 benchmark results are displayed in Figure 8.This BEM method contains both the axial momentum correction and the angular momentum correction in its original inviscid form (S1). Error levels are unchanged in the swirl-free case of λ = ∞, but have diminished further for the finite λ cases, especially for λ = 10.For the cases λ = 5 and λ = 2 the improved agreement with the AD reference is confined to the mid and outer part of the rotor, but the inboard (center) part of the rotor still has sizeable error levels.Note on this inboard part that the BEM and BEM-E methods from Figures 5 and 6 underpredict the axial velocity u d , whereas the BEM-E-S1 method overpredicts by roughly the same magnitude.The constant spanwise loading from tip at 100% R to root at 0.5% R (see Appendix A for details) causes formation of a strong center vortex characterized by large azimuthal induced velocities and strong radial gradients for pressure and velocity along the entire infinite length of the wake.This center-wake vortex is perfectly sustained in the BEM-E-S1 analytical formulation but will in the AD computations loose strength as it propagates downstream due to numerical diffusion of radial pressure and velocity gradients.The BEM-E-S1 method's overprediction of inboard axial velocities is attributed to this.
The center-wake diffusion issue is not only an AD CFD model numerical effect.In real life any strong vortex will decay over time.Vortex decay is the result of laminar and turbulent (eddy) viscosity.It is therefore relevant to investigate if the ideal inviscid center-vortex which approaches singularity towards the center-axis can be somehow relaxed, such that the act of viscosity and vortex decay is mimicked.With this objective in mind the BEM-E-S2 method was formulated.The underlying assumption is that angular wake momentum is still conserved downstream of the disk, but instead of imposing conservation for each annular element, the exchange of azimuthal stress between adjacent annuli is included.Angular wake momentum is then conserved in a lumped sense for all the annuli.In the far-wake limit the swirl will therefore be of solid-body type and no longer predominantly irrotational as before.This is modeled by replacing Equation ( 5) by Equations ( 12) and ( 13) in the set of governing BEM equations, as earlier stated.The BEM-E-S2 benchmark results are displayed in Figure 9.The mid and outboard parts of the rotor are only marginally impacted by the changed formulation of the far-wake swirl.On the inboard (center) part of the rotor, the BEM-E-S2 method underpredicts the axial velocity u d compared to the AD reference for all finite values of tip-speed ratio λ.The wake swirl solid-body assumption causes transfer of angular momentum from the wake center annuli to the outer annuli, so the swirl-induced wake suction of the BEM-E-S2 method is weaker.This explains the u d underprediction at the rotor center.The center-wake diffusion issue is not only an AD CFD model numerical effect.In real life any strong vortex will decay over time.Vortex decay is the result of laminar and turbulent (eddy) viscosity.It is therefore relevant to investigate if the ideal inviscid center-vortex which approaches singularity towards the center-axis can be somehow relaxed, such that the act of viscosity and vortex decay is mimicked.With this objective in mind the BEM-E-S2 method was formulated.The underlying assumption is that angular wake momentum is still conserved downstream of the disk, but instead of imposing conservation for each annular element, the exchange of azimuthal stress between adjacent annuli is included.Angular wake momentum is then conserved in a lumped sense for all the annuli.In the far-wake limit the swirl will therefore be of solid-body type and no longer predominantly irrotational as before.This is modeled by replacing Equation ( 5) by Equations ( 12) and ( 13) in the set of governing BEM equations, as earlier stated.The BEM-E-S2 benchmark results are displayed in Figure 9.The mid and outboard parts of the rotor are only marginally impacted by the changed formulation of the far-wake swirl.On the inboard (center) part of the rotor, the BEM-E-S2 method underpredicts the axial velocity   compared to the AD reference for all finite values of tip-speed ratio .The wake swirl solid-body assumption causes transfer of angular momentum from the wake center annuli to the outer annuli, so the swirl-induced wake suction of the BEM-E-S2 method is weaker.This explains the   underprediction at the rotor center.Therefore, the BEM-E-S1 method with ideal inviscid far wake swirl, and the BEM-E-S2 method with smeared out viscous far wake swirl, respectively overpredicts and underpredicts the rotor center flow-rate.A logical final approach is therefore to combine S1 and S2 such that the rotor center flowrate is neither overpredicted nor underpredicted.In the proposed BEM-E-S12 hybrid method Equation ( 5) is only partially replaced by Equations ( 12) and ( 13), meaning that the stress-induced Therefore, the BEM-E-S1 method with ideal inviscid far wake swirl, and the BEM-E-S2 method with smeared out viscous far wake swirl, respectively overpredicts and underpredicts the rotor center flow-rate.A logical final approach is therefore to combine S1 and S2 such that the rotor center flow-rate is neither overpredicted nor underpredicted.In the proposed BEM-E-S12 hybrid method Equation ( 5) is only partially replaced by Equations ( 12) and ( 13), meaning that the stress-induced cancellation of azimuthal shear in the far-wake leading to solid-body rotation is only partial.The expression for ∆p w (r) in Equation ( 11) is then a weighted average: U θw1 (r) is computed as in BEM-E-S1 from Equation ( 5), and U θw2 (r) is computed as in BEM-E-S2 from ( 12) and (13).w is a weight factor.The BEM-E-S12 benchmark results are displayed in Figure 10.The weight factor is w = 0.6.This value leads to the best obtainable agreement between the AD reference and the new BEM method with wake expansion and wake swirl corrections.Compared to the classical BEM discrepancies in Figure 6, error reductions of roughly one order of magnitude are obtained across all investigated values for the tip-speed ratio λ.Apart from the weight factor w in BEM-E-S12 no heuristic fitting is involved in the improved accuracy obtained from the axial and angular momentum correction BEM models.A physical alternative to the BEM-E-S12 hybrid method could be to introduce an azimuthal shear stress term to the angular momentum conservation Equation (5).A careful formulation of the shear stress term would be needed to minimize the created interdependency between adjacent annular elements.Recent investigation by van Kuik [24] on circumventing an apparent power paradox for very small tip-speed ratios is relevant here.The infinite power paradox was resolved by introducing a small amount of viscosity to the core part of the wake center vortex.
Energies 2018, 11, x FOR PEER REVIEW 16 of 28 1 () is computed as in BEM-E-S1 from Equation ( 5), and  2 () is computed as in BEM-E-S2 from Equations ( 12) and (13). is a weight factor.The BEM-E-S12 benchmark results are displayed in Figure 10.The weight factor is  = 0.6 .This value leads to the best obtainable agreement between the AD reference and the new BEM method with wake expansion and wake swirl corrections.Compared to the classical BEM discrepancies in Figure 6, error reductions of roughly one order of magnitude are obtained across all investigated values for the tip-speed ratio .Apart from the weight factor  in BEM-E-S12 no heuristic fitting is involved in the improved accuracy obtained from the axial and angular momentum correction BEM models.A physical alternative to the BEM-E-S12 hybrid method could be to introduce an azimuthal shear stress term to the angular momentum conservation Equation (5).A careful formulation of the shear stress term would be needed to minimize the created interdependency between adjacent annular elements.Recent investigation by van Kuik [24] on circumventing an apparent power paradox for very small tip-speed ratios is relevant here.The infinite power paradox was resolved by introducing a small amount of viscosity to the core part of the wake center vortex.Matching of the radial distribution of axial flow-rate through the rotor is achieved, so a similar match of the radial distribution of produced power should be expected.The BEM and AD power coefficient distributions are displayed in Figure 11.The constant load levels from ID # 3, 7, 11, and 17 (see Table 1) respectively for the tip-speed ratios ∞, 10, 5, and 2 are used, which correspond to nearly power-optimal conditions.The power coefficient is defined as standard   () =  ( 12   () ∞ 3 ) ⁄ , and is evaluated for each discretized annular element along the rotor radial dimension.The specific Matching of the radial distribution of axial flow-rate through the rotor is achieved, so a similar match of the radial distribution of produced power should be expected.The BEM and AD power coefficient distributions are displayed in Figure 11.The constant load levels from ID # 3, 7, 11, and 17 (see Table 1) respectively for the tip-speed ratios ∞, 10, 5, and 2 are used, which correspond to nearly power-optimal conditions.The power coefficient is defined as standard C p (r) = P/ 1  2 A d (r)ρU 3 ∞ , and is evaluated for each discretized annular element along the rotor radial dimension.The specific expression for C p (r) for the AD reference is given in Appendix A. The specific expression for C p (r) for the BEM methods is: Figure 11 shows coinciding C p -distributions at λ values ∞, 10, and 5 for the AD reference and all the BEM methods containing angular momentum conservation and wake expansion correction, i.e., BEM-E, BEM-E-S1, BEM-E-S2, and BEM-E-S12.Especially the reference match of the BEM-E method is interesting, since it does not contain the far wake swirl-induced pressure drop ∆p w .This indicates that the angular momentum conservation, Equation ( 5), plays a more important role than ∆p w in computing the flowfield correctly.In the challenging case of λ = 2 the small differences in C p between the corrected BEM methods become visible, and it is seen that inclusion of the far wake pressure drop leads to better agreement with the AD reference for the C p -distribution.Also comment on the obvious result that the classical BEM method exhibits a clear discrepancy from all the other methods across the entire radial range and for all tip-speed ratios.

Impact on Modern Wind Turbine Blades
The corrected BEM models remove up to 5% error in flowrate through the rotor computed by the classical BEM method at normal operation thrust levels, and even more at the root and tip regions.5% flowrate error corresponds to approximately 0.5 degrees error in predicted angle-of-attack on the blade.The performance impact from a systematic 0.5 degrees AoA error is case dependent.The power and load consequence will depend on the designed stall margin of the blade layout, as well as the lift overshoot corresponding to the stall margin along the span of the blade.Two errors occur from not applying the corrected models.Firstly, the classical BEM's 3-5% underprediction of flowrate through the inner and mid part of the rotor at typical thrust levels leads to 0.3-0.5 degrees underprediction of the local AoA.This AoA misprediction will be well within the typical stall margin, but during gusts and periods with profound atmospheric shear, the onset of local blade stall will occur slightly sooner than expected.Secondly, the classical BEM's 10-20% and even higher overprediction of flowrate through the outermost and tip part of the rotor at typical thrust levels leads to several degrees overprediction of the local AoA.This will cause a lower power takeout than expected.On the other hand, the tip-loss effect, that reduces the power takeout potential towards the tip, will alleviate the power drop consequence from overpredicting AoA.More studies are needed to quantify the resulting impact in the tip region.Overall, a C p−rotor increase of about 0.005 is gained at a tip-speed ratio of 10 as seen in Table 2.This corresponds to roughly 1% more predicted power.Note, however, that the eventual extra power obtainable on site from a blade optimized with the axial and angular momentum correction models applied, has not been quantified, and should not be confused with the 1% extra predicted power.A specific study is needed for this, a study that addresses impact on both power and loads.The obtained results in the present investigation with increased power potential on the inboard and mid parts in expense off decreased power in the tip region, show that the spanwise bulk of power production moves slightly inboard.This could be beneficial for blade loading, both FLS and ULS.

Summary and Perspectives
Two BEM correction models for axial and angular momentum conservation have been derived and quantified.The corrections contain three parts:

Axial momentum
Axial momentum correction through annular element forcing coefficient, ∆C ts .Far-wake swirl-induced pressure drop, ∆p w .

• Angular momentum
Angular momentum conservation enforcement, Equation ( 5) These corrections contain no empiricism.In an alternative version (S2) of the far-wake swirl-induced pressure drop, solid-body far-wake rotation was assumed, and slightly improved agreement with the actuator disk reference results was achieved in the hybrid wake swirl correction model (S12), where the far-wake pressure drop is computed from a weighted average of the inviscid swirl S1 model, and the viscosity impacted solid-body swirl S2 model.
All three correction parts mentioned above imply that the annular independency assumption must be abandoned.This is the primary reason for the classical BEM not to include these corrections.Classical BEM goes back to before the advent of modern computation power.The introduced annular element dependency is still rather weak, and the proposed BEM methods (BEM-E, BEM-E-S1, BEM-E-S2, BEM-E-S12) still converge easily, and still compute extremely fast compared to CFD-based alternatives.The convergence of BEM-E-S1 method with fully inviscid swirl cannot maintain full blade loading all the way to the rotor center due to the strong center vorticity, and blade loading must therefore be shed at a radial position of a few percent relative radius.This poses no conflict with modern blades where the innermost part of the blade is a regular cylinder carrying no bound vorticity.The BEM-E and BEM-E-S2 methods do not have this issue.
The introduction of annular element dependency makes extra iterations necessary for static BEM analysis.For transient BEM analysis, as used in numerous aero-elastic WTG simulators such as FLEX, HAWC, Bladed, etc., each BEM iteration can be performed at a new timestep.This approach is already adopted for the iteration of the transient wake response from dynamic inflow conditions.Typical timestep size is 0.01-0.02s which is 1-3 orders of magnitude less than the time constants associated with the blade and wake aerodynamics, so iterating the BEM during time-stepping instead of introducing new BEM iterations at each time-step would be physically sound and align well with existing methodology.
The two correction models for axial and angular momentum are "nice to have" modeling features, but not "need to have".The present advanced and mature wind turbines demonstrate this by their mere existence.Yet, the corrections circumvent systematic errors in the classical BEM and can be implemented without any adverse effect or incompatibility issue with the existing BEM or the various existing BEM add-on corrections.Agreement is now obtained between actuator disk results and the corrected BEM in a physically correct non-empirical way.The possible WTG end product value from improved modeling accuracy is to be revealed through specific blade design studies.
In the present investigation annular element interexchange of pressure force axial momentum has been computed and mapped through axi-symmetric AD CFD simulations to quantify directly the final missing term in the axial momentum conservation equation (along with the far-wake pressure drop).In a broader view, this methodology of quantifying the non-analytic pressure gradient induced annular element forcing part of the BEM through brute-force CFD is generic and can be applied to bridge other remaining BEM inaccuracies.In the example of a rotor not strictly confined to a plane, e.g., due to coning, pre-deflection, and load-driven deflection, the forcing from the blade on the passing stream will also have a radial projection, but otherwise, the governing BEM Equations ( 1)-(9) will remain the same, and accurate quantification of the aerodynamic blockage impact from coning etc. can be found from mapping ∆C ts for, say, a 1-parameter family of deflection levels.No such correction for out-of-plane blades currently exist.The extent to which a wind turbine rotor is out-of-plane increases with blade length and flexibility.A blade out-of-plane BEM correction is therefore increasingly relevant as the wind turbines continue to grow in size.

Conclusions
Inclusion of the missing momentum terms and strict enforcement of momentum conservation in the momentum theory contained in the BEM model has revealed that:

•
The corrected BEM model matches AD computations very well.

•
The standard BEM model underpredicts the local power and flowrate through the inner part of the disk by 10-20% at normal operation thrust levels.

•
The standard BEM model underpredicts the local power and flowrate through the mid part of the disk by up to 5% at normal operation thrust levels.

•
The standard BEM model overpredicts the local power and flowrate through the tip of the disk by 10-20% at normal operation thrust levels.

•
The power underprediction and overprediction by the standard BEM on the inner-mid and outer part of the blade respectively cancel out for infinite tip-speed ratios.A small overall power gain of 1% is observed for the corrected BEM methods for λ = 10, increasing gradually to 8% for λ = 2.
Specific blade studies are needed to quantify the net power/load impact of the corrected BEM in the blade tip region in combination with a tip loss model, as well as overall performance impacts on modern wind turbine rotors with high tip speed ratios around 10.
As a final remark re-state that the improved modelling accuracy, resulting from strictly enforcing axial and angular momentum conservation, concerns the axi-symmetric static BEM part of the complete BEM model (see Figure 1).Other sources of error from other parts of the BEM model, e.g., the several add-ons, might yield similar or even larger errors than those corrected for in this paper.As such, these momentum conservation corrections are just one step towards further improved BEM. favored because of higher execution speed.The lack of higher-order precision was compensated for by employing high mesh resolution.Table A1 shows the change in rotor power efficiency C p−rotor when higher order discretization schemes are used on the same mesh (configurations 2-3) and when combining higher order discretization with increased mesh resolution (configurations 4).C p−rotor decreases marginally with up to 0.1% when increasing discretization order and mesh resolution.Satisfactory discretization independency is shown within 3 significant digits.The tendency of slightly decreasing C p−rotor is due to improved capture of the slip discontinuity between wake and surrounding flow with improved discretization accuracy and resolution.The inevitable introduction of even the slightest numerical diffusion along the slip-line will cause a slight transfer of energy from the high-momentum surrounding flow to the low-momentum wake, thereby artificially increasing the potential for power extraction, although marginally.This is also evident from Figure A1, where the radial distribution of the discretization-induced change in C p reaches a maximum towards the disk tip, indicating that the tip area is particularly sensitive to discretization and numerical diffusion.The simulated domain for all AD computation is shown in Figure A2.The actuator disk is adjacent to the symmetry axis at = 0 and extents perpendicularly to the tip at r = R.The inflow Dirichlet condition on axial velocity is located 50R upstream, and the outflow Neumann condition on pressure is located 50R downstream.The radial farfield condition on pressure is also of Neumann type and is on average located 50R away from the symmetry axis.The disk thickness is 0.02R.
The streamlines emanating upstream and downstream from the disk in the step 1 solution are used to define the annuli subdomains, see Figure A3.The step 3 mesh contains the streamlines as internal mesh lines and is therefore slightly different for each of the 19 configurations.A mesh example is shown in Figure A4.All annuli borderlines (streamlines) continue 10-20 R upstream and downstream of the disk.The step 3 meshes are fully unstructured and contain 500,000-600,000 elements.A step 3 solution is shown in Figure A5.Notable flow features include the pressure drop across the disk, the wake expansion, the swirl-induced innerwake speed-up and associated sustained low pressure.These axisymmetric 2D domain AD computations converge after approximately 150 pseudo-timesteps with residual reductions of minimum 4-5 orders of magnitude for high-load configurations and machine precision convergence, 10-12 orders of magnitude, for the lower load cases.The step 3 mesh contains the streamlines as internal mesh lines and is therefore slightly different for each of the 19 configurations.A mesh example is shown in Figure A4.All annuli borderlines (streamlines) continue 10-20 R upstream and downstream of the disk.The step 3 meshes are fully unstructured and contain 500,000-600,000 elements.A step 3 solution is shown in Figure A5.Notable flow features include the pressure drop across the disk, the wake expansion, the swirl-induced inner-wake speed-up and associated sustained low pressure.These axisymmetric 2D domain AD computations converge after approximately 150 pseudo-timesteps with residual reductions of minimum 4-5 orders of magnitude for high-load configurations and machine precision convergence, 10-12 orders of magnitude, for the lower load cases.In the final post-processing step 4 the ∆ distribution along the radial extension of the disk is computed for each of the 19 AD configurations, see Equation (17).
Other post-processing steps include computation of the thrust and power coefficients, see [19] for derivation details.  is the disk subdomain.

Figure 1 .
Figure 1.Constituents of the BEM model.The shaded box indicates the part of the BEM model which is modelled and corrected in this work.The four red boxes indicate the chain of dependency for the momentum conservation corrections.

Figure 1 .
Figure 1.Constituents of the BEM model.The shaded box indicates the part of the BEM model which is modelled and corrected in this work.The four red boxes indicate the chain of dependency for the momentum conservation corrections.

Figure 2 .
Figure 2. (Left): Cross-section schematic of domain, flow quantities, and streamlines.The flow enters the control volume at the far upstream inlet, passes the disk, and exits the control volume at the far downstream outlet.Four annular elements bounded by five stream tubes (incl.the center-axis) are shown.(Right): Schematic of forces from rotor blade on flow.

Figure 2 .
Figure 2. (Left): Cross-section schematic of domain, flow quantities, and streamlines.The flow enters the control volume at the far upstream inlet, passes the disk, and exits the control volume at the far downstream outlet.Four annular elements bounded by five stream tubes (incl.the center-axis) are shown.(Right): Schematic of forces from rotor blade on flow.

28 Figure 3 .
Figure 3. (Left): An annular element bounded by four surfaces: Inlet, outlet, outer surface _ , inner surface _ .(Right): Cross-section axi-symmetric CFD domain part surrounding the disk.Black streamlines show the annular element radial expansion (one is highlighted black), colored lines are isobar contours.Note that due to the wake (streamline) expansion and the pressure gradients in the disk's vicinity, the pressure integrals in Equation (16) will be non-zero.

Figure 3 .
Figure 3. (Left): An annular element bounded by four surfaces: Inlet, outlet, outer surface _out, inner surface s_in.(Right): Cross-section axi-symmetric CFD domain part surrounding the disk.Black streamlines show the annular element radial expansion (one is highlighted black), colored lines are isobar contours.Note that due to the wake (streamline) expansion and the pressure gradients in the disk's vicinity, the pressure integrals in Equation (16) will be non-zero.

Figure 5 .
Figure5.AD results: Radial distribution of normalized axial flow velocity through the rotor,   , for different levels of rotor loading (color code, see Table1) and tip-speed ratio .

Figure 5 .
Figure 5. AD results: Radial distribution of normalized axial flow velocity through the rotor, u d , for different levels of rotor loading (color code, see Table1) and tip-speed ratio λ.

Figure 6 .
Figure 6.BEM results compared to AD.

Figure 6 .
Figure 6.BEM results compared to AD.

Figure 7 .
Figure 7. BEM-E results compared to AD.

Figure 7 .
Figure 7. BEM-E results compared to AD.

Figure 7 .
Figure 7. BEM-E results compared to AD.

Energies 2018 ,
11, x FOR PEER REVIEW 18 of 28 underprediction is balanced out by the outer (tip) part overprediction, and it means that the Betz value remains unaffected.

Figure 11 .
Figure 11.Power coefficient comparisons at optimal loading for the AD and BEM methods.Figure 11.Power coefficient comparisons at optimal loading for the AD and BEM methods.

Figure 11 .
Figure 11.Power coefficient comparisons at optimal loading for the AD and BEM methods.Figure 11.Power coefficient comparisons at optimal loading for the AD and BEM methods.

Figure A1 .
Figure A1.Relative difference in C p for different discretizations.

Figure A3 .
Figure A3.Step 1 solution streamlines defining a set of annular elements.

Figure A3 .
Figure A3.Step 1 solution streamlines defining a set of annular elements.

Table 1 .
Overview of AD computations.

Table 2 .
C p−rotor comparisons for the BEM methods.