https://doi.org/10.5194/wes-8-1853-2023
https://doi.org/10.5194/wes-8-1853-2023
Research article |  | 13 Dec 2023

# An analytical linear two-dimensional actuator disc model and comparisons with computational fluid dynamics (CFD) simulations

Abstract

The continuous up-scaling of wind turbines enabled by more lightweight and flexible blades in combination with coning has challenged the assumptions of a plane disc in the commonly used blade element momentum (BEM)-type aerodynamic codes for the design and analysis of wind turbines. The objective with the present work is thus to take a step back relative to the integral 1-dimensional (1-D) momentum theory solution in the BEM model in order to study the actuator disc (AD) flow in more detail.

We present an analytical, linear solution for a two-dimensional (2-D) AD flow with one equation for the axial velocity and one for the lateral velocity, respectively. Although it is a 2-D model, we show in the paper that there is a good correlation with axis-symmetric and three-dimensional (3-D) computational fluid dynamics (CFD) simulations on a circular disc. The 2-D model has thus the potential to form the basis for a simple and consistent rotor induction model.

For a constant loading, the axial velocity distribution at the disc is uniform as in the case of the classical momentum theory for an AD. However, an important observation of the simulated flow field is that immediately downstream of the disc the axial velocity profiles change rapidly to a shape with increased induction towards the edges of the disc and less induction on the central part. This is typically what is seen at the disc in full non-linear CFD AD simulations, which is what we compare with in the paper.

By a simple coordinate rotation the analytical solution is extended to a yawed disc with constant loading. Again, a comparison with CFD, now with a 3-D simulation on a circular disc in yaw, confirms a good performance of the analytical 2-D model for this more complicated flow.

Finally, a further extension of the model to simulate a coned disc is obtained using a simple superposition of the solution of two yawed discs with opposite yaw angles and positioned so the two discs just touch each other. Now the validation of the model is performed with results from axis-symmetric CFD simulations of an AD with a coning of both 20 and 20. In particular, for the disc coned in the downwind direction there is a very good correlation between the simulated normal velocity to the disc, whereas some deviations are seen for the upwind coning.

The promising correlation of the results for the 2-D model in comparison with 3-D simulations of a circular disc with CFD for complicated inflow like what occurs at yaw and coning indicates that the 2-D model could form the basis for a new, consistent rotor induction model. The model should be applied along diagonal lines on a rotor and coupled to an angular momentum model. This application is sketched in the outlook and is a subject for future research.

1 Introduction

The most striking characteristic of the modern, industrial wind turbine development over the last 4–5 decades is the up-scaling. In the beginning of the 1980s there was already an industrial series production of turbines with a rated power of 50–100 kW and rotor diameters of 15–20 m. Many were installed in big wind farms like Palm Springs and Tehachapi in the United States. Today turbines of the maximum size of around 15 MW with rotor diameters of 230–240 m have been introduced in the market, which gives an up-scaling factor of around 15 on rotor diameter over this time span.

It is clear that this up-scaling correspondingly increased the requirements for the aerodynamic and aeroelastic design tools like the modelling of unsteady inflow over the rotor disc, wind shear and wind veer, yawed flow, and wake flow . The origin of the aerodynamic modelling of the common aeroelastic model tools today – like FLEX5 (Flex4) (Øye1996), FAST , BLADED (Bossanyi2003), GAST , Cp-Lambda , FOCUS (WMC2019) and HAWC2 – is the blade element momentum (BEM) theory based on the propeller theory by which originally assumed homogeneous and stationary flow. Various modifications have been implemented in the different codes to adapt to the conditions mentioned above. Most codes, e.g., now use a more local representation of the induction over the rotor plane than on a ring element as originally proposed to better account for varying inflow over the rotor disc from wind shear and turbulence.

An excellent overview of the capability of most of the above-mentioned codes has recently been presented by . This was done by both bench-marking and validation of experimental data as well as of higher-fidelity model data. One conclusion is that the BEM-type codes still show higher deviations for yawed inflow than higher-fidelity codes like computational fluid dynamics (CFD) and vortex-type codes. For yawed flow cases the variability between BEM-type codes is double the variability between CFD codes. This illustrates the lack of a unified implementation of the sub-models in BEM-type codes, e.g., for the modelling of yaw.

The objective with the present work is thus to take a step back relative to the integral momentum theory solution for an actuator disc, which is the core of the BEM model. We will study the actuator disc (AD) flow in more detail to provide a better background for understanding the shortcomings and approximations in the BEM-implemented momentum theory and deviations from higher-fidelity modelling of AD flow such as (1) the velocity distribution at the disc for constant loading, (2) independence of annular stream tubes, (3) flow field for a yawed AD and (4) flow field for a coned AD.

With respect to the axial velocity distribution for a constant loaded AD where the integral momentum theory predicts (assumes) a uniform, constant profile, a number of CFD-based AD simulations have shown higher induction towards the edges of the disc and less induction on the central part . Also the shortcomings of the momentum theory to predict the flow through a coned rotor have been illustrated by CFD AD simulations of and , where considerable deviations from the constant axial velocity profile assumed by the momentum theory are found.

Studies within these four areas have been carried out by many researchers in the past, and it seems that the modelling can be grouped into three approaches: (1) vortex models, (2) combined analytical and numerical models, and (3) Euler-equation-based models which the present model belongs to.

The approach based on vortex models is by far the most common way of studying the AD flow field and not least the use of a vortex cylinder (VC) model with a none-expanding wake, infinite number of blades, infinite tip speed ratio and uniform loading. Some recent studies in relation to the above-mentioned research areas using this model are (relations between momentum theory and vortex theory), (modelling of non-planar rotors), (yawed disc), and (independence of annular stream tubes).

In the category of studies combining analytical and numerical solutions presented a modified BEM for coned rotors. The improvement in the BEM model comprises a proper consideration of the relative placement of the wake and includes the radial-induced velocity. A comparison of the improved BEM with the CFD results of and shows reasonable correlation as concerns the slope of the profiles on the inner part of the disc, whereas the non-linear outboard part required a further modelling of the expansion of the stream tubes.

Another study by shows the non-uniformity of the axial velocity profile caused by the pressure acting at the stream tube annuli, thus making them not independent from each other in contradiction to the independence found by . Finally in this group of models the work by can be mentioned, who used conservation of axial and angular momentum and the calculus of variations to derive the optimal performance of horizontal-axis wind turbines over a wide range of tip speed ratios. In the two-dimensional (2-D) model to be presented here the analytical solution is limited to a constant loading.

In the third model group, the Euler equations and the continuity equation together with specified body forces are the starting point for modelling the AD flow. It seems to be a unique approach, and only the so-called asymptotic acceleration potential method by has some resemblance to the method used in the present research.

It might be considered surprising to study a 2-D AD model for the representation of rotor induction. However, it is expected that the main difference between an axis-symmetric and a 2-D model will be related to the analogy between radial velocity and lateral velocity for the axis-symmetric and 2-D model, respectively, where the last mentioned is expected to have the highest velocities.

The model and the background for its derivation are presented in Sect. 2 followed by flow characteristics of the AD model in axial flow in Sect. 3, including comparisons with axis-symmetric CFD simulations. Section 4 extends the model to yawed inflow and the validation with three-dimensional (3-D) CFD AD simulations. Section 5 presents a further extension for modelling a coned disc and comparisons with axis-symmetric CFD AD simulation of a disc with both 20 and −20 coning. Outlook and conclusions are given in Sects. 6 and 7, respectively.

2 The analytical 2-D AD model

The analytical linear solution (2Dl) of the 2-D AD model was first presented by during a work on the development of the actuator cylinder (AC) model for modelling the aerodynamics of vertical-axis wind turbines (VAWTs). The basic idea behind the AC model was the need to extend the well-known actuator disc concept to an actuator surface (AS) of arbitrary shape where the AS coincides with the swept surface of the considered turbine type (Madsen1983, 1985, 1988). For a straight-bladed VAWT rotor this means a cylindrical shape which thus can be modelled with the AC model.

However, this extension of the AD concept to an AS raised the question of how to compute the flow field through and around an AS. The well-known momentum theory solution in an integral form for an AD cannot just be transferred for use on an AS of arbitrary shape, which would violate the assumptions in the AD momentum model.

## 2.1 Derivation of the descriptive system of equations

We have chosen an approach for the derivation of the equations that follows a general method of for the analysis of flow fields with a main velocity V under the influence of external body forces f. used the method to develop a general theory for wings of finite span and considered in most cases the induced velocities from the action of external forces as perturbations to the main flow field with the velocity V.

presented an analysis of the flow around “the ideal propeller” described as a thin disc with axial forces acting over the disc, which is based on the method of and thus also forms the basis for the present work. However, studied the axis-symmetric flow problem, where we here will stick to a 2-D flow. The reason to focus on the 2-D flow case is that, as we will see later, a complete and simple analytical solution can be derived, which gives valuable and relevant insight into different flow features, whereas for the analysis of the axis-symmetric case had to introduce approximations or assumptions, e.g., that the pressure solution found along the centre line of the disc was assumed to be valid at other radial positions. As we also will see later, the 2-D AD flow has characteristics resembling the 3-D AD flow and thus the potential to form the basis for a wind turbine rotor induction model. Furthermore the simple approach to model a coned rotor with two yawed disc is a big advantage of the 2-D model.

It should be noted that the full 3-D set of equations for the flow around an arbitrary 3-D AS has been presented by and .

### 2.1.1 The full non-linear 2-D solution

The approach for derivation of the descriptive system of equations is based on the Euler equations on differential form and the equation of continuity with the external body forces f specified on the AD. Assuming steady flow the 2-D Euler equations take the following form:

$\begin{array}{}\text{(1)}& {v}_{x}\frac{\partial {v}_{x}}{\partial x}+{v}_{y}\frac{\partial {v}_{x}}{\partial y}=-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial x}+\frac{\mathrm{1}}{\mathit{\rho }}{f}_{x},\text{(2)}& {v}_{x}\frac{\partial {v}_{y}}{\partial x}+{v}_{y}\frac{\partial {v}_{y}}{\partial y}=-\frac{\mathrm{1}}{\mathit{\rho }}\frac{\partial p}{\partial y}+\frac{\mathrm{1}}{\mathit{\rho }}{f}_{y},\end{array}$

and the equation of continuity is as follows:

$\begin{array}{}\text{(3)}& \frac{\partial {v}_{x}}{\partial x}+\frac{\partial {v}_{y}}{\partial y}=\mathrm{0}.\end{array}$

Following the approach by we rewrite the velocity components as

$\begin{array}{}\text{(4)}& {v}_{x}=V+{w}_{x},\text{(5)}& {v}_{y}={w}_{y}.\end{array}$

Here V is the velocity of the main, uniform flow field; wx and wy the induced velocities by the action of the body forces; and vx and vy the final velocity components.

We then non-dimensionalize the velocity components with V, the pressure p with V2ρ, and volume forces fx and fy with V3ρ. However, although we now have non-dimensionalized the equations, the notation of velocity, body force components and pressure will not be changed.

Inserting the new velocity component notation and the non-dimensionalization into Eqs. (1), (2) and (3) and rearranging, we get

$\begin{array}{}\text{(6)}& \frac{\partial {w}_{x}}{\partial x}=-\frac{\partial p}{\partial x}+{f}_{x}-\left({w}_{x}\frac{\partial {w}_{x}}{\partial x}+{w}_{y}\frac{\partial {w}_{x}}{\partial y}\right),\text{(7)}& \frac{\partial {w}_{y}}{\partial x}=-\frac{\partial p}{\partial y}+{f}_{y}-\left({w}_{x}\frac{\partial {w}_{y}}{\partial x}+{w}_{y}\frac{\partial {w}_{y}}{\partial y}\right),\end{array}$

and the continuity equation

$\begin{array}{}\text{(8)}& \frac{\partial {w}_{x}}{\partial x}+\frac{\partial {w}_{y}}{\partial y}=\mathrm{0}.\end{array}$

Following the terminology of , the terms in the brackets on the right-hand side of Eqs. (6) and (7) can be interpreted as induced or second-order forces gx and gy, defined as

$\begin{array}{}\text{(9)}& {g}_{x}=-\left({w}_{x}\frac{\partial {w}_{x}}{\partial x}+{w}_{y}\frac{\partial {w}_{x}}{\partial y}\right),\text{(10)}& {g}_{y}=-\left({w}_{x}\frac{\partial {w}_{y}}{\partial x}+{w}_{y}\frac{\partial {w}_{y}}{\partial y}\right).\end{array}$

Inserting Eqs. (9) and (10) into Eqs. (6) and (7), we get

$\begin{array}{}\text{(11)}& \frac{\partial {w}_{x}}{\partial x}=-\frac{\partial p}{\partial x}+{f}_{x}+{g}_{x},\text{(12)}& \frac{\partial {w}_{y}}{\partial x}=-\frac{\partial p}{\partial y}+{f}_{y}+{g}_{y}.\end{array}$

Now we take the divergence of Eqs. (11) and (12):

$\begin{array}{}\text{(13)}& \frac{{\partial }^{\mathrm{2}}{w}_{x}}{\partial {x}^{\mathrm{2}}}=-\frac{{\partial }^{\mathrm{2}}p}{\partial {x}^{\mathrm{2}}}+\frac{\partial {f}_{x}}{\partial x}+\frac{\partial {g}_{x}}{\partial x},\text{(14)}& \frac{{\partial }^{\mathrm{2}}{w}_{y}}{\partial x\partial y}=-\frac{{\partial }^{\mathrm{2}}p}{\partial x\partial y}+\frac{\partial {f}_{y}}{\partial y}+\frac{\partial {g}_{y}}{\partial y}.\end{array}$

Differentiating the equation of continuity Eq. (8) with respect to x and inserting Eqs. (13) and (14), we get

$\begin{array}{}\text{(15)}& \frac{{\partial }^{\mathrm{2}}p}{\partial {x}^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}p}{\partial {y}^{\mathrm{2}}}=\left(\frac{\partial {f}_{x}}{\partial x}+\frac{\partial {f}_{y}}{\partial y}\right)+\left(\frac{\partial {g}_{x}}{\partial x}+\frac{\partial {g}_{y}}{\partial y}\right).\end{array}$

This is a type of Poisson equation, and with the boundary conditions that the body forces are zero at infinity the solution takes the following form:

$\begin{array}{}\text{(16)}& p\left(f\right)=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\iint \phantom{\rule{0.33em}{0ex}}\frac{{f}_{x}\left(x-\mathit{\xi }\right)+{f}_{y}\left(y-\mathit{\eta }\right)}{\left(x-\mathit{\xi }{\right)}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\xi }\mathrm{d}\mathit{\eta },\end{array}$

and the pressure part from the g terms takes the form

$\begin{array}{}\text{(17)}& p\left(g\right)=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\iint \phantom{\rule{0.33em}{0ex}}\frac{{g}_{x}\left(x-\mathit{\xi }\right)+{g}_{y}\left(y-\mathit{\eta }\right)}{\left(x-\mathit{\xi }{\right)}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\xi }\mathrm{d}\mathit{\eta }.\end{array}$

The pressure for the full solution can then be written as

$\begin{array}{}\text{(18)}& p=p\left(f\right)+p\left(g\right).\end{array}$

This shows that the full solution for the pressure is derived as the sum of two parts: one part, the linear part, being a function of the prescribed forces f on the disc, and the non-linear part, which is a function of the g forces called second-order or induced forces.

When the pressure is derived, the velocity components can be determined by integration of Eqs. (11) and (12).

We will now in the remaining part of the paper focus on the linear solution of the AD flow as this part is much simpler than the full solution. This is because the f forces are prescribed and only different from zero on the disc, whereas the g forces are different from zero in the whole wake region. See for more details.

Furthermore, it is important to notice that the linear solution is still a valid solution satisfying the equations, however, for a new flow problem where also the g forces, derived on the basis of the linear solution, have to be considered as external body forces. The evaluation and interpretation of the linear solution with this in mind is valuable for understanding the flow characteristics of the linear AD flow solution.

This insight into the linear solution is also the basis for the procedure for deriving the full non-linear solution in an iterative manner (see ). In a first iteration the derived g forces from the linear solution are inserted as volume forces but with opposite sign. Then a new solution, which now is a function of both the f and g forces, is derived comprising the solution of the Poisson equation over the area where the g forces are acting. This iterative solution procedure continues until the solutions do not change from iteration to iteration. The technique has been used for the solution of the full non-linear AC flow model (Madsen1983).

Figure 1Notation for the 2-D actuator disc.

## 2.2 The 2Dl solution for the 2-D AD

In the linear solution the flow field is only acted upon by the body forces f on the infinitely thin actuator disc as the induced body forces g are assumed to be zero (see Fig. 1 for notation). The pressure field is then derived from Eq. (15), which now has the following form:

$\begin{array}{}\text{(19)}& \frac{{\partial }^{\mathrm{2}}p}{\partial {x}^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}p}{\partial {y}^{\mathrm{2}}}=\mathrm{0}.\end{array}$

This is the Laplace equation. The pressure p will thus be a potential function outside the disc showing a jump Δp over the disc. Following the derivation by the potential function for the pressure can be obtained by covering the surface with doublets with strength Δp. The pressure jump Δp is linked to the body forces fx as

$\begin{array}{}\text{(20)}& \mathrm{\Delta }p\left(\mathit{\eta }\right)=\underset{S}{\int }\mathrm{d}\mathit{\eta }\underset{\mathit{ϵ}\to \mathrm{0}}{lim}\underset{-\mathit{ϵ}}{\overset{\mathit{ϵ}}{\int }}{f}_{x}\mathrm{d}\mathit{\xi }.\end{array}$

As mentioned above derived the solution for an axis-symmetric disc, but we will derive the pressure potential for a 2-D disc or in fact an actuator strip with a width of 2b and extending from z to z and with doublets with the intensity Δp distributed over the surface S. For any point not situated on the disc the potential will, according to , be given by

$\begin{array}{}\text{(21)}& p=\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\underset{S}{\int }\mathrm{\Delta }p\frac{\partial }{\partial n}\left(\frac{\mathrm{1}}{\mathit{\omega }}\right)\mathrm{d}S,\end{array}$

where n denotes the positive normal to the disc and ω the distance from any point on the disc to the point A where the potential is calculated. The details of the derivation can be found in Appendix A.

The resulting pressure function has been derived to

$\begin{array}{}\text{(22)}& p=-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{\Delta }p\frac{x}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta }.\end{array}$

where x and y have been non-dimensionalized with the half width b of the disc.

The induced axial velocity wx can now be found by integrating Eq. (11) from −∞ to x:

$\begin{array}{}\text{(23)}& \begin{array}{rl}{w}_{x}& =-p+\underset{\mathrm{\infty }}{\overset{x}{\int }}{f}_{x}\mathrm{d}\mathit{\xi }=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{\Delta }p\frac{x}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta }\\ & +\mathrm{\Delta }p\phantom{\rule{0.25em}{0ex}}\left(\text{last term only in the wake},\phantom{\rule{0.25em}{0ex}}x\phantom{\rule{0.125em}{0ex}}>\phantom{\rule{0.125em}{0ex}}\mathrm{0}\right).\end{array}\end{array}$

The wy component is determined through Eq. (12) with fy=0. It means first differentiating p with respect to y and then integrating $\partial p/\partial y$ from −∞ to x. The derivation found in Appendix B leads to the following equation for wy:

$\begin{array}{}\text{(24)}& {w}_{y}=\underset{\mathrm{\infty }}{\overset{x}{\int }}\frac{\partial p}{\partial y}\mathrm{d}\mathit{\xi }=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{\Delta }p\frac{y-\mathit{\eta }}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta }.\end{array}$

With a constant loading Δp on the disc, the pressure and velocity components can be derived by integrating Eqs. (22), (23) and (24) over the width of the disc. The derivation of the pressure function shown in Appendix C leads to

$\begin{array}{}\text{(25)}& p=-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\mathrm{\Delta }p\left[\mathrm{arctan}\frac{\left(\mathrm{1}-y\right)}{x}+\mathrm{arctan}\frac{\left(\mathrm{1}+y\right)}{x}\right],\end{array}$

and then inserting the pressure equation into Eq. (23), we get

$\begin{array}{}\text{(26)}& \begin{array}{rl}& {w}_{x}=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\mathrm{\Delta }p\left[\mathrm{arctan}\frac{\left(\mathrm{1}-y\right)}{x}+\mathrm{arctan}\frac{\left(\mathrm{1}+y\right)}{x}\right]+\mathrm{\Delta }p\\ & \left(\text{last term only in the wake},x\phantom{\rule{0.125em}{0ex}}>\phantom{\rule{0.125em}{0ex}}\mathrm{0}\right).\end{array}\end{array}$

Likewise the integration of wy for constant loading, shown in Appendix D, yields

$\begin{array}{}\text{(27)}& {w}_{y}=\frac{\mathrm{\Delta }p}{\mathrm{4}\mathit{\pi }}\mathrm{ln}\left(\frac{{x}^{\mathrm{2}}+\left(y+\mathrm{1}{\right)}^{\mathrm{2}}}{{x}^{\mathrm{2}}+\left(y-\mathrm{1}{\right)}^{\mathrm{2}}}\right).\end{array}$

Finally, according to Eqs. (4) and (5) the vx and vy velocity components are derived as

$\begin{array}{}\text{(28)}& {v}_{x}=\mathrm{1}+{w}_{x},\text{(29)}& {v}_{y}={w}_{y}.\end{array}$

In this section we will study the flow characteristics through the 2-D actuator disc by numerical exploration based on the analytical linear solution, 2Dl, for a constant, uniform loading of the disc.

First we relate the pressure jump, Δp, to the thrust, T, and thrust coefficient, CT:

$\begin{array}{}\text{(30)}& T=\mathrm{2}b\mathrm{\Delta }p,\end{array}$

where 2b is the width of the actuator disc.

From the definition of the thrust coefficient, CT, we get

$\begin{array}{}\text{(31)}& {C}_{\mathrm{T}}=\frac{T}{\mathrm{2}b\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{V}^{\mathrm{2}}}.\end{array}$

Combining Eqs. (30) and (31) and keeping in mind that Δp is non-dimensionalized with ρV2, we get

$\begin{array}{}\text{(32)}& {C}_{\mathrm{T}}=\mathrm{2}\mathrm{\Delta }p.\end{array}$

## 3.1 Velocity profiles vx and vy

For a loading corresponding to a thrust coefficient CT=0.4 the vx and vy velocity profiles as a function of y are shown for different upstream and downstream positions in Fig. 2. The vx velocity profiles are seen to develop from a Gaussian-type profile upstream of the disc, to be sharper when approaching the disc, and finally at the disc to be completely rectangular with free-stream velocity just outside the disc. Accelerated flow above free-stream velocity is then seen to occur just downstream of the disc and peaking before x=2. Further downstream in the wake the velocity profiles become almost rectangular at x=10 and with a deficit within the projected area of the disc, which means a non-expanding wake.

Figure 2(a) The axial velocity component vx is shown as a function of the lateral coordinate y for different streamwise positions x from upstream of the disc to far downstream in the wake for a thrust coefficient of 0.4. The velocity profiles for the momentum solution are also shown at the disc and in the far wake. (b) The lateral velocity component vy is shown for the same positions.

The vy profiles are seen to be anti-symmetric around the x axis as expected and have a strong peak at the edge of the disc. Downstream the profiles decay fast.

The momentum solution for the vx component, also shown in Fig. 2, is seen to be lower than the actual linear solution, which will be explored next.

## 3.2 Scaling CT to fit the 2Dl solution to momentum theory

As seen above the 2Dl solution does not fit exactly with the classical momentum theory solution for an actuator disc. The computed induction is lower, and it is an important shortcoming for the wider use of the model as a rotor induction model.

To investigate that further we derive the relation between the thrust coefficient CTl of the 2Dl model and the induction coefficient, al, for the 2Dl solution presented above.

The axial velocity vxd at the disc is derived from Eq. (26) by $x\to -\mathrm{0}$,

$\begin{array}{}\text{(33)}& {w}_{xd}=\frac{\mathrm{\Delta }p}{\mathrm{2}\mathit{\pi }}\left(\frac{\mathit{\pi }}{\mathrm{2}}+\frac{\mathit{\pi }}{\mathrm{2}}\right)=\frac{\mathrm{\Delta }p}{\mathrm{2}},\end{array}$

and from the definition of a,

$\begin{array}{}\text{(34)}& {a}_{l}=\mathrm{1}-{v}_{xd}=\mathrm{1}-\left(\mathrm{1}-{w}_{xd}\right),\text{(35)}& {a}_{l}={w}_{xd}=\frac{\mathrm{\Delta }p}{\mathrm{2}}=\frac{{C}_{\mathrm{Tl}}}{\mathrm{4}},\end{array}$

or

$\begin{array}{}\text{(36)}& {C}_{\mathrm{Tl}}=\mathrm{4}{a}_{l}.\end{array}$

The classical momentum theory derived as an axis-symmetric solution for a circular actuator disc gives the following relation between the thrust coefficient, CTm, and the induction, am (the subscript m denotes momentum theory):

$\begin{array}{}\text{(37)}& {C}_{\mathrm{Tm}}=\mathrm{4}{a}_{\mathrm{m}}\left(\mathrm{1}-{a}_{\mathrm{m}}\right).\end{array}$

Comparing the solution from the two models, the momentum theory gives a velocity at the disc, vxm=0.887, for a CT of 0.4, while the 2Dl model gives vx=0.9 for a similar CT.

The reduced induction by the 2Dl model is due to the action of the induced body forces defined in Eqs. (9) and (10) that accelerate the flow in the wake region and also through the disc and upstream of the disc through the pressure field from the g forces.

Without a detailed derivation of the g forces, we can conclude that there will be gx forces along the lines $y=±\mathrm{1}$ in the wake region, which is responsible for the acceleration of the flow through the disc, requiring increased external volume forces, fx (increased CTl), in order to correlate with the velocity at the disc from the momentum solution. A check of the x-momentum balance of a control volume around the disc at some distance will show that the balance is only fulfilled if the g forces are considered as external body forces. This has been demonstrated for the AC flow models based on the same solution approach (Madsen1982).

The analogy to use a fixed wake vortex model for the same flow problem would be a vortex sheet along the lines $y=±\mathrm{1}$. This vortex sheet will require external body forces to be fixed (non-expanding wake) because the flow in the vicinity of the disc is not aligned with the vortex sheet and thus induces forces on the sheet.

In order to force the 2Dl model to give the same induction in the rotor plane and in the far-field wake as the momentum theory, we can simply scale the thrust coefficient based on Eqs. (36) and (37):

$\begin{array}{}\text{(38)}& \frac{{C}_{\mathrm{Tl}}}{{C}_{\mathrm{Tm}}}=\frac{\mathrm{4}{a}_{l}}{\mathrm{4}{a}_{\mathrm{m}}\left(\mathrm{1}-{a}_{\mathrm{m}}\right)},\end{array}$

and enforcing al to be equal to am, we get

$\begin{array}{}\text{(39)}& \begin{array}{rl}{C}_{\mathrm{Tl}}& ={C}_{{\mathrm{T}}_{\mathrm{m}}}\frac{\mathrm{4}{a}_{m}}{\mathrm{4}{a}_{\mathrm{m}}\left(\mathrm{1}-{a}_{\mathrm{m}}\right)}={C}_{{\mathrm{T}}_{\mathrm{m}}}\frac{\mathrm{1}}{\mathrm{1}-{a}_{\mathrm{m}}}\\ & ={C}_{\mathrm{Tm}}\frac{\mathrm{2}}{\mathrm{1}+\sqrt{\mathrm{1}-{C}_{{\mathrm{T}}_{\mathrm{m}}}}}.\end{array}\end{array}$

(further down we will denote CTm as CT).

For example for a thrust coefficient of 0.89 we shall increase CTl to 1.34 to match the momentum solution as shown in the left graph in Fig. 3. As mentioned above this considerable increase in CT is necessary due to the strong impact of the g forces on the wake flow, which is seen clearly in the right graph in Fig. 3. Here the vx velocity profiles together with the enclosing streamlines of the disc show that part of the flow passing through the disc and thus de-accelerated due to the f forces at the disc plane is accelerated up to free-stream velocity at a short distance downstream in the wake region. Although the g forces are only acting from the disc plane and downstream along the lines $y=±\mathrm{1}$, the pressure field from these forces will also influence the flow upstream of the disc.

Figure 3(a) The axial velocity component vx is shown as a function of the lateral coordinate y for different streamwise positions x from upstream of the disc to far downstream in the wake for a thrust coefficient of 1.34 in the 2Dl model and compared with the momentum solution for a CT of 0.89. (b) The enclosing streamlines of the disc for the same loading are shown together with vx profiles at different streamwise positions. The streamlines are derived on the basis of the local flow angle $\mathrm{arctan}\left(\frac{{v}_{y}}{{v}_{x}}\right)$ with starting points at the edge of the disc.

It can finally be noticed that the same scaling procedure for correcting a linear solution is used in the actuator cylinder model for the VAWT's implemented in HAWC2. However, in this case the correction is applied directly as the scaling of the velocity components .

## 3.3 Superposition of solutions for more discs

As we have a linear model type we can easily study the interaction of more discs with, e.g., different loadings or locations to each other. Later we will show that the superposition is one way to model a coned disc by the superposition of two yawed discs.

Let a new disc 2 have the origin xD2,yD2 and a loading ΔpD2 in the original coordinate system as shown in Fig. 1. We then get the contributions vxD2 and vyD2 to the velocity field by inserting them into Eqs. (26) and (27):

$\begin{array}{}\text{(40)}& \begin{array}{rl}{w}_{x\mathrm{D}\mathrm{2}}& =\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\mathrm{\Delta }{p}_{\mathrm{D}\mathrm{2}}\left[\mathrm{arctan}\frac{\mathrm{1}-\left(y-{y}_{\mathrm{D}\mathrm{2}}\right)}{x-{x}_{\mathrm{D}\mathrm{2}}}\\ & +\mathrm{arctan}\frac{\mathrm{1}+\left(y-{y}_{\mathrm{D}\mathrm{2}}\right)}{x-{x}_{\mathrm{D}\mathrm{2}}}\right]+\mathrm{\Delta }{p}_{\mathrm{D}\mathrm{2}}\\ & \text{(last term only in the wake of disc 2)},\end{array}\text{(41)}& {w}_{yD\mathrm{2}}=\frac{\mathrm{\Delta }{p}_{\mathrm{D}\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\mathrm{ln}\left(\frac{\left(x-{x}_{\mathrm{D}\mathrm{2}}{\right)}^{\mathrm{2}}+\left(\left(y-{y}_{\mathrm{D}\mathrm{2}}\right)+\mathrm{1}{\right)}^{\mathrm{2}}}{\left(x-{x}_{\mathrm{D}\mathrm{2}}{\right)}^{\mathrm{2}}+\left(\left(y-{y}_{\mathrm{D}\mathrm{2}}\right)-\mathrm{1}{\right)}^{\mathrm{2}}}\right).\end{array}$

However, the set-up does not work for overlapping discs.

Here the superposition will first be used to study an important characteristic of the model for discs at the same streamwise position with different loadings as in the example in the left part of Fig. 4. The discs are positioned so close that they form one disc with double width, 4b, and with a thrust coefficient of 0.89 and 0.45 on each half. Upstream and downstream of the discs there is a clear interaction between the vx velocity profiles of the discs, but exactly in the plane of the discs there is no interaction of the axial velocity component, although the thrust coefficient jumps from 0.89 on the lower part to half, 0.445, on the upper part.

Figure 4In (a) we show the flow characteristics of two discs placed side by side so they just touch. The lower disc has the double loading of the upper disc. Enclosing streamlines for the two discs are shown together with vx velocity profiles at different axial positions. (b) The same type of data is shown, however, for the same loading on the two discs and with disc 2 shifted a distance of half the disc width in both x and y directions.

This can be used to guide the implementation of the BEM model in aerodynamic and aeroelastic codes as proposed by using a VC model for a similar case. The result justifies a local, point-wise evaluation of the momentum balance, e.g., used in the implementation of the BEM model in the HAWC2 code . Finally, it should be noticed that there is an interaction of the lateral velocity component in the disc plane.

When the two discs do not have the same streamwise position, we see an interaction of the vx profiles, as in the case for the example in the right part of Fig. 4. A slightly increased induction can be seen on the upper part of disc 1. Furthermore, an accelerated axial flow velocity between the two discs is noted.

## 3.4 An important characteristic of the velocity profiles of the 2Dl model at a small distance downstream of the disc and comparison with CFD simulations

In Fig. 3 we saw that at the disc and far downstream in the wake the vx profiles were rectangular, but in between some bending or curved profiles could be seen, e.g., at x=1 and x=2. We will now go much closer to the disc where we find the development of the profiles as shown in the left graph in Fig. 5. It is clearly seen that the vx velocity profiles immediately develop from a uniform distribution to shapes with stronger deceleration close to $y=±\mathrm{1}$. The same is typically seen in full, non-linear simulations of AD flow .

Figure 5In (a) we show the axial velocity vx profiles at positions close behind the disc as well as further downstream for a thrust coefficient of 0.89. For comparison the momentum solution for vx at the disc and in the far wake is also included. (b) The vx velocity profile at the disc from an axis-symmetric CFD AD simulation is compared with the 2Dl profile at a downstream distance of 0.08 scaled by a factor of 1.05 to account for the decrease in velocity downstream behind the disc and to give the best fit to the CFD results extracted at the disc.

We use here the results of from an axis-symmetric simulation of an AD. Comparing the axial velocity profile at the disc from the CFD simulation with profiles from the 2Dl model extracted at different downstream positions, we find the best correlation in shape using the 2Dl profiles at x=0.08 and scaled by a factor of 1.05, as shown in the right graph of Fig. 5.

Figure 6(a) A zoom of the right graph in Fig. 5. (b) 2Dl vy velocity profiles are compared with the radial velocity profile from an axis-symmetric simulation. At the same downstream position x=0.08 as for the vx component, the best fit for vy is obtained by a scaling factor of 0.67.

The reason for this development of the vx profiles close behind the disc is probably that in this region the g forces along the lines $y=±\mathrm{1}$ are acting both upstream and downstream of the axial position where we extract the velocity profiles. Therefore, there might be some local cancelling of their influence so that the deficit is closer to the full solution. The left graph of Fig. 6 shows a very good correlation with the CFD solution when applying the above scaling. It is expected that such scaling is necessary as we have computed the velocity profile a short distance downstream of the disc where the axial velocity decays rapidly. For a wider application of the 2Dl model it is important that the scaling factor is constant. Therefore, we will use the same scaling in the remaining flow examples in the paper.

Finally, the vy profiles are compared with the profiles of the radial velocity component, vr, from CFD simulations in the right graph of Fig. 6. When vy is computed at the same position, x=0.08, as for vx the best fit is obtained by scaling down by a factor of 0.67. The need for a strong down-scaling is expected as a 2-D solution will have a much stronger lateral velocity component than the radial velocity in an axis-symmetric solution.

4 The 2Dl solution for yawed flow and comparisons with 3-D CFD AD computations with constant loading

We continue now to derive the analytical solution for a 2-D AD in yawed flow (see Fig. 7). This is an important flow situation for a turbine which causes increased loading and reduced power. However, yaw is now frequently studied and used for wake steering in wind farms to increase the power of downstream turbines by upstream wake deflection.

Figure 7Notation for the 2-D actuator disc in yaw.

The objective with the present model derivation is the further study of the capability of the simple model to handle complex flow like a yawed AD. A further objective is, as the next step (next section), to combine two yawed discs with opposite yaw angles to model a coned rotor.

The starting point is the flow solution for wx and wy in Eqs. (26) and (27) in the non-yawed x,y coordinate system which then is rotated by the yaw angle, Θ, into the new system, ${x}^{\prime },{y}^{\prime }$ (see Fig. 7). The velocity components, wx and wy, in the original system, x,y, are then derived as a function of ${x}^{\prime },{y}^{\prime }$.

The following coordinate transformation is introduced:

$\begin{array}{}\text{(42)}& {x}^{\prime }=\left(x\mathrm{cos}\mathrm{\Theta }-y\mathrm{sin}\mathrm{\Theta }\right),\text{(43)}& {y}^{\prime }=\left(x\mathrm{sin}\mathrm{\Theta }+y\mathrm{cos}\mathrm{\Theta }\right).\end{array}$

The coordinate transformations are inserted into Eqs. (26) and (27):

and

$\begin{array}{}\text{(45)}& {w}_{y}=\frac{\mathrm{\Delta }p}{\mathrm{4}\mathit{\pi }}\mathrm{ln}\left(\frac{\begin{array}{c}\left(x\mathrm{cos}\mathrm{\Theta }-y\mathrm{sin}\mathrm{\Theta }{\right)}^{\mathrm{2}}+\\ \left(\left(x\mathrm{sin}\mathrm{\Theta }+y\mathrm{cos}\mathrm{\Theta }\right)+\mathrm{1}{\right)}^{\mathrm{2}}\end{array}}{\begin{array}{c}\left(x\mathrm{cos}\mathrm{\Theta }-y\mathrm{sin}\mathrm{\Theta }{\right)}^{\mathrm{2}}+\\ \left(\left(x\mathrm{sin}\mathrm{\Theta }+y\mathrm{cos}\mathrm{\Theta }\right)-\mathrm{1}{\right)}^{\mathrm{2}}\end{array}}\right).\end{array}$

The total velocity components, vx and vy, for the yawed disc in the x,y coordinate system are then derived as

$\begin{array}{}\text{(46)}& \begin{array}{rl}& {v}_{x}=\mathrm{1}+{w}_{x}-\mathrm{\Delta }p\phantom{\rule{0.25em}{0ex}}\left(\text{last term only in the wake of the}\\ & \text{yawed disc:}\phantom{\rule{0.25em}{0ex}}{x}^{\prime }>\mathrm{0}\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}y>-\mathrm{cos}\mathrm{\Theta }\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}y<\mathrm{cos}\mathrm{\Theta }\right),\\ & {v}_{y}={w}_{y}.\end{array}\end{array}$

## 4.1 Streamlines and velocity profiles for a yawed disc of 30, 45 and 60∘ in comparison with 3-D CFD computations

We will now explore the flow field of a yawed disc of 30, 45 and 60 by presenting enclosing stream lines of the disc and axial velocity profiles from upstream to downstream positions (left graph in Figs. 8, 9 and 10). Further, a comparison with 3-D CFD AD simulations by and is shown in the graphs to the right.

Figure 8(a) The enclosing streamlines of the yawed disc plus the streamline through the centre of the disc for a yaw of 30. Further, vx profiles are shown at three x positions. We show also the line yd along the disc at a distance of Δx=0.08, which is the line along which the normal velocity to the disc is shown in (b) in comparison with 3-D CFD simulation of the AD.

Figure 9Similar graphs as in Fig. 8 above but for a yaw of 45.

We compare the normal velocity to the disc vn along a line yd at a distance of Δx=0.08 (see the streamline plots where the yd line is included), which was the distance we found above to give the best correlation with CFD for aligned flow. The normal velocity vn to the disc is important as it is the velocity component in a rotor simulation with a major influence on the angle of attack (AoA).

We also show vn along other lines parallel to the disc at other distances Δx, where one is at the disc and one is downstream at Δx=0.2. This is to illustrate the variation in the velocity profiles close behind the disc.

Starting with the evaluation of the results for the yaw case of 30 in Fig. 8 we find overall a very good correlation for the vn profiles with considerable higher induction on the downstream edge of the disc. Also the non-linear shape in this region fits well with the CFD results. Similar non-linear velocity profiles are seen in the CFD simulations by and by the skewed cylindrical vortex wake model by .

Typically in engineering yaw modelling in aerodynamic and aeroelastic codes this velocity distribution is approximated with a linear correlation using the model as, e.g., in the implementation of the HAWC2 code . We further notice the acceleration of vn on the upstream edge of the disc and outside, which is stronger than in the CFD results. This can mainly be ascribed to the difference between a 2-D and 3-D flow solution with higher lateral velocities in the 2-D flow.

Continuing with the 45 yaw case in Fig. 9 there is still overall a good correlation between the 2Dl model results and the CFD computations. As above the correlation is almost perfect on the downstream half of the disc, whereas the deviations now are slightly increased on the other half where the 2Dl model computes less induction than the CFD results. Again this is probably due to the 2-D modelling of the flow that results in a stronger acceleration of the lateral flow, in particular outside the disc but now also influencing part of the flow through the disc close to the edge of the disc.

Finally in the 60 yaw case in Fig. 10 we now see major deviations, although the slope of the vn curves is comparable. Overall the 2Dl model computes less induction over the whole disc, and this might be due to the scaling of CT based on aligned flow where the wake with the g forces now is offset considerably from one side of the disc to the other.

Figure 10Similar graphs as in Fig. 8 above but for a yaw of 60.

Based on the above comparison cases with 3-D CFD results we can conclude that the 2Dl model simulates yawed flow cases up to 3040 with quite good accuracy and includes details, like the bending of the vn curves, not modelled by common engineering yaw models.

5 Modelling of coned rotors with the 2Dl model

Modelling of coned rotors has been of great importance for decades in the wind research community as turbine concept studies often have included downwind coned rotors. This is due to potential advantages by such concepts by alleviating the flapwise blade root moment due to the deloading by centrifugal forces on a coned rotor or blade. In fact the CFD computations on coned rotors , to be used here for comparisons, were carried out historically to provide input for the design of a small two-bladed rotor that could be coned up to 90 for surviving storm situations.

Recent concept studies that include downwind coned rotors are, e.g., , and . These studies still rely on BEM-type simulation tools which do not model any impact on induction from the coning, which is in contrast to results with higher-fidelity models like the CFD AD simulations presented here.

## 5.1 Modelling a coned disc by superposition of two yawed discs

The modelling of the coned rotor is based on the superposition of two yawed discs with an opposite yaw angle of ±30, as in the example shown in the left graph of Fig. 11. In this case there is still a distance between the edges of the discs, and we can notice an accelerated flow between the discs.

Figure 11In (a) we show two separated actuator discs with a yaw angle of 30 and −30. In (b) they are moved towards each other so the edges just coincide and form one coned disc of 30 coning.

However, when we now move the discs so close that the two edges just coincide, as shown in the right graph of Fig. 11, we see one velocity deficit behind the two yawed discs and we can call it one coned disc with a cone angle of 30 and the double size of the two yawed discs.

## 5.2 Comparison of the 2Dl model with an axis-symmetric CFD simulation for a coned disc of 20 and −20∘

Next we compare the 2Dl solution for a coned rotor with axis-symmetric CFD simulations, as shown in the right graph of Fig. 12. As for the yawed disc we compare vn as this is the most important velocity component in a rotor simulation with direct influence on angle of attack.

Figure 12In (a) we show the enclosing streamlines of a coned disc with 20 coning and a CT=0.89. In the same graph we also show the yd lines parallel with the two disc parts at a distance of $\mathrm{\Delta }x=\mathrm{2}×\mathrm{0.08}$, which are the lines along which we extract the velocities in the 2Dl model for comparison with the CFD results in (b). However, the CFD velocities were extracted at the disc.

To be consistent with the distance behind the disc where we extract the 2Dl velocities, as tuned for aligned flow, we use the same relative distance behind the disc, which means a double absolute distance $\mathrm{\Delta }x=\mathrm{2}\cdot \mathrm{0.08}$ as the coned disc has double the size of the aligned and yawed disc.

Comparing now the 2Dl simulations with the CFD results we see a very close correlation and just a small underestimation of vn when using the data extracted at $\mathrm{\Delta }x=\mathrm{2}\cdot \mathrm{0.08}$. An almost complete coincidence is seen for the data extracted slightly closer to the disc at $\mathrm{\Delta }x=\mathrm{2}\cdot \mathrm{0.04}$.

In Fig. 13 we compare next the 2Dl simulations with CFD results for a cone angle of −20. First of all we notice that the shape of the curves are quite different from the velocity profiles for downwind coning as shown above as there is a decrease in the velocity close to the edge. However, the 2Dl model predicts this quite well. The difference occurs when we apply the scaling of the velocities as it shifts the profile slightly towards higher velocities than the CFD results show.

Figure 13Same type of data as in Fig. 12 but for a cone angle of −20.

## 5.3 On the complexity of the flow solution for a coned disc

As mentioned in the Introduction (Sect. 1), has proposed a correction to BEM for improving the modelling of a coned rotor. To illustrate with the present model how complex the induction is in a coned disc and thus how challenging a single-point correction of BEM like the one of can be, we show in Fig. 14 the vx and vy contributions from both sides of the coned disc which then sum up to the total component. Disc 1 is the lower part of the downwind coned rotor, and in the left graph in Fig. 14 the green line is thus the vx induction from disc 1 itself, whereas the red line is the contribution from the other half, disc 2 of the coned rotor. This contribution is seen to be accelerated flow towards the centre of the coned disc. Likewise we show the same for the vy component in the graph to the right in Fig. 14. It is seen that the total component is about 0.4 at the tip and therefore contributes significantly to vn. This illustrates the importance of including the radial component for the derivation of the induction in a coned rotor.

Figure 14In (a) we show the vx components contributing to the total induction on the lower part, i.e. disc 1. The same for vy in (b).

## 5.4 Flow solution for a range of cone angles

With the 2Dl model we have simulated the normal velocity, vn, for a range of cone angles from 60 to −60 for a thrust coefficient of 0.89.

Figure 15Simulated vn distributions with the 2Dl model for a range of positive cone angles in (a) and for negative cone angles (b).

In Fig. 15, in the graph to the left with the vn result for positive cone angles, the shape of the profiles becomes more slender with an increasing difference in velocity at the centre and at the edge. At low cone angles the increase in the normal velocity at the central part of the disc above the velocity for a plane disc can be noticed.

For the negative cone angles in the right graph of Fig. 15 the characteristic sweep on the outer part of the velocity profile towards higher induction can be seen up to a cone angle of about −40.

## 5.5 Coned disc in yaw

Finally, we have simulated a coned disc of 30 in 20 yawed flow to explore the capability of the 2Dl model to simulate such complex inflow. In the right graph in Fig. 16 we show the enclosing streamlines of the disc. A considerable downward deflection of the wake flow is seen. It means that a bigger part of the flow passing through the lower part of the disc is deflected outside the projected area of the disc and thus accelerated up to the free-stream velocity due to the induced g forces.

In the right graph of Fig. 16 the normal velocity to the disc along the line yd parallel to the disc at a distance of $\mathrm{\Delta }x=\mathrm{2}×\mathrm{0.08}$ shows as expected a jump passing through the centre of the disc. Further it can be seen that the induction on the two sides of the yawed disc are quite different with a stronger variation from edge to centre on the downwind part.

Figure 16In (a) we show the enclosing streamlines of a 30 coned disc with a yaw angle of 20 together with axial velocity profiles at different streamwise positions. In (b) the normal velocity to the disc along the coned and yawed disc is shown in comparison with the velocities for an aligned coned disc.

6 Outlook

The promising results of the 2Dl model based on comparisons with 3-D and axis-symmetric CFD simulations indicate that although the model is 2-D, it could form the basis for a consistent, general rotor induction model for wind turbine rotors with one equation for vx and one for vy including modelling of yawed and coned disc flow. This simplicity would also be of big advantage for the model to be used in rotor optimizations in the derivation of analytical gradients.

One prerequisite for the model to be used for a rotor induction model is a coupling to an angular momentum model for tangential induction. Following the model approach by this can be done through the following equation derived from angular momentum balance over the actuator disc combined with a blade element analysis:

$\begin{array}{}\text{(47)}& {v}_{\mathrm{t}}=\frac{{V}_{\mathrm{0}}{C}_{Q}}{\mathrm{2}\left(\mathrm{1}-a\right)},\end{array}$

where vt is the tangential-induced velocity and CQ is the non-dimensional tangential load coefficient defined as

$\begin{array}{}\text{(48)}& {C}_{Q}=\frac{{V}_{\mathrm{r}}^{\mathrm{2}}{C}_{x}c{N}_{\mathrm{B}}}{{V}_{\mathrm{0}}^{\mathrm{2}}\mathrm{2}\mathit{\pi }r}\end{array}$

and

$\begin{array}{}\text{(49)}& {C}_{x}={C}_{l}\mathrm{sin}\mathit{\phi }-{C}_{\mathrm{d}}\mathrm{cos}\mathit{\phi }.\end{array}$

Vr is the relative velocity to the blade section, c local chord, NB number of blades, φ the local flow angle from the rotor plane to Vr, Cl and Cd the local airfoil coefficients, and a the axial induction factor which is computed with the 2Dl model.

The specific implementation of the tangential induction model can be discussed, e.g., if it is based on local values or probably best as the average CQ and a over a ring element.

Likewise the input to the 2Dl model in the form of CT has to be derived from the blade element analysis and takes a similar form as for CQ :

$\begin{array}{}\text{(50)}& {C}_{\mathrm{T}}=\frac{{V}_{\mathrm{r}}^{\mathrm{2}}{C}_{y}c{N}_{\mathrm{B}}}{{V}_{\mathrm{0}}^{\mathrm{2}}\mathrm{2}\mathit{\pi }r},\end{array}$

with

$\begin{array}{}\text{(51)}& {C}_{y}={C}_{l}\mathrm{cos}\mathit{\phi }+{C}_{\mathrm{d}}\mathrm{sin}\mathit{\phi }.\end{array}$

Further on the implementation of the 2Dl model in a rotor analysis, the induction should be computed along diagonal lines of the rotor using the local yaw and tilt angle and cone angle in the 2Dl model. Such an implementation would fit well to the polar grid implementation of BEM in HAWC2 as described in . The approach could be the computation of the induction along 12–16 diagonals at equally spaced azimuth positions which would give a good resolution of azimuthal variations due to non-uniform inflow and yaw and tilt angles.

It would also be necessary to relax the constant loading assumed for the integration across the disc of Eqs. (23) and (24). However, these integrals can be integrated numerically without slowing down the computations too much or as in the implementation of the actuator cylinder model where the loading is assumed to be piecewise constant which allows the integration of influence coefficients only once in a time simulation.

The capability of the 2Dl model to simulate the flow through a coned rotor in yaw by using two yawed discs of different yaw angles with opposite sign is a big advantage of the model. This allows also for a straightforward simulation of a coned disc in yaw by yawing the one disc more than the other.

The superposition of several differently yawed discs with constant loading could also be used to approximate an arbitrary rotor shape and loading by a segmented diagonal line of discs just coinciding at the edges. This has some resemblance to the approach by where a staggering of VC models shifted in the streamwise direction is used to model non-planar rotors.

A final comment on the 2Dl model is that the evaluation of the velocities at a distance 5 %–10 % behind the disc means that the induction even for a plane rotor is not completely local as there would be velocity contributions from nearby points along the diagonal line. This is considered as an advantage as test cases with high loading in the aeroelastic code HAWC2 using the local induction approach indicate that some averaging of the local CT should be used.

7 Conclusions

We have presented an analytical linear solution for a 2-D AD flow that is derived from the Euler equations and the equation of continuity. The induction at the disc is uniform for a constant loading like in the momentum theory, but we have shown that immediately behind the disc the axial velocity profiles are non-linear and match by simple scaling a full non-linear axis-symmetric flow solution for an AD.

By a simple coordinate rotation the analytical solution was extended to a yawed disc. This flow solution matches well a full 3-D CFD simulation of a circular AD when comparing the normal velocity to the disc in the plane with skewed inflow up to 30–40.

The superposition of more discs enabled studying turbine interaction and is also used to form a coned disc by positioning two yawed discs with opposite yaw angles so the edges of the discs just coincide. An excellent agreement for a downwind coned disc of 20 was found when comparing the normal velocity to the disc with an axis-symmetric CFD simulation of a coned circular disc. For a disc with upwind coning of −20 the typical shape of the velocity profile with increased induction on the most outboard part matched well the CFD simulations, but some offset of the velocity profile was found.

The promising results of comparing the 2Dl model with 3-D and axis-symmetric CFD simulations on a circular AD have indicated that the 2Dl model could form the basis for a consistent, simple induction model for aerodynamic and aeroelastic wind turbine rotor simulations. First of all the model should for such an application be combined with an angular momentum model which has been sketched in the paper. Then the computation of induction should be used along diagonal lines of the rotor in 12–16 different azimuth positions where the skewed inflow angle from yaw error or tilt is input to the model. Finally, an extension of the present implementation to simulate a variable loading across the disc could be done by numerical integration of the equations for vx and vy or by superposition of several yawed discs, each with a constant loading and yaw angle forming a segmented line approximating the real diagonal line of the rotor shape.

Appendix A: Derivation of the pressure function

We show here the derivation of pressure potential for a 2-D disc or in fact an actuator strip extending from z to z and with doublets with the intensity Δp distributed over the surface S. For any point not situated on the disc the potential will according to be given by

$\begin{array}{}\text{(A1)}& p=\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\underset{S}{\int }\mathrm{\Delta }p\frac{\partial }{\partial n}\left(\frac{\mathrm{1}}{\mathit{\omega }}\right)\mathrm{d}S,\end{array}$

where n denotes the positive normal to the disc and ω the distance from any point on the disc to the point A where the potential is calculated. ω will then be

$\begin{array}{}\text{(A2)}& \mathit{\omega }=\sqrt{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}},\end{array}$

and inserting into Eq. (21), we get

$\begin{array}{}\text{(A3)}& \begin{array}{rl}p& =\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\underset{-b}{\overset{{}^{b}}{\int }}\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{\partial }{\partial x}\left(\frac{\mathrm{1}}{\sqrt{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}}}\right)\mathrm{d}\mathit{\zeta }\\ & =\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\underset{-b}{\overset{{}^{b}}{\int }}\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{-\mathrm{0.5}\cdot \mathrm{2}x}{\left({x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}+{\mathit{\zeta }}^{\mathrm{2}}{\right)}^{\mathrm{3}/\mathrm{2}}}\mathrm{d}\mathit{\zeta }.\end{array}\end{array}$

The integral is of the type (207) in the handbook of :

$\begin{array}{}\text{(A4)}& \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{\mathrm{d}t}{T\sqrt{T}}=\left[\frac{\mathrm{2}\left(\mathrm{2}ct+b\right)}{q\sqrt{T}}{|}_{-\mathrm{\infty }}^{\mathrm{\infty }}\right],\end{array}$

where $T=a+bt+c{t}^{\mathrm{2}}$ and $q=\mathrm{4}ac-{b}^{\mathrm{2}}$. Now,

$\begin{array}{ll}& a={x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}},\\ & b=\mathrm{0},\\ & c=\mathrm{1},\\ & q=\mathrm{4}\left({x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}\right)\right),\\ & T={x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}+{t}^{\mathrm{2}}.\end{array}$

We integrate Eq. (A3):

$\begin{array}{}\text{(A5)}& p=-\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\underset{-b}{\overset{{}^{b}}{\int }}\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\phantom{\rule{0.33em}{0ex}}\left[x\frac{\mathrm{2}\left(\mathrm{2}ct+b\right)}{q\sqrt{T}}{|}_{-\mathrm{\infty }}^{\mathrm{\infty }}\right],\text{(A6)}& p=-\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\underset{-b}{\overset{{}^{b}}{\int }}\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\phantom{\rule{0.33em}{0ex}}\left[x\frac{\mathrm{4}t}{\begin{array}{c}\mathrm{4}\left[{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}\right]\\ \sqrt{\left(x-\mathit{\xi }{\right)}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}+{t}^{\mathrm{2}}}\end{array}}{|}_{-\mathrm{\infty }}^{\mathrm{\infty }}\right],\text{(A7)}& \begin{array}{rl}p=& -\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\underset{-b}{\overset{{}^{b}}{\int }}\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\phantom{\rule{0.33em}{0ex}}\left[\frac{x}{\left[{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}\right]}\\ & +\frac{x}{\left[\left(x-\mathit{\xi }{\right)}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}\right]}\right],\end{array}\text{(A8)}& p=-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-b}{\overset{{}^{b}}{\int }}\mathrm{\Delta }p\frac{x}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta }.\end{array}$
Appendix B: Derivation of the velocity component wy

The starting point is Eq. (12) with the volume forces ${f}_{y},{g}_{y}=\mathrm{0}$:

$\begin{array}{}\text{(B1)}& \frac{\partial {w}_{y}}{\partial x}=-\frac{\partial p}{\partial y}\end{array}$

and thus

$\begin{array}{}\text{(B2)}& {w}_{y}=\underset{-\mathrm{\infty }}{\overset{x}{\int }}\frac{\partial {w}_{y}}{\partial x}\mathrm{d}x=-\underset{-\mathrm{\infty }}{\overset{x}{\int }}\frac{\partial {p}_{y}}{\partial y}\mathrm{d}x.\end{array}$

We differentiate the pressure Eq. (A8) with respect to y:

$\begin{array}{}\text{(B3)}& \frac{\partial p}{\partial y}=-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{\Delta }p\frac{-\mathrm{2}\left(y-\mathit{\eta }\right)x}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta }.\end{array}$

Inserting Eq. (B3) into Eq. (B2) gives

$\begin{array}{}\text{(B4)}& {w}_{y}=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{2}\left(y-\mathit{\eta }\right)\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\underset{-\mathrm{\infty }}{\overset{x}{\int }}\frac{x}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}x.\end{array}$

For the integration of the last integral in Eq. (B4) we note the integral is of the type (113) in the handbook of :

$\begin{array}{}\text{(B5)}& \int \frac{t}{T}=\frac{\mathrm{2}a+bt}{qT}-\frac{b}{q}\int \frac{\mathrm{d}t}{T},\end{array}$

where

$T=a+bt+c{t}^{\mathrm{2}}.$

We have

$\begin{array}{ll}& T={t}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}},\\ & a=\left(y-\mathit{\eta }{\right)}^{\mathrm{2}},\\ & b=\mathrm{0},\\ & c=\mathrm{1},\\ & q=\mathrm{4}ac-{b}^{\mathrm{2}}=\mathrm{4}\left(y-\mathit{\eta }{\right)}^{\mathrm{2}},\end{array}$
$\begin{array}{}\text{(B6)}& \begin{array}{rl}{w}_{y}& =\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{2}\left(y-\mathit{\eta }\right)\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\cdot \phantom{\rule{0.33em}{0ex}}\left[\frac{\mathrm{2}a\right)}{qT}\right]{|}_{-\mathrm{\infty }}^{x}\\ & =\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{2}\left(y-\mathit{\eta }\right)\mathrm{\Delta }p\mathrm{d}\mathit{\eta }\\ & \cdot \frac{\mathrm{2}\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}{\left(\mathrm{4}\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}\right)\left({x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}\right)},\end{array}\text{(B7)}& {w}_{y}=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{\Delta }p\frac{\left(y-\mathit{\eta }\right)}{\left({x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}\right)}\mathrm{d}\mathit{\eta }.\end{array}$

The pressure function has in Eq. (A8) been derived to

$\begin{array}{}\text{(C1)}& p=-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{\Delta }p\frac{x}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta }.\end{array}$

When the loading Δp is constant, we can derive the analytical solution for the pressure.

The integral is of the type (107) in the handbook of :

$\begin{array}{}\text{(C2)}& \int \frac{\mathrm{d}t}{T}=\frac{\mathrm{2}}{\sqrt{q}}\mathrm{arctan}\frac{\mathrm{2}ct+b}{\sqrt{q}},\end{array}$

where

$T=a+bt+c{t}^{\mathrm{2}}.$

We have

$\begin{array}{ll}& T={t}^{\mathrm{2}}-\mathrm{2}yt+\left({x}^{\mathrm{2}}+{y}^{\mathrm{2}}\right),\\ & a={x}^{\mathrm{2}}+{y}^{\mathrm{2}},\\ & b=-\mathrm{2}y,\\ & c=\mathrm{1},\\ & q=\mathrm{4}ac-{b}^{\mathrm{2}}=\mathrm{4}\left({x}^{\mathrm{2}}+{y}^{\mathrm{2}}\right)-\left(-\mathrm{2}y{\right)}^{\mathrm{2}},\\ & \sqrt{q}=\mathrm{2}x,\end{array}$

$\begin{array}{}\text{(C3)}& \begin{array}{rl}p& =-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}x\mathrm{\Delta }p\phantom{\rule{0.33em}{0ex}}\left[\frac{\mathrm{2}}{\mathrm{2}x}\mathrm{arctan}\left(\frac{\mathrm{2}t-\mathrm{2}y\right)}{\mathrm{2}x}\right)\right]{|}_{-\mathrm{1}}^{\mathrm{1}}\\ & =-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}x\mathrm{\Delta }p\phantom{\rule{0.33em}{0ex}}\left[\frac{\mathrm{2}}{\mathrm{2}x}\mathrm{arctan}\left(\frac{\mathrm{2}-\mathrm{2}y}{\mathrm{2}x}\right)\\ & -\mathrm{arctan}\left(\frac{-\mathrm{2}-\mathrm{2}y}{\mathrm{2}x}\right)\right],\end{array}\text{(C4)}& p=-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\mathrm{\Delta }p\left[\mathrm{arctan}\left(\frac{\mathrm{1}-y}{x}\right)+\mathrm{arctan}\left(\frac{\mathrm{1}+y}{x}\right)\right].\end{array}$

wy has in Eq. (24) been derived to

$\begin{array}{}\text{(D1)}& {w}_{y}=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{-\mathrm{1}}{\overset{{}^{\mathrm{1}}}{\int }}\mathrm{\Delta }p\frac{y-\mathit{\eta }}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta }.\end{array}$

The integrals are of the type (112) in the handbook of ,

$\begin{array}{}\text{(D2)}& -\int \frac{\mathrm{d}t\cdot t}{T}=-\frac{\mathrm{1}}{\mathrm{2}c}\mathrm{ln}T+\frac{b}{\mathrm{2}c}\int \frac{\mathrm{d}t}{T},\end{array}$

and (107)

$\begin{array}{}\text{(D3)}& y\int \frac{\mathrm{d}t}{T},\end{array}$

where

$T=a+bt+c{t}^{\mathrm{2}}$

and

$\begin{array}{ll}& T={t}^{\mathrm{2}}-\mathrm{2}yt+\left({x}^{\mathrm{2}}+{y}^{\mathrm{2}}\right),\\ & a={x}^{\mathrm{2}}+{y}^{\mathrm{2}},\\ & b=-\mathrm{2}y,\\ & c=\mathrm{1},\\ & q=\mathrm{4}ac-{b}^{\mathrm{2}}=\mathrm{4}\left({x}^{\mathrm{2}}+{y}^{\mathrm{2}}\right)-\left(-\mathrm{2}y{\right)}^{\mathrm{2}}=\mathrm{4}{x}^{\mathrm{2}},\\ & \sqrt{q}=\mathrm{2}x.\end{array}$

We first notice that the last integral in Eq. (B6) cancels out the integral in Eq. (D3), as we have

$\begin{array}{}\text{(D4)}& \frac{b}{\mathrm{2}c}\int \frac{\mathrm{d}t}{T}=\frac{-\mathrm{2}y}{\mathrm{2}}\int \frac{\mathrm{d}t}{T}=-y\int \frac{\mathrm{d}t}{T}.\end{array}$

For the other integral in Eq. (B6), when inserting into Eq. (D1), we have

$\begin{array}{}\text{(D5)}& {w}_{y}=-\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\mathrm{\Delta }p\phantom{\rule{0.33em}{0ex}}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathrm{ln}\left({t}^{\mathrm{2}}-\mathrm{2}yt+\left({x}^{\mathrm{2}}+{y}^{\mathrm{2}}\right)\right)\right]{|}_{-\mathrm{1}}^{\mathrm{1}}=>\phantom{\rule{0.25em}{0ex}},\text{(D6)}& \begin{array}{rl}{w}_{y}& =\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\mathrm{\Delta }p\left[\mathrm{ln}\left(\mathrm{1}-\mathrm{2}y+\left({x}^{\mathrm{2}}+{y}^{\mathrm{2}}\right)\right)\\ & -\mathrm{ln}\left(\mathrm{1}+\mathrm{2}y+\left({x}^{\mathrm{2}}+{y}^{\mathrm{2}}\right)\right)\right]=>\phantom{\rule{0.25em}{0ex}},\end{array}\text{(D7)}& {w}_{y}=-\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\mathrm{\Delta }p\left[\mathrm{ln}\left(\frac{{x}^{\mathrm{2}}+\left(y+\mathit{\eta }{\right)}^{\mathrm{2}}}{{x}^{\mathrm{2}}+\left(y-\mathit{\eta }{\right)}^{\mathrm{2}}}\right)\right].\end{array}$
Appendix E

Table E1List of symbols.

Code availability

The results from the 2Dl model have been simulated with different Python codes which are not well suited for sharing due to many temporary parts for output and other testing.

Data availability

The CFD simulations used for comparisons can be shared by contacting the author.

Competing interests

The author has declared that there are no competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

I thank Kenneth Thomsen, head of division “Wind Turbine Design” in the department “Wind and Energy Systems” at DTU, and Christian Back, head of section “ARD” in the department “Wind and Energy Systems” at DTU, for providing internal funding for this work.

I would also like to thank Flemming Rasmussen in section “ARD” in the department “Wind and Energy Systems” at DTU for carrying out a pre-review.

Likewise the pre-review of the paper by Mads Holst Aagaard Madsen, postdoc at the Technical University of Munich, is acknowledged.

Finally, I thank the reviewer David Wood and an anonymous reviewer for their valuable input, as well as Jaime Liew at MIT for drawing attention to a few typos in the equations.

Financial support

This research has been supported by the department “Wind and Energy Systems” at the Technical University of Denmark.

Review statement

This paper was edited by Alessandro Bianchini and reviewed by David Wood and one anonymous referee.

References

Boorsma, K., Schepers, G., Aagard Madsen, H., Pirrung, G., Sørensen, N., Bangga, G., Imiela, M., Grinderslev, C., Meyer Forsting, A., Shen, W. Z., Croce, A., Cacciola, S., Schaffarczyk, A. P., Lobo, B., Blondel, F., Gilbert, P., Boisard, R., Höning, L., Greco, L., Testa, C., Branlard, E., Jonkman, J., and Vijayakumar, G.: Progress in the validation of rotor aerodynamic codes using field data, Wind Energ. Sci., 8, 211–230, https://doi.org/10.5194/wes-8-211-2023, 2023. a

Bortolotti, P., Kapila, A., and Bottasso, C. L.: Comparison between upwind and downwind designs of a 10 MW wind turbine rotor, Wind Energ. Sci., 4, 115–125, https://doi.org/10.5194/wes-4-115-2019, 2019. a

Bortolotti, P., Ivanov, H., Johnson, N., Barter, G. E., Veers, P., and Namura, N.: Challenges, opportunities, and a research roadmap for downwind wind turbines, Wind Energy, 25, 354–367, https://doi.org/10.1002/we.2676, 2022. a

Bossanyi, E.: Bladed 4.15 Documentation, https://dnvgldocs.azureedge.net/Bladed 4.15 Theory Manual/index.html (last access: 11 December 2023), 2003. a

Bottasso, C., Croce, A., Savini, B., Sirchi, W., and Trainelli, L.: Aero-servo-elastic modeling and control of wind turbines using finite-element multibody procedures, in: Multibody System Dynamics 16, 291–308, https://doi.org/10.1007/s11044-006-9027-1, 2006. a

Branlard, E. and Gaunaa, M.: Cylindrical vortex wake model: right cylinder, Wind Energy, 18, 1973–1987, https://doi.org/10.1002/we.1800, 2015. a, b

Crawford, C.: Re-examining 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., Troldborg, N., and Branlard, E.: A simple vortex model applied to an idealized rotor in sheared inflow, Wind Energ. Sci., 8, 503–513, https://doi.org/10.5194/wes-8-503-2023, 2023. a, b, c

Glauert, H.: Airplane Propellers, in: Division L in Aerodynamic Theory, vol. IV, edited by: Durand, W. F., 169–360, Springer, https://doi.org/10.1007/978-3-642-91487-4_3, 1935. a, b

Jonkman, J., Hayman, G., Jonkman, B., and Damiani, R.: AeroDyn v15 User’s Guide and Theory Manual, NREL, https://www.nrel.gov/wind/nwtc/assets/pdfs/aerodyn-manual.pdf (last access: 8 December 2023), 2016. a

Koning, C.: Influence of the Propeller on Other Parts of the Airplane Structure, Aerodynamic Theory, 361–430, https://doi.org/10.1007/978-3-642-91487-4_4, 1935. a, b, c, d, e, f, g

Larsen, T. and Hansen, A.: How 2 HAWC2, the user's manual, Denmark, Forskningscenter Risoe, Risoe-R-1597, ISBN 978-87-550-3583-6, 2007. a

Li, A., Gaunaa, M., Pirrung, G. R., and Horcas, S. G.: A computationally efficient engineering aerodynamic model for non-planar wind turbine rotors, Wind Energ. Sci., 7, 75–104, https://doi.org/10.5194/wes-7-75-2022, 2022. a, b

Madsen, H. A.: The actuator cylinder :, Institute of Industrial Constructions and Energy Technology, Aalborg University Centre, ISBN 87-87681-04-8, 1982. a

Madsen, H. A.: On the ideal and real energy conversion in a straight bladed vertical axis wind turbine, PhD thesis, Aalborg University, ISBN 87-87681-06-04, 1983. a, b, c

Madsen, H. A.: Actuator Cylinder: A Flow Model for Vertical Axis Wind Turbine, Proceedings of the BWEA Wind Energy Conference (british Wind Energy Association), 147–154, ISBN 0-85298-576-2, 1985. a, b

Madsen, H. A.: Application of Actuator Surface Theory on Wind Turbines, in: IEA R&D WECS: Joint action on aerodynamics of wind turbines, https://www.researchgate.net/publication/326032585_Application_of_Actuator_Surface_Theory_on_Wind_Turbines (lst access: 8 December 2023), 1988. a, b

Madsen, H. A.: Forskning i aeroelasticitet EFP-2000, Chap. 6 Yawberegning i HAWC, 41–50, ISBN 87-550-2891-8, 2001.  a

Madsen, H. A.: A CFD analysis of the actuator disc flow compared with momentum theory results, 109–124, https://www.researchgate.net/publication/281319996_A_CFD_analysis_of_the_actuator_disc_flow_compared_with_momentum_theory_results (last assess: 8 December 2023), 1997. a, b

Madsen, H. A. and Rasmussen, F.: The influence on energy conversion and induction from large blade deflections, Wind Energy for the Next Millennium, Proceedings, 138–141, 1999. a, b, c, d

Madsen, H. A., 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, b

Madsen, H. A., Larsen, T. J., Schmidt Paulsen, U., and Vita, L.: Implementation of the Actuator Cylinder Flow Model in the HAWC2 code for Aeroelastic Simulations on Vertical Axis Wind Turbines, Proceedings of 51st Aiaa Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 12 pp., ISBN 978162410181, 2013. a, b

Madsen, H. A., 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/wes-5-1-2020, 2020. a, b, c, d, e, f

Mikkelsen, R., Sørensen, J. N., and Shen, W. Z.: Modelling and analysis of the flow field around a coned rotor, Wind Energy, 4.3 121–135, https://doi.org/10.1002/we.50, 2001. a, b, c, d

Øye, S.: FLEX4 simulation of wind turbine dynamics, Proceedings of the 28th IEA Meeting of Experts Concerning State of the Art of Aeroelastic Codes for Wind Turbine Calculations, Lyngby, Denmark, edited by: Maribo Pedersen, B., https://www.osti.gov/etdeweb/servlets/purl/283712 (last access: 11 December 2023), 1996. a

Pao, L. Y., Zalkind, D. S., Griffith, D. T., Chetan, M., Selig, M. S., Ananda, G. K., Bay, C. J., Stehly, T., and Loth, E.: Control co-design of 13 MW downwind two-bladed rotors to achieve 25 % reduction in levelized cost of wind energy, Annu. Rev. Control, 51, 331–343, https://doi.org/10.1016/j.arcontrol.2021.02.001, 2021.  a

Riziotis, V. and Voustinas, S.: Gast: A General Aerodynamic and Structural Prediction Tool for Wind Turbines, in: Proceedings of the EWEC'97, Dublin Castle, Ireland, October 1997, ISBN 9780953392209, 1997. a

Rumble, J. R. (Ed.): CRC Handbook of Chemistry and Physics, 104th Edition (Internet Version 2023), CRC Press/Taylor & Francis, Boca Raton, FL, 2023. a, b, c, d

Sørensen, J. N. and Kock, C. W.: A model for unsteady rotor aerodynamics, J. Wind Eng. Ind. Aerod., 58, 259–275, https://doi.org/10.1016/0167-6105(95)00027-5, 1995. a, b

Troldborg, N., Sørensen, J. N., Mikkelsen, R., and Sørensen, N.: A simple atmospheric boundary layer model applied to large eddy simulations of wind turbine wakes, Wind Energ., 17, 657–669, https://doi.org/10.1002/we.1608, 2014. a, b

van Bussel, G.: The use of the asymptotic acceleration potential method for horizontal axis wind turnime rotor aerodynamics, J. Wind Eng. Ind. Aerod., 39, 161–172, https://doi.org/10.1016/0167-6105(92)90542-I, 1992. a

van Kuik, G. A. and Lignarolo, L. E.: Potential flow solutions for energy extracting actuator disc flows, Wind Energy, 19, 1391–1406, https://doi.org/10.1002/we.1902, 2016. a

Von Kármán, T. and Burgers, J. M.: Aerodynamic Theory: General aerodynamic theory: Perfect fluids, 2, J. Springer, https://books.google.dk/books?id=RC0HaAEACAAJ (last access: 11 December 2023), 1935. a, b, c, d, e, f, g

Vølund, P. and Rasmussen, F.: Designverifikation for fleksibel to-bladet mølle, https://backend.orbit.dtu.dk/ws/portalfiles/portal/7729706/ris_r_1073.pdf (last access: 11 December 2023), 1999. a

WMC: FOCUS6 – The integrated wind turbine design suite, Knowledge Centre WMC, https://dokumen.tips/documents/focus6-the-integrated-wind-turbine-design-suite-wmceu-on-and-offshore-wind.html?page=1 (last access: 8 December 2023), 2019. a

Wood, D. H. and Hammam, M. M.: Optimal performance of actuator disc models for horizontal-axis turbines, Frontiers in Energy Research, 10, 971177, https://doi.org/10.3389/fenrg.2022.971177, 2022. a