the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A computationally efficient engineering aerodynamic model for nonplanar wind turbine rotors
Mac Gaunaa
Georg Raimund Pirrung
Sergio González Horcas
In the present work, a computationally efficient engineering model for the aerodynamic load calculation of nonplanar wind turbine rotors is proposed. The method is based on the vortex cylinder model and can be used in two ways: either used as a correction to the currently widely used blade element momentum (BEM) method or used as the main model, replacing the BEM method in the engineering modeling complex. The proposed method needs the same order of computational effort as the ordinary BEM method, which makes it ideal for timedomain aeroservoelastic simulations. The results from the proposed method are compared with results from two higherfidelity aerodynamic models: a liftingline method and a Navier–Stokes solver. For planar rotors, the aerodynamic loads are identical to the current BEM model when the drag force is excluded during the calculation of the induced velocities. For nonplanar rotors, the influence of the blade outofplane shape, measured by the difference of the load between the nonplanar rotor and the planar rotor, is in very good agreement with higherfidelity models. Meanwhile, the existing BEM methods, even with a correction of radial induction included, show relatively large deviations from the higherfidelity method results.
Please read the corrigendum first before continuing.

Notice on corrigendum
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.

Article
(6745 KB)

The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(6745 KB)  Fulltext XML
 Corrigendum
 BibTeX
 EndNote
The blade element momentum (BEM) method has long been dominant in the lowfidelity aerodynamic modeling of horizontalaxis wind turbines. Until now, it is the main working horse for wind turbine aeroservoelastic simulations and is widely used in the wind turbine design and optimization framework. There are many explicit and implicit assumptions in the BEM method. The BEM method explicitly assumes that uniform inflow is applied to the rotor that is operating at a high tipspeed ratio and the stream tubes are independent of each other. The model also implicitly assumes a planar rotor with straight blades and using quasisteady aerodynamics. There has been extensive work on the modifications and corrections to the BEM method, such as dynamic stall model (Leishman and Nguyen, 1990; Hansen et al., 2004; Larsen et al., 2007), dynamic inflow model (Schepers and Snel, 1995; Yu et al., 2019), polargrid based unsteady BEM (Madsen et al., 2020a), modeling of turbulent inflow (Mann, 1994), highthrust correction (Spera, 1994; Madsen et al., 2020a; Burton et al., 2021) and corrections for operation in yawed conditions (Leishman, 2005).
The results from the BEM method generally show surprisingly good agreement with higherfidelity models, at least on the integral level. However, due to the progress of wind turbine technology, modern multimegawatt designs are generally more flexible than the stiff machines of the 1980s. It implies that modern wind turbine blades typically have more prebend, larger cone angle and larger deformations. The influence of blade outofplane shapes on the aerodynamics is then more pronounced and can not simply be neglected. Some new developments have even more pronounced outofplane shapes. For example, a downwind wind turbine designed for lowwind conditions could have large cone and prebend and possibly dramatic outofplane deformations (Madsen et al., 2020b). In addition, some wind turbines are equipped with winglets to reduce the drag force and also the noise. Modern wind turbines are generally designed using mainly the BEMbased codes. Higherfidelity tools such as liftingline method (LL) or fully resolved Navier–Stokes solvers (often referred to as computational fluid dynamics, CFD) are mostly used for comparison or for very specific load cases due to the high computational effort. However, when the blades have large outofplane shapes due to prebend, deformation or cone, the results from these BEM codes will have relatively large differences compared to the results from higherfidelity tools (Madsen and Rasmussen, 1999). This is because the influence of blade outofplane shapes is not correctly captured by these BEMbased codes. As a result, a design from an optimization tool using a BEMbased aerodynamics module could be far from the actual optimal solution. There may even exist aeroelastic instabilities that are not correctly captured by the BEMbased tools.
On the other hand, rotorresolved CFD and liftingline method are computationally too expensive for extensive use in current and nearfuture design optimization processes. Navier–Stokes solvers with fully resolved rotor geometry have no difficulties predicting nonplanar rotor effects. But for the liftingline method, care should be taken on the influence of curved bound vortex on itself (Li et al., 2020), the correct directions of applying the lift and drag force, and also the possible noncirculatory lift during the aerodynamic load calculation, to correctly predict the effects. Therefore, a lowfidelity model that could capture the most important features of the aerodynamics of nonplanar rotors, while maintaining approximately the same level of computational effort as the current BEM methods, would be of great value to both the scientific and the commercial wind turbine communities. Both for the design optimization and for the aeroelastic simulations.
In order to correctly account for the outofplane shapes of the wind turbine blades in the lowfidelity model, the physics behind the problem should be analyzed, and then the most important aspects should be proactively modeled while less important features can be neglected. In the present work, the force on the nonplanar rotor is firstly analyzed in a physically consistent way using the Kutta–Joukowski theorem. The conclusion from the analysis is that the streamwiseshifted starting position of the trailed vorticity, due to the nonplanar bound vortex surface swept by the blades, will influence both axial and radial induction, and both have direct influences on the aerodynamic loads. Therefore, we consider the vortex cylinder model (Branlard and Gaunaa, 2015a) has the potential to capture these most important features in the aerodynamics of the nonplanar rotors. There has been previous work by Crawford (2006) using a vortex cylinder model for the aerodynamic calculation of coned rotors. In that work, the same idea of using the axial and radial induction from the vortex cylinder model is proposed, and the closed equations for the model are given in the framework of momentum theory. Some comparisons with actuator disc results for uniformly loaded cases are shown for the planar rotor. And for the coned rotors, the axial induction and aerodynamic loads are compared with the results from actuator disc simulations performed by Mikkelsen (2004). The good agreement shows the potential of the vortex cylinder model for the nonplanar rotor. However, the work is limited to the momentum theory framework while the equivalence of the vortex cylinder model and the momentum theory for planar rotors is not highlighted. Furthermore, other important features of the vortex cylinder model for nonplanar rotors, such as the similarity to the planar rotors or the impact of unsteady airfoil aerodynamics on the steadystate results, are not described. Nevertheless, the pioneering work by Crawford (2006) inspires the authors and is a good starting point for the current study, which also builds on previous efforts on superposition of vortex cylinders (Branlard and Gaunaa, 2015a).
In the present work, a detailed analysis of the vortex cylinder model for nonplanar rotors will be performed. A method based on the vortex cylinder model for the aerodynamic load calculation of such nonplanar rotors is then proposed. The description of the implementation of the proposed method is in the framework of the HAWC2 code (Larsen and Hansen, 2007). Some details of the implementation may be different compared to other BEMbased aeroelastic codes. In the present work, only the outofplane shape of the blade is considered, which means the blade is assumed to have no inplane sweep. The engineering aerodynamic model for the blades with only inplane shapes is described by Li et al. (2021). The structure of the work is as follows: the Kutta–Joukowski analysis of the nonplanar rotor is performed in Sect. 2. Then, the vortex cylinder model and its relationship with the momentum theory are briefly introduced in Sect. 3. Some important aspects for the implementation of the BEM method and the proposed method for nonplanar rotors are described in Sect. 4. The coupling of the blade element theory with the vortex cylinder model, including details on the tiploss correction and a summary of the algorithm, is described in Sect. 5. The lowfidelity and higherfidelity aerodynamic models for comparisons are described in Sect. 6. The setup of the test cases is described, and the results from the lowfidelity models are compared with the results from higherfidelity models in Sect. 7. Finally, the conclusions are drawn and the future work is summarized in Sect. 8.
For the planar rotor with straight blades, the Kutta–Joukowski analysis was previously used to derive the similarity between the superposition of the vortex cylinders and the BEM method by Branlard and Gaunaa (2015a). In this section, the influence of the blade outofplane shapes on the aerodynamics is investigated using the Kutta–Joukowski theorem (Okulov et al., 2015). The blade is assumed only possible to have outofplane shapes (dihedral or cone) but no inplane shapes (blade sweep) in the following analysis.
The coordinate system is defined as follows and is illustrated in Fig. 1. The x axis is the axial direction, which is positive in the incoming wind direction. The rotation vector of the rotor is in the positive x direction. The y axis is the radial direction, which is positive in the direction of increasing radius of the blade. The z axis is normal to both the x axis and the y axis, and its direction is defined so that a righthanded system is found. The z axis is defined as the tangential direction. The airfoils are aligned perpendicular to the main axis of the halfchord line.
The local dihedral angle is defined to be positive when the blade is tilting upwind and can be calculated using the blade mainaxis geometry.
The analysis is applied to a nonplanar rotor with given blade bound circulation and induced velocity at each blade section. For section i of a blade with radius of r_{i}, the bound circulation strength is ${\mathrm{\Gamma}}_{i}^{\mathrm{B}}$ and the local dihedral angle is κ_{i}. The axial, tangential and radial induced velocity are u_{a,i}, u_{t,i} and u_{r,i}.
The relative velocity experienced by the blade section is V_{rel,i}. The bound circulation of a blade section ${\mathbf{\Gamma}}_{i}^{\mathrm{B}}$ is tangent to the local blade section.
With the Kutta–Joukowski theorem in threedimensional vector form, the lift force on the blade is obtained.
In Eq. (3), the force f_{i} is corresponding to lifting force per unit curved blade length. The lifting force per unit radius, corresponding to what is used in momentum theory analysis, is ${\mathit{f}}_{i}^{*}={\mathit{f}}_{i}\frac{\mathrm{d}s}{\mathrm{d}r}$, as also shown by Madsen et al. (2020a). The $\frac{\mathrm{d}s}{\mathrm{d}r}$ term is representing the ratio of the local change of curved blade length and local change of radius. For blades with only outofplane shapes, it is equal to $\mathrm{1}/\mathrm{cos}{\mathit{\kappa}}_{i}$.
The force ${\mathit{f}}_{i}^{*}$ is divided into two parts: a part with and a part without the direct contribution of the local dihedral angle κ_{i}.
For a nonplanar rotor with upwind direction dihedral (κ_{i}>0), such as prebend or upwind cone, some conclusions can be obtained according to Eq. (5). Comparing to the corresponding planar rotor (as shown later on in Fig. 4), the nonplanar rotor will have outboard radial force. Furthermore, there will be additional tangential driving force due to the radial induction, which for the wind turbine case is positive. So, to get the correct tangential load distribution and consequently the aerodynamic power, it is not only necessary to correctly model the influence of the nonplanarity of the rotor on the axial and tangential induced velocity, but also on the radial induced velocity. However, the radial induced velocity is not available from the momentum theory.
There are two tracks to modify the BEM method to model the nonplanar effects. The first track is based on the previous work on the development of the radial induction correction for the application in the BEM method, derived based on an analytical 2D actuator disc/strip model combined with an engineering fit of the numerical actuator disc simulations (Madsen, 1997; Madsen et al., 2010):
where C_{T,av}(r) is the averaged thrust coefficient as function of the radial position and is defined as
However, the radial induction and the axial induction correspond to planar rotors.
Apart from the momentum theory, which effectively only applies to a planar rotor disc, it is possible to calculate the induced velocity, including the radial component, from analytical equations at each blade section with the superposition of vortex cylinders. So, the second track is based on the vortex cylinder model where the assumption of the planar rotor in the previous work by Branlard and Gaunaa (2015a) is relaxed. This approach inherently includes the effect of axial displacement of the cylindrical wake of the nonplanar rotor in a physically consistent manner. The model even has the potential to completely replace the momentum theory in the BEM method and will be described in the following sections. In the present work, both methods will be tested numerically in Sect. 7.
The vortex cylinder model is a simplified representation of the vortex system of a horizontalaxis wind turbine rotor. The model consists of superposition of bound vortex discs, nonexpanding vortex cylinders with both tangential and longitudinal vorticity, and root vortices (Branlard and Gaunaa, 2015a). An illustration of different components in the vortex cylinder model is shown in Fig. 2.
This model can be considered as the special case of the Joukowski rotor model for the limiting case of the number of blades tending to infinity. It has been shown by Branlard and Gaunaa (2015a) that for a planar rotor, the induced velocities from the version of the vortex cylinder model where the wake rotation effect is neglected are identical to the induced velocities obtained from the momentum theory when the drag force component is not included in the force balancing from which the induced velocities are calculated. It was also argued by Branlard (2017) that the correct way of implementing the blade element momentum (BEM) method should exclude the drag during the calculation of the induced velocities and include the drag in the aerodynamic load calculation afterwards.
One advantage of the vortex cylinder model over the momentum theory is that the induced velocity at any arbitrary point in the flow field is known. In contrast, from the momentum theory, only the axial and tangential velocity at the rotor disc and at infinitely far upstream and downstream of the rotor plane are known. This advantage has been used in the application of the vortex cylinder model in the calculation of the induction zone of a wind turbine (Branlard and Meyer Forsting, 2015) and the wind farm blockage effects (Branlard and Meyer Forsting, 2020). Another advantage of the vortex cylinder model over the momentum theory is that it does not require the assumptions of a planar rotor and the flow with constant speed being perpendicular to it. Other applications of the vortex cylinder model include modeling of wind turbines in yaw (Branlard and Gaunaa, 2016) and modeling of the dynamic inflow effects (Yu et al., 2019). The results from all the aforementioned applications compare well with higherfidelity tools, indicating that the main mechanisms are captured using this framework.
When applying the vortex cylinder method to a nonplanar rotor, the starting position of the cylindrical vortex sheets follows the curved bound vortex surface and will be displaced upstream or downstream compared to the case of a planar rotor. The induced velocity on the nonplanar rotor surface will therefore be different from the induced velocity of a planar rotor. This effect can be modeled by the superposition of the vortex cylinders according to the curved bound vortex surface swept by the blades that have outofplane shapes. However, the possibility of using the vortex cylinder model for the nonplanar rotor is not well recognized and is thus not widely utilized. The work of Crawford (2006) is on this topic, but the system closure for the nonplanar rotor is different compared to the current study. In the current implementation, the system closure is determined in the far wake, thus assuming the same method of system closure as for the nonplanar rotors. Then, the equations of the inductions of the nonplanar rotor are in concise forms following this assumption. With the current system closure, which is an important assumption in this work, there are clear physical connections between the vortex cylinder model of a nonplanar rotor and a planar rotor and subsequently the connection to the momentum theory. In the following content of this work, zero yaw error, no rotor tilt and uniform inflow are assumed. The vortex cylinder is then a right cylinder as opposed to an oblique cylinder used in yawed flow analysis (Branlard and Gaunaa, 2016).
3.1 The right vortex cylinder
The equations of the inductions of a right vortex cylinder have been derived in detail by Branlard and Gaunaa (2015b). The most important equations and conclusions are summarized in this section. A cylindrical vortex sheet can be decomposed into tangential and longitudinal vorticity components. The strength of the tangential vorticity on the vortex cylinder is the ratio of the total vorticity strength to the helical pitch h. The total vorticity strength of the vortex cylinder ΔΓ_{tot} is equal to the trailed vorticity strength of all blades.
The tangential vorticity contributes to both axial and radial induced velocities. For the vortex cylinder with radius R at axial position x=0, the axial and radial induced velocity at the calculation point with radius of r and axial position of x are shown to be in the form of complete elliptic integrals (Branlard and Gaunaa, 2015b).
where
and K(m), E(m) and Π(n,m) are the complete elliptic integral of the first, second and third kind.
The other components of the vortex cylinder, which are the bound vortex disc, the longitudinal vorticity and the root vortex line, only have contribution to the tangential velocity. The tangential induced velocity of the entire flow field is derived by Branlard and Gaunaa (2015b) as follows:
3.2 Superposition of vortex cylinders for planar rotors
Consider the superposition of the Joukowski rotor model to achieve radially varying bound circulation. There will be helical trailed vorticities emanated along each blade with the strength equal to the derivative of the bound circulation strength with respect to the radius. In Joukowski's rotor model, it is assumed that the radial distribution of the bound circulation of all blades is the same. Then consider the corresponding vortex cylinder model that is the limiting case of Joukowski's rotor model, where the number of blades tends to infinity. It is consisted of a superposition of cylindrical vortex sheets with both tangential and longitudinal vorticity. Details of the superposition of vortex cylinders have been described by Branlard and Gaunaa (2015a). The most important aspects are summarized in this section.
With the superposition of the vortex cylinders, the bound circulation is assumed to be piecewise constant along the blade. The blade is discretized radially into n sections, and there will be one calculation point for each section. Consequently, there will be (n+1) trailing points corresponding to (n+1) vortex cylinders. For the innermost part of the rotor, which is usually the rotor hub and is defined from the center of rotation to the beginning of the first blade section, the bound circulation strength is zero since there are no blades. For the ease of notation and calculation, two ghost sections with the index of 0 and n+1 with zero circulation strength are introduced. For the system, the number of unknown variables is n, which is equal to the number of sections. A sketch of the superposition of the cylindrical vortex system is shown in Fig. 3.
The closure of the system determines the tangential vorticity strength of each vortex cylinder. The closure of the system is determined at the far wake, and the cylindrical vortex sheet is assumed to convect at a constant speed equal to the mean of the two farwake velocities surrounding the vortex sheet (Branlard and Gaunaa, 2015a). It can be shown that the system closure can be performed in the form of helical pitch h(r) or the annulus axial induction factor a(r) at the rotor disc, and the two formulations are equivalent to each other. The system closure in the form of the annulus axial induction factor will be used in this work and will be briefly described.
With the system closure, the axial induction factor of section i is calculated as follows:
The effective thrust coefficient ${C}_{\mathrm{T},\mathrm{eff},i}$ is equal to the thrust coefficient from the Kutta–Joukowski analysis ${C}_{\mathrm{T},\mathrm{KJ},i}$ minus the contribution of wake rotation ${C}_{\mathrm{T},\mathrm{rot},i}$.
where
The wake rotation effect increases toward the rotor rotational axis and decreases when the tipspeed ratio increases. For typical modern wind turbine designs the effect of this term is rather small and may be neglected.
In Eqs. (15) and (16), k_{s,i} is the total nondimensional bound circulation of section i, and λ_{r} is the local speed ratio for the position with radius r. The tangential induction factor ${a}_{i}^{\prime}$ is defined as follows and can be calculated according to Eq. (12) with the condition of x=0.
Considering Eqs. (15) and (16), when the wake rotation effect is included, the aerodynamic loading of a section is dependent on all sections that are further outboard compared to it. The system could be solved from outboard to inboard using Eq. (13), and the system closure is completed when the annulus axial induction factor of all blade sections is calculated. Then, the tangential vorticity of the vortex cylinder that is just inside section i is obtained from the annulus axial induction factor of this section and the neighboring section inside. The equation can be derived from Eq. (9) using the condition of x=0.
3.3 Highthrust correction
In the vortex cylinder model, the relationship between the axial induction factor and the effective thrust coefficient is in the same form as in the momentum theory: ${C}_{\mathrm{T},\mathrm{eff},i}=\mathrm{4}{a}_{i}(\mathrm{1}{a}_{i})$, by inversing Eq. (13). However, when the thrust coefficient (C_{T}) is high, especially when C_{T}>1, the momentum theory breaks down and the vortex cylinder model gives unphysical results. Then, corrections should be made for these highthrust conditions. Different highthrust corrections are available, such as the linear extrapolation by Spera (1994) and the polynomial function of a and C_{T} by Madsen et al. (2020a). To be consistent with the BEM module implemented in the HAWC2 code (Larsen and Hansen, 2007), the polynomial function of a and C_{T} by Madsen et al. (2020a) is chosen. Then, in the system closure, the equation of axial induction factor from the thrust coefficient in Eq. (13) should be replaced by
where the coefficients k_{1} … k_{3} are defined: k_{1}=0.2460, k_{2}=0.0586 and k_{3}=0.0883.
3.4 System closure of nonplanar rotors
The system closure of the vortex cylinder model for planar rotors has been described in Sect. 3.2. It is assumed that each vortex cylinder convects with a constant velocity that is determined in the far wake. With this method of system closure, the relationship between the vortex cylinder model and the momentum theory is revealed by Branlard and Gaunaa (2015b), and the results are generally in good agreement with higherfidelity models. The reason is probably that the error introduced when assuming the convective velocity is constant balances the error introduced when assuming the nonexpanding wake as shown for the uniformly loaded disc (Øye, 1990; Madsen et al., 2007). In the present work, we assume the same method of system closure as for the planar rotor case: the closure is determined at the far wake and can be used for nonplanar rotors with moderate outofplane shapes. With this assumption, the equations are in concise forms and can clearly show the connections between the model for a nonplanar rotor and a planar rotor. However, this assumption does not necessarily hold for extreme cases. To investigate this, a numerical test of blades with relatively large prebend and cone angle will be shown in Sect. 7.
For illustration, we show the superposition of the vortex cylinders for a nonplanar rotor and the corresponding planar rotor with the same radial discretization in Fig. 4. For each section, the corresponding planar rotor has the same total bound circulation (of all blades) ${\mathrm{\Gamma}}_{i}^{\mathrm{pl}}$ as that of the nonplanar rotor ${\mathrm{\Gamma}}_{i}^{\mathrm{np}}$.
3.4.1 Similarity of thrust coefficient
The first similarity is the calculation of the thrust coefficient of the nonplanar rotor and the corresponding planar rotor. For the nonplanar rotor, the Kutta–Joukowski thrust coefficient of section i is calculated from the x component of the force in Eq. (5), obtained using the Kutta–Joukowski analysis.
For the planar rotor, the thrust coefficient is also obtained using Eq. (5):
Since ${\mathrm{\Gamma}}_{i}^{\mathrm{np}}={\mathrm{\Gamma}}_{i}^{\mathrm{pl}}$, and assuming the tangential induction factor a^{′} of the nonplanar rotor and the corresponding planar rotor are the same, which is proven analytically in Sect. 3.4.3, then the Kutta–Joukowski thrust coefficients of the nonplanar rotor and the planar rotor are identical. In addition, the contribution of the wake rotation to the thrust coefficient will also be identical. So, with the given bound circulation distribution of the nonplanar rotor, the thrust coefficient distribution can be directly calculated as if the rotor is planar. This is true regardless of whether the wake rotation effect is included or excluded.
3.4.2 Similarity of tangential and longitudinal vorticity in the wake
The second similarity is that the two vortex wake systems have the same tangential and longitudinal vorticity strength distribution. For the two rotors with the same radial distribution of bound circulation, it can be easily shown that the trailed vorticity strengths between each section are the same.
According to the assumption, the closure of the superposition of the vortex cylinders is determined at the far wake (infinitely far downstream). Therefore, there is no influence of the changed starting position of the vortex cylinders. As a result, both the tangential and longitudinal vorticity of the nonplanar rotor wake is the same as that of the corresponding planar rotor that has the same bound circulation distribution.
According to the description in Sect. 3.4.1, the annulus axial induction of the corresponding planar rotor can be calculated from the thrust coefficient of the nonplanar rotor since the thrust coefficients of the two rotors are identical. Then, since the nonplanar rotor and the corresponding planar rotor have the same tangential vorticity strength, the tangential vorticity strength of the nonplanar rotor can be calculated from the annulus axial induction factor of the corresponding planar rotor using Eq. (20). This means that with a given bound circulation distribution, the tangential vorticity of the nonplanar rotor can also be calculated as if the rotor is planar. The same argument can also be made for the longitudinal vorticity in the wake.
3.4.3 Similarity of tangential induction
It is assumed that the radial distribution of the tangential induction factor of the nonplanar rotor is the same as the corresponding planar rotor in Sect. 3.4.1; this will be analytically proven in this section. Firstly, consider the superposition of the planar vortex cylinders consisting of bound vortex discs as well as tangential and longitudinal trailed vorticities, as shown in Fig. 4. The total strength of the root vortex is zero. Since the tangential vorticities have no contribution to the tangential induction, only the bound vortex discs and the longitudinal vorticity are considered here. For a better illustration, the side view of a part of the axisymmetric vortex system that was illustrated in Fig. 4 is shown in Fig. 5.
Consider the point × in Fig. 5 with radial position of r_{i} that is between R_{i} and R_{i+1} and is just inside the vortex cylinder ($x={\mathrm{0}}^{+}$). The circular contour line C_{1} is perpendicular to the flow, passes point × with radius r_{i} and is centered on the line of r=0. With the axisymmetry of the flow, the tangential velocity has the same value all along the circular contour.
Recall the definition of circulation:
The relationship between the velocity along the contour line C_{1}, which is the tangential velocity at r_{i} and the net circulation through the contour C_{1}, is obtained using the definition of circulation in Eq. (25).
So, the tangential velocity at point × that is just inside the vortex cylinder is
For the point + in Fig. 5 with radial position of r_{i} and just outside the vortex cylinder ($x={\mathrm{0}}^{}$), consider the circular contour line C_{2} with radius r_{i} that is perpendicular to the flow, passing point + and is centered on the line of r=0. Similarly, using the definition of circulation in Eq. (25), the tangential velocity at point + is zero since there is no net circulation passing through the contour C_{2}.
There is a jump of the tangential velocity when the flow passes through the bound vorticity disc. The tangential velocity at the disc should be the mean value of the tangential velocities at the two sides of the disc.
The same result was obtained for the planar rotor case by Branlard and Gaunaa (2015b) by evaluation of the contribution of each component of the cylindrical vortex system to the tangential velocity. As for the planar rotor, the side view of a part of the vortex system of the nonplanar rotor that was illustrated in Fig. 4 is shown in Fig. 6.
Similar as for the planar rotor case, consider the point × in Fig. 6 with radial position of r_{i} that is between R_{i} and R_{i+1} and is just inside the vortex cylinder. With the axisymmetry of the flow and the definition of the circulation in Eq. (25), the tangential velocity is derived to be identical to the planar rotor case in Eq. (27). It is important to point out that the net circulation through the contour (C_{1} or C_{2}) is not influenced by the path of the circulation on either side of the contour.
For the point + that is just outside the vortex cylinder, the tangential velocity is derived to be zero and is identical to the planar rotor case in Eq. (28). Then, with the same argument as for the planar rotor, the tangential velocity at the curved bound vortex surface of the nonplanar rotor should be the mean value of the tangential velocities at the two sides and is in identical form as that for the planar rotor case in Eq. (29).
3.5 Relationship between vortex cylinder model and momentum theory
The relationship between the vortex cylinder model and the momentum theory will be described separately for the planar rotor and the nonplanar rotor.
3.5.1 Planar rotor
For a planar rotor at a high tipspeed ratio, when excluding the contribution of drag to the momentum balancing for determining the induced velocities, the converged results from the momentum theory are equal to those from the vortex cylinder model (Branlard and Gaunaa, 2015a). As the tipspeed ratio decreases, results from the vortex cylinder model and basic 2D momentum theory start differing, especially toward the rotor axis. As also shown in Branlard and Gaunaa (2015a), this difference stems from the pressure drop caused by centrifugal forces due to wake rotation, which is not included in the classic 2D momentum framework that the BEM method is built on. The contribution of wake rotation to the thrust coefficient derived from the vortex cylinder model in Eq. (16) can be applied to the momentum theory as a modification to account for this effect (Branlard and Gaunaa, 2015a). When the wake rotation effect is included in the momentum theory, the annulus axial induced velocity and the tangential induced velocity from the momentum theory and the vortex cylinder model are identical. The radial induced velocity is available from the vortex cylinder model but is not modeled in the momentum theory. However, for straight blades forming the planar rotor, the radial induced velocity has no effect on the convergence calculation or the aerodynamic load calculation. This is because the radial velocity has no contribution to the projection of the velocity into the 2D airfoil section for the straight, nonswept blades that are perpendicular to the rotor axis. These relationships can be written in the following mathematical form, where the subscript “VC” represents the vortex cylinder model, the subscript “MT” represents the momentum theory and the superscript “pl” represents the planar rotor.
3.5.2 Nonplanar rotor
For a nonplanar rotor, the inductions from the momentum theory with wake rotation effect included are equivalent to the inductions from the vortex cylinder model for the corresponding planar rotor and excluding the radial induced velocity. This means the momentum theory equivalently assumes the rotor is planar when calculating inductions. Then, the vortex cylinder model for the nonplanar rotor is equivalent to the momentum theory with the following corrections: for the annulus axial induced velocity, the correction is the difference of the results of the nonplanar rotor and the corresponding planar rotor from the vortex cylinder model. The tangential induced velocity from both methods are the same, as described in Sect. 3.4.3. The radial induction is not available from the momentum theory, so the correction should be the complete radial induction of the nonplanar rotor from the vortex cylinder model.
These relationships between the momentum theory and the vortex cylinder model for the nonplanar rotor are summarized in the equations as follows, where the superscript “np” represents the nonplanar rotor.
where the corrections are
Some important aspects of the implementation of the lowfidelity models that use blade element theory and rely on the 2D airfoil data are briefly discussed. They are important for the load calculation and to get good agreement with the higherfidelity models.
4.1 Impact of unsteady airfoil aerodynamics on steady state
For the blade with outofplane shapes, it is necessary to include the unsteady airfoil aerodynamics model (usually referred to as the dynamic stall model), even for the steadystate simulation. Otherwise, the results of the tangential forces will have a visible error. The reason originates from the conclusions of unsteady 2D aerodynamics: the correct circulatory lift can be obtained if the magnitude of the effective angle of attack is determined at the threequarterchord point, but the direction of it should be determined by the velocity at the quarterchord point (Gaunaa, 2002; Bergami and Gaunaa, 2012). For the lowfidelity model using blade element theory, the angle of attack is usually determined at only one calculation point per section. For instance, for the BEM module in the HAWC2 code, it corresponds to the threequarterchord point. The direction of the circulatory lift is then equivalently transformed to perpendicular to the velocity at the quarterchord point by including an additional torsion rate drag that is proportional to the circulatory lift. In addition, the noncirculatory part of the lift force should be correctly included. For the blade with outofplane shapes, even when the rotational speed is constant, the midchord point will experience a component of acceleration that is perpendicular to the chord due to the projection of centrifugal acceleration. In addition, the angular velocity vector will also have a projection in the 2D section that will result in an effective pitching motion of the airfoil section. The contribution of both the midchord acceleration and the torsion rate should be included when calculating the aerodynamic force. For details, see the reports by Hansen et al. (2004) and Pirrung and Gaunaa (2018).
4.2 Curved blade length projection correction
For blades with inplane or outofplane shapes, the curved blade length in an elementary annulus (ds) is different from the change of the radius (dr). Then, for the momentum analysis in a stream tube, it is necessary to multiply the local thrust and torque coefficient with the term $\frac{\mathrm{d}s}{\mathrm{d}r}$ to account for this difference (Madsen et al., 2020a). For the Kutta–Joukowski analysis of the nonplanar rotor in Eq. (5), the term $\frac{\mathrm{d}s}{\mathrm{d}r}$ is already included in the equation of force. So, it is not necessary to include this term again in the thrust coefficient in Eq. (22) during the convergence calculation using the vortex cylinder model.
The BEM method is the blade element theory coupled with the momentum theory. Similarly, the vortex cylinder model should be coupled with the blade element theory for the aerodynamic load calculation of a rotor with finite number of blades. The link between the blade element theory and the vortex cylinder method is the relationship of the blade bound circulation and the trailed vorticity strength of the vortex cylinder. The trailed vorticity strength is calculated from the total bound circulation of the two neighboring sections using Eq. (24). The total bound circulation is calculated from the blade bound circulation and assuming all blades have the same bound circulation strength.
The blade bound circulation is calculated from the circulatory part of the lift coefficient ${C}_{\mathrm{L}}^{\mathrm{C}}$. For quasisteady simulations, the circular part of the lift coefficient can be replaced by the quasisteady lift coefficient ${C}_{\mathrm{L}}^{\mathrm{QS}}$.
In this section, the coupling of the vortex cylinder model and the blade element theory will be firstly described in Sect. 5.1. The application of the tiploss factor for the nonplanar rotor is important and is described in Sect. 5.2. Finally, the implementation of the proposed vortex cylinder model is summarized in Sect. 5.3.
5.1 Vortex cylinder model as a correction to BEM or as the full model
According to the relationship between the momentum theory and the vortex cylinder model described in Sect. 3.5, there are two possible methods of using the vortex cylinder model when coupling with the blade element theory. Firstly, the vortex cylinder model can be used as a modification to the existing BEM model, which is named and labeled as the BEMVC model. Otherwise, the vortex cylinder model can be considered as the complete replacement of the momentum theory in the BEM method. With the blade element theory coupled with the vortex cylinder model, we have the blade element vortex cylinder (BEVC) model, which does not include any momentum theory results. The reader may argue that the BEVC model should completely replace the BEM model for the prediction of the aerodynamic loads in an existing aeroelastic code. As a first step, the authors recommend the use of the vortex cylinder model as a correction to the BEM method (BEMVC). This is because the framework of the BEM method with many submodels and corrections has been implemented in the aeroelastic codes and has been intensively tested.
5.2 Tiploss correction
Prandtl's tiploss factor is commonly applied to the BEM method to account for the difference between a finite number of blades and the assumption of an infinite number of blades in the momentum theory (Glauert, 1935; Sørensen, 2015). Similarly, the vortex cylinder model also assumes infinite number of blades. Also, considering the relationship between the momentum theory and the vortex cylinder model discussed in the previous sections, a tiploss correction should be applied to the vortex cylinder model. The tiploss factor F presented by Glauert (1935) was implemented in the BEM module in the HAWC2 code (Madsen et al., 2020a):
where φ is the inflow angle.
The tiploss correction is applied by scaling the thrust coefficient with the inverse of the tiploss factor when calculating the blade axial induction.
where the subscript “B” indicates the induction at the blade.
5.2.1 Tiploss for nonplanar rotor
Care should be taken when applying the tiploss correction to the nonplanar rotor. The first aspect is the angle to use when calculating the tiploss factor in Eq. (35). The inflow angle in the rotor coordinate system, which is the flow angle seen by the rotor plane, is usually used (Madsen et al., 2020a). Another possible choice is the flow angle in the sectional coordinate system, which is the flow angle seen by the 2D airfoil section. For planar rotors, it is not necessary to distinguish between them because they are identical. Since the tiploss factor is developed for planar rotors, it is not possible to analytically show which one is better than the other when applied to nonplanar rotors. In a preliminary numerical investigation that is not reported in the present work, it was discovered that results when using the sectional flow angle to calculate the tiploss factor are in slightly better agreement with the higherfidelity models compared to using the inflow angle. As a result, the sectional flow angle is recommended when calculating the tiploss factor and is used in the BEMVC method when calculating the results in Sect. 7.
The second aspect is that the tiploss factor is only to model the amplified axial induction at the blade compared to the annulusaveraged axial induction, and it should not directly change the trailed (tangential or longitudinal) vorticity strength of the vortex cylinders. Recall the similarity of the vortex cylinder model for the nonplanar rotor and the corresponding planar rotor with the same circulation distribution described in Sect. 3.4. The correct implementation of the tiploss correction in the vortex cylinder model for the nonplanar rotor could be considered as a twostep approach.
In the first step, the axial induction factor at the blade of the corresponding planar rotor with the tiploss correction is calculated using Eq. (36). The second step is to calculate the difference between the annulus axial induction of the nonplanar rotor and the planar rotor using the vortex cylinder model. From the effective thrust coefficient, the annulusaveraged axial induction factor of the planar rotor is calculated.
where the subscript “∞” represent infinite number of blades or annulusaveraged value.
The tangential vorticity of the vortex cylinder is calculated from the annulus axial induction of the planar rotor ${a}_{\mathrm{\infty}}^{\mathrm{pl}}$ using Eq. (20) and is duplicated here with the update notation.
Then, the annulus axial induction factor of the nonplanar rotor ${a}_{\mathrm{\infty}}^{\mathrm{np}}$ can be calculated using Eq. (9). Finally, the axial induction at the blade section i of the nonplanar rotor is then equal to the sum of the blade axial induction of the planar rotor and the difference of the annulus axial induction of the nonplanar rotor and the planar rotor.
The tiploss factor is only applied to the axial induction but not applied to the tangential or radial induction, which is following the application of the tiploss correction in the BEM module in the HAWC2 code.
5.2.2 Erroneous implementation
If the model is used without clearly distinguishing between the axial induction on the blade and the annulusaveraged axial induction, the resulting system closure could be wrong. If using the blade axial induction factor ${a}_{\mathrm{B}}^{\mathrm{pl}}$ instead of the annulus axial induction factor ${a}_{\mathrm{\infty}}^{\mathrm{pl}}$ to calculate tangential vorticity in Eq. (38), the tangential vorticity will be directly scaled by the tiploss factor, which is unphysical. Then, the annulus axial induction and the radial induction calculated using Eqs. (9) and (10) will be directly scaled by the tiploss factor due to the wrong tangential vorticity.
For the planar rotor with straight blades, the calculated aerodynamic loads on the blade using the erroneous method will still be correct. The tangential vorticity from the erroneous method is wrong and the radial induction calculated using Eq. (10) will then be wrong. However, the blade axial induced velocity and the tangential induced velocity are correctly calculated. In addition, the radial induction has no contribution to the aerodynamic loads because it has no contribution to the flow seen by the 2D section when the blade is straight and the rotor is planar, as described in Sect. 3.5.1.
5.2.3 Other implementation of tiploss correction
The tiploss correction used in the present work is scaling the thrust coefficient when calculating the blade axial induction and is actually only applied to the planar part of the axial induction. There are other definitions of the tiploss factor, such as the ratio of the blade axial induction and the annulusaveraged axial induction (Branlard and Gaunaa, 2014). Then, it is possible to directly utilize the tiploss factor as
However, consider that the tiploss factor is originally developed for planar rotors. As in Prandtl's simple model of system of material sheets, the flow will go around the vortex disc edges. Also for the modern definition of the tiploss factors (Branlard and Gaunaa, 2014), a planar rotor disc is always assumed. Then, the method of directly applying Prandtl's tiploss factor in Eq. (35) to the nonplanar axial induction as in Eq. (40) is then without a clear physical background and thus not recommended.
5.3 Algorithm of proposed vortex cylinder models
As has been described previously in this section, the proposed vortex cylinder model can be used in two ways: either used as a correction to the BEM method (BEMVC) or solely used and coupled with the blade element theory (BEVC). Details of the implementation of both methods have been described previously in this work and are summarized in Algorithm 1.
The higherfidelity models for the comparison are the Navier–Stokes solver EllipSys3D (Michelsen, 1992, 1994; Sørensen, 1995) and the liftingline module in the aerodynamic solver MIRAS (RamosGarcía et al., 2016), both developed at the Technical University of Denmark (DTU). The lowerfidelity aerodynamic models used for comparison are the BEM method, the BEM method with radial induction correction by Madsen et al. (2020a) (BEMu_{r}) and the proposed BEMVC method. Details of these model setups are given in this section.
6.1 Navier–Stokes solver
The pressurebased incompressible threedimensional solver EllipSys3D was used to solve the Reynoldsaveraged Navier–Stokes equations, using a finitevolume discretization. An inlet/outlet strategy was followed for the boundary conditions of the outer limit of the CFD domain. The flow was assumed to be fully turbulent, and the k−ω SST model (Menter, 1994) was employed. These higherfidelity simulations are labeled in the present work as CFD.
Several rotorresolved meshes were built. They were generated in two consecutive steps, which were fully scripted in order to ensure a similar resulting grid quality. First, a structured mesh of the blade surface was generated with the openly available Parametric Geometry Library (PGL) tool (Zahle, 2019). A total of 128 cells were used in the spanwise direction, and the chordwise direction was discretized with 256 cells. Secondly, the surface mesh was radially extruded with the hyperbolic mesh generator HypGrid (Sørensen, 1998) to create a volume grid. A total of 256 cells were used in this process, and the resulting outer domain was located at approximately 11 rotor diameters. A boundary layer clustering was taken into account, with an imposed first cell height of $\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ m, in order to target y^{+} values lower than unity. The resulting volume meshes accounted for a total of 14.2 million cells. The grid topology is illustrated in Fig. 7, through the particular case of the baseline straight blade.
While a steady solver was used, unsteady separation is expected near the root of the wind turbine blade in operation. To mitigate the effects that this can have on the conclusions of the present work, all the CFD quantities were averaged for the last 350 iterations.
6.2 Liftingline solver
The liftingline module in the aerodynamic solver MIRAS (RamosGarcía et al., 2016) is implemented as a timemarching approach and uses the 2D airfoil data. This work uses a modified version of the liftingline module that includes the influence of the curved bound vortex on the induced velocity as described by Li et al. (2020), which is labeled as LLmod in that work. The bound vorticity is located at the quarterchord line, and the calculation points are placed on the threequarterchord line. The influence of the curved bound vortex is modeled by adding the difference of the induced velocity due to the 3D bound vorticity and an imaginary 2D bound vorticity (infinitely long line vortex) evaluated at the threequarterchord point to the induced velocity of the blade section (Li et al., 2020). The curved bound vortex influence is assumed to be constant along the chord. The angle of attack to determine the circulatory lift and drag coefficient is the angle of attack at the threequarterchord point. The flow environment to determine lift and drag, especially the direction of the lift and drag force, is at the quarterchord point. In addition, as the discussion in Sect. 4.1, the noncirculatory lift is included in the liftingline model when calculating the aerodynamic forces. For the setup in this study, each blade is discretized into 50 sections with cosine spacing. Each simulation is calculated for 20 000 time steps and each step correspond to 1.5^{∘} of azimuthal angle, resulting in a total of 83.3 revolutions. The airfoil data are from 2D fully turbulent CFD results (Bortolotti et al., 2019). The vortex core size in the calculations is 0.1 % of the local chord length. The first row of trailed vorticities begins from the trailing edge of the blade.
6.3 Lowfidelity models
Three lowfidelity aerodynamic models are used for the comparison. The first one is the BEM method implemented in the HAWC2 code version 12.8 (Larsen and Hansen, 2007). The second one is the BEM method with radial induction correction in Eq. (6) (BEMu_{r}) in the same version of HAWC2 code and is described by Madsen et al. (2020a). The third one is the BEMVC method proposed in this work and is implemented in a test version of the HAWC2 code based on version 12.8. As has been discussed in Sect. 5.1, the proposed BEMVC method utilizes the vortex cylinder model as a correction to the existing BEM method. The results should be identical to the BEVC model when the drag is excluded in the momentum balancing. So, the results from the BEVC model, which solely uses the vortex cylinder model and does not directly use the momentum theory, are not shown. For the lowfidelity models in the HAWC2 code, each time step corresponds to 0.01 s and each simulation is calculated for 700 s to get the steadystate value. Each blade is discretized radially into 80 sections. The airfoil data are also from 2D fully turbulent CFD results and are identical to those used in the liftingline method. As has been described in Sect. 5.2, the flow angle seen by the airfoil section is used to calculate the tiploss factor for the BEMVC method. For the BEM method and BEMu_{r} method, the original implementation of the tiploss factor in the HAWC2 code using the inflow angle is applied. Since all three lowfidelity models are implemented in the HAWC2 code, it is then guaranteed the transformation of the velocity and force between different coordinate systems during the computation is consistent.
In this section, the distributed aerodynamic load in the axial and tangential direction, as well as the integrated loads of aerodynamic thrust and power from different lowfidelity models, is compared with results from higherfidelity models. The axial and tangential loads are defined to be positive when aligned with x and z coordinate respectively as defined in Fig. 1. The higherfidelity models are the Navier–Stokes solver (CFD) and the liftingline method (LL) as described in Sect. 6. By comparing the results from the three lowfidelity models with higherfidelity models, it will be highlighted to which extent the influence of the nonplanar rotor geometry can be correctly modeled by each of the lowerfidelity models.
7.1 Test cases
There are five different wind turbine blades used for the comparison; all of them are based on the IEA10.0198 10 MW reference wind turbine (RWT) (Bortolotti et al., 2019). The baseline straight blade is modified by aligning the halfchord line to a straight main axis. For the upwind dihedral blades, the main axes that determine the planforms are obtained from modified Bézier curves which are parameterized with dihedral ratio ${\stackrel{\mathrm{\u203e}}{r}}_{s}$, dihedral magnitude Δd and tip dihedral angle Λ_{tip}, which are illustrated in Fig. 8. In addition, some cases with large cone angles are applied to these dihedral blades to exploit the range of capability of the models. The radius of the unconed rotor is 99 m, of which the hub radius is 2.8 m. The parameterization of the dihedral blades is very similar to that for the previous study of swept blades (Li et al., 2018). The dihedral blades W1 to W4 correspond to Blade1 to Blade4 in the previous study but with outofplane shapes (dihedral) instead of inplane shapes (sweep). The blades are assumed to be stiff, which means the effect of elastic deformation is not included. In addition, the pitch angle is zero for all test cases.
The parameters of these upwind dihedral blades are summarized in Table 1. The main axes of these dihedral blades are illustrated in Fig. 9. The purpose of having dihedral blades with different dihedral magnitude and different tip dihedral angle is to represent different possible shapes of the dihedral blades.
The airfoils are aligned perpendicular to the curved main axis, which is the halfchord line. The chord and twist distribution of the dihedral blades remain unchanged compared to the baseline straight blade. The radius of the dihedral blades is identical to that of the baseline straight blade, but the curved blade length is increased due to the dihedral. For the simulations in this section, the uniform inflow of 8 m s^{−1} with no yaw error is applied to the rotor with a constant rotational speed of 0.855 rad s^{−1}. For the unconed rotors, the tipspeed ratio is 10.58. At this operational condition, the thrust coefficient of the unconed rotor with baseline straight blades is 0.90 and the rotor power coefficient is 0.46, as predicted using the BEM method. At a radius of 70 m, the angle of attack predicted by the BEM method is 5.76^{∘}.
7.2 The distributed load
For the test cases described in Sect. 7.1, the distributed aerodynamic loads calculated from different aerodynamic models are summarized and are compared in this section. For the calculation of the aerodynamic loads, both lift and drag force are included. In this study, the focus is on the influence of the blade dihedral on the loads. The nearroot region (i.e., up to an approximate radius of 20 m) experienced flow separation in the CFD solution, and it is not the focus of this study.
7.2.1 Baseline blade with zero cone
Firstly, the steadystate results of the baseline straight blade without cone calculated from different models are compared and plotted in Fig. 10. Please note that the distributed loads plotted from all models correspond to aerodynamic force per unit radius. The three BEM methods give identical results as expected. The results from the higherfidelity models are similar to the results from the BEM models.
7.2.2 Dihedral blade with zero cone
The steadystate results of the different upwind dihedral blades with zero cone angle are calculated with different aerodynamic models. The axial load and tangential load of the dihedral blade W1 is shown in Fig. 11. The distributed loads correspond to aerodynamic force per unit radius; the curved blade length projection correction described in Sect. 4.2 is applied.
Comparing the distributed load of the baseline straight blade and the dihedral blade W1, it is difficult to draw conclusions for the axial load or the tangential load, because no clear trends can be seen. In order to clearly show the influence of the blade dihedral on the loads predicted by different aerodynamic models, the difference of the loads of the dihedral blade W1 with respect to the baseline straight blade is shown in Fig. 12. This is done by directly subtracting the load of the baseline straight blade from the load of the dihedral blade for the same radius. Consistent with this remark, the following comparisons throughout this work only include the difference of the distributed loads. The absolute loads are not shown since it is more difficult to draw conclusions from them.
It can be seen that for the difference of both the axial and tangential load, both higherfidelity models (CFD and LL) predict a fairly similar pattern of spanwise load redistribution. For the spanwise location that is further inboard compared to where the blade starts to become dihedral, the axial and tangential loads of the upwind dihedral blade are lower compared to the baseline straight blade. When moving from the spanwise location where the blade starts to become dihedral towards halfway until the blade tip, both axial and tangential loads are also lower compared to the baseline. When moving further towards the tip, both axial and tangential loads are then increased compared to the baseline until the blade tip.
A similar pattern of spanwise load redistribution was also observed for swept blades in the previous works (Li et al., 2018, 2020, 2021). However, the redistribution of the loads only takes place where the blade is swept. The loads of the swept blade and the straight blade are almost identical for the inboard part of the blade, where the blade is still straight. Instead, for the dihedral blade, the influence of the blade dihedral has a pronounced influence on the inboard part of the blade that is still straight. This means the blade dihedral at the outboard part of the blade has an influence throughout the spanwise locations of the blade instead of only on the part of the blade that has a dihedral shape. This could be explained by the outofplane blade shape moving the starting position of the trailed vortex system in the axial direction. That, in turn, could effectively move the inner parts of the rotor further into or out of the induction field created by the vortex sheets trailed from the outer sections. By contrast, the influence of blade sweep on the trailed vortex is only on the azimuthal starting position of the trailed vortex. If we consider the trailed vortex as a frozen helical wake, sweeping the blade would only result in an azimuthal twisting of the helical wake. Therefore, the shape of the vortex wake would almost remain unchanged, so that the influence on the induction would not be global.
For both the ordinary BEM method and the BEM method with radial induction correction (BEMu_{r}), the influence of the blade dihedral is not correctly predicted. For the inboard part of the blade that has no dihedral, both methods predict zero offset of loads. The BEM method predicts lowered axial load for the entire portion of the blade that has a dihedral shape. For the tangential load, the BEM method predicts negligible differences compared to the baseline straight blade. The performance of the BEM method is as expected because of the assumption of radial independence in the stream tube theory and the changed streamwise starting position of the trailed vorticity due to blade dihedral is not modeled. The BEMu_{r} method is able to predict the increase in the load near the tip, for both axial and tangential loads. However, the decrease in the loads near where the blade starts to become dihedral and also further inboard is not predicted, because for the BEMu_{r} method, the radial dependency is modeled to a limited extent and is only on the radial induction but not on the axial induction (Madsen et al., 2020a). In addition, the axial and radial induction from the BEMu_{r} method corresponds to a planar rotor.
In comparison, the proposed BEMVC method correctly predicts the pattern and the magnitude of the load redistribution for both axial and tangential loads. The decrease in the loads further inboard compared to where the blade starts to become dihedral is also well predicted. It should be highlighted that the spanwise location of the crossing of the zero load difference is also in good agreement with the results from higherfidelity models (CFD and LL). The largest error with the BEMVC method is mainly for the tipmost part: the increase in the load is overpredicted for both axial and tangential loads. This could be due to the use of Prandtl's tiploss correction in the model. This model is based on a wake shape corresponding to that from a planar rotor with straight blades. The authors believe that a more advanced aerodynamic model that can replace the current tiploss correction, if coupled with the proposed vortex cylinder model, could have better agreement with higherfidelity models. An example is the nearwake model (Madsen and Rasmussen, 2004; Pirrung et al., 2017), which approximately models the near wake as the helical trailed vorticity and is currently coupled with a farwake model that is based on the momentum theory.
It should be mentioned that the difference between the higherfidelity models (CFD and LL) and the BEMVC method for the baseline straight blade is of similar magnitude as the influence due to blade dihedral. However, since the model predicts the sensitivity of changes in dihedral relatively well, it is favorable to be used for parameter studies and to be eventually integrated in a multifidelity aerodynamic optimization framework. For example, in order to design a rotor with dihedral blades, higherfidelity models could be used for the initial design of a straight blade. Then, the proposed vortex cylinder model could be used to explore the sensitivity of different dihedral parameters on the aerodynamic loads, with a relatively low computational effort.
The results of the other three upwind dihedral blades are shown in Appendix B1. The same conclusions as that for the blade W1 also hold for these results.
7.2.3 Upwind cone
To exploit the range of validity of the proposed method, a large upwind cone of 15^{∘} is applied to the baseline straight blade as well as the blades with upwind dihedral (W1 to W4). For the coned cases, the radius of the rotor will decrease compared to the radius of the rotor having the same blades but with zero cone. For better comparison, the abscissa in the figures correspond to the radius of the blade without cone, and the loads are defined as force per unit radius. The factor of $\frac{\mathrm{d}s}{\mathrm{d}r}$, which now equals to the secant of the sum of the cone angle (θ_{c}, positive when cone upwind) and the dihedral angle κ, is multiplied by the loads of the coned blades.
For the baseline straight blade with 15^{∘} of cone upwind, the difference of the loads compared to the straight blade without cone is plotted in Fig. 13. For the upwind dihedral blades with 15^{∘} of cone further upwind, the difference of the loads compared to the baseline blade with the same upwind cone angle is calculated. The results of the upwind coned dihedral blade W1 is shown in Fig. 14. The results of the other dihedral blades with further upwind cone are summarized in Appendix B2.
For the straight blade with large cone in Fig. 13, the results from the proposed method are in good agreement with the results from higherfidelity models (CFD and LL). The proposed method predicts the same trends as the higherfidelity models, while the BEM method and BEMu_{r} method predict different trends. For the upwind dihedral blades with large cone angle in Fig. 14 and in Appendix B2, the results from the proposed method are also in better agreement with the higherfidelity model compared to the BEM method and BEMu_{r} method. However, the results from the proposed method have some differences compared to the higherfidelity model probably due to the limitation of the current tiploss correction, especially near the blade tip that has large blade dihedral.
7.2.4 Downwind cone
To further exploit the range of validity of the proposed method, a large downwind cone angle of 15^{∘} is applied to the baseline straight blade as well as the blades with upwind dihedral (W1 to W4). As has been discussed, the radius of the coned rotor will change compared to the radius of the rotor with the same blades but without cone. For better comparison, the abscissa of the figures corresponds to the radius of the blade without cone, and the loads are again defined as force per unit radius. The term $\frac{\mathrm{d}s}{\mathrm{d}r}$ calculated from Eq. (41) should be multiplied by the loads of the coned blades. For the downwindconed straight blade, the difference of the loads compared to the baseline straight blade without cone is plotted in Fig. 15.
For the upwind dihedral blades with 15^{∘} of cone downwind on top of it, the difference of the loads compared to the baseline blade with the same cone angle downwind is calculated. The results of the coned dihedral blade W1 are shown in Fig. 16. The results of the other upwind dihedral blades with downwind cone are summarized in Appendix B3.
It can be seen that for both the straight blade and the upwind dihedral blade with large downwind cone, the results from the proposed method (BEMVC) are in good agreement with the results from the higherfidelity models (CFD and LL). On the other hand, the BEM method and the BEM method with radial induction correction (BEMu_{r}) are not able to correctly predict the difference of the loads.
7.3 Integrated aerodynamic loads
The integrated aerodynamic loads, which are the aerodynamic power and thrust from different models, are compared in this section. Please note that when comparing the integrated aerodynamic loads, errors in the distributed loads may cancel out. So, it is important to bear in mind that the performance of the different aerodynamic models is not fully represented by their abilities to predict the total aerodynamic power or thrust of the rotor. The aerodynamic force (per unit length of radius) F^{*} on each blade section is composed of the axial force ${F}_{\mathrm{a}}^{*}$, the tangential force ${F}_{\mathrm{t}}^{*}$ and the radial force ${F}_{\mathrm{r}}^{*}$. They are defined to be positive when aligned with x, z and y coordinate respectively as defined in Fig. 1. For the calculation of the aerodynamic load, both lift and drag force are included. The position of applying the force on the blade section is p. For simplicity, we use the halfchord point coordinate as p. This means we neglect the distance between the halfchord point and the quarterchord point and also the contribution of the twist to the vector p. In addition, the contribution of the sectional airfoil aerodynamic moment (calculated from C_{m}) to the aerodynamic momentum of the rotor is also neglected. Then, the distributed aerodynamic moment from each blade section is
The x component of the aerodynamic moment M is the contribution to aerodynamic torque. Then, the aerodynamic power of the rotor is the integrated contribution of M_{x} of all N_{B} blades at the rotational speed of Ω:
The aerodynamic thrust of the rotor is the total contribution of the axial force of all N_{B} blades:
The aerodynamic power and thrust of the rotor with baseline straight blades without cone are summarized in Table 2. For the rotors with dihedral blades without cone, it is difficult to directly draw conclusions from the absolute value of power and thrust. To better illustrate and compare the integral effects of the rotor dihedral represented by the aerodynamic power and thrust predicted using different methods, the relative difference of the aerodynamic power and thrust with respect to the baseline rotor from each method are calculated and are summarized in the bar plots in Figs. 17 and 18.
For the aerodynamic power, the relative change predicted by LL is underestimated compared to the prediction by CFD, but the results are showing similar trends. One of the reasons could be the use of the 2D airfoil data in the liftingline method. The ordinary BEM method predicts almost no influence of blade dihedral on power, except for W1, which predicts the correct direction in which the power increases but underestimates the magnitude. The BEMu_{r} method predicts the same direction as CFD and LL in which the power of the dihedral rotors are increased compared to the baseline. However, the magnitude of the increment is overestimated by a factor of approximately 2. Comparing to the other two BEM methods, the relative increment of power predicted by the proposed BEMVC method is in better agreement with the predictions by CFD and LL.
For the aerodynamic thrust, the magnitude of the relative decrement is overestimated by approximately 20 % by LL compared to CFD. The BEM method predicts the correct trend that the thrust decreases but the magnitude is underestimated compared to LL and CFD. The BEMu_{r} method predicts a small increase in aerodynamic thrust instead of a decrease in thrust as predicted by the higherfidelity models. The relative change of the thrust predicted by the proposed BEMVC method is very similar to the results predicted by LL, and the relative difference of the predicted relative change is less than 20 %.
In summary, the proposed BEMVC model is in better agreement with higherfidelity models when predicting the integrated aerodynamic power and thrust of the dihedral rotor, compared to the ordinary BEM method. The BEM method with radial induction correction (BEMu_{r}) predicts worse compared to the BEM method for both aerodynamic power and thrust.
7.4 Low loading cases
The results shown up to this point in this work all correspond to fairly high thrust coefficients around approximately 0.9. In this section, the results for operational conditions corresponding to lower thrust coefficients are shown. For simplicity, three of the operational conditions defined in the IEA Wind TCP Task 37 report (Bortolotti et al., 2019) that anticipated lower thrust coefficients are used. The operational conditions of these lower loading cases are summarized in Table 3, by means of wind speed, tipspeed ratio and pitch angle θ_{p}. For each case, the rotational speed is 0.909 rad s^{−1}, and uniform inflow with no yaw error is applied to the rotor. For each operational condition, the baseline straight blade is pitched with θ_{p}, and the thrust coefficient predicted by BEM is also included in Table 3. In order to allow the comparability of the results of this straight blade with the blade accounting for a dihedral angle, the effects of pitching were introduced as an additional twist instead of a rotation around the pitching axis. The motivation for this was related to the introduction of inplane geometry components by pitching the blade. The twisted blade is named W1U^{*}, where ^{*} is the wind speed.
For the results shown in the previous sections, the liftingline (LL) results were in good agreement with the CFD results. Therefore, only the LL method is used to generate the higherfidelity results for the comparison in this section.
First, the differences of the axial load and tangential load of the dihedral blade W1U12 and the pitched straight blade are shown in Fig. 19. As for the high loading case in Fig. 12, both the ordinary BEM method and the BEM method with radial induction correction (BEMu_{r}) are not able to correctly predict the influence of the blade dihedral. In comparison, the proposed BEMVC method correctly predicts the shape and magnitude of the load redistribution for both axial and tangential loads. The decrease in the load further inboard compared to where the blade starts to become dihedral is also well predicted. As for the high loading case, the largest difference between the LL and the BEMVC results is mainly for the tipmost part. As previously mentioned, this could be due to the use of Prandtl's tiploss correction in the model.
The results for the dihedral blades W1U15 and W1U20, at wind speeds of 15 and 20 m s^{−1} respectively, are shown in Figs. 20 and 21. The same conclusion can be made from these lower loading results. However, with the decrease in the thrust coefficient, the difference between the prediction by the BEM method and the LL method decreases. This is because the relative importance of the induced velocities decreases in these cases, making the performance more directly dictated by the velocity components from the inflow and the rotor rotation and less from the induced velocities.
7.5 Computational effort
The computational efforts to obtain the steadystate results used in the present work, measured in CPU time, are summarized in this section. The CFD computations using EllipSys3D were performed on DTU's highperformance computing (HPC) cluster Jess, in which each node has 20 cores running at 2.8 GHz. All the CFD simulations of the present work required a wall clock time of approximately 3.5 h when using 216 cores. The liftingline (LL) computations using the MIRAS code were performed on the Sophia HPC cluster, in which each node has 32 cores running at 2.9 GHz. Each of the LL simulations in the present work required a wall clock time of approximately 100 h when using 32 cores. Note that the computational time for the LL method in the MIRAS code in this study is relatively high, because the settings were chosen to achieve the highest possible fidelity irrespective of the computational cost. Therefore the computational effort for the MIRAS calculations in this work is not indicative of the performance for normal use of the tool. Settings that increased the computational effort in this work are small time steps, not using farwake cutoff, etc. The computational time is expected to be largely decreased if efforts are dedicated to improving the simulation setup. However, this is beyond the scope of the present work.
The computations using the HAWC2 code were performed on a single core of a 2018 workstation at 4.8 GHz. The simulations were performed with structural properties included and with large stiffness to approximate stiff structures. The simulations were run for 600 s in the simulation time to reach steady state. The simulations required a wall clock time of approximately 600 and 650 s for the BEM method and the BEMVC method, respectively. For a standalone version of the BEM method or the BEMVC/BEVC method, one steadystate computation can be done in much less than 1 s using a single CPU core.
A new computationally efficient method for the aerodynamic load calculation of nonplanar rotors is described. The method is based on the vortex cylinder model and can be used in two ways: either as a correction to the currently widely used blade element momentum (BEM) method or as the main model, replacing the BEM method in the engineering modeling complex. For uniform inflow that is perpendicular to the rotor plane, the influence of the blade outofplane shapes on the distributed aerodynamic loads, measured by the difference of the loads between the nonplanar rotor and the planar rotor, is shown to be in good agreement with higherfidelity models. The predicted distributed and integrated aerodynamic loads are in better agreement with higherfidelity models than the baseline BEM method and also a BEM method with a radial induction correction. While the present work focused on stiff geometries, the developed framework would be able to handle outofplane deflections during aeroelastic simulations accounting for blade elasticity, without any loss of generality. The new model is approximately as numerically efficient as ordinary BEMbased models, which makes it favorable for aeroservoelastic simulation as well as design optimization of horizontalaxis wind turbines whose blades have outofplane shapes. Therefore, the authors recommend the use of the proposed model as a correction to the existing BEM codes.
For the future work on the model applications, it would be interesting to use both the standard BEM method and the proposed method for the aerodynamic or aeroelastic design of a nonplanar rotor under the same constraints. Higherfidelity models, such as CFD or liftingline method, could be used for the benchmark of the different designs, as done in the present work. The method is also favorable for integration in a multifidelity aerodynamic design framework. There are also several ways in which future work could improve the model. Firstly, it would be favorable to have modifications to the existing Prandtl tiploss correction. For example, it is possible to use the distance between the tip vortex and the calculation point when calculating the correction for a nonplanar rotor, instead of using the radial distance as currently implemented in the model. Secondly, it would be beneficial to further develop the model for the application of blades with both inplane and outofplane shapes. One possible track of the development is to couple the vortex cylinder model and the nearwake model, which approximately models the near wake as helical trailed vorticities and is currently coupled with a farwake BEM method. Thirdly, it would be interesting to investigate the unsteady effects of the nonplanar rotor, such as aerodynamic damping and dynamic inflow effect. Fourthly, it would be beneficial to further develop the vortex cylinder model for the application of nonplanar rotors in yawed flow. Finally, further development of the model focusing on analytical gradients would be favorable for application in a gradientbased wind turbine design optimization framework.
a  axial induction 
a^{′}  tangential induction 
C_{L}  lift coefficient 
C_{D}  drag coefficient 
C_{m}  moment coefficient 
C_{Q}  torque coefficient 
C_{T}  thrust coefficient 
C_{T,eff}  effective thrust coefficient 
C_{T,KJ}  Kutta–Joukowski thrust coefficient 
C_{T,rot}  thrust coefficient due to wake rotation 
C_{T,av}  averaged coefficient for the radial induction function 
Δd  dihedral magnitude 
f  lift force vector on the blade with the definition of force per unit length of curved blade length 
f^{*}  lift force vector on the blade with the definition of force per unit radius 
F^{*}  aerodynamic force vector on the blade with the definition of force per unit radius 
F  tiploss factor 
h  helical pitch 
k  factor for the calculation of the elliptic integral 
k_{1}, k_{2}, k_{3}  factors for the relationship between axial induction and thrust coefficient 
k_{s}  normalized sectional circulation of the vortex cylinder 
M  the aerodynamic moment vector 
N_{B}  number of blades 
P  aerodynamic power of the rotor 
r  radius of the calculation point 
R  radius of the vortex cylinder 
R_{tot}  radius of the rotor 
T  aerodynamic thrust of the rotor 
u_{a}  axial induced velocity 
u_{t}  tangential induced velocity 
u_{r}  radial induced velocity 
Δu_{a}  the correction to the axial induced velocity 
Δu_{r}  the correction to the radial induced velocity 
U_{0}  wind speed 
V_{rel}  relative velocity 
x  axial position of the calculation point with respect to the vortex cylinder 
Greek letters  
Γ  bound vorticity strength of the vortex cylinder, equal to the bound vorticity strength of all blades 
Γ^{B}  blade bound vorticity strength 
Γ_{root}  root vortex 
ΔΓ  trailed vorticity strength of the vortex cylinder 
γ_{b}  radial bound vorticity strength 
γ_{l}  longitudinal vorticity strength of the vortex cylinder 
γ_{t}  tangential vorticity strength of the vortex cylinder 
Λ_{tip}  tip dihedral angle 
φ  inflow angle 
ρ  density of air 
θ_{c}  cone angle 
κ  dihedral angle 
Ω  rotor speed 
λ_{r}  speed ratio at radius r 
Subscripts  
a  in the axial direction 
t  in the tangential direction 
r  in the radial direction 
i  at blade section i 
B  the value at the blade 
∞  the annulusaveraged value 
eff  effective value 
tot  the total value 
MT  from the momentum theory 
BEM  from the BEM method 
VC  from the vortex cylinder model 
Superscripts  
B  the value at the blade 
np  the nonplanar rotor 
pl  the planar rotor 
QS  quasisteady 
C  circulatory part 
B1 Zero cone angle
The difference of the loads of the dihedral blades (W2 to W4) with zero cone compared to the baseline straight blade without cone.
B2 Upwind cone
The difference of the loads of the dihedral blades (W2 to W4) with 15^{∘} of cone upwind compared to the baseline straight blade with the same upwind cone.
The code is not provided because it is part of the inhouse tool HAWC2. However, the reader is able to recover the code using the algorithm provided in this paper.
The airfoil data are from 2D fully turbulent CFD results (https://www.osti.gov/biblio/1529216ieawind tcptasksystemsengineeringwindenergywp2referencewindturbines; Bortolotti et al., 2019).
AL conducted the study as part of his PhD research. The idea of the proposed model originated from MG. The proposed model was jointly developed by MG, AL and GRP. The similarities of the superposition of the vortex cylinder model of nonplanar rotor and planar rotor were described by AL and MG. The application of the tiploss correction as well as highthrust correction for the nonplanar rotor was described by AL and MG. The implementation of the proposed model in HAWC2 code and the computations using the HAWC2 code were performed by AL with contribution from GRP. The CFD method was introduced by SGH, and the CFD results were computed by SGH. The postprocessing of the CFD results was performed by SGH with contribution from AL. The liftingline results were computed by AL, and the postprocessing was performed by AL. All authors jointly drew the conclusions of the work and contributed to writing this paper.
DTU Wind Energy develops and distributes HAWC2 on commercial and academic terms.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank our colleague Néstor Ramos García in DTU Wind Energy for the suggestions in the liftingline simulation using MIRAS, a vortex code mainly developed by him.
This research has been supported by the Smart Tip project, funded by Innovationsfonden (grant no. 704600023B).
This paper was edited by Alessandro Bianchini and reviewed by two anonymous referees.
Bergami, L. and Gaunaa, M.: ATEFlap Aerodynamic Model, a dynamic stall model including the effects of trailing edge flap deflection, Danmarks Tekniske Universitet, Risø Nationallaboratoriet for Bæredygtig Energi, 52 pp., ISBN 8755039340, 9788755039346, 2012. a
Bortolotti, P., Tarrés, H. C., Dykes, K., Merz, K., Sethuraman, L., Verelst, D., and Zahle, F.: Systems Engineering in Wind Energy – WP2.1 Reference Wind Turbines, Tech. rep., NREL – National Renewable Energy Laboratory [data set], available at: https://www.osti.gov/biblio/1529216ieawindtcptasksystemsengineeringwindenergy wp2referencewindturbines (last access: 2 April 2021), 2019. a, b, c, d
Branlard, E.: Wind Turbine Aerodynamics and VorticityBased Methods, Springer, https://doi.org/10.1007/9783319551647, 2017. a
Branlard, E. and Gaunaa, M.: Development of new tiploss corrections based on vortex theory and vortex methods, J. Phys.: Conf. Ser., 555, 012012, https://doi.org/10.1088/17426596/555/1/012012, 2014. a, b
Branlard, E. and Gaunaa, M.: Superposition of vortex cylinders for steady and unsteady simulation of rotors of finite tipspeed ratio, Wind Energy, 19, 1307–1323, https://doi.org/10.1002/we.1899, 2015a. a, b, c, d, e, f, g, h, i, j, k
Branlard, E. and Gaunaa, M.: Cylindrical vortex wake model: right cylinder, Wind Energy, 18, 1973–1987, https://doi.org/10.1002/we.1800, 2015b. a, b, c, d, e
Branlard, E. and Gaunaa, M.: Cylindrical vortex wake model: skewed cylinder, application to yawed or tilted rotors, Wind Energy, 19, 345–358, https://doi.org/10.1002/we.1838, 2016. a, b
Branlard, E. and Meyer Forsting, A.: Using a cylindrical vortex model to assess the induction zone infront of aligned and yawed rotors, in: Proceedings of EWEA Offshore 2015 Conference, European Wind Energy Association (EWEA), EWEA Offshore 2015 Conference, 10–12 March 2015, Copenhagen, Denmark, 2015. a
Branlard, E. and Meyer Forsting, A.: Assessing the blockage effect of wind turbines and wind farms using an analytical vortex model, Wind Energy, 23, 2068–2086, https://doi.org/10.1002/we.2546, 2020. a
Burton, T., Jenkins, N., Sharpe, D., Bossanyi, E., and Graham, M.: Wind energy handbook, Wiley, ISBN 9781119451167, 2021. a
Crawford, C.: Reexamining the precepts of the blade element momentum theory for coning rotors, Wind Energy, 9, 457–478, https://doi.org/10.1002/we.197, 2006. a, b, c
Gaunaa, M.: Unsteady Aerodynamic Forces on NACA 0015 Airfoil in Harmonic Translatory Motion, PhD thesis, available at: https://orbit.dtu.dk/en/publications/unsteadyaerodynamicforcesonnaca0015airfoilinharmonictran (last access: 18 January 2022), 2002. a
Glauert, H.: Airplane Propellers, in: Division L in Aerodynamic Theory, vol. IV, edited by: Durand, W. F., Springer, 169–360, https://doi.org/10.1007/9783642914874_3, 1935. a, b
Hansen, M. H., Gaunaa, M., and Madsen, H. A.: A BeddoesLeishman type dynamic stall model in statespace and indicial formulations, RisøR1354, Roskilde, Denmark, 2004. a, b
Larsen, J., Nielsen, S., and Krenk, S.: Dynamic stall model for wind turbine airfoils, J. Fluid. Struct., 23, 959–982, https://doi.org/10.1016/j.jfluidstructs.2007.02.005, 2007. a
Larsen, T. and Hansen, A.: How 2 HAWC2, the user's manual, Risø National Laboratory, 71 pp., ISBN 9788755035836, 8755035833, available at: https://orbit.dtu.dk/en/publications/how2hawc2theusersmanual (last access: 18 January 2022), 2007. a, b, c
Leishman, J. G.: Principles of helicopter aerodynamics, Cambridge University Press, ISBN 0521858607, 2005. a
Leishman, J. G. and Nguyen, K. Q.: Statespace representation of unsteady airfoil behavior, AIAA J., 28, 836–844, https://doi.org/10.2514/3.25127, 1990. a
Li, A., Pirrung, G., Madsen, H. A., Gaunaa, M., and Zahle, F.: Fast trailed and bound vorticity modeling of swept wind turbine blades, J. Phys.: Conf. Ser., 1037, 062012, https://doi.org/10.1088/17426596/1037/6/062012, 2018. a, b, c
Li, A., Gaunaa, M., Pirrung, G. R., RamosGarcía, N., and Horcas, S. G.: The influence of the bound vortex on the aerodynamics of curved wind turbine blades, J. Phys.: Conf. Ser., 1618, 052038, https://doi.org/10.1088/17426596/1618/5/052038, 2020. a, b, c, d
Li, A., Pirrung, G. R., Gaunaa, M., Madsen, H. A., and Horcas, S. G.: A computationally efficient engineering aerodynamic model for swept wind turbine blades, Wind Energ. Sci. Discuss. [preprint], https://doi.org/10.5194/wes202196, in review, 2021. a, b
Madsen, H. and Rasmussen, F.: The influence on energy conversion and induction from large blade deflections, in: Wind energy for the next millennium, Proceedings, edited by: Petersen, E. L., Hjuler Jensen, P., Rave, K., Helm, P., and Ehmann, H., James and James Science Publishers, 1999 European Wind Energy Conference and Exhibition, 1–5 March 1999, Nice, France, 138–141, ISBN 190291600X, 1999. a
Madsen, H. Å.: A CFD analysis of the actuator disc flow compared with momentum theory results, in: IEA Joint action. Aerodynamics of wind turbines, edited by: Pedersen, M. B., Technical University of Denmark, Department of Fluid Mechanics, in: 10th IEA Meeting on Aerodynamics, IEA, 16–17 December 1996, Edinburgh, UK, 109–124, 1997. a
Madsen, H. Å. and Rasmussen, F.: A near wake model for trailing vorticity compared with the blade element momentum theory, Wind Energy, 7, 325–341, https://doi.org/10.1002/we.131, 2004. a
Madsen, H. Å., Mikkelsen, R., Øye, S., Bak, C., and Johansen, J.: A Detailed investigation of the Blade Element Momentum (BEM) model based on analytical and numerical results and proposal for modifications of the BEM model, J. Phys.: Conf. Ser., 75, 012016, https://doi.org/10.1088/17426596/75/1/012016, 2007. a
Madsen, H. Å., Bak, C., Døssing, M., Mikkelsen, R. F., and Øye, S.: Validation and modification of the Blade Element Momentum theory based on comparisons with actuator disc simulations, Wind Energy, 13, 373–389, https://doi.org/10.1002/we.359, 2010. a
Madsen, H. Å., Larsen, T. J., Pirrung, G. R., Li, A., and Zahle, F.: Implementation of the blade element momentum model on a polar grid and its aeroelastic load impact, Wind Energ. Sci., 5, 1–27, https://doi.org/10.5194/wes512020, 2020a. a, b, c, d, e, f, g, h, i, j, k
Madsen, H. Å., Zahle, F., Meng, F., Barlas, T., Rasmussen, F., and Rudolf, R. T.: Initial performance and load analysis of the LowWind turbine in comparison with a conventional turbine, J. Phys.: Conf. Ser., 1618, 032011, https://doi.org/10.1088/17426596/1618/3/032011, 2020b. a
Mann, J.: The spatial structure of neutral atmospheric surfacelayer turbulence, J. Fluid Mech., 273, 141–168, 1994. a
Menter, F. R.: Twoequation eddyviscosity turbulence models for engineering applications, AIAA J., 32, 1598–1605, 1994. a
Michelsen, J. A.: Basis3D – a Platform for Development of Multiblock PDE Solvers, Tech. Rep. AFM 9205, Technical University of Denmark, 1992. a
Michelsen, J. A.: Block structured Multigrid solution of 2D and 3D elliptic PDE's, Tech. Rep. AFM 9406, Technical University of Denmark, 1994. a
Mikkelsen, R.: Actuator Disc Methods Applied to Wind Turbines, PhD thesis, Technical University of Denmark, available at: https://orbit.dtu.dk/en/publications/actuatordiscmethodsappliedtowindturbines (last access: 18 January 2022), 2004. a
Okulov, V., Sørensen, J., and Wood, D.: The rotor theories by Professor Joukowsky: Vortex theories, Prog. Aerosp. Sci., 73, 19–46, https://doi.org/10.1016/j.paerosci.2014.10.002, 2015. a
Øye, S.: A Simple Vortex Model, in: Proceedings of the third IEA Symposium on the Aerodynamics of Wind Turbines, ETSU, Harwell, UK, 4.1–4.15, 1990. a
Pirrung, G., Riziotis, V., Madsen, H., Hansen, M., and Kim, T.: Comparison of a coupled near and farwake model with a freewake vortex code, Wind Energ. Sci., 2, 15–33, https://doi.org/10.5194/wes2152017, 2017. a
Pirrung, G. R. and Gaunaa, M.: Dynamic stall model modifications to improve the modeling of vertical axis wind turbines, DTU Wind Energy E0171, DTU, Roskilde, Denmark, 2018. a
RamosGarcía, N., Sørensen, J., and Shen, W.: Threedimensional viscousinviscid coupling method for wind turbine computations, Wind Energy, 19, 67–93, 2016. a, b
Schepers, J. G. and Snel, H.: Joint investigation of dynamic inflow effects and implementation of an engineering method., Tech. rep. ECNC94107, ECN, available at: https://publications.ecn.nl/E/1995/ECNC94107 (last access: 18 January 2022), 1995. a
Sørensen, J. N.: General Momentum Theory for Horizontal Axis Wind Turbines, in: vol. 4, Springer, ISBN 9783319221144, 2015. a
Sørensen, N. N.: General Purpose Flow Solver Applied to Flow over Hills, RisøR827(EN), Risø National Laboratory, Roskilde, Denmark, 1995. a
Sørensen, N. N.: HypGrid2D a 2D Mesh Generator, RisøR1035(EN), Risø National Laboratory, Roskilde, Denmark, 1998. a
Spera, D. A.: Wind turbine technology: fundamental concepts of wind turbine engineering, ASME Press, ISBN 13: 9780791812051, ISBN 10: 0791812057, 1994. a, b
Yu, W., Tavernier, D., Ferreira, C., van Kuik, G. A. M., and Schepers, G.: New dynamicinflow engineering models based on linear and nonlinear actuator disc vortex models, Wind Energy, 22, 1433–1450, https://doi.org/10.1002/we.2380, 2019. a, b
Zahle, F.: Parametric Geometry Library (PGL), Tech. rep., DTU Wind Energy, available at: https://gitlab.windenergy.dtu.dk/frza/PGL (last access: 23 November 2020), 2019. a
 Abstract
 Introduction
 Kutta–Joukowski analysis
 Vortex cylinder model
 Some important aspects in models using blade element theory
 Blade element theory with vortex cylinder model
 The models for comparison
 Results
 Conclusions and future work
 Appendix A: Nomenclature
 Appendix B: Results of the distributed load
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(6745 KB)  Fulltext XML
 Abstract
 Introduction
 Kutta–Joukowski analysis
 Vortex cylinder model
 Some important aspects in models using blade element theory
 Blade element theory with vortex cylinder model
 The models for comparison
 Results
 Conclusions and future work
 Appendix A: Nomenclature
 Appendix B: Results of the distributed load
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References