Development and Validation of an Advanced Actuator Line Model for Wind Turbines

As wind turbine technology proceeds towards the development of more advanced and complex machines, modelling tools with fidelity higher than the ubiquitous Blade Element Momentum (BEM) method are needed. Among them, the Actuator Line Method (ALM) stands out in terms of accuracy and computational cost. Moving from this background, an advanced ALM method has been developed within the commercial solver CONVERGE®. As elements of novelty, this tool features a Lagrangian method for sampling the local inflow velocity and a piece-wise smearing function for the force projection process. Various sub-models for both Horizontal Axis Wind Turbines (HAWTs) (e.g. the Shen tip loss correction) and Vertical Axis Wind Turbines (VAWTs) (e.g. the MIT dynamic stall model) has also been included. Aim of the research is to address the new challenges posed by modern machines. HAWTs are in fact getting larger and larger, shifting the research focus on the interaction of increasingly deformable blades with the atmosphere at the microand mesoscale level. VAWTs on the other hand, whose popularity has arisen in the last years, thanks to their advantages in non-conventional applications, e.g. floating offshore installations, are extremely complex machines to study, due to their inherently unsteady aerodynamics. The approach has been validated on selected test cases, i.e. the DTU 10MW turbine and a real 2-blade H-rotor, for which both high-fidelity CFD and experimental data are available.


Introduction
As wind turbine technology proceeds towards the development of more advanced and complex machines, modelling tools with fidelity higher than the ubiquitous Blade Element Momentum (BEM) method are needed in an increasing number of applications. Inter alia, hybrid simulation techniques like the Actuator Line Method (ALM) [1] are being increasingly used, thanks to the much better description of the flow field gained from Computational Fluid Dynamics (CFD). The Actuator Line Method is a hybrid approach developed by Shen and Sørensen [1], which combines a lumped parameters representation of the blade-flow interaction with an unsteady Navier-Stokes equations solver, used to solve the flow field across the turbine rotor and in the wake. In further detail, each blade is replaced by a dynamically equivalent actuator line, introducing aerodynamic forces into the computational grid based on its current position and local flow conditions. Thanks to this strategy, this method does not require the resolution of the boundary layer encompassing the blade surfaces, with a corresponding strong decrease in the required computational cost.
Moving from this background, this study presents an advanced ALM method, that has been developed within the commercial solver CONVERGE ® . This simulation environment, ubiquitously known for its potential in case of combustion processes and confined flows, indeed provides very interesting features that make it attractive also in turbomachines operating in external flows. In particular, the software allows for autonomous meshing, i.e. it automatically creates the mesh at runtime, dynamically adapts the mesh throughout the simulation. This is a key factor to reduce the operator-dependence of the results and to speedup the preparation of each run. Compared to other existing tools of the same kind, the code also features a novel method for sampling the local inflow velocity, based on a Lagrangian time average, as well as an improved smearing function for the insertion of the computed aerodynamic forces into the computational domain (Regularization Kernel) [2]. Within the framework of an ongoing collaboration with academia, important sub-models, both for Horizontal-Axis Wind Turbines (HAWTs) (e.g. the Shen tip loss correction [3]) and Vertical Axis Wind Turbines (VAWTs) (e.g. a customized version of the MIT dynamic stall model [4], specifically tailored for use in this type of rotors) have been implemented. Moreover, a method of interpolating lift and drag coefficients based on the blade airfoils relative thickness has proven to be effective in reproducing blade forces calculated by a full blade-resolved CFD simulation.
Based on these characteristics, the developed software is supposed to be able to address the numerous challenges that arise when simulating the latest generation of wind turbines. In the case of HAWTs for example, the growing size of rotors is shifting the research focus on understanding their behaviour when operating in complex motion and inflow conditions, mainly deriving from the interaction of the increasingly deformable blades with the turbulent structures generated in the atmosphere at the micro-and mesoscale level [5]. On top of that, the increment in wind farm extension and clustering introduces new levels of interaction of a turbine with each other and with the surrounding environment. VAWTs on the other hand, for which renewed interest has arisen in the last years thanks to their advantages in nonconventional applications, e.g. deep-sea, floating offshore installations [6], are intrinsically complex machines to be studied, due to their inherently unsteady rotor aerodynamics.
In order to verify such capability, an extensive validation of the code has been also carried out on selected test cases. In particular, the HAWT version of the code was benchmarked on the ubiquitous DTU 10MW turbine, while for the VAWT one a real 2-blade H-rotor, for which both high-fidelity CFD and experimental data were available, was chosen. Upon comparison of the ALM data with the benchmark ones, the developed model has proven to be a robust and reliable tool in terms of both performance and flow field description.

Actuator Line Model
In order to assess the capability of the proposed approach to simulate different turbine architectures, two test cases, both characterized by the availability of extensive documentation and benchmark data, were selected. In the present work, the Actuator Line Method was implemented in the CONVERGE 3.0 solver via a dedicated User Defined Function (UDF). The corresponding workflow, whose core formulation is the same for both HAWTs and VAWTs, can be synthetically described as follows. For each timestep: (a) the Aerodynamic Center (AC) of each blade element, i.e. the corresponding airfoil quarter-chord point, is located on the grid based on the turbine characteristics and operating condition; (b) the local undisturbed velocity V is extracted from the resolved CFD flow field and is used to evaluate the blade angle of attack α and chord-based Reynolds number Re; (c) the airfoil lift (CL) and drag (CD) coefficients are obtained via linear interpolation of the available tabulated aerodynamic data; (d) the blade element forces are inserted into the computational domain as volume momentum sources; (e) the modified unsteady Navier-Stokes set of equations is solved.
As in every ALM code, the steps (b) the insertion of aerodynamic forces into the computational domain and (d) the sampling of the local inflow velocity from the resolved flow field of the procedure above described are particularly critical for the stability and accuracy of the simulation.

Regularization Kernel
Computed volume forces are commonly plugged into the computational domain after having been smeared out by means of an ad hoc function (Regularization Kernel), as in Figure 1. As a matter of fact, their direct insertion into the mesh cell intersected by the actuator line would cause a singularity in the solution, with the corresponding instability and numerical oscillations issues [7]. As highlighted in many studies [8,9], the proper selection and calibration of such function is pivotal for the robustness and accuracy of the ALM approach.
In the present work in particular, a novel three-dimensional piece-wise δ kernel distribution, introduced by one of the authors in a recent publication [2], was adopted. Upon comparison of this new formulation with the Gaussian kernel commonly used in the majority of ALM tools, as in Figure 2, numerous advantages can be highlighted: • higher stability: at given kernel width β, the proposed piece-wise function is less concentrated in correspondence of the actuator line, with a corresponding improvement in the simulation stability. • reduced truncation error: as demonstrated by Xie [2], the δ function is fully conservative for every projection width, in contrast with the Gaussian one, which needs large β values to reduce the truncation error to acceptable values.

Fig. 1.
Comparison of the δ kernel adopted for the study with the standard Gaussian one.

Velocity Sampling
Similarly to the kernel function, the sampling of the blade inflow velocity exploits a novel approach recently developed by one of the authors [2], which is based on the Langrangian tracking of the fluid particles interacting with the turbine rotor, as schematically illustrated in Figure 2. More in detail, the instantaneous velocity u is averaged over time along the path outlined by each particle: where z(t') is the trajectory at earlier time t' (with t'< t). The time constant T from Eq. 1 defines the memory length of the averaging process and it is scaled on the norm of the strain rate tensor Sij via a user-defined constant C: This technique presents numerous advantages with respect to the other sampling approaches currently available [12], in particular: • no dependence on the sampling point position: since each fluid particle interacting with the rotor is tracked over time, there is no need to define a priori the velocity sampling position • higher flexibility: as the integration interval T scales with the local turbulence intensity, this approach is able to maintain the large-scale fluctuations characterizing the upstream flow, while filtering the instantaneous local effects associated to the presence of the blade forces. The weight of the latter is in this case given by the constant C • no dependence on the grid structure For the present work, C=4 was used for the DTU, C=1 for the H-rotor simulation as a result of a dedicate sensitivity analysis [11]. In order to assess the capability of the proposed approach to simulate different turbine architectures, two test cases, both characterized by the availability of extensive documentation and benchmark data, were selected. For HAWTs, the present study focused on the ubiquitous DTU 10 MW RWT, originally designed by Bak et al. [13] at Denmark Technical University (DTU) with the aim of creating a benchmark configuration for the next generation of large scale off-shore machines. As shown in Figure 3a, the single rotor blade is characterized by a consistent chord reduction along the span, resulting in a local Reynolds variation from Re=6x10 6 at the root to Re=12x10 6 at the tip for the design condition. In virtue of the lightweight design of the rotor, the FFA-W3-xxx airfoils series was selected, corresponding to a relative thickness range between 21.1% and 36.0% (excluding the cylindrical section from hub to root). One particular feature of this blade is the presence of wedge-shaped Gurney flaps on the first 40% of the span, where they are employed as a passive device to alleviate the negative effects of root stall. The selected test case was simulated at different wind speeds V0 Î [4, 6, 8, 12, 16] m/s below and above the rated one (V0, rated=11.4 m/s), at a rotational speed ω = 9.6 rpm.
For VAWTs on the other hand, a real two-blade, H-Shaped Darrieus turbine, tested in different studies [14,15] in the Open Jet Facility (OJF) at the Delft University of Technology, was chosen. Dimensions of the rotor and the blades, which are equipped with a NACA0018 airfoil, are shown in Figure 3b. For this machine, both three-dimensional flow fields [15] and azimuthal blade forces profiles [14] are available from ad hoc post-processing of measured Particle Image Velocimetry (PIV) data. Two functioning points have been considered for the experimental campaign, corresponding to the turbine unstable (TSR=2, V0=10.2 m/s), i.e. characterized by dynamic stall over a major portion of one rotor revolution, and stable (TSR=4.5, V0=9.1 m/s) operation, respectively, in which on the other hand the stall phenomenon is almost absent.

HAWT Model
For the 10MW machine considered in the present study, input polar data came from the CFD investigation performed by DTU on the FFA-W3-xxx airfoils series used to design the turbine blade [13]. The corresponding two-dimensional RANS simulations were carried out for Re Î [6 x 10 6 , 12 x 10 6 ] and an AoA Î [-32°, 32°] . For the modelling of turbulence, the 4-equations γ-Reθ transition model was employed. The obtained data was then extended to the whole post-stall region (AoA=±180°) based on the flat plate assumption. In order to account for the inherently three-dimensional effects occurring at the root and tip of the blade, where the local flow behavior strongly deviates from the two-dimensional one assumed for the Blade Element approach used for the construction of the ALM method, ad hoc corrections were employed in the proposed simulation tool.
At the root, stall delay effects were embedded into the simulation by correcting a priori the input tabulated polar data with the model developed by Bak et al. [16]. In that region in fact, due to the centrifugal pumping exerted by the blade rotation on the local boundary layer, stall resistance is higher than the one predicted by static two-dimensional tests.
At the tip on the other hand, the formation and shedding of coherent vortex-like structures at the extremity of the blade, in virtue of the local pressure gap between pressure and suction sides, leads to a notable load reduction along the blade span. Among the different corrections developed over the years to account for this phenomenon, the one from Shen et al. [3] was selected for the present work. The latter, based on the classical approach by Prandtl, corrects the computed normal and tangential blade force coefficients via the factor F1: with Nb number of blades, r local radius and φ incidence angle. The factors G and n are computed as semi-empirical functions of Nb, the tip-speed ratio, and local chord radial gradient.
It must be noted how, since the ALM already resolves the full three dimensional flow field around the actuator line, the use of a tip loss correction might seem redundant. In many studies [17][18][19] however, it has been observed how the base ALM formulation alone is unable to account for this effect, probably because of the use of the regularization kernel to insert the volume forces into the domain. Therefore, its adoption is justified in this work.

VAWT Model
For the small VAWT tested in the present study, pre-stall aerodynamic coefficients were calculated with the viscid panel method software XFoil [20]. The pre-stall curves were then extended over the deep stall region, i.e. up to AoA=±180°, via the Viterna-Corrigan model [21], as recommended by Marten et al. [22] for VAWT applications. The corresponding dataset is reported in Figure 4 for the lowest and highest Re considered for the analysis.
It must be also noted how, in order to account for the flow curvature phenomena in the simulation campaign, the tabulated data were derived for a virtually cambered airfoil, obtained by proper conformal mapping from the original NACA0018 mounted on the test turbine [23]. In several studies in fact, it has been observed how the cycloidal motion of VAWT blades with respect to the turbine axis of rotation is responsible for the deformation of streamlines around them, increasing their effective degree of curvature with respect to the freestream flow. This reflects on the aerodynamic performance of the blade, which behaves aerodynamically (i.e. in terms of lift, drag and moment coefficients) like a virtually transformed equivalent airfoil, with a camber line defined by its arc of rotation and oriented outwards with respect to the turbine axis, as apparent in Figure 4.  The presence of a proper dynamic stall model, capable of accounting for the oscillating nature of VAWT loads during unstable operation, is pivotal for an accurate simulation of VAWTs, especially at medium-low TSRs [25,26]. In the present work, the modified MIT model developed by Noll and Ham [4] was selected, since it is able to provide a fairly accurate description of the Leading Edge Vortex (LEV) shedding cycle, which is directly related to the airfoil loads in the deep stall region, at a reduced implementation effort. The latter is represented schematically in Figure 5.
In this case, the empirical constant γ, which regulates the sensitivity of the dynamic stall angle αDS to the oscillation reduced frequency, was set to 1.0 rad, as in [24]. For further details about the implementation of the model in an ALM code, please refer to [24].

Numerical Set-up
All three-dimensional simulations were carried out with an unsteady Reynolds-Averaged Navier-Stokes (URANS) approach via the commercial solver CONVERGE®, in the 3.0 release. The PISO algorithm was employed for the resolution of the momentum equations, in a pressure-based formulation. The selection of the turbulence model varied with the simulated machine: k-ε RNG model, which is the default one in the CONVERGE® solver, for HAWTs, k-w SST for VAWTs, according to the guidelines of Balduzzi et al. [27]. For detailed information about the settings for adopted numerical schemes, convergence definition, maximum inner iterations and residuals, please refer to [11,28]. An open-field computational domain was employed, characterized by the far-field boundary conditions commonly adopted for the simulation of wind turbines: freestream velocity V0, turbulence intensity and length scale at the inlet, atmospheric pressure at the outlet, symmetry on the lateral patches. The corresponding shape and dimensions changed according to simulated machine, as shown in Figure 6. For the DTU 10MW, a cylindrical domain was selected, with a length of 12D and a diameter of 10D [11]. The latter also features the surface of the turbine nacelle, which was treated as a no-slip wall. For the Darrieus rotor instead, a more traditional 35D x 20D x 20D rectangular box was used [28].
The discretization of both computational domains was performed exploiting the automatic meshing functionality offered by the CONVERGE® software, which generates a structured, hexahedral mesh, whose cell density is automatically adjusted in each region based on the prescribed refinement level. The resulting grids count 8'748'324 elements for the HAWT and 9'578'382 for the VAWT. An example can be seen in Figure 6 for the DTU 10 MW rotor.

Results
In order to gain a better understanding of the code prediction capability, the results from ALM simulations were compared with the benchmark data available for each class of machine. The analysis focused in particular on the blade loading profiles, in virtue of their relevance for the assessment of rotor both aerodynamic and structural performance. Figure 7 reports the comparison in terms of normal (FN) and tangential (FT) force profiles along the blade span between ALM and the blade-resolved CFD data from Bak et al. [13], which has been taken as a reference for this work.

HAWT
In general terms, the proposed methodology is able to provide a good estimation of the blade loads along the whole span, especially for wind speeds below the rated one (Vrated=11.4 m/s). Only discrepancies can be observed for the tangential load FT at the root (0.05 < r/R < 0.25) and tip (0.95 < r/R < 1.0) of the blade, where the inherently three-dimensional nature of the local flow makes accurate predictions challenging for approaches such as the ALM, which are based on the Blade Element assumption, i.e. that each blade section behaves as a two-dimensional airfoil. In spite of the implemented corrections in fact, the ALM still tends to slightly underestimate the stall delay effect at the root of the blade, whilst, on the contrary, overestimating the blade load reduction in the upper part of the blade due to tip vortexes. Raising the inflow wind speed beyond the rated one, it can be observed how the discrepancy between the ALM prediction and the CFD reference data tends to increase, especially in the mid part of the blade (0.4 < r/R < 0.8). Such phenomenon might be related from a physical point of view to a slight overestimation of the rotor blockage, i.e. the slowing down of the flow across the machine, by the ALM. It must be also noted that there is a small uncertainty about the blade pitch value used for the CFD data, since it is not clearly specified in the corresponding report. This also reflects in the underestimation of the power produced by the rotor beyond the rated point, as shown in Figure 7. Figure 8 shows the comparison in terms of blade loading at the rotor midspan between the ALM, the experimental and the two-dimensional blade-resolved CFD data from [29] over one rotor revolution under unstable operation, i.e. for TSR=2. It must be noted that all forces are here normalized over the freestream momentum qꝏ=0.5ρV0 2 , scaled over the turbine radius R.

VAWT
In average terms, i.e. not considering the oscillations observed in the CFD data due to the shedding of the LEV, the proposed approach is able to provide a reasonable estimation of the blade loads along one rotor revolution. The introduction of the MIT model is in fact related to a more physical representation of the lift and drag behavior under severe dynamic stall conditions. In particular, the model is able to capture not only the stall delay observed in the CFD data for increasing AoA (pitch-up) and corresponding to the state 2 of the MIT formulation, but also the drag overshoot associated to the LEV shedding and downstream convection event (state 3). In the upwind region (0° < θ < 180°) of the rotor revolution, this reflects on a quite accurate estimation of the blade loads trends up to their peak values (θ=60°). For decreasing angle of attack (pitch-down) nonetheless, several anomalies in the predicted tangential force profile can be observed, especially for 125° < θ < 180°, where it goes to zero in a discontinuous way. In the authors' view, this issue might be related to the specific formulation of the adopted MIT model, in particular of the state modelling the flow re-attachment phase (state 4). Since the drag coefficient is semi-empirically computed as CD=CL× tan(α) in fact, its combination with the lift contribution leads inevitably to the nullification of the tangential force value. This is even more evident in the downwind region (180° < θ < 360°), where the zero-torque region extends to ca. 80 degrees of rotation.
Increasing the turbine rotational speed, the effects of dynamic stall on the blade loading profile progressively become negligible, until a fully stable operation is achieved. In this conditions, as clearly visible from Figure 9, the blade flow remains attached for the entire duration of a rotor revolution. In this case, the ALM is able to provide a quite good estimation of the blade normal and tangential forces azimuthal trends, although it tends on one hand to slightly underestimate the magnitude of the upwind torque peak, on the other hand to overestimate the normal load in the downwind region. This is a well-known issue of the ALM, which currently does not properly capture the blockage, i.e. slowing down of the flow due to energy extraction, between the upwind and downwind sectors of the rotor. In the present study, an advanced ALM method for both HAWTs and VAWTs has been developed within the commercial solver CONVERGE®, featuring as novelty a Lagrangian method for sampling the local inflow velocity and a piece-wise smearing function for the force projection process. Various sub-models for both HAWTs (e.g. the Shen tip loss correction) and VAWTs (e.g. the MIT dynamic stall model) has also been included.
From the validation process, it has emerged a good capability of the code to predict the performance, in terms of blade normal and tangential loads, of a large-scale horizontal-axis turbine, i.e. the DTU 10MW turbine, across the entire range operation of the machine. Agreement with reference blade-resolved CFD data is particularly good for wind speeds below the rated one, while a slight discrepancy is observed when the rotor is subjected to pitch regulation. In these conditions, the ALM tends in fact to underestimate blade performance.
When applied to a small 2-blade H-rotor, the proposed tool has been able to provide a quite good estimation of the blade normal and tangential forces trends over one rotor revolution for the turbine stable operation, i.e. TSR=4.5. Lowering the rotational speed down to TSR=2, a relevant accuracy degradation is observed, due to the increasing weight of the dynamic stall phenomenon. The adoption of the MIT dynamic stall model, although notably improving the description of the lift and drag hysteresis loops, introduces some anomalies in the predicted torque history, due to flaws intrinsically related to its formulation.