Lifting Line Free Vortex Wake

In QBlade the aerodynamic forces acting on a rotor can be modeled using the Lifting Line Free Vortex Wake method (LLFVW). Similar to the Blade Element Momentum Method, in the LLFVW model the blade forces are calculated using two dimensional sectional airfoil polar data. The main difference is that the rotor wake, shed from the blades, is explicitly resolved. This is a large improvement over the commonly used Blade Element Momentum Method style approaches, which necessitate the introduction of a large number of empirical corrections into the simulated system. Modeling the wake dynamics explicitly avoids the dependency on such correction models and often lead to more physically sound results. Simulation results are improved especially in cases where the assumptions of the Blade Element Momentum Method are violated. These include unsteady operation, large blade deformations and high tip speed ratios where the turbulent wake state is approached. Such conditions become more and more prevalent with the ongoing trend towards larger rotor sizes and offshore floating wind turbines.

Overview of LLFVW Theory

Blade and wake elements

Fig. 5 Basic elements of the blade and wake model inside the LLFVW algorithm.

The rotor is represented by a lifting line, located at the quarter chord position of the 2D airfoil sections (see Fig. 5). Each blade panel is represented by a vortex ring which consists of four straight vortex filaments. The circulation of the bound vortex lines, forming the lifting line, is calculated from the relative inflow velocity and the lift and drag coefficients that are obtained from tabulated airfoil data. The sectional circulation \(\partial\Gamma\) is calculated according to the Kutta-Joukowski theorem:

\[\begin{align} \partial F_L(\alpha) = \rho V_{rel} \times \partial\Gamma , \end{align}\]

where \(\partial F_L\) is the sectional lift force and \(\rho\) is the fluid density. The relative velocity \(V_{rel}\) is obtained from a simple vector addition of the free stream velocity \(V_\infty\), the blade motion \(V_{mot}\) and the induced velocity \(V_{ind}\), which is calculated from the contribution of all vortex elements on the blade and in the wake through the Biot-Savart equation:

\[\begin{align} V_{ind} = -\frac{1}{4\pi}\int\Gamma\frac{\vec{r}\times d\vec{l}}{r^3} . \end{align}\]
Aerodynamic Flowchart

Fig. 6 Flowchart for a single timestep of the aerodynamic calculations in QBlade

At the beginning of each time step, the circulation distribution along the blade is calculated. This is carried out with an iterative procedure which ensures that the forces predicted by the Kutta-Joukoswki theorem and the blade element theorem coincide. During the iteration only the bound vorticity distribution is updated, while the induction of the wake elements on the blade is only evaluated once. After convergence is obtained, the rotor rotation is advanced for a single time step. All free wake vortex elements are convected with the local inflow and local induced velocity. After the wake convection step, new vortex elements are created between the trailing edge of each blade panel and the last row of wake vortices that were convected away from the trailing edge. As a last step, the circulation is computed and assigned to the new released vortex lines through the Kutta condition:

\[\begin{align} \Gamma_{trail} = \frac{\partial{\Gamma_{bound}}}{\partial x}\Delta x . \end{align}\]
\[\begin{align} \Gamma_{shed} = \frac{\partial{\Gamma_{bound}}}{\partial t}\Delta t . \end{align}\]

Wake Lattice and Connectivity

Wake Lattice

Fig. 7 Visualization of the wake lattice structure with wake nodes and filaments.

Fig. 7 shows the wake lattice structure. Shed- and trailing vortices are interconnected via common vortex nodes. During the free wake convection step the evolution of the wake is evaluated by advancing the positions of the vortex nodes in time. Each newly created vortex node is attached to at least one shed and one trailing vortex filament, thus the total number of vortex nodes is approximately half the number of vortex filaments. Consequently, the Biot-Savart equation has to be evaluated around:

\[\begin{align} N_{nodes}\cdot N_{vortices} \approx \frac{N^2_{vortices}}{2} \end{align}\]

times for a fully populated (assuming that no vortex elements have been removed) infinite wake lattice. Compared to a vortex particle discretization, where no inter-connectivity exists, this means a reduction in computational cost by a factor of 2, due to the inter-connectivity of the wake lattice. To facilitate strategies that reduce the number of free vortices within the wake a method to remove individual vortices from the wake mesh has been implemented whereby vortex filaments are detached from their corresponding nodes. A check is performed during every step of the simulation that removes isolated vortex nodes which are not attached to any vortex filament. The more vortices have been removed from the wake lattice, the lower the aforementioned leverage of the interconnections.

Vortex Core Desingularization


Fig. 8 Velocity distribution around the vortex core.

The Biot-Savart equation exhibits a singularity at the core where \(\vec{r}=0\) (Fig. 8). To prevent this singularity from affecting the stability of the simulation and also to model the viscous core of the bound and free vortices more accurately, a model for a viscous vortex core was implemented. Many different models that describe the tangential velocity distribution around the core exist, such as the Rankine, Lamb-Oseen or Ramasay and Leishman models (see Hommes et al.1). In QBlade a simple cut-off radius is used, which is added to the denominator of this equation in the form of \(r_{c}^2\), and ensures that the induced velocity smoothly approaches zero in the vicinity of the core. This is a computationally efficient implementation as the viscous core modeling is directly implemented in the calculation of the induced velocity. For other vortex models a viscous parameter needs to be evaluated from the relative vortex positions in addition to the Biot-Savart equation. This has a severe effect on the simulation performance, as the evaluation of the viscous parameter is carried \(N^2_{vortices}/2\) times per time step. When shed from the trailing edge of the blade, a vortex is release with an initial core-size \(r_c\) (a value of around 10% of local chord is proposed from experience). The core-size is updated every time step according to:

\[\begin{align} r_c = r_0+\sqrt{\frac{4a\delta_v \nu \Delta t}{1+\epsilon}} . \end{align}\]

where \(a = 1.25643\) is a constant, \(\delta_v\) is the turbulent viscosity coefficient (a value depending on rotor size, see Sant2), \(\nu\) is the kinematic viscosity and \(\Delta t\) the time step size. The strain rate of the vortex filament is computed as:

\[\begin{align} \epsilon = \frac{\Delta l}{l} . \end{align}\]

The desingularized Biot-Savart equation then becomes:

\[\begin{align} V_{ind} = -\frac{1}{4\pi}\int\Gamma\frac{\vec{r}\times\partial\vec{l}}{r^3+r_c^2} . \end{align}\]

T. Hommes, J. Bosschers, and H. W.M. Hoeijmakers. Evaluation of the radial pressure distribution of vortex models and comparison with experimental data. Journal of Physics: Conference Series, 656(1):7–11, 2015. doi:10.1088/1742-6596/656/1/012182.


Tonio Sant. Improving BEM ‐ based Aerodynamic Models in Wind Turbine Design Codes Improving BEM ‐ based Aerodynamic Models Tonio Sant Tonio Sant Improving BEM-based Aerodynamic Models in Wind Turbine Design Codes. TU Delft, 2007. ISBN 978‐99932‐0‐483‐1.