Modeling the influence of streamwise flow field acceleration on the aerodynamic performance of an actuator disk

Wind turbines operating in complex terrain can be subject to a background flow field that varies in the streamwise direction. This variation can influence the induced velocity at the turbine and thus also the power performance. In the present work, a simple model is derived for the situation of an actuator disk (AD) operating in a background flow field featuring a constant streamwise velocity gradient. Reynolds-averaged Navier–Stokes (RANS) simulations of this scenario are performed, showing that a positive acceleration yields a reduction of induction, while a negative acceleration leads to an increase in induction. The new model accurately captures this behavior and significantly reduces the prediction error compared to classical momentum theory, where the effect of the background flow acceleration is disregarded. The model indicates that the maximum power coefficient and the corresponding values of the optimal induction and thrust coefficient depend on the flow acceleration. This possibly impacts optimal operational strategies under such conditions as well as wind turbine and wind farm design considerations, which often rely on the assumption of a streamwise uniform flow field.

1 Introduction

Classical momentum theory as developed by Rankine (1865), Froude (1878), and Froude (1889) yields insights into the relation between thrust and induction of an idealized wind turbine operating in a uniform flow field. Also, the theoretical limit of aerodynamic performance, referred to as the Betz–Joukowsky limit, can be derived from it (Betz, 1920; Joukowsky, 1920; Okulov and Van Kuik, 2012). The aerodynamic performance referred to here is the aerodynamic power, which is the product of thrust and axial velocity through the turbine plane. Today, momentum theory plays a fundamental role in the design, power performance predictions, and aeroelastic analysis of wind turbines. However, due to its simplifying assumptions, which do not necessarily apply to the actual operating conditions of modern wind turbines, various modifications to momentum theory have been proposed. When turbines operate in the atmospheric boundary layer, they are usually subject to wind shear, i.e., a variation of streamwise velocity over height. Various models exist to account for this. Wagner et al. (2011) suggested using a reference velocity for performance predictions based on the kinetic energy flux through the turbine disk plane and validated this model with measurements of a full-scale turbine. Later, Chamorro and Arndt (2013) derived an analytical correction of the Betz limit by including the effect of shear upstream and downstream of the wind turbine in momentum theory. Because shear would also lead to a variation of thrust force over the disk, Draper et al. (2016) included this additional effect in their model. Based on a vortex theory framework, Gaunaa et al. (2023) showed that momentum theory should be applied locally within momentum-based design and analysis tools to properly account for the effect of shear. The problem of yaw misalignment and its effect on induction was approached in the early days by Glauert (1926) and more recently by Heck et al. (2023) and then by Tamaro et al. (2024), who incorporate both shear and yaw misalignment in their model. Mikkelsen and Sørensen (2002) developed a correction to account for the effect of wind tunnel blockage on the aerodynamic performance of turbines in wind tunnels. The potential of using a diffuser to improve the aerodynamic turbine performance was studied by, e.g., Jamieson (2009). Sørensen (2016) discusses the theory of both wind tunnel blockage and diffuser modeling in detail.

What has not yet been considered in modeling is the effect of a flow that is generally inhomogeneous in the streamwise direction. This can be the case in complex terrain or in dense wind farms where obstacles or turbines interact with the atmospheric boundary layer. However, various authors have shown that these variations in the streamwise flow field development can have a significant impact on the aerodynamic performance of a wind turbine, suggesting that it could also be beneficial to approach this problem from an analytical perspective. The previous literature has mainly focused on describing the phenomenon, while few studies have attempted to model it. Troldborg et al. (2022) drew attention to this phenomenon in a large eddy simulation (LES) study, in which they investigated the power performance of a full-scale wind turbine placed on top of the first of two ridges, which were part of the Perdigão field campaign. In the study, they varied the surface roughness and showed that this affects both the separation and the flow development behind the ridge: for low surface roughness, the flow stayed attached, while it detached for high values. However, the undisturbed wind speed at the turbine position when it was not operating remained nearly the same in all cases. Despite this, once the turbine was operating, the predicted mean power performance varied significantly between the different cases of surface roughness. This indicated that the way the wake is advected downstream behind the turbine also affects the performance. Zengler et al. (2024) performed idealized Reynolds-averaged Navier–Stokes (RANS) simulations of an actuator disk (AD) located on top of a quasi-two-dimensional Gaussian hill. Simulations were conducted for both different surface roughnesses and hill widths. The thrust force was distributed evenly over the AD, and no turbine controller was involved, allowing a clearer picture of how exactly the performance of the AD is influenced by the flow field. Similar to Troldborg et al. (2022), they found that for low roughness, the flow remained attached to the surface behind the ridge, leading to a deceleration and eventually to a higher induction at a given thrust coefficient, which is equivalent to a lower power performance. Here, the thrust coefficient was based on the total thrust due to the AD acting on the flow and the rotor-averaged wind speed at the position of the turbine when it was not operating. This approach is also adopted in the present work. As the surface roughness increased, the situation became more complicated, especially when the undisturbed background flow detached from the surface already at the top of the ridge, as was the case when the hill width was very small; in this case, a deceleration was still noticeable at hub height. However, the induction would be lower compared to the flat reference case, which is equivalent to an increase in power performance. Dar et al. (2023) performed wind tunnel experiments of a more idealized scenario to study how wakes behave under streamwise pressure gradients. The pressure gradient was imposed by placing the turbine model on a ramp in the tunnel, causing the flow to accelerate or decelerate, depending on whether the ramp was tilted upward or downward. In addition, the power coefficient was reported as showing an increase for an accelerating background flow and a decrease for a decelerating background flow. Revaz and Porté-Agel (2024) performed an LES of a wind turbine on the ridge of a quasi-two-dimensional hill parameterized by a trigonometric function and varied geometrical features of the setup such as hill height, hill width, plateau length, and turbine hub height. In all cases, the background flow field decelerated behind the ridge, and they reported a trend of decreasing power performance with increasing deceleration. An outlier from this trend was identified for a steep hill and a strong separation forming behind it, which is consistent with the previously mentioned findings by Zengler et al. (2024). Further works in which the influence of streamwise flow variations on power performance is reported exist (Cai et al., 2021; Mishra et al., 2024). However, it is often the case that the effect is not investigated in isolation, making it difficult to draw conclusions. For example, the presence of a controller raises the question of how this reacts in an accelerating flow field. In addition, it is often unclear what the reference wind speed is at the turbine position, which is essential to investigate the bias in power performance due to the development of the flow field downstream. Currently, most research shows that flow deceleration behind the turbine can lead to a decrease in performance, while flow acceleration results in an increase in performance.

Efforts to include the effect of acceleration in turbine modeling have been mainly focused on the wake behavior under such circumstances (Shamsoddin and Porté-Agel, 2018; Dar and Porté-Agel, 2022; Dar et al., 2023). Only Cai et al. (2021) attempted to also model the power performance in the presence of pressure gradients, which is, in a sense, equivalent to modeling flow acceleration. For this purpose, they used a wake model developed for pressure gradients in combination with a linearized flow solver and a control volume analysis. Their approach is based on the linear perturbation theory for the flow over hills (Jackson and Hunt, 1975; Hunt et al., 1988) and solves a Laplace equation for the disturbance in the velocity field. Although this model can predict both wake profiles and power performance, which are in sound agreement with measurements, their approach comes with an additional computational burden compared to fully analytical models and also cannot be directly incorporated into aeroelastic design because it does not yield local information about the flow state in the turbine plane. Apart from these efforts, no other work, to the authors' knowledge, has aimed at developing an analytical model to account for the effect of streamwise flow acceleration on aerodynamic performance. The present work is an attempt to model the effect of a constant streamwise flow acceleration on the induction and thus on the aerodynamic performance of a wind turbine.

The remainder of this article is structured as follows. First, a simplified model is derived to account for the effect of streamwise flow acceleration in momentum theory. Second, the model is validated with RANS simulations of an AD exposed to various constant velocity gradient flows. Last, implications and modeling details are discussed.

Figure 1 Sketch of the problem and notation for theoretical analysis.

2 Model derivation

In the following, a model is derived that relates the induced velocity in the AD plane with the thrust coefficient of the AD in an accelerating flow field. It should be stressed that the derivation is not without contradictions from a strict, physical perspective, which will be discussed in more detail in Sect. 5.1. However, the final model proves to yield comparably accurate predictions, and the derivation itself allows for a clearer understanding of which effects the model captures. The theory development is based on the idea of modifying momentum theory to account for the fact that the underlying background flow field changes between the inlet and outlet of the streamtube surrounding the AD. The derivation is very similar to classical momentum theory and mainly differs when it comes to the treatment of the pressure forces acting on the streamtube and when one needs to relate the undisturbed velocities at the streamtube inlet, the AD position, and the streamtube outlet with each other. The latter step is not necessary in momentum theory, as the undisturbed background flow field does not change between these positions.

As in momentum theory, a steady, incompressible, and inviscid flow past an AD is considered. Figure 1 shows a sketch of the problem and its notation, which is based on the following logic: U and u are the undisturbed and disturbed velocities, respectively. A subscript _∞ refers to quantities in the undisturbed flow, while the subscript ₀ denotes quantities ahead of the disk plane, _R denotes quantities in the disk plane, and ₁ denotes the ultimate wake. In addition, the subscript _u denotes the case of a uniform flow field without a velocity gradient. The undisturbed flow field is subject to a constant velocity gradient, i.e.,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\text{d}{U}_{\mathrm{\infty }}}{\text{d}x}}=C.\end{array}$$

Because an infinite axial extent of the gradient would eventually lead to negative velocities upstream or downstream of the AD, the extent of the streamtube considered in Fig. 1 is bounded to regions of positive velocities. The limiting case of negative velocities in the wake is discussed in Appendix B. For now, the model derivation is continued, keeping this detail in mind.

By virtue of the Bernoulli equation, the pressure jump Δp across the AD fulfills

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }p={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{U}_{\mathrm{\infty },\mathrm{0}}^{\mathrm{2}}+p_{\mathrm{\infty },\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{u}_{\mathrm{1}}^{\mathrm{2}}-p_{\mathrm{1}}.\end{array}$$

Further, the undisturbed flow quantities ahead and behind the disk are related by the Bernoulli equation as well by

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{U}_{\mathrm{\infty },\mathrm{0}}^{\mathrm{2}}+p_{\mathrm{\infty },\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{U}_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}+p_{\mathrm{\infty },\mathrm{1}},\end{array}$$

which yields for the pressure jump

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{\Delta }p={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{U}_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}+p_{\mathrm{\infty },\mathrm{1}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{u}_{\mathrm{1}}^{\mathrm{2}}-p_{\mathrm{1}}.\end{array}$$

By assuming, as in classical momentum theory (Rankine, 1865; Froude, 1878, 1889), that the pressure difference between the wake and the free stream is zero (i.e., $p_{\mathrm{1}}=p_{\mathrm{\infty },\mathrm{1}}$), the pressure drop over the AD is seen to represent the ultimate difference in kinetic energy between free stream and wake:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{\Delta }p={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }\left(U_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}-u_{\mathrm{1}}^{\mathrm{2}}\right).\end{array}$$

Next, momentum balance on the streamtube sketched in Fig. 1 is considered, yielding

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{\rho }{A}_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{2}}-\mathit{\rho }{A}_{\mathrm{0}}{U}_{\mathrm{\infty },\mathrm{0}}^{\mathrm{2}}=-T+F_{\mathrm{p}},\end{array}$$

with the thrust force T related to the pressure jump as

$$\begin{array}{}\text{(7)}& {\displaystyle }T=\mathrm{\Delta }p{A}_{\mathrm{R}}\end{array}$$

and the axial contribution of the pressure acting on the streamtube, F_p. The inlet and outlet areas of the streamtube A₀ and A₁ are related to the rotor area A_R by conservation of mass:

$$\begin{array}{}\text{(8)}& {\displaystyle }{U}_{\mathrm{\infty },\mathrm{0}}{A}_{\mathrm{0}}=u_R{A}_{\mathrm{R}}=u_{\mathrm{1}}{A}_{\mathrm{1}}.\end{array}$$

In contrast to classical momentum theory, the pressure term F_p is not zero, as the pressure in the undisturbed flow is not constant in the axial direction. Here, it is approximated to be

$$\begin{array}{}\text{(9)}& {\displaystyle }{F}_{\mathrm{p}}\approx \left(p_{\mathrm{\infty },\mathrm{0}}-p_{\mathrm{1}}\right)A_{\mathrm{R}}.\end{array}$$

Combining the equation for momentum balance Eq. (6) with the mass balance Eq. (8) yields

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathit{\rho }{A}_{\mathrm{R}}{u}_R{u}_{\mathrm{1}}-\mathit{\rho }{A}_{\mathrm{R}}{u}_R{U}_{\mathrm{\infty },\mathrm{0}}=-T+F_P.\end{array}$$

Further, the thrust force T is replaced by the product of the pressure jump and disk area as defined in Eq. (7), and the term for the pressure force (Eq. 9) is introduced into the momentum balance as well:

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathit{\rho }{A}_{\mathrm{R}}{u}_R{u}_{\mathrm{1}}-\mathit{\rho }{A}_{\mathrm{R}}{u}_R{U}_{\mathrm{\infty },\mathrm{0}}=-\mathrm{\Delta }p{A}_{\mathrm{R}}+\left(p_{\mathrm{\infty },\mathrm{0}}-p_{\mathrm{1}}\right)A_{\mathrm{R}}.\end{array}$$

Expressing Δp via Eq. (5) and $p_{\mathrm{\infty },\mathrm{0}}-p_{\mathrm{1}}$ via Eq. (3) yields

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathit{\rho }{A}_{\mathrm{R}}{u}_R{u}_{\mathrm{1}}-\mathit{\rho }{A}_{\mathrm{R}}{u}_R{U}_{\mathrm{\infty },\mathrm{0}}=& -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }\left(U_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}-u_{\mathrm{1}}^{\mathrm{2}}\right)A_{\mathrm{R}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }\left(U_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}-U_{\mathrm{\infty },\mathrm{0}}^{\mathrm{2}}\right)A_{\mathrm{R}}.\end{array}\end{array}$$

This equation can be solved for the rotor velocity, leading to

$$\begin{array}{}\text{(13)}& {\displaystyle }{u}_R={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(U_{\mathrm{\infty },\mathrm{0}}+u_{\mathrm{1}}\right).\end{array}$$

This result is similar to the classical result by Froude (1889), which says that in non-accelerating flows, the velocity in the rotor plane is the average between the undisturbed, upstream velocity and the wake velocity. In the case of an accelerating flow, the velocity varies in the streamwise direction. Thus, it needs to be decided where to evaluate U_∞ to obtain U_∞,0 but also U_∞,1, which is necessary to obtain the wake velocity u₁ based on Eq. (5). For this purpose, a length scale L is introduced, which relates the undisturbed velocity in the disk plane with the undisturbed velocities upstream and downstream of the disk due to the constant gradient condition by

$$\begin{array}{}\text{(14)}& {\displaystyle }{U}_{\mathrm{\infty },\mathrm{0}}=U_{\mathrm{\infty },R}-L{\displaystyle \frac{\mathrm{d}{U}_{\mathrm{\infty }}}{\mathrm{d}x}},\text{(15)}& {\displaystyle }{U}_{\mathrm{\infty },\mathrm{1}}=U_{\mathrm{\infty },R}+L{\displaystyle \frac{\mathrm{d}{U}_{\mathrm{\infty }}}{\mathrm{d}x}}.\end{array}$$

These relations close the presented system of equations for a given length scale. Lastly, the induction factor and the thrust coefficient are introduced as

$$\begin{array}{}\text{(16)}& {\displaystyle }a& {\displaystyle }=\mathrm{1}-{\displaystyle \frac{u_R}{U_{\mathrm{\infty },R}}},\text{(17)}& {\displaystyle }{C}_{\mathrm{T}}& {\displaystyle }={\displaystyle \frac{T}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{A}_{\mathrm{R}}{U}_{\mathrm{\infty },R}^{\mathrm{2}}}}.\end{array}$$

Combining Eq. (17) with Eqs. (7), (5), and (13) leads to

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{C}_{\mathrm{T}}& ={\displaystyle \frac{\mathrm{\Delta }p{A}_{\mathrm{R}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{A}_{\mathrm{R}}{U}_{\mathrm{\infty },R}^{\mathrm{2}}}}={\displaystyle \frac{U_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}-u_{\mathrm{1}}^{\mathrm{2}}}{U_{\mathrm{\infty },R}^{\mathrm{2}}}}\\ & ={\displaystyle \frac{U_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{u}_R^{\mathrm{2}}+\mathrm{4}{u}_R{U}_{\mathrm{\infty },\mathrm{0}}-U_{\mathrm{\infty },\mathrm{0}}^{\mathrm{2}}}{U_{\mathrm{\infty },R}^{\mathrm{2}}}}.\end{array}\end{array}$$

The free stream velocities ahead and behind the AD are replaced by the relations in Eqs. (14) and (15), and after some rearranging, one obtains

$$\begin{array}{}\text{(19)}& {\displaystyle }{C}_{\mathrm{T}}={\displaystyle \frac{\mathrm{4}{u}_R\left(U_{\mathrm{\infty },R}-u_R\right)+\mathrm{4}\left(U_{\mathrm{\infty },R}-u_R\right)L\frac{\mathrm{d}{U}_{\mathrm{\infty }}}{\mathrm{d}x}}{U_{\mathrm{\infty },R}^{\mathrm{2}}}}.\end{array}$$

Applying the definition of the induction Eq. (16) and introducing the non-dimensional length scale $l=L/D$ and the non-dimensional velocity gradient $\mathit{\beta }=\frac{D}{U_{\mathrm{\infty },R}}\frac{\mathrm{d}{U}_{\mathrm{\infty }}}{\mathrm{d}x}$, with D denoting the AD diameter, yields the final form of the model:

$$\begin{array}{}\text{(20)}& {\displaystyle }{C}_{\mathrm{T}}=\mathrm{4}a(\mathrm{1}-a)+\mathrm{4}al\mathit{\beta }.\end{array}$$

From this equation, it can be seen that, for an accelerating flow field, the relation for the thrust in a uniform flow field is modified by an additive term proportional to the free stream velocity gradient and the induction. The term lβ is equivalent to a relative change in velocity over the distance l, with l itself being unknown. Crespo et al. (1999) stated that the length of the initial wake expansion region, thus the region where wake and free stream pressure equalize, is effectively approximately one diameter long, a value which was also adopted by Dar and Porté-Agel (2022). Because the present derivation also assumes that the pressures equalize at the non-dimensional position l, a reasonable first estimate of this value corresponds to one diameter, l=1. The central assumptions to obtain this result are symmetry with respect to the axial extent of the control volume and the approximation of the pressure contribution to the momentum balance. Because the first component on the right-hand side of Eq. (20) simply is the thrust-induction relation in a uniform flow, a reasonable generalization might be that arbitrary thrust-induction curves obtained in a uniform flow field can be corrected for acceleration by simply adding the last term on the right-hand side of Eq. (20). Thus,

$$\begin{array}{}\text{(21)}& {\displaystyle }{C}_{\mathrm{T}}\left(a\right)=C_{T,u}\left(a\right)+\mathrm{4}al\mathit{\beta },\end{array}$$

with the thrust-induction relation C_T,u(a) for a uniform background flow field. Equation (21) is the new model, which will be validated in Sect. 4.

Figure 2 Boundaries of the computational domains considered, simulating a constant velocity gradient. In the lateral direction, the domain width is 11.390 D. The gray area indicates a region of mostly cubic cells of a side length of approximately 0.047 D. Dashed lines indicate two domain configurations that are not considered within the model validation. The color-coding for the different domains is used consistently in all subsequent figures.

Figure 3 Computational grid for a deceleration, with only every eighth grid line shown. (a) Grid in the x–z plane and (b–d) grid in the y–z plane at x equal to (b) −20 D, (c) 0 D, and (d) 50 D,

Figure 4 Comparison of Eulerian and RANS simulations for the baseline uniform case with theoretical results from momentum theory.

3 Simulation setup

Three-dimensional RANS simulations of an AD operating in a flow field featuring a constant velocity gradient in the streamwise direction are performed in order to validate the simple model later in Sect. 4. The thrust-induction curves for a wide range of positive and negative accelerations are evaluated. In the following, the simulation setup is presented in detail.

Figure 5 Simulation results with (a) the undisturbed streamwise velocity along the centerline of the AD, (b) the streamwise gradient, and (c) the predicted thrust-induction curves.

3.1 Grid design

For every considered velocity gradient, a separate computational grid is generated. In all grids, the Cartesian coordinate system is centered around the AD, with x, y, and z denoting the streamwise, lateral, and vertical directions, respectively. The cross-sectional area A_S(x) of each grid is defined from mass conservation by the constant velocity gradient condition, i.e.,

$$\begin{array}{}\text{(22)}& {\displaystyle }{A}_S\left(x\right)={\displaystyle \frac{U_R{A}_{S,R}}{U_R+x\frac{\mathrm{d}{U}_{\mathrm{\infty }}}{\mathrm{d}x}}},\end{array}$$

with U_R and A_S,R being the undisturbed velocity and the cross-sectional area in the AD plane. These quantities are held constant for all grids. The latter is quadratic with a side length of 11.390 D, resulting in a blockage ratio of 0.605 %¹. Based on Eq. (22), the upper and lower boundaries of the grids are modified vertically for the different velocity gradients. The lateral boundaries are not modified, resulting in a quasi-two-dimensional domain. The diffuser region has a length of 34.170 D, with the AD being placed at the end of its first third. Inlet and outlet regions are appended to the diffuser with a length of 11.390 and 34.170 D, respectively, and the transitions between these regions are smoothened by a Laplacian smoothing algorithm. The vertical boundaries of the domains are shown in Fig. 2. In the case of negative velocity gradients, it has been found that a recirculation zone can develop in the wake. In the case of the strongest negative velocity gradient simulated (later identified as the case where $\mathit{\beta }=-\mathrm{0.029}$), a recirculation zone forms 8 D behind the disk at C_T=0.6 and moves as close as 3 D towards the disk at $C_{\mathrm{T}}=\frac{\mathrm{8}}{\mathrm{9}}$. In some more moderate scenarios of negative acceleration, the recirculation does not seem to affect the induction significantly, as will be shown later in Sect. 4.1. Because recirculation is not taken into account within momentum theory or in the presented model, the cases that seem to affect the induction noticeably are disregarded in the model validation process and marked by dashed lines in Fig. 2 and all subsequent figures. A more thorough discussion on this can be found in Appendix B. Returning to the domain design, the generated grids are curvilinear and consist of $\mathrm{320}\times \mathrm{128}\times \mathrm{128}$ cells in the x, y, and z direction, respectively. Refinement regions exist at the transitions between the diffuser and the inlet and outlet regions. Further, from a distance of 1.708 D ahead of the AD to 5.125 D behind it, a refinement region with a lateral and vertical side length of approximately 2.278 D exists. In this region, cells are nearly cubic, with a resolution of approximately 21 cells per diameter. In Fig. 3, a computational grid for the case of a deceleration (later identified as $\mathit{\beta }=-\mathrm{0.015}$) is shown.

3.2 Actuator disk model

The AD is simulated with a uniform distribution of thrust force per disk area. As described by Zengler et al. (2024), the outer 25 % of the radial boundaries are linearly smeared out to improve solution convergence. The induction is evaluated as the area weighted mean induction over the part of the disk, which is not affected by the smearing. The actuator shape method is applied to project the force onto the computational grid (Réthoré et al., 2014; Troldborg et al., 2015).

3.3 Boundary conditions, RANS closure, and solver

Inlet and outlet conditions are applied at the streamwise boundaries. The inlet velocity is specified such that U_R is reached in the AD plane. The chosen U_R corresponds to a diameter-based Reynolds number of 3.798 million. Inlet turbulence properties correspond to a turbulence intensity of 1 %, with the specific dissipation rate ω set to 1 million. The lateral boundaries are periodic, and slip conditions are applied at the upper and lower boundaries. As RANS closure, the two-equation k–ω–SST model without modification of the standard model coefficients is applied (Menter, 1994). The incompressible, three-dimensional, Navier–Stokes solver EllipSys3D (Sørensen, 1995; Michelsen, 1992) is employed. The used solution algorithm is a SIMPLE-like procedure based on the non-relaxed momentum equations (Sørensen, 2018). Rhie–Chow interpolation is included to avoid decoupling of the pressure and velocity fields (Rhie and Chow, 1983), and the convective terms are discretized with the QUICK scheme (Leonard, 1979).

3.4 Verification and sensitivity of simulation setup

In Fig. 4, the thrust-induction curves obtained using the RANS and Euler equations for a uniform background flow are compared with momentum theory. It is seen that the Eulerian and RANS simulations do not perfectly agree with momentum theory. This discrepancy can be attributed to the numerical discretization, as previously reported (Hodgson et al., 2021; Mikkelsen, 2004). If better agreement with the momentum theory results is desired, it is necessary to increase the grid resolution or to probe the induction slightly behind the position of the AD as done in a previous study (Zengler et al., 2024). However, within this work, no correction is performed; thus, the curve representing the RANS results is used as the zero-acceleration/uniform curve in the following.

In Appendix A, a detailed sensitivity study of the simulation setup is performed. These results can be summarized as follows. A mesh sensitivity analysis shows that the change in induction due to halving the cell length in all directions is less than 0.3 %. Domain blockage and radial smearing influence the results on the order of 1 %, while the effect of the diffuser length on the results is negligible.

Figure 6 Model predictions for l=1: (a) collapsed induction curves for C_T−4alβ and (b) MAE for induction predictions without and with correction for acceleration.

4 Model validation

First, the simulation results from the previously introduced RANS simulations are presented for the various velocity gradients considered. Subsequently, the model predictions for these cases are analyzed in detail.

4.1 Simulation results

The simulation results are presented in Fig. 5. The color-coding is consistent with that in Fig. 2 and the subsequent Figs. 6 and 8, where each color represents one specific velocity gradient. In panels (a) and (b) of Fig. 5, the normalized streamwise velocities U_∞ along the centerline and the respective normalized gradients for the different domains presented previously are shown. In all cases, a nearly constant streamwise gradient is obtained ahead and behind the AD. Panel (c) shows the predicted thrust-induction curves obtained by applying a certain thrust coefficient and evaluating the corresponding induction. As the thrust coefficient increases, the differences between the individual curves grow. In cases of a positive acceleration, the induction is reduced, while it is increased in cases of a negative acceleration. For the two cases with the strongest negative acceleration – dashed lines in Fig. 5 and not considered for model validation – a more complex induction behavior can be observed for thrust coefficients higher than 0.5. This can be related to the development of a recirculation zone in the wake and is discussed in Appendix B in more detail.

4.2 Model predictions

In the following, it is investigated to what extent the model represented by Eq. (21) is capable of predicting the trends observed in Fig. 5. The non-dimensional velocity gradient β is evaluated in the undisturbed flow field at the position of the AD in the center, thus in Fig. 5b at $x/D=\text{0}$. As stated before, the length of the initial wake expansion region is approximately one diameter long (Crespo et al., 1999).

Model predictions are presented in Fig. 6. If the model works as expected, solving Eq. (21) for $C_{T,u}\left(a\right)=C_{\mathrm{T}}\left(a\right)-\mathrm{4}al\mathit{\beta }$ and inserting the respective simulation results into the right-hand side of this formula should lead to a collapse of all simulation results on the same curve. This is presented in panel (a) of Fig. 6. Indeed, it can be seen that all curves previously shown in Fig. 5c tend to collapse on the zero-acceleration curve. However, at thrust coefficients higher than 0.75, increasing discrepancies between the single curves can be observed.

To quantify the model accuracy, the mean absolute error (MAE) is evaluated as

$$\begin{array}{}\text{(23)}& \text{MAE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n|a-a_{\mathrm{pred}}|}{n}}\end{array}$$

over all thrust coefficients and all n=11 simulations. While a is the induction evaluated from a simulation at a specific acceleration at a given C_T, a_pred is the predicted induction for this acceleration either based on the thrust-induction curve obtained from the uniform case or from the thrust-induction curve corrected by the acceleration term. Because the power coefficient is calculated as $C_{\mathrm{P}}=C_{\mathrm{T}}(\mathrm{1}-a)$, the MAE directly translates to a variation in C_P for a given C_T. Results are presented in Fig. 6b. If no correction is considered, the mean error increases with an increasing thrust coefficient up to a value of around 0.024. Correcting the predictions for flow acceleration yields a reduction of error for all thrust coefficients. Specifically, for C_T<0.5, the error is reduced by 54 % on average, while for C_T≥0.5, the error is reduced by over 78 % on average. The maximum error in this case at high thrust is approximately 0.007. A significant increase in the mean error is noticeable above a thrust coefficient of 0.75. At this point, it is important to note two things. First, the MAE was only calculated for the presented simulations. It does not yield insights into specific cases but rather gives an estimate for the potential error reduction. Second, l=1 was set without further investigation. Based on the presented simulations, it is possible to extract the length scale, which would lead to an exact collapse of thrust-induction curves obtained in a non-uniform flow on the uniform flow curve. This will be discussed in Sect. 5.3.

5 Discussion

The above results demonstrate that the presented model captures the essential physical mechanisms influencing the induction of an AD subject to a flow featuring a constant velocity gradient. The model structure indicates that in an accelerating flow field, the axial induction behaves approximately as if the uniform case was subject to an additional thrust coefficient proportional to −4alβ.

5.1 Modeling assumptions

Within the derivation of the model, it was assumed that the pressure force acting on the streamtube can be easily quantified as the product of the disk area and the pressure difference between upstream and downstream pressure. This assumption eventually led to Eq. (13), which states that the rotor velocity depends on the undisturbed velocity upstream and the wake velocity downstream. With both being dependent on the axial distance, the model results would depend on their position of evaluation. To overcome this problem, the length scale l, which closes the system of equations, was introduced. For nonzero values of l, model predictions depend on this free parameter. If all elements of the model were strictly correct, the arbitrary choices of the position of the inlet and outlet of the control volume would not affect the model results, as long as they were chosen far enough away from the disk to justify the assumption of a negligible pressure disturbance from the disk, $p_{\mathrm{1}}=p_{\mathrm{\infty },\mathrm{1}}$, used in the derivation. Further, rearranging Eq. (13) and applying Eq. (14) indicates that the wake velocity varies as follows:

$$\begin{array}{}\text{(24)}& {\displaystyle }{u}_{\mathrm{1}}=\mathrm{2}{u}_R-U_{\mathrm{\infty },R}+L{\displaystyle \frac{\mathrm{d}{U}_{\mathrm{\infty }}}{\mathrm{d}x}}.\end{array}$$

For arbitrary values of L, equivalent to l, the energy equation yields

$$\begin{array}{}\text{(25)}& {\displaystyle }{u}_{\mathrm{1}}=\sqrt{U_{\mathrm{\infty },\mathrm{1}}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{\Delta }p}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }}}}.\end{array}$$

According to Eq. (24), the wake velocity varies linearly with L. Meanwhile, Eq. (25), which is based on first principles, shows a nonlinear dependency of L. These observations above reveal the effects of the inconsistencies in the assumption about the pressure force on which the model is built on.

For the used simulation setup, it was inherently assumed that the extent of the gradient is sufficiently large, so variations of the gradient further upstream or downstream do not affect the induction. In other words, the control volume of the model, with its axial extent defined by l, always lies completely within the constant velocity gradient region. In the sensitivity analysis in Appendix A, it can be seen that shortening the diffuser length by 33 % has an insignificant effect on the resulting induction. For a future work, it would be interesting to investigate how short the extent of the velocity gradient can be without affecting the induction. This could shed light on the length scale l.

The influence of the pressure force F_p in Eq. (9) was assumed to be approximately the net force acting on the AD due to the pressure difference upstream and downstream. Its exact calculation would require the evaluation of the integral

$$\begin{array}{}\text{(26)}& {\displaystyle }{F}_{\mathrm{p}}=-\underset{CV}{\oint }p\mathit{d}\mathit{A}\cdot {\mathit{e}}_x\end{array}$$

over the streamtube with the infinitesimal surface area dA and the normal vector in the axial direction e_x. To solve this, the exact development of the flow field within the streamtube would need to be known, which on the other hand requires the evaluation of this integral.

Despite the inconsistencies of the simple model pointed out in the discussion above, the model has shown to capture the relevant physics that an accelerating flow field has on the induction of an AD. This justifies further investigations using the model in the next section.

Figure 7 Dependence of thrust and power coefficients C_T and C_P on induction based on Eq. (20) for various β and l=1. Solid lines represent the uniform background flow curve.

5.2 Influence on power performance

Because aerodynamic power can generally be calculated as the product of thrust and velocity, the modification of the relationship between these two also yields a different power performance of wind turbines. Based on the general relation between the thrust and power coefficient

$$\begin{array}{}\text{(27)}& {\displaystyle }{C}_{\mathrm{P}}=C_{\mathrm{T}}(\mathrm{1}-a),\end{array}$$

their curves with respect to induction are shown in Fig. 7 for several velocity gradients. It can be readily seen that a negative acceleration leads to a general degradation of power performance, while a positive acceleration leads to an improvement. This is consistent with previous findings (Troldborg et al., 2022; Dar et al., 2023; Zengler et al., 2024; Revaz and Porté-Agel, 2024) and stresses the insight that, especially in complex terrain, the power performance gain due to terrain speed up ahead of a turbine on a hill might be reduced due to flow deceleration behind the hill.

An obvious question is whether the optimal operational point and thus the optimal thrust and hub height velocity are altered as well. This can be answered within the scope of the simple model by solving for the optimal induction based on Eqs. (20) and (27), which, for a constant l, yields

$$\begin{array}{}\text{(28)}& {\displaystyle }{a}_{\mathrm{opt}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}l\mathit{\beta }-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\sqrt{\mathrm{1}+l\mathit{\beta }+l^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}}.\end{array}$$

In the case of a vanishing acceleration, $a_{\mathrm{opt}}=\frac{\mathrm{1}}{\mathrm{3}}$, which is the classical result from momentum theory. However, for a non-vanishing acceleration, the optimal induction changes with the acceleration and the length scale. This is also shown in Fig. 7. In fact, the optimum induction increases in the case of a positive acceleration and decreases in the case of a negative acceleration. With respect to the previously mentioned situation of a wind turbine operating on a hill, this also means that, in addition to the physical limitation of maximum extractable energy, optimal turbine operation is also altered, which, depending on the turbine controller, may lead to performance degeneration relative to the highest reachable performance.

Figure 8 Extracted length scale l from simulated thrust-induction curves.

5.3 On the length scale

Up to now, it was assumed that the length scale l has a constant value, which was adopted from the literature (Crespo et al., 1999; Dar and Porté-Agel, 2022). However, it is possible to determine it for the presented simulations by solving Eq. (21) for l. In this context, it is important to mention that l does not necessarily need to be interpreted as a physical parameter but can be seen as a free parameter that determines how strong the influence of acceleration on the induction is. Calculating l based on Eq. (21) therefore quantifies how far the model deviates from the simulated physics. This is shown in Fig. 8, where l is plotted as a function of the thrust coefficient. It can be seen that l depends on both the thrust coefficient and the velocity gradient, indicating that both are not sufficiently represented in the present model formulation. It is remarkable that two separate branches of l develop, depending on whether the acceleration is negative or positive. Both branches grow towards larger values of l as the thrust coefficient increases, indicating that within the simplified modeling framework, the influence of flow field acceleration on aerodynamic performance increases as the thrust coefficient increases. At low thrust – ignoring the two separate branches – l is nearly independent of the acceleration. Very close to zero thrust, the general behavior seems to be highly nonlinear, and the two branches tend to diverge from each other, with the thrust coefficient approaching zero. The explanation for this behavior is not known to the authors. Note that this branching effect is not directly noticeable on the a−C_T curves in Fig. 5. It should be emphasized that the length scale estimation close to zero thrust is sensitive to the grid design and numerical errors due to the small numbers involved.

The previously employed value of a reference length of one diameter (corresponding to l=1) is in the range of values l takes within the full range of thrust coefficients and velocity gradients. Considering the prediction results presented earlier, this value can be seen as a compromise between low and high thrust cases.

5.4 Limitations and future work

As very idealized simulations were performed, a natural limitation is that the model still needs to prove its general applicability for more complex flow cases. As such, the atmospheric flow over complex terrain or the flow within a dense wind farm comes to mind, as mentioned earlier. However, in these scenarios, it is likely that the streamwise velocity gradient is not constant. One way to apply the presented model in these cases would be to devise a method to obtain an effective value for the combined quantity lβ from the more generally varying undisturbed background velocity.

As seen in Fig. 7, the optimal point of operation changes when a turbine is subjected to an accelerating flow field. This gives rise to the following question: can existing control strategies cope with this, or would it result in the turbine operating in a suboptimal state with respect to both loads and power performance? This question should be addressed in future works.

Lastly, only scenarios of a comparably low induction (a<0.5) were investigated numerically. It is of interest how the model performs in high-induction scenarios.

6 Conclusions

Streamwise acceleration of the background flow field is usually not considered for wind turbine performance predictions or the design of controllers. In this work, it was shown that this effect can have a significant impact on the relation between induction and thrust and therefore on aerodynamic performance predictions. A simple model, which extends classical momentum theory to the case of a constantly accelerating streamwise flow field, was presented and validated with RANS simulations. The model reduces the average prediction error in comparison with classical momentum theory for thrust coefficients lower than 0.5 by 54 % and for thrust coefficients higher than 0.5 by more than 78 %. The model indicates that the induced velocity for a given thrust depends on the background flow acceleration. The axially induced velocity decreases when the background flow has a positive streamwise acceleration and increases for a negative streamwise acceleration. The model furthermore shows that the maximum power coefficient as well as the corresponding values of optimal induction and thrust coefficients depend on the flow acceleration. This implies that wind turbines on sites with streamwise velocity gradients may not operate optimally under such conditions. Nonetheless, the idealized conditions considered within this work do not represent realistic operating conditions of wind turbines, indicating a need for further research in order to check the applicability of the presented model in practice.

Appendix A: Sensitivity analysis

To ensure that the presented conclusions are independent of the computational settings and geometrical choices in the RANS simulations, the sensitivity of the results with respect to cell size dx, blockage ratio b, diffuser length L_Diff, and radial smearing factor α at a thrust coefficient of $\frac{\mathrm{8}}{\mathrm{9}}$ is analyzed. All results are summarized in Table A1. All sensitivity investigations were carried out in a non-accelerating flow (β=0.0), except in the diffuser length study, which assumed $\mathit{\beta }=-\mathrm{0.004}$.

Table A1 Sensitivity of the simulation results to cell size, blockage ratio, diffuser length, and radial smearing. Bold values are used within this study.

The influence of the cell size is studied by performing two additional simulations, one at a grid that has only half the number of cells in each dimension ($\mathrm{160}\times \mathrm{64}\times \mathrm{64}$) and another that contains twice the number of cells in each dimension ($\mathrm{640}\times \mathrm{256}\times \mathrm{256}$), thus effectively halving or doubling the cells per diameter in the disk region. The first row in Table A1 shows that doubling the cell resolution in the disk region leads to a change in induction of 0.229 %. Based on the grid convergence index, the relative error with respect to the exact solution is estimated to be 0.266 %, with an estimated order of convergence of p=2.831 (Roache, 1994, 1997).

To investigate the blockage effect – and thus the influence of the flow confinement due to the limited domain size – a disk with half the radius of the original disk in a domain with twice the cell resolution is simulated, therefore keeping the cells per diameter constant in the disk region. The blockage ratio is defined as $b=\frac{A}{A_{S,R}}$, with the disk area A and the cross-sectional area in the disk plane A_S,R. The results in Table A1 indicate that reducing the blockage ratio by a factor of 4 leads to an increase in induction of 1.378 %.

Within the derived model, it is assumed that the velocity in the background has a constant axial gradient. In the simulations, it is necessary to limit the axial extent of the acceleration region. To investigate if this has an influence on the simulation results, two additional simulations with a decreasing diffuser length are performed. The reference diffuser is 34.170 D long, while the two shorter ones are 28.475 and 22.780 D long, respectively. The AD is always positioned 11.390 D behind the beginning of the diffuser. In all cases, the mildest negative velocity gradient is simulated. The change in induction, as reported in Table A1, is, for both of the short diffusers, less than 1 %, indicating a negligible influence on the results. One might notice that the variation in induction is lower for the shortest diffuser compared to the intermediate one. A possible explanation for this is that the grid in both cases is slightly different, mainly because the transition regions between the diffuser and the inlet and outlet need additional cell refinement, leading to the grid itself influencing the results here as well.

To improve solution convergence, a radial smearing factor α of the thrust force as described by Zengler et al. (2024) is applied in the AD model. It has a range between 0 and 1, where 0 means that no smearing is applied while 1 results in smearing along the entire radius. In addition to the employed value of 0.25, simulations with 0.15 and 0.35 are conducted. Table A1 shows that the expected variation in induction is on the order of 1 %. It must be noted that changing the smearing factor eventually leads to a variation of the blockage ratio and also of the cell resolution in the disk area, an effect which is not considered here.

Overall, it can be concluded that the results are nearly independent of the cell size and diffuser length, while domain blockage and disk smearing influence the evaluated induction on the order of 1 %. Further, because the simulations results are mainly compared among each other within this study, systematic errors in the simulation setup cancel each other out in these cases.

Appendix B: Wake breakdown

For the cases of highest negative acceleration, it can be observed that a recirculation region develops in the wake at certain thrust coefficients. This is shown in Fig. B1, where the axial velocity along the centerline at $C_{\mathrm{T}}=\frac{\mathrm{8}}{\mathrm{9}}$ is shown. At a certain axial distance, the velocity becomes negative. This is referred to as breakdown of the wake. This observation can be roughly explained by an energy balance between the wake and the surrounding flow. Equation (5) shows that for a fixed pressure jump across the disk, the energy difference between the free stream and wake velocities is constant. When the free stream velocity decreases, the wake velocity also needs to decrease. Because kinetic energy per mass scales with the velocity squared, the fixed energy difference between free stream and wake results in a larger decrease in the lower wake velocity than the free stream. At a certain point, this velocity becomes zero, leading to a breakdown of the wake and the development of a recirculation zone. From Eq. (5), one can therefore derive a criterion for this situation by setting u₁=0. This result can be expressed in terms of C_T as

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{U_{\mathrm{\infty }}\left(x\right)}{U_{\mathrm{\infty },R}}}<\sqrt{C_{\mathrm{T}}}.\end{array}$$

When the normalized free stream velocity reaches a value smaller than the square root of the thrust coefficient, a breakdown of the wake can be expected. At $C_{\mathrm{T}}=\frac{\mathrm{8}}{\mathrm{9}}$, this is the situation for $U_{\mathrm{\infty }}\left(x\right)<\mathrm{0.94}{U}_{\mathrm{\infty },R}$. Comparing this value with Fig. 5, it would be expected that in the two strong deceleration cases, the velocity becomes negative around 1.94 and 2.62 D behind the disk. Simulation results show that in both cases, this happens further downstream than predicted. This finding may be partly explained by the entrainment of momentum from the free stream into the wake. For future modeling, one could consider a high-thrust correction for the wake breakdown due to deceleration.

Figure B1 Normalized centerline velocity for $C_{\mathrm{T}}=\frac{\mathrm{8}}{\mathrm{9}}$ for different background velocity gradients. Gray area indicates negative velocities.