Unsteady BEM

The unsteady BEM accounts for the dynamic variations of the relative velocity at the blade element, which can be produced by
  • a rotating rotor with yaw/tilt/cone angles or in a non-spatially uniform wind (for example with shear)
  • the controller pitching the blades or accelerating/decelerating the rotor
  • the vibration of the blades or the motion of the tower
  • turbulent wind
In Ashes, the dynamic variations of the induced velocity are computed with the Dynamic Wake Model suggested by Øye (1991a). This model consists of a filter on the induced velocity, which calculates the induced velocity at a time step 
 based on the wind conditions at system position at
 and the induced velovity at the previous time step

Note: at t = 0, since no previous time step is available, the induced velocity is calculated with the Steady BEM algorithm even if the Unsteady BEM scheme is selected. The ensuing time steps are computed with the Unsteady BEM

The computation of the induced velocity at a time step 
is done performing the following steps:

1 Calculate the relative velocity and the angle of attack

The relative velocity is taken as the sum of the wind velocity, the rotational velocity of the element and the induced velocity of the previous time step
, noted 
. The angle of attack is then calculated following step 2 of the Steady BEM algorithm, as illustrated in the figure below

2 Compute the aerodynamic coefficients accounting for dynamic stall

First, the Reynolds number is calculated and the aerodynamic coefficients are looked up from the Airfoil polar file, as in step 3 of the Steady BEM algorithm.
Then, the lift coefficient CL is adjusted to account for dynamic stall effects following the work by Øye (1991a).

Note: the dynamic stall model is not relevant for cylinders

In this work, the dynamic stall is modelled through a so called separation function 
such as 
$$C_L(\alpha) = fC_{L,inv}(\alpha) + (1-f)C_{L,fs}(\alpha)$$
 is the lift coefficient for inviscid flow without separation and 
 is the lift coefficient for fully separated flow. In order to solve this equation, we first apply it to static conditions, such as
$$C_L^{st}(\alpha) = f^{st}C_{L,inv}^{st}(\alpha) + (1-f^{st})C_{L,fs}^{st}(\alpha)$$
where the superscript 'st' refers to the static solution.
corresponds to the value looked up in the polar and 
 is obtained by extrapolating the linear region of the lift curve, as illustrated below. Note that by definition, 

As suggested by Hansen et al. (2004f)
 is given by
Note that the theoretical upper limit for 
 is 1. If a higher value is obtained, then 
 is taken. Once 
 is calculated, 
 is computed by assuming that it comes back to the static value following
This forumla can be analytically integrated to give
$$f(t)=f^{st}(t)+f(t-\Delta t)-f^{st}(t)\exp{\frac{-\Delta t}{\tau}}$$
$$\Delta t$$
 is the time step of the simulation and 
$$f(t-\Delta t)$$
 corresponds to the separation function at the previous time step. 
is a time constant taken as
$$\tau = \frac{4c}{|W|}$$
 is the chordlength and 
 is the relative velocity.
 has been determined, the next step to solve equation (1) is to calculate 
. By reorganizing equation (2), we can express 
, then 
 is computed as 
$$C_{L,fs}(\alpha) = C_L^{st}(\alpha)/2$$
 have been calculated, they are insterted in equation (1) to determine the lift coefficient.

3 Calculate the quasi-static induced velocity

Once the aerodynamic coefficients have been calculated and adjusted for dynamic stall, the lift and drag forces are computed as in step 4 of the Steady BEM algorithm, and a quasi-static induced velocity 
is computed from those forces as in step 5 of the Steady BEM.
If requested in the Aerodynamics tab of the Analysis parameters window, this steps will include tip and hub loss correction based on Prandlt's factor and Glauert correction.

4  Compute the dynamic induced velocity

Based on the quasi-static induced velocity, a dynamic induced velocity is calculated following the equations presented in Hansen (2008d). Two time constants are defined such as 


 is the axial induction factor based on the induced velocity at the previous time step and the current wind speed, 
 is the rotor radius, 
 is the magnitude of the wind velocity and 
 is the distance from the current element to the axis of rotation of the rotor.
Note: if 
, the value 
 is used instead

The dynamic induced velocity at the time 
is then given by 

 is an intermediate value of 
 dependent on the quasi-static induced velocity through the following equation:


Note: according to Martin Hansen (author of the book on which these equations are based), the value of
"is empirically calibrated to fit measurements, and maybe this value should be different for the very large WTs being build now." (Jan-2024)

In practice, this equations are solved assuming that their right-hand sides are constant and taking the values at time step 
. The right-hand side of the previous equation is given by
$$H=W_{qs,n}+k\cdot\tau_1\frac{W_{qs,n}-W_{qs,n-1}}{\Delta t}$$

$$\Delta t=t_{n}-t_{n-1}$$
. Then, 
 is obained analytically through
$$W_{int,n}=H+(W_{int,n}-H)\exp{\left(-\frac{\Delta t}{\tau_1}\right)}$$

 is then obtained analytically through
$$W_{i,n}=W_{int,n}+(W_{i,n-1}-W_{int,n})\cdot\exp{\left(-\frac{\Delta t}{\tau_2}\right)}$$


5 Apply the Yaw/Tilt factor

The Yaw/tilt model as described in step 7 of the Steady BEM algorithm is then applied to the dynamic induced velocity.

6 Compute the aerodynamic loads

Once the induced velocity is known, the aerodynamic loads are calculated as in step 4 of the Steady BEM algorithm assuming the aerodynamic coefficients calculated in step 2