Articles | Volume 7, issue 6
Research article
15 Dec 2022
Research article |  | 15 Dec 2022

Comparison of large eddy simulations against measurements from the Lillgrund offshore wind farm

Ishaan Sood, Elliot Simon, Athanasios Vitsas, Bart Blockmans, Gunner C. Larsen, and Johan Meyers

Numerical simulation tools such as large eddy simulations (LESs) have been extensively used in recent years to simulate and analyze turbine–wake interactions within large wind farms. However, to ensure the reliability of the performance and accuracy of such numerical solvers, validation against field measurements is essential. To this end, a measurement campaign is carried out at the Lillgrund offshore wind farm to gather data for the validation of an in-house LES solver. Flow field data are collected from the farm using three long-range WindScanners, along with turbine performance and load measurements from individual turbines. Turbulent inflow conditions are reconstructed from an existing precursor database using a scaling-and-shifting approach in an optimization framework, proposed so that the generated inflow statistics match the measurements. Thus, five different simulation cases are setup, corresponding to five different inflow conditions at the Lillgrund wind farm. Operation of the 48 Siemens 2.3 MW turbines from the Lillgrund wind farm is parameterized in the flow domain using an aeroelastic actuator sector model (AASM). Time-series turbine performance metrics from the simulated cases are compared against field measurements to evaluate the accuracy of the optimization framework, turbine model, and flow solver. In general, results from the numerical solver exhibited a good comparison in terms of the trends in power production, turbine loading, and wake recovery. For four out of the five simulated cases, the total wind farm power error was found to be below 5 %. However, when comparing individual turbine power production, statistical significant errors were observed for 16 % to 84 % of the turbines across the simulated cases, with larger errors being associated with wind directions resulting in configurations with aligned turbines. While the compared flapwise loads in general show a reasonable agreement, errors greater than 100 % were also present in some cases. Larger errors in the wake recovery in the far wake region behind the lidar installed turbines were also observed. An analysis of the observed errors reveals the need for an improved controller implementation, improvement in representing meso-scale effects, and possibly a finer simulation grid for capturing the smaller scales of wake turbulence.

1 Introduction

Recent years have seen the emergence of wind farm simulation tools that cover the whole chain from flow-coupled aeroelastic models to power-grid models. The complexity of these models ranges from analytical tools, which simplify wake expansion and merging, to complex computational fluid dynamics (CFD) solvers, which represent the turbines and their influence on the surrounding flow field. Amongst all of these numerical tools, large eddy simulations (LESs) have been extensively used in recent years to resolve detailed representation of the turbulent flow in and around large wind farms (Calaf et al.2010; Lu and Porté-Agel2011; Archer et al.2013; Martínez-Tossas et al.2015; Ghaisas and Archer2016; Munters and Meyers2018; Lin and Porté-Agel2019). This increased detail in simulating the physics governing wind farm flows has facilitated the study of wind farm aerodynamics and enabled the analysis of phenomena like turbine–wake interactions, gusts, atmospheric stratification, and the effect of wind farms on local wind climate. A comprehensive review of LES for wind farm simulations addressing the aforementioned effects can be found in the references (Mehta et al.2014; Porté-agel et al.2020). Additionally, LES has also been used to investigate and develop coordinated wind farm control strategies, which could provide the benefits of power maximization, asset life extension, and grid frequency regulation, thus improving the performance and capabilities of wind farms (Goit and Meyers2015; Yılmaz and Meyers2018; Bossanyi2018; Boersma et al.2019; Frederik et al.2020). However, to give credibility to these studies it is essential to validate numerical solvers against reliable measurement data. While wind-tunnel experiments provide a useful avenue for testing, their accuracy in representing full-scale wind farms is limited due to the size and measurement constraints of wind tunnels (Bastankhah and Porté-Agel2017). Therefore, proper validation of wind farm numerical models requires accurate reference data in the form of detailed flow field and performance measurements from existing wind farms. To this end, a measurement campaign was carried out at the Lillgrund wind farm, located 10 km off the coast of southern Sweden, as part of the Horizon 2020 TotalControl project. The measurement campaign made use of three long-range lidars, which measure the inflow conditions for the farm while also resolving the flow field in a part of the Lillgrund wind farm for wake measurements. The flow field data were supplemented by simultaneous power and structural-load measurements from individual wind turbines. The combination of lidar inflow field data, turbine performance data, wake data, and loading data provide a unique data set for the validation of coupled flow and aeroelastic solvers.

In this work, SP-Wind, an in-house aeroelastic LES solver, which has previously been used extensively for wind farm modeling and control optimization, is used to simulate the operation of the Lillgrund wind farm during the measurement campaign. While previous wind farm validation studies have been carried out using LES (Wu and Porté-Agel2011, 2013; Nilsson et al.2014; Wu and Porté-Agel2015; Draper et al.2016; Simisiroglou et al.2018), this work differs on three fronts. First, the atmospheric conditions at the Lillgrund site are recreated in the numerical domain by analyzing inflow lidar measurements. Second, instead of initializing the flow field from scratch, as is the convention for LES precursor simulations (Stevens et al.2014), the current work proposes a framework for reusing a previously generated precursor flow database for matching the conditions during the measurement campaign, substantially reducing the associated computational costs and time for LES wind farm validation studies. Third, the current study utilizes a novel aeroelastic actuator sector model (AASM) for parameterizing the turbine forces in the numerical domain (Vitsas and Meyers2016). Compared to other actuator models such as the actuator disk model (ADM) and the actuator line model (ALM), the AASM has the advantage of accurately representing rotating turbine blades while allowing for coarser time steps through spatial and temporal filtering and decoupling of the LES time step and the time step of the flexible multibody model that is a part of the AASM. The present article is organized as follows: Sect. 2 details the specifics of the measurement campaign and available data sets. Section 3 presents the specifications of the numerical solver used in this study, while Sect. 4 outlines the optimization framework developed to recreate the atmospheric conditions at Lillgrund using a previously generated precursor data set. A comparison of individual turbine performance results from the numerical solver against field measurements and a wake deficit analysis are presented in Sect. 5. Finally, Sect. 6 outlines a summary of the validation and the challenges associated with LES validation studies of wind farms.

2 Lillgrund offshore measurement campaign

2.1 Lillgrund offshore wind farm

The Lillgrund offshore wind farm is located approximately 10 km off the coast of southern Sweden, just south of the Öresund Bridge, where average wind speeds are overall close to 8.5 ms−1 (Sebastiani et al.2021). The wind farm contains 48 wind turbines (Siemens SWT-2.3-93) with a total capacity of 110 megawatts (MW). The farm's turbines have a rotor diameter of 93 m, hub height of 65 m, and a tip height of 115 m. The farm is known to suffer from performance losses, as the turbines originally intended for the farm were replaced by larger models leading to a tighter layout when normalized by turbine diameter (Simisiroglou et al.2018).

2.2 Lidar measurements

During the measurement campaign from September 2019 to February 2020, three long-range WindScanners, i.e., pulsed scanning Doppler wind lidars (Vasiljević et al.2016), were installed on the Lillgrund wind turbine transition pieces and used to measure the flow field both upstream and within the farm, in the layout shown in Fig. 1.

Figure 1Top-down view of the Lillgrund site and scan areas for inflow and wake measuring lidars installed on turbines B08 “Vara”, A07 “Levante”, and C07 “Sterenn”. Coordinates are relative to the turbine B08 position, spaced in multiples of five rotor diameters. Vara performs PPI sector scans which are reconstructed into wind vectors along the magenta inflow centerline. Sterenn and Levante perform coordinated dual-Doppler complex trajectory scans which provide time and space synchronized measurements at three heights along the wake transect indicated with black crosses. The overlapping area between the areas scanned in red by Levante and blue by Sterenn are shaded in green. The positions within the green area which lie off the synchronized wake transect line are resolved through time averaging and applying dual-Doppler retrieval to the 10 min averaged wind field.


2.2.1 Inflow lidar

The inflow measuring system on turbine B08, “Vara”, performed repeating plan position indicator (PPI) sector scans with a constant elevation angle of 8 and azimuth sweep of 60. The center line of the arc scan was intended to lie parallel with the B-row of turbines; however no hard targets were visible from its install location to allow for a precise alignment. Instead, the system was coarsely aligned during installation and, subsequently, the static misalignment was determined and later corrected for using the turbine's calibrated nacelle direction signal. The wake-scan positions (red and blue areas) were determined through an optimization where the lidar positions and points along the wake transect were first defined, and a trajectory for each lidar was generated which matched the position and timing of both lidars, along with kinematic controls to control the scan speed, motor acceleration, etc. The scanned areas are mostly symmetric but not exactly due to the actual installed positions and wake transect point locations chosen.

The inflow lidar data were processed firstly by removing periods of low-signal-quality, i.e periods with carrier-to-noise ratio (CNR) below 26 dB. Partial scans with a low proportion (< 80 %) of valid radial speed values were also filtered out. The remaining valid data were reconstructed into two-component horizontal wind vectors using the integrating velocity–azimuth process (iVAP) method (Liang2007). Lastly, the direction misalignment due to imperfect installation of the system was corrected by determining the static offset between the nacelle direction measurements from turbine B08 and lidar reconstructed wind direction. The static direction misalignment was found to be 20.47. This method does not account for the fact that the turbine yaw controller has a delay in wind direction tracking due to actuator limits. Taken on average over the entire campaign however, it provides the best estimation possible when a hard target alignment is not feasible to perform. Note that this correction has no effect on the wind speed result from the lidar measurements, only the resolved wind direction.

A time series of processed wind speed and direction inflow measurements corresponding to the turbines' hub height (65 m) is shown in Fig. 2, with the equivalent wind rose shown in Fig. 3. This location corresponds to the lidar range gate at a distance of 430 m, given the lidar's inclined beam and origin at 5 m above mean-sea-level. Range gates between 70 and 1490 m were sampled in steps of 20 m.

Figure 2Wind-speed (a) and direction (b) measurements at hub height of 65 m processed from the inflow lidar mounted on turbine B08. Missing data correspond to equipment downtime and filtering of low signal quality periods.


Figure 3Wind rose depicting hub height wind speed and direction from the inflow measuring lidar.


2.2.2 Wake scanning lidars

Two additional scanning lidars identical to the inflow measuring system were installed within the wind farm to measure wake effects and intra-farm flows. The two wake scanning lidars designated as “Levante” on A07 and “Sterenn” on C07 performed coordinated dual-Doppler complex trajectory scans within the green overlapping area shown in Fig. 1. The points placed along the wake transect path indicated with black crosses were repeated at three heights (6, 70, 130 m a.s.l.) and were time and space synchronized using the DTU WindScanner software. Overlapping scan areas which lie away from the wake transect path are not time–space synchronized but have been averaged over a 10 min period to produce a 3D wind field. From this, a horizontal slice has been taken to obtain the 2D wind field at constant hub height. The dual-Doppler wind retrieval method used in this study follows the process outlined in Simon and Courtney (2016).

2.3 SCADA data

Data from the wind farm monitoring system were provided by Vattenfall. This included the following channels from all 48 turbines in the wind farm: active power, blade pitch angles (A/B/C), nacelle direction, rotor speed, and wind speed. The raw streaming data format did not maintain a specific time resolution, where deviations in the sampling rate existed between channels and over time. The data set was processed into a constant sampling rate of 0.5 Hz time steps across all channels. This period was chosen as it was the highest time resolution present in the raw data files. The periods and channels with lower time resolutions were upsampled (interpolated) to match. The raw nacelle direction signals were found to have differing offsets and were then calibrated on a per-turbine basis and corrected in the final data set. Given the fact that met mast data were not available, the calibration was performed using sets of turbines each consisting of a front turbine and the corresponding closest downstream turbine. The inflow direction generating the largest wake loss at the downstream turbine location corresponds to the direction from the front turbine to the downstream turbine, whereby the offset of the measured nacelle direction at the front turbine follows directly. To cover the complete wind rose, four direction sectors of 90 were defined, and within each of these the nacelle direction calibration was performed using 2–3 sets of turbines.

2.4 Load data

Six wind turbines (designations B06, B07, B08, C08, D07, and D08) have been outfitted in the past with load measuring equipment and could be used during the TotalControl measurement campaign. These data are available starting November 2019 at 10 Hz sampling frequency. Strain gauge measurements are available both at the blade root (installed 1.5 m from blade root) and tower base (installed 8.52 m from tower base). The tower base position includes two sensors oriented 90 apart. The sensors were installed and calibrated by Siemens Gamesa.

3 Numerical solver

3.1 Fluid solver

SP-Wind is a wind farm large eddy simulation code built on a high-order flow solver developed over the last 15 years at KU Leuven (Calaf et al.2010; Allaerts and Meyers2015; Munters and Meyers2018). The three-dimensional, unsteady, and spatially filtered Navier–Stokes momentum and temperature equations


are solved. In these equations, ũ=[u1̃,u2̃,u3̃] is the filtered velocity field. Further, θ̃ is the filtered potential temperature field, and θ0 is the background adiabatic base state. The pressure gradient is split into a mean background pressure gradient p driving the mean flow and a fluctuating component p̃. The very high Reynolds numbers in the atmospheric boundary-layer flow combined with typical spatial resolutions in LES justify the omission of resolved effects of viscous momentum transfer and diffusive heat transfer. Instead, these are represented by modeling the subgrid-scale stress tensor τs and the subgrid-scale heat flux qs originating from spatially filtering the original governing equations (Allaerts and Meyers2015). Coriolis effects are included through the angular velocity vector ω=Ωsin ϕ, where Ω is the earth's rotation and ϕ is the latitude of the wind farm. Thermal buoyancy is represented by g(θ̃-θ0)/θ0, with g the gravitational acceleration, θ̃ the filtered potential temperature, and θ0 a reference temperature. The effect of surface friction is included using a standard wall-stress model, corresponding to a logarithmic velocity profile with a roughness length z0 (Bou-Zeid et al.2005) to represent marine boundary layers, without the need to resolve wave effects on the mesh. Finally, Ftot represents body forces added to the flow and consists of F and FR, which represent forces exerted by the wind turbines on the flow and a fringe-region approach to introduce inflow boundary conditions in the domain respectively.

Spatial discretization is performed in the horizontal and span-wise directions by using pseudo-spectral schemes, while vertical fourth-order energy-conservative finite differences are used in the vertical direction. The equations are marched in time using an explicit fourth-order Runge–Kutta scheme, and grid partitioning is achieved through a scalable pencil decomposition approach. Subgrid-scale stresses are modeled with a standard Smagorinsky model with Mason and Thomson wall damping (Allaerts and Meyers2015).

3.2 Structural solver

Deformation of the turbine blade and tower is employed by a finite-element floating frame of reference formulation (Shabana2013). Each element is described by reference coordinates which specify its position and orientation and elastic coordinates that define its deformation with respect to the body coordinate system. Bryant angles are used to describe the orientation of the rotor's body reference frame; however, only the rotation of the turbine rotor is assumed to contribute to its dynamic behavior. Deformations along the tilting, yawing, and pre-coned axis are taken into account quasi-statically. The governing equation for the system can be written as (Shabana2013)


where M, C, and K are the mass, damping, and stiffness matrices, respectively, computed using the structural specifications of the Siemens SWT-2.3-93 turbines. The vector q represents the generalized coordinates, while q˙ and q¨ represent their first and second time derivatives. Φq and λ are the constrain Jacobian matrix and Lagrange multipliers, respectively, and Qg represents the gravitational loads acting on the rotor and tower elements. The vector Qa=FArtrFAtwr contains the aerodynamic loads evaluated at the rotor and tower nodes, as described in Sect. 3.3. Finally, Qv is composed of the Coriolis and gyroscopic loads (Shabana2013). Further details regarding the coordinate system used and the derivation of the equations of motion is given in Appendix A.

To solve the equations of motion, first an eigenvalue problem is solved without damping and external loading to extract the mode shapes and natural frequencies of the structure. Then, the order of the system is reduced by a common modal transformation technique. Hence, the rotor blades are represented by six modes (two flap wise, two edgewise, one torsional, and one axial) and the tower is represented by four modes (two side to side, two fore–aft). Finally, the reduced order system is integrated in time by using the generalized-α method, with a time step of 0.01 s and spectral radius of 0.9 for low numerical damping (Arnold et al.2007).

3.3 Turbine model

As it is computationally prohibitive to fully resolve a wind turbine structure and its forces, actuator methods have been extensively used in research to parameterize wind turbine forces onto the flow grid. The most widely used amongst these models is the actuator disk model (ADM) and the actuator line model (ALM) (Churchfield et al.2017). The ADM is a simple representation of a wind turbine in the flow domain, in which a disk is used to parameterize the turbine forces on the flow. Though simple, the ADM does not accurately represent wind turbine operation due to a lack of discrete rotating turbine blades, and hence often requires additional tuning to improve performance. Comparatively, the ALM provides a more accurate parametrization by representing the turbine blades as actuator lines, with each point exerting a force on the flow based on its local inflow velocity and airfoil distribution. However, the ALM suffers from a limitation that the movement of the actuator line tip over a time step is limited to the size of a cell, thus requiring very fine time steps and hence high computational costs. To overcome this limitation, an actuator sector model (ASM) was developed, which swept the rotor forces across a sector area and hence allowed for coarser time steps. The ASM therefore represents an intermediary between the ALM and ADM, with the ability to resolve structures in the near wake with greater detail than the ADM due to the presence of rotating actuator lines while avoiding the high computational cost associated with the fine time steps of the ALM (Storey et al.2015). The ASM was later extended to an aeroelastic actuator sector model (AASM) by incorporating two-way fluid structure interaction (FSI) coupling to account for structural deformations (Vitsas and Meyers2016).

In this work, the Siemens 2.3 MW turbines are modeled by using the AASM coupled with the nonlinear flexible multibody dynamics model described in the previous section. Since the LES computations are more intensive than integrating the structural equations, a sub-cycling process is employed, for which the aeroelastic coupling scheme is shown in Fig. 4. The relative velocity Vrel is evaluated at each airfoil element along the blade based on the induced velocity field u^x at the airfoil’s deflected position and on the blade’s out-of-plane and in-plane motion represented by q˙OoP and q˙iP, respectively. The relative velocity Vrel also includes the effect of rotor angles (yaw, tilt, and precone) by using rotation matrices to transform the incoming flow field. The flow angle ϕ is then determined from the involved velocity triangle, which comprises the pitch angle β, the torsional deflection τ, and the angle of attack α. The lift L, drag D, and pitching moment M at each airfoil section is then determined using 2D airfoil look-up tables and used to evaluate the aerodynamic forces FAm comprising of the normal, tangential, and span-wise forces FNm, FTm, and FSm. This is done at every sub-cycle m and serves as an input to the equations of motion given by Eq. (3), which are subsequently solved at each sub-step. The forces are then spatially and temporally filtered to obtain the body forces F, which serve as an input to the flow solver in Eqs. (1). Further details of the coupling of the turbine model, flow solver, and filtering are given in Appendix B.

Figure 4The AASM comprises the following steps: (a) evaluation of the angle of attack from the airfoil’s cross-section velocity triangle and blade’s motion, (b) the local cross-section forces are computed from 2D airfoil data, (c) the blades sweep a sector area using a sub-cycling scheme, and (d) the sector area forces are time-filtered. Panel (d) denotes that more weight is given to the forces in the end position. Figure taken from Vitsas and Meyers (2016).

The pitch and rotational speed of all the turbines are controlled using an implementation of the DTU wind energy controller (Hansen and Henriksen2013). A comparison of the simulated power output and thrust in SP-Wind using AASM under the influence of a range of uniform velocities against reference data for the Siemens 2.3 MW turbines is shown in Fig. 5. Slight differences can be seen in both the simulated power and thrust, which can be attributed to the coarse grid resolution across the turbine blades in SP-Wind, which was restricted due to the resolution of the precursor simulations used as they were originally developed for a bigger 10 MW turbine and a larger wind farm. Using the same grid resolution for the smaller 2.3 MW turbines at Lillgrund leads to differences in induction at the rotor plane.

Figure 5Comparison of simulated power and thrust coefficient in SP-Wind using the AASM against reference data (Göçmen and Giebel2016).

Figure 6Flow profiles for PDBL cases. (a) Mean velocity. Dashed lines indicate log-law profiles. (b) Resolved Reynolds shear stress and turbulent kinetic energy. Markers are plotted every 10 data points for clarity.


Figure 7Flow profiles for CNBL cases. Top left: horizontal velocity. Top right: total (resolved + subgrid) shear stress. Bottom left: wind veer. Bottom right: potential temperature. Markers are plotted every 10 data points for clarity. Figure taken from Sood et al. (2020a).

4 Recreating the inflow conditions at Lillgrund

4.1 Precursor database

The turbulent inflow conditions for wind farm inflow are obtained from the publicly available precursor data from the TotalControl Flow database (Munters et al.2019a, b, c, d). The precursor data contain unsteady three-dimensional flow data of an unperturbed atmospheric boundary layer (i.e., without the influence of turbines). The database comprises two pressure-driven boundary layers (PDBLs) and three conventionally neutral boundary layers (CNBLs), spanning different surface roughness lengths and boundary layer heights. The TotalControl flow fields are initialized using mean velocity profiles upon which random divergence-free perturbations are added. These initial conditions are then advanced in time for 20 physical hours, so that the influence of the unphysical perturbations has disappeared and the flow has reached a fully turbulent and statistically stationary state. The PDBL is a simple representation of a neutral atmospheric boundary layer which ignores the effects of rotation and thermal stratification, while the CNBL provides a more realistic representation by including these effects. Both these boundary layer types have been extensively used in wind farm LES. Specifications of the five different boundary layers are given in Table 1, while their flow profiles are shown in Figs. 6 and 7. The database has previously been used in a study to determine the effect of CNBL height on wind farm performance (Sood et al.2020b). All the boundary layers have been initialized using a friction velocity of 0.28 ms−1, which is a typical value for marine boundary layers (Brost et al.1982). For the PDBL, this requires a driving pressure gradient of p/ρ=-5.2267×10-5ms-2 according to the equation u*=-Hρp1/2. Further information regarding the initialization of the boundary layers and the precursor database is available in the provided references (Anderson et al.2020; Allaerts and Meyers2015). A stream-wise slab of the velocity and temperature field was stored to disk when running the precursor, and is later introduced in the wind farm domain by means of body forces in a so-called fringe region (Stevens et al.2014; Munters et al.2016). To match the inflow conditions measured by the lidar measurement campaign, the data from the precursor data set can be transformed to different flow conditions by re-scaling and shifting the flow variables. Using friction velocity u* for velocity scaling, different wind speeds can be attained by re-scaling the entire flow field by a different target friction velocity u*t. This is possible for offshore wind farms, as the solution is scale invariant at high Reynolds numbers. Additionally, in line with the classical outer layer similarity hypothesis (Townsend1976), for offshore atmospheric boundary layers at high Reynolds numbers the roughness elements are much smaller than the boundary-layer height, and hence the roughness acts merely to increase surface stress without any structural changes in the flow (Castro2007; Jiménez2004). The effect of a different target roughness lengths z0t can thus be imposed by applying an offset on the mean flow in line with the difference in surface roughness. Hence, denoting the imposed reference friction velocity and roughness length in the current cases by u*r and z0r respectively, the flow can be re-scaled and shifted as

(5) u t ( x , t ) = u * t u r ( x , t ) u * r + 1 κ ln z 0 r z 0 t e 1 ,

where κ is the Von Karman constant and e1 is a unit vector in the positive x direction, symbolizing that the flow shifting is applied to the velocity u1̃ in the stream-wise direction. The reference velocity field which is transformed, ur(x,t), consists of the stream-wise and span-wise velocities at time instant t. The same friction velocity and surface roughness is used for both the velocity components. Scaling the velocity also leads to scaling of the timescales, according to the following equation:

(6) Δ t t = Δ t u * t u * r .

For the CNBL data sets, a similarity parameter h*=|fc|h/u* relates the height of the atmospheric boundary layer with the Rossby–Montgomery scale u*/|fc| (Arya1975, 1978). Furthermore, an empirical relation for the height of the boundary layer with the strength of the capping inversion is given as

(7) h = A θ 0 g Δ θ u * 2 ,

where A is an empirical parameter with the value 500 and Δθ is the strength of the capping inversion (Csanady1974). Based on the above relations, re-scaling via a new u* implies a change in latitude and θ0/Δθ, if h is to remain constant.

Table 1Specifications of the TotalControl flow database (Anderson et al.2020; Munters et al.2019a, b, c, d).

Download Print Version | Download XLSX

4.2 Optimization framework

While SP-Wind does support changing wind directions during a simulation run (Munters et al.2016), each simulation is restricted to a single wind direction θ and a time frame of 75 min to limit computational costs. Thus, the available inflow data from the measurement campaign were divided into numerous 75 min overlapping time windows. We define P as one of the LES data sets from the five available precursor cases listed in Table 1 and Q as a data set consisting of a 75 min time block extracted from the available lidar database, with both P and Q containing only the stream-wise and span-wise velocity time series as the lidar measurements did not include vertical velocity measurements. Additionally, the velocities in the LES data sets are extracted from the dimensional TotalControl inflow database domain at the equivalent range gate locations as the lidar field measurements. Rows in the two data sets represent range gate locations from the measurement campaign, while columns represent time-series data. While various methods exist to compare the similarity between multi-dimensional time-series data (Salarpour and Khotanlou2018), we limit our analysis to two simple distance metrics to determine the differences between the lidar measurements and the available LES precursors. These metrics are defined as


The first metric d1 is a Euclidean norm which provides a measure of the differences between the time-averaged profiles at range gate locations with the over bar denoting 75 min time averages. The second metric d2, a Frobenius norm, is a measure of the differences between the co-variances of the two data sets, accounting for spatial correlations and also providing a measure for the differences in turbulence intensities across the rotor area. Without the inclusion of the covariance distance, it was observed that while the mean velocity profiles in obtained solutions matched well, they had large differences in the turbulence intensities. The metrics are also assigned weights according to a vector w as shown in Fig. 8a, giving highest preference to the range gates spanning the rotor area 20 % around the hub height, followed by the remaining rotor area and finally the rest of the vertical domain. As discussed in Sect. 4.1, the available LES data sets can be modified for each case P by using the scaling-and-shifting parameters u*t/u*r and z0t/z0r, henceforth collectively referred to as the transformation vector ζ=[u*t,z0t], to obtain a new LES realization denoted by Pζ. Therefore, the distance metrics d1(Pζ,Q) and d2(Pζ,Q) can be determined over the entire 6 month measurement campaign between each LES data set P and a lidar data set Q, for different combinations of the transformation parameters ζ. The sum of both the metrics can be used as a measure of similarity between the LES and lidar data, with lower values indicating greater similarity. Thus for each data set P, a minimization problem can be defined to determine the transformation vector ζ and the time window Q from the measurement campaign which returns the least distance between the LES data and lidar measurements, indicating highest similarity. The cost function of the optimization problem can be defined as

(10) min ζ , Q d 1 ( P , Q , ζ ) + d 2 ( P , Q , ζ )

Figure 8(a) Height varying weights used for the distance metrics in the minimization problem. (b) Locations of 72 range gates, with range gate 1 being the lowest and closest to the turbine B08 and 72 being the highest and furthest away.


and is solved using the SLSQP solver from the SciPy Python package (Virtanen et al.2020). As the distances are calculated for both the stream-wise and span-wise velocity components, the wind direction error is implicitly accounted for and minimized. In its current form, both the distance metrics are given equal weights in the minimization problem, however future work could explore a Pareto front for determining optimal weights for the two metrics. After sweeping through the entire measurement campaign, five unique time windows of 75 min length each, corresponding to five different LES flow realizations which best matched the lidar data are obtained. The first three matches are obtained by transforming the PDk TotalControl LES data set, while the fourth and fifth matches are obtained from the CNk4 and CNk8 data set, respectively, without any transformations. While additional matches were identified, which also contained flow realizations obtained from the PDkhi and CNk2 LES data sets, we restrict further analysis to the best five cases with highest similarity due to computational limitations. A comparison of the mean vertical profiles at range gate location for these cases are shown in Fig. 9 and their specifications are outlined in Table 2. A comparison of the average turbulence intensity across the range gates spanning the rotor area is also presented in Table 3.

Figure 9Comparison of vertical mean total velocity profile (a) and wind veer (b) between the selected LES and corresponding lidar data. Shaded area represents 95 % confidence intervals on the mean evaluated using the block bootstrap method.


Table 2Specifications of selected validation cases.

Download Print Version | Download XLSX

Table 3Turbulence intensity (TI) comparison of validation cases, averaged over the rotor area. Error estimates represent 95 % confidence interval and are calculated using the block bootstrap method with block length of 600 and 1000 bootstrap iterations.

Download Print Version | Download XLSX

From Fig. 9, it can be seen that majority of the selected LES cases match the lidar measurements within error limits across the rotor area. As the sample size of comparison is limited, bootstrapping is used to determine the measure of accuracy of the computed means. Since a traditional bootstrap approach of randomly re-sampling the original time-series data is inappropriate for time series with intrinsic correlation, the moving block bootstrap method is utilized (Kunsch1989). The length of individual blocks in the moving bootstrap method was set to 10 min as a compromise between having enough bootstrap blocks from the measured time-series data and keeping the block lengths large enough to ensure that individual blocks can be assumed to be independent from each other. Through a sensitivity study, it was determined that 1000 bootstrap iterations were sufficient to obtain converged uncertainty estimates for the mean. Larger deviations in total mean wind speed can be observed at heights above the rotor tip, which can be attributed to the preference given to the hub height and rotor disk area through weights in the optimization problem. The PDBL simulations have the largest error when comparing the mean veer between the LES and lidar measurements as, by definition, the PDBL simulations have zero veer and are hence incapable of representing a veered flow. While the CNBL simulations do include veer, the TotalControl precursor database was not designed to cover large veer conditions, thus still leading to errors when compared to the veer in measurement data. Nevertheless, the absolute error never exceeds more than 7 over the rotor area for any of the cases at a given range gate. Good comparison is also seen between the measured and simulated turbulence intensity, with the maximum error never exceeding 1 %.

4.3 Sources of error

The developed optimization framework is therefore able to successfully recreate the inflow conditions at the Lillgrund wind farm for the chosen time periods without significant errors in the mean inflow velocity across the rotor area, as evident by Fig. 9. However, a few shortcomings of the optimization framework must be addressed. Atmospheric stability and stratification can play a significant role in the performance of a wind farm by affecting wake recovery and ambient turbulent intensity (Magnusson and Smedman1994). Additionally, due to the close land proximity of the Lillgrund wind farm, the atmospheric conditions at the farm can be significantly affected by a land–sea system as a combination of marine and land boundary layers. The current methodology does not take the meso-scale effects at the Lillgrund site into account, instead relying on a database comprising of marine pressure driven and neutral boundary layers to match the inflow conditions. This is due to the lack of atmospheric stability and temperature measurements at the Lillgrund site during the measurement campaign. If the data were available, the inflow database used could have been modified before the transformation procedure to consist of precursor boundary layers which best represent the stability conditions at the wind farm site. In its current form, the considered cases may be biased as they are chosen simply based on minimum error in the distance metrics and do not account for stability or land–sea effects. Additionally, without the inclusion of appropriate weights or additional distance metrics to penalize wind veer error, the current procedure may have been biased to overemphasize mean speed error as compared to wind veer error. The error in the lidar and LES velocity profiles above the rotor area could also affect the vertical transport into the wind turbine array atmospheric boundary layer (WTABL), which may affect the overall wind farm performance (Calaf et al.2010).

5 Large eddy simulations of the Lillgrund wind farm

5.1 Numerical setup

The simulation domain in SP-Wind has a size of 16×16×1.5km3 in the stream-wise, span-wise, and vertical directions, respectively. The grid resolution is 13.33×13.33×6.66m3, resulting in a computational grid of 1200×1200×225=324×106 grid points. The choice of domain size is restricted by the one used for the precursor data sets, which was initially designed for simulations with the much larger TotalControl reference wind farm to avoid blockage effects (Sood et al.2020b). Wind farm simulations in SP-Wind are performed in a sequence of steps via the concurrent precursor framework (Munters et al.2016). First, the inflows from the previously generated TotalControl precursor database are made to advance in time in a domain without wind turbines, called the precursor domain. Concurrently, the flow is transformed using the identified transformation parameters through the optimization framework to obtain the five cases identified in Table 2 and fed into a second domain which contains wind turbines as body forces FR through a fringe region. For each of the five cases, the Lillgrund wind farm is rotated to simulate different wind directions. The flow is allowed to pass through the wind farm for 1800 time steps to account for start-up transients, after which data are collected for evaluating the performance of the farm for 9000 time steps. While the original precursor inflow database have a LES time step ΔtLES=0.5 s, the time step of the wind farm simulations is altered due to the scaling of the velocity field, as per Eq.  (6). The multibody aeroelastic computations are performed with a smaller time step of ΔtMB for a finer resolution of the structural loading. The general domain and time parameters of the simulations are summarized in Table 4.

Table 4Summary of the general domain parameters.

Download Print Version | Download XLSX

5.2 Time-averaged flow fields

Time-averaged hub height flow fields for all the selected validation cases are shown in Fig. 10. It can be seen that out of all the cases CNk41 has the highest inflow velocity and PDk3 has the lowest, in accordance with results shown in Fig. 9. The different wind directions spanning the five cases lead to different operation states of the same turbines within the Lillgrund wind farm, due to changes in available hub height wind speed but also due to individual turbines operating in a waked or un-waked condition as per the orientation of the upstream turbines. For instance, the cases PDk1 and CNk81 with wind directions 119 and 222, respectively, have a larger number of turbines operating under fully waked conditions compared to the other three cases. This leads to a data set which allows us to evaluate the performance of Lillgrund turbines when subjected to varying operating conditions.

Table 5Comparison of percentage of turbines with statistically significant errors based on 95 % confidence intervals between LES simulations and field measurements.

Download Print Version | Download XLSX

Figure 10Time-averaged stream-wise hub height velocity for all selected cases. Different wind directions are realized by rotating the entire wind farm in the simulation domain and feeding the inflow velocity in the horizontal-x direction.


5.3 Performance comparison

Comparison of the mean wind farm power output obtained from SP-Wind against the field measurements from Lillgrund is presented in Fig. 11. From Fig. 11b, it can be seen that the total farm power production in LES for three out of five cases are in good agreement with the measurement data, with the errors within 95 % confidence intervals of each other. While significant errors can be seen for the cases PDk3 and CNk41, the relative power error for the CNBL case is about 4.2 %, which is still quite low. The distribution of individual turbine power production for the case PDk1 is shown in Fig. 12, including uncertainty estimates again determined by the moving bootstrap method. From Figs. 11a and 12 it can be seen that individual turbine power trends across the farm are captured well, where similar trends are observed for power production peaks and valleys for un-waked and waked turbines, indicating that on average the wind direction in the LES cases captures the real world field conditions during the time windows. However, when comparing individual power production on a turbine level, statistically significant errors are evident for 58 % of the turbines for the PDk1 where error bars for the 95 % confidence intervals do not overlap when comparing LES predictions and field measurements. The same analysis is performed for all the simulated cases, and the proportion of turbines with statistically significant errors are presented in Table 5. In general, it was observed that the turbines at the front of the farm did not have statistically significant errors; however, turbines operating within the wakes of an upstream turbine had higher errors which were significant for wind directions resulting in an aligned configuration.

Figure 11Comparison of LES time-averaged power output of individual turbines (a) and total wind farm power output (b) against field SCADA measurements. Error bars represent 95 percent confidence intervals on the mean and are computed using the block–bootstrap method and time windows of 600 s and 1000 bootstrap samples.


Figure 12Bar plot of turbine power production for the case PDk1. Error bars represent 95 % confidence intervals on the mean and are computed using the block–bootstrap method and time windows of 600 s and 1000 bootstrap samples.


Figure 13Effect of wind direction on relative power error between LES calculations and field SCADA measurements for the five simulated cases.


The effect of the wind direction on the individual turbine power errors is presented in Fig. 13. It can be observed that for all the cases, the relative power error is generally low for all the turbines at the front of the farm facing the incoming flow field, further indicating that the replicated flow fields match the inflow conditions at the Lillgrund site. However, turbines within the farm and operating in strong wakes exhibit significantly higher errors, which gradually appear to reduce towards to back of the farm. This behavior is not visible in the case PDk2 and CNk4, with most of the turbines operating in an un-waked state due to the wind direction. This effect can also be witnessed when observing the proportion of turbines with significant errors, as presented in Table 5. The highest errors can be seen for the PDk3 case for almost all the turbines operating inside the farm, with some of the power errors exceeding 100 %. It was observed that while in SP-Wind all the turbines for this case were operational and producing power, most of the turbines in the field were either shut or producing negligible power. To further investigate the discrepancies in individual turbine power output across the five cases, we investigate the performance of the implemented controller in the wind farm. Figure 14 shows a comparison between the rotor speed and pitch actuation in SP-Wind against field measurements for all the turbines. It can be seen that while the performance of the pitch actuation and rotor speed controller show a good comparison in general across the five cases, larger errors can be seen for a few turbines in terms of both pitch and rotational speed. In particular, pitch measurements from the wind farm exhibit non-zero pitch angles for the majority of the turbines, even though they are operating in below-rated conditions and should traditionally have zero pitch, as is the case for the turbines operating in the LES simulations. The highest pitch errors are exhibited by the case PDk3, indicating that instead of operating in the Region 2 control regime as expected, the turbines within the farm are operating in a start-up region with significant pitching action at low wind speeds. This demonstrates the differences between the actual field turbine controller and the one implemented in SP-Wind, as detailed information about the controller and all of its modes of operation was not available. Larger pitch angles would lead to reduced power production for the turbines, which can be observed for the PDk3 case with the maximum amount of turbines with non-zero pitch, leading to the larger errors in power production as observed in Figs. 11 and 13. Additional controller tuning was attempted by adjusting the gains of the standard torque proportional controller based on the power and rotor speed measurements from the SCADA data. While a reduction in overall farm power error was observed, it was below 1 %, which is not significant. Another source of error worth investigating is the yaw misalignment between the simulated cases and the field measurements, which is presented in Fig. 15. While all the turbines in the simulation domain of SP-Wind are statically aligned with the mean wind direction without the use of a local yaw angle controller, it can be observed that this was not the case for the field turbines, which dynamically change their orientation in accordance with local wind direction changes within the farm. In fact, the local flow angle for the turbines operating in the LES cases were found to be negligible, indicating that even if a yaw angle controller was incorporated in SP-Wind, its effect would not have been significant. An analysis of turbine orientation error across the wind farm layout is shown in Fig. 16 for cases PDk1 and PDk3, as they exhibited the highest relative error of the five cases. It can be observed that for both the cases, which also have a similar mean wind direction, there is a gradient in the turbine orientation error which increases towards the bottom of the farm. This indicates a gradient in the wind direction across the farm potentially due to meso-scale effects, which the turbine yaw controllers are reacting to. As this gradient is not represented in the LES, we end up with higher turbine orientation errors for a sub-set of the wind turbines.

Figure 14Errors between LES time-averaged collective pitch (a) and rotational speed (b) of individual turbines and field SCADA measurements. Missing data points correspond to missing data.


Figure 15Errors between LES time-averaged wind direction at individual turbines and field SCADA measurements. All turbines in the LES domain face the wind and have zero yaw misalignment. Missing data points correspond to missing data.


Figure 16Effect of wind direction on relative turbine orientation error between LES calculations and field SCADA measurements for the cases PDk1 and PDk3. Missing data points correspond to missing SCADA data.


From the previous analysis, it can be concluded that multiple factors contribute to the errors in performance. It was observed that turbines operating in a strongly waked state exhibit the highest power errors across all the cases. This indicates that the wakes in the LES environment do not recover at the same rate as they do in the field and is further confirmed by a wake analysis in the following section. As previously remarked in Sect. 4.3, this could be due to the fact that the optimization framework used to match the inflow conditions at the Lillgrund site did not account for atmospheric stability and meso-scale effects. The larger errors in the mean velocity profile above the rotor area could also be affecting the wake recovery due to vertical transport, and the grid resolution may also be a contributing factor. An analysis of the turbine orientation error across the farm indicated gradients in wind direction in the field which were not recreated in the simulation environment, which could lead to differences in turbine–wake overlaps and hence power production errors across the farm. Finally, differences in controller operation in different regimes due to a lack of detailed information about the field controllers, in particular excessive pitching at low wind speeds, can have a significant effect on the observed errors.

Comparison of the mean flap-wise blade root bending moments for 5 of the 48 turbines from the Lillgrund wind farm for which loading data were available is shown in Figs. 17 and 18 for the PDBL and CNBL simulations, respectively. To determine the effect of fatigue, we use the damage equivalent loads (DELs) to compare the load histories of the same turbines across the LES and field measurement data. DEL is computed using the Palmgren–Miner rule and the Wöhler equation to account for accumulating fatigue damage caused to the wind turbine components by the fluctuating structural loads (Sutherland1999). The loads time series are counted and binned into individual cycles using the rainflow counting algorithm (Socie and Downing1982), and for the wind turbine blades the components follow the Wöhler’s curve with a slope coefficient equal to 10 (Freebury and Musial2000). The moving block bootstrap methodology is again utilized for evaluating the mean flap-wise moment and the corresponding DEL, with block lengths of 10 min and 1000 bootstrap iterations. Fatigue analysis was not conducted for the CNBL simulations as the data were not logged consistently for the turbines in the time periods of these simulations, making the data unfit for rainflow analysis. Both the average blade root flap-wise moments and the corresponding DELs exhibit a good comparison in the trends of the PDBL simulations in Fig. 17. For the PDk2 case, while the moments of the un-waked upstream turbines B08, C08, and D08 show lower errors, significant errors are observed for the turbine B06 which is operating in a fully waked state behind turbine C08. For turbine B07, which is operating in a partially waked state within the farm, while the differences in the flapwise moments is not significant, large errors are visible when comparing the flapwise DEL. A similar observation can be made for the CNBL simulations in Fig. 18. Turbines B08, C08, and D08 are operating in free-stream conditions for both the CNBL simulations, hence exhibiting better comparison against field measurements when compared to turbines B06 and B07 which are operating in fully waked conditions. A possible explanation for the larger errors observed in the loads of the waked turbines could be the relatively coarse grid resolution utilized in SP-Wind when looking at the number of cells across the rotor diameter, such that not all relevant turbulent structures upstream of turbine wakes that contribute to loads are captured. In the current work, the use of a finer grid for the simulations was not feasible, as the computational cost at the current grid resolutions for the wind farm simulations is already quite high. Additionally, the larger power error for the waked turbines as observed in Fig. 13 would result in differences in blade loading.

Figure 17Comparison of LES time-averaged blade root flap-wise moments (a) and DEL (b) against field measurements. Error bars represent 95 % confidence intervals on the mean. Results are normalized by maximum SCADA data for each case.


Figure 18Comparison of LES time-averaged blade root flap-wise moments for (a) CNk41 (b) CNk81 against field measurements. Error bars represent 95 % confidence intervals on the mean. Results are normalized by maximum SCADA data for each case.


5.4 Wake analysis

Due to equipment failure, data from the wake measuring lidars were unfortunately available for only two cases, PDk3 and CNk41, from the five selected validation time periods from the measurement campaign. Comparing wake recovery for the B06 turbine in the CNk41 case in Fig. 19, we see good agreement in wake location and recovery downstream from the turbines. It can be observed that while SP-Wind provides an accurate representation of the near-wake region behind turbine B06, the far-wake region characterized by a downstream distance greater than four rotor diameters exhibits higher errors. This can further be observed from the stronger wakes in SP-Wind in the far-wake region for turbine C06 in Fig. 20, leading to a lower inflow velocity at turbine C05 and hence the underestimation of power production for turbine C05 as observed in Fig. 11. Observing the wakes from the B05 turbine from PDk3 in Fig. 21, the effect of incorrectly representing the turbine orientation can be seen. While partial waked conditions are observed in SP-Wind for turbines A05, B05, and C05, the lidar measurements exhibit a higher wake overlap. This is observable through the wake deficit analysis shown in Fig. 21, where while at first for the LES case the wake behind the turbine B05 appears to be aligned with field measurements, larger deviations are visible further downstream from the turbine where the wake shifts to the left when observed from an upstream point of view. This leads to a larger inflow velocity across the downstream rotor at turbine C05 in LES as compared to lidar measurements and hence the greater reported power production in simulations as compared to SCADA data.

Figure 19Comparison of LES time-averaged velocity (b) against lidar wake measurements (a) for case CNk41 and turbine B06. (c) Comparison of velocity wake deficit at downstream locations from turbine B06. Line at y/D=0 is perpendicular to the turbine orientation from SCADA data. Shaded area represents 95 % confidence interval on the mean.


Figure 20Comparison of far-wake velocity wake deficit at downstream locations from turbine C06 for the case CNk41. Line at y/D=0 is perpendicular to the turbine orientation from SCADA data. Shaded area represents 95 % confidence interval on the mean.


Figure 21Comparison of LES time-averaged velocity (b) against lidar wake measurements (a) for case PDk3 and turbine B05. (c) Comparison of velocity wake deficit at downstream locations from turbine B05. Lidar data at x/D=4.0 were not available. Line at y/D=0 is perpendicular to the turbine orientation from SCADA data. Shaded area represents 95 % confidence interval on the mean.


6 Discussion, conclusions and future work

In this work, a validation study was conducted to compare SP-Wind, a high-fidelity large eddy simulation solver, against field measurements obtained from the Lillgrund offshore wind farm near the coast of Sweden. To recreate the atmospheric conditions at the Lillgrund site, a framework was developed to create the inflow conditions for the wind farm in the numerical domain by reusing a precursor database through scaling and shifting of the velocity. Thus, the cost intensive step of developing multiple precursor simulations for different atmospheric conditions spanning the duration of the measurement campaign was eliminated. Five time periods from the measurement campaign were selected for simulations in the LES environment. Upon comparison against field measurements, results from SP-Wind show good comparison in terms of the recreated inflow velocity field, power produced by individual turbines as well as the entire wind farm, mean and fatigue blade-root loading, and individual turbine wakes and recovery. However, limitations of the flow solver were exhibited in certain instances. Higher errors were observed in the performance of turbines operating in a strongly waked state. These errors were evident in both the turbine power and flapwise load measurements.The discrepancy in the performance of waked turbines could be attributed to the differences in atmospheric stability between the used precursor boundary layers and the atmospheric boundary layers at the Lillgrund site, or larger errors in the simulated profiles above the rotor area which could affect vertical transport and wake recovery behind turbines. Another source of error could be that the grid resolution in the numerical domain may not be fine enough to capture enough of the relevant smaller turbulent structures in turbine wakes. In the present study, the grid resolution was limited by the precursor database and computational costs. Additionally, controller mismatch due to lack of information of the field controllers also lead to discrepancies in the produced power. Particularly, higher errors were evident in a case which was operating closer to the cut-in wind speed of the turbines located at the Lillgrund site. Nevertheless, the results from the validation study are promising, proving the capability of a high-fidelity numerical solver to represent on-field conditions and performance output of a large wind farm.

The analysis in this work thus highlights certain areas in wind farm LES which still require further research to make accurate performance predictions for offshore wind farms. In particular, the lack of prior knowledge of meso-scale effects, field turbine control logic, and insufficiently fine grid resolution are highlighted as significant factors which can lead to errors in predictions. Therefore, future work should focus on improving on these effects, in particular for situations with aligned wind turbines. Further investigation into the role of atmospheric stability and meso-scale effects at the Lillgrund site should also be investigated to determine its role in wake recovery. This can be achieved by expanding the precursor database with finer resolution data sets, and data sets covering more atmospheric conditions. A closer co-operation with the turbine manufacturers to obtain a detailed controller description for the low wind speed zones and local yaw control dynamics would also be beneficial to improve modeling performance. Having exhibited the capability of the numerical solver in representing normal wind farm operation, validation studies could also be conducted to evaluate the effect of coordinated wind farm control strategies, such as wake steering and induction control, to improve wind farm performance.

Appendix A: Multibody model

A1 Coordinate system

The rotor of the wind turbine is considered as one single body with each blade modeled using a number of interconnected beam elements, as can be seen in Fig. A1. The origin of the body reference X1iX2iX3i is located at the center of the turbine rotor and coincides with the root nodes of the finite element beam representations of the respective blades. The X1i axis is directed along the length of the first blade, while the axis X3i is directed along the axis of rotation of the rotor in the upwind sense. The location of the origin of the body reference with respect to the global coordinate system X1X2X3 (not shown in Fig. 1) is denoted by the vector of Cartesian coordinates Ri, while the orientation of the body reference w.r.t. the global coordinate system is denoted by the vector of rotational coordinates θi. Only rotation along the turbine's main axis is assumed to contribute dynamically to the behavior of the turbine, while tilting and yawing motions are taken into account quasi-statically. Finally, the elastic deformations of the rotor are described by the vector of elastic coordinates qfi, which describe the displacements of the finite element nodes and the local derivatives thereof. Summarizing, the configuration of the turbine rotor can be described by the following vector of generalized coordinates:

(A1) q i = R i T θ i T q f i T T .

Figure A1Finite element beam representation of the wind turbine blades.


The global position of an arbitrary point on the jth beam element of the ith body of the multibody system can be written as

(A2) r i j = R i + A i u i j ,

where uij is the displacement vector of the ijth element and Ai is the transformation matrix of body i, which defines the orientation of the body reference with respect to the global reference. The yaw (ϕ), tilt (ψ), rotation (θ), and precone angles (γ) are used as rotational coordinates about their respective axes, and the resulting transformation matrix cased by these rotations is given by

(A3) A i = R y - 1 R t - 1 R r - 1 R p - 1 ,

where Ry,Rt,Rr, and Rp are the 3-D yaw, tilt, roll, and precone rotation matrices.

A2 Energy of the turbine rotor

The global velocity of a selected point can be determined by differentiating rij with respect to time to obtain

(A4) r ˙ i j = R ˙ i + A ˙ i u i j + A i u ˙ i j ,

which can be simplified by taking into account the assumption of quasi-static variation of yaw, tilt, and precone (ϕ˙,ψ˙,γ˙=0). Thus, the kinetic energy Tij of element ij in the finite element representation of the rotor can be obtained by the following formula:

(A5) T i j = 1 2 V i j ρ i j r i j ˙ T r ˙ i j d V i j ,

where Vij and ρij are the mass density and volume of the ijth element, respectively. The total kinetic energy of the turbine rotor can be determined by summing up the kinetic energies of all the elements. Using the floating frame of reference approach (Shabana2013), the expression for the elastic potential energy takes a very simple form, since only the elastic (i.e., not the rigid body) displacement of the body contributes to the elastic potential energy. Consequently, the potential energy of the ijth element is given by

(A6) Π i j = 1 2 q f i T K ff i j q f i ,

where Kffij is the element stiffness matrix expressed with respect to the body reference frame.

A3 Equations of motion

The equations of motion of the turbine rotor are developed from Lagrange’s equations for constrained systems. These equations can be written in the following form:

(A7) d d t T i q ˙ i T - T i q i T + Π i q i T + C q i λ T = Q i ,

where Ti and Πi are the kinetic and potential energy of the ith body, qi is the vector of generalized coordinates of the ith body, and Qi is the vector of generalized forces associated with the coordinates of the ith body. Furthermore, λ is the vector of Lagrange multipliers and is the Cqi constraint Jacobian matrix, defined as

(A8) C q i = C q i = C 1 q 1 i C 1 q 2 i C 1 q n i C 2 q 1 i C 2 q 2 i C 2 q n i C n c q 1 i C n c q 2 i C n c q n i ,

where C=C(q,t)=(C1C2Cnc)T is the vector of linearly independent constraint functions that satisfy the holonomic constraint equations of the multibody system

(A9) C ( q , t ) = 0 .

After evaluating and expanding the partial derivatives in equation A7 in terms of the mass and stiffness elements of the rotor structure, we can obtain the equation of motion of the multibody structure as

(A10) M i q i ¨ + K i q i + C q T λ = Q e i + Q v i ,

where Qei is the vector of generalized external forces containing the aerodynamic and gravity forces at each body element and Qvi is the quadratic velocity vector of body i, as defined by

(A11) Q v i = 1 2 q i q ˙ i T M i q ˙ i T - M ˙ i q ˙ i .
Appendix B: Coupling of turbine model with flow solver

During each LES time step, the blades sweep a sector area where the loads and their dynamic response are evaluated in a two-way fluid–structure interaction (FSI) manner. Loads acting on the turbine and tower structure lead to deformations, which are evaluated using Eq. (3), and the subsequent loads are then computed on the structure's deformed positions, before being added to the flow Eq. (1) as body force terms F. Before being added to the flow equations, the body forces of are processed by spatial and time filtering, as detailed in the following subsections.

B1 Spatial filter of rotor-swept forces

First, the unsteady forces F(u^,q) are smeared out in the surrounding LES mesh nodes by taking their convolution with a Gaussian kernel Gn, resulting in the spatially filtered forces F^(x):

(B1) F ^ ( x ) = n = 1 N t j = 1 N b = 3 0 R F ( u ^ , q , r ) G n ( x - r e j ) d r ,

where Nt is the number of turbines, Nb the number of blades, and ||x-rej|| is the Euclidean distance between the LES grid point and the deflected actuator line point, accounting for structural deformations.

B2 Time filter of rotor-swept forces

The body forces F of the Navier–Stokes equations are then calculated by time-filtering the Gaussian-filtered forces F^ through a first-order low-pass filter, which gives more weight to the last sub-iterations (see Fig. 4d). The time-filter is given as follows:

(B2) d F i d t = 1 τ f F ^ - F i ,

where Fi is the time-filtered force at Runge–Kutta stage i, F^ are the spatially filtered forces through the sub-cycles, and τf is the time filter constant. The filter constant τf defines the effective sector angle, which is chosen to be equal to the LES time step. Equation (B2) is integrated during the sub-iterations of the multibody solver using an implicit Euler scheme. The choice of using a first-order filter is justified by its simplicity; however, future work could focus on studying the influence of higher order filters.

Code and data availability

The SP-Wind flow solver is a proprietary software of KU Leuven. The TotalControl inflow database used in this work to recreate the inflow conditions at the Lillgrund wind farm is publicly available on Zenodo (, Munters et al.2019a;, Munters et al.2019b;, Munters et al.2019c;, Munters et al.2019d). Time-averaged LES inflow, SCADA and inflow lidar measurements are available online at (Sood2022). The rest of the raw data of the simulation results can be obtained by contacting the corresponding author.

Author contributions

IS and JM jointly set up the concept and objectives of the current work. IS performed code implementations in the LES code SP-Wind and carried out the wind farm simulations. AV and BB developed and implemented the multibody framework and turbine model. ES conducted the measurement campaign and data analysis under the guidance of GCL. IS selected the validation cases with support from JM. IS performed the data analysis and data visualization. IS and JM wrote the paper with contributions and revisions from ES, AV, and GCL.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Wind Energy Science. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government department EWI. The authors also acknowledge DTU, Siemens Gamesa, and Vattenfall for providing the lidar, loads, and SCADA data, respectively.

Financial support

This research has been supported by the European Commission, Horizon 2020 Framework Programme (TotalControl (grant no. 727680)).

Review statement

This paper was edited by Julie Lundquist and reviewed by two anonymous referees.


Allaerts, D. and Meyers, J.: Large eddy simulation of a large wind-turbine array in a conventionally neutral atmospheric boundary layer, Phys. Fluids, 27, 065108,, 2015. a, b, c, d

Anderson, S. J., Meyers, J., Sood, I., and Troldborg, N.: TotalControl D 1.04 Flow Database for reference wind farms, (last access: 9 December 2021), 2020. a, b

Archer, C. L., Mirzaeisefat, S., and Lee, S.: Quantifying the sensitivity of wind farm performance to array layout options using large-eddy simulation, Geophys. Res. Lett., 40, 4963–4970,, 2013. a

Arnold, M., Brüls, O., Arnold, M., Brüls, O., Arnold, M., and Brüls, O.: Convergence of the generalized-α scheme for constrained mechanical systems, Multibody System Dynamics, Springer Verlag, 85, 187–202,, 2007. a

Arya, S.: Comparative Effects of Stability, Baroclinity and the Scale-Height Ratio on Drag Laws for the Atmospheric Boundary Layer, J. Atmos. Sci., 35, 40–46,<0040:CEOSBA>2.0.CO;2, 1978. a

Arya, S. P. S.: Comments on “Similarity Theory for the Planetary Boundary Layer of Time-Dependent Height”, J. Atmos. Sci., 32, 839–840,<0839:COTFTP>2.0.CO;2, 1975. a

Bastankhah, M. and Porté-Agel, F.: A new miniaturewind turbine for wind tunnel experiments. Part I: Design and performance, Energies, 10, 923,, 2017. a

Boersma, S., Doekemeijer, B. M., Siniscalchi-Minna, S., and van Wingerden, J. W.: A constrained wind farm controller providing secondary frequency regulation: An LES study, Renew. Energ., 134, 639–652,, 2019. a

Bossanyi, E.: Combining induction control and wake steering for wind farm energy and fatigue loads optimisation, J. Phys.-Conf. Ser., 1037, 032011,, 2018. a

Bou-Zeid, E., Meneveau, C., and Parlange, M.: A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows, Phys. Fluids, 17, 1–18,, 2005. a

Brost, R., Lenschow, D., and Wyngaard, J.: Marine Stratocumulus Layers. Part 1: Mean Conditions., J. Atmos. Sci., 39, 800–817,<0800:MSLPMC>2.0.CO;2, 1982. a

Calaf, M., Meneveau, C., and Meyers, J.: Large eddy simulation study of fully developed wind-turbine array boundary layers, Phys. Fluids, 22, 015110,, 2010. a, b, c

Castro, I. P.: Rough-wall boundary layers: Mean flow universality, J. Fluid Mech., 585, 469–485,, 2007. a

Churchfield, M. J., Schreck, S., Martínez-Tossas, L. A., Meneveau, C., and Spalart, P. R.: An advanced actuator line method for wind energy applications and beyond, 35th Wind Energy Symposium, Grapevine, Texas, 9–13 January 2017,, 2017. a

Csanady, G. T.: Equilibrium theory of the planetary boundary layer with an inversion lid, Bound.-Lay. Meteorol., 6, 63–79,, 1974. a

Draper, M., Guggeri, A., and Usera, G.: Validation of the Actuator Line Model with coarse resolution in atmospheric sheared and turbulent inflow, J. Phys.-Conf. Ser., 753, 082007,, 2016. a

Frederik, J. A., Doekemeijer, B. M., Mulders, S. P., and van Wingerden, J. W.: The helix approach: Using dynamic individual pitch control to enhance wake mixing in wind farms, Wind Energy, 23, 1739–1751,, 2020. a

Freebury, G. and Musial, W.: Determining equivalent damage loading for full-scale wind turbine blade fatigue tests, 2000 ASME Wind Energy Symposium, Reno, NV, USA, 10–13 January 2000, 287–297,, 2000. a

Ghaisas, N. S. and Archer, C. L.: Geometry-Based Models for Studying the Effects of Wind Farm Layout, J. Atmos. Ocean. Tech., 33, 481–501,, 2016. a

Göçmen, T. and Giebel, G.: Estimation of turbulence intensity using rotor effective wind speed in Lillgrund and Horns Rev-I offshore wind farms, Renew. Energ., 99, 524–532,, 2016. a

Goit, J. P. and Meyers, J.: Optimal control of energy extraction in wind-farm boundary layers, J. Fluid Mech., 768, 5–50,, 2015. a

Hansen, M. H. and Henriksen, L. C.: Basic DTU Wind Energy controller, DTU Wind Energy, (last access: 1 December 2021), 2013. a

Jiménez, J.: Turbulent flows over rough walls, Annu. Rev. Fluid Mech., 36, 173–196,, 2004. a

Kunsch, H. R.: The Jackknife and the Bootstrap for General Stationary Observations, Ann. Stat., 17, 1217–1241,, 1989. a

Liang, X.: An Integrating Velocity–Azimuth Process Single-Doppler Radar Wind Retrieval Method, J. Atmos. Ocean. Tech., 24, 658–665,, 2007. a

Lin, M. and Porté-Agel, F.: Large-eddy simulation of yawedwind-turbine wakes: comparisons withwind tunnel measurements and analyticalwake models, Energies, 12, 1–18,, 2019. a

Lu, H. and Porté-Agel, F.: Large-eddy simulation of a very large wind farm in a stable atmospheric boundary layer, Phys. Fluids, 23, 065101,, 2011. a

Magnusson, M. and Smedman, A.-S.: Influence of Atmospheric Stability on Wind Turbine Wakes, Wind Eng., 18, 139–152, (last access: 1 December 2021), 1994. a

Martínez-Tossas, L. A., Churchfield, M. J., and Leonardi, S.: Large eddy simulations of the flow past wind turbines: actuator line and disk modeling, Wind Energy, 18, 1047–1060,, 2015. a

Mehta, D., van Zuijlen, A. H., Koren, B., Holierhoek, J. G., and Bijl, H.: Large Eddy Simulation of wind farm aerodynamics: A review, J. Wind Eng. Ind. Aerod., 133, 1–17,, 2014. a

Munters, W. and Meyers, J.: Dynamic strategies for yaw and induction control of wind farms based on large-eddy simulation and optimization, Energies, 11, 177,, 2018. a, b

Munters, W., Meneveau, C., and Meyers, J.: Turbulent Inflow Precursor Method with Time-Varying Direction for Large-Eddy Simulations and Applications to Wind Farms, Bound.-Lay. Meteorol., 159, 305–328,, 2016. a, b, c

Munters, W., Sood, I., and Meyers, J.: Precursor dataset PDk, Zenodo [data set],, 2019a. a, b, c

Munters, W., Sood, I., and Meyers, J.: Precursor dataset PDkhi, Zenodo [data set],, 2019b. a, b, c

Munters, W., Sood, I., and Meyers, J.: Precursor dataset CNk2, Zenodo [data set],, 2019c. a, b, c

Munters, W., Sood, I., and Meyers, J.: Precursor dataset CNk4, Zenodo [data set],, 2019d. a, b, c

Nilsson, K., Ivanell, S., Hansen, K. S., Mikkelsen, R., Sørensen, J. N., Breton, S.-P., and Henningson, D.: Large-eddy simulations of the Lillgrund wind farm, Wind Energy, 18, 449–467,, 2014. a

Porté-agel, F., Bastankhah, M., and Shamsoddin, S.: Wind-Turbine and Wind-Farm Flows: A Review, Bound.-Lay. Meteorol., 44, 1573–1472,, 2020. a

Salarpour, A. and Khotanlou, H.: An empirical comparison of distance measures for multivariate time series clustering, International Journal of Engineering, Transactions B: Applications, 31, 250–262, 2018. a

Sebastiani, A., Castellani, F., Crasto, G., and Segalini, A.: Data analysis and simulation of the Lillgrund wind farm, Wind Energy, 24, 634–648,, 2021. a

Shabana, A. A.: Dynamics of Multibody Systems, Cambridge University Press, 4 edn.,, 2013. a, b, c, d

Simisiroglou, N., Polatidis, H., and Ivanell, S.: Wind farm power production assessment: a comparative analysis of two actuator disc methods and two analytical wake models, Wind Energ. Sci. Discuss. [preprint],, 2018. a, b

Simon, E. and Courtney, M.: A Comparison of sector-scan and dual Doppler wind measurements at Høvsøre Test Station – one lidar or two?, Report, (last access: 30 November 2021), 2016. a

Socie, D. and Downing, S.: Simple Rainflow Counting Algorithms, Int. J. Fatigue, 4, 31–40,, 1982. a

Sood, I.: Comparison of Large Eddy Simulations against measurements from the Lillgrund offshore wind farm – Manuscript data, Zenodo [data set],, 2022. a

Sood, I., Meyers, J., and Lanzilao, L.: TotalControl D 1.8 Coupling of Gaussian wake merging to background ABL model, (last access: 1 December 2021), 2020a. a

Sood, I., Munters, W., and Meyers, J.: Effect of conventionally neutral boundary layer height on turbine performance and wake mixing in offshore windfarms, J. Phys.-Conf. Ser., 1618, 062049,, 2020b. a, b

Stevens, R. J., Graham, J., and Meneveau, C.: A concurrent precursor inflow method for Large Eddy Simulations and applications to finite length wind farms, Renew. Energ., 68, 46–50,, 2014. a, b

Storey, R. C., Norris, S. E., and Cater, J. E.: An actuator sector method for efficient transient wind turbine simulation, Wind Energy, 18, 699–711,, 2015. a

Sutherland: Fatigue analysis of wind turbines. Technical report, Sandia National Laboratories, Tech. rep.,, 1999. a

Townsend: The Structure of Turbulent Shear Flow. Cambridge University Press, ISBN 9780521298193, 1976. a

Vasiljević, N., Lea, G., Courtney, M., Cariou, J.-P., Mann, J., and Mikkelsen, T.: Long-Range WindScanner System, Remote Sensing, 8, 896,, 2016. a

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods, 17, 261–272,, 2020. a

Vitsas, A. and Meyers, J.: Multiscale aeroelastic simulations of large wind farms in the atmospheric boundary layer, J. Phys.-Conf. Ser., 753, 082020,, 2016. a, b, c

Wu, Y. T. and Porté-Agel, F.: Large-Eddy Simulation of Wind-Turbine Wakes: Evaluation of Turbine Parametrisations, Bound.-Lay. Meteorol., 138, 345–366,, 2011. a

Wu, Y. T. and Porté-Agel, F.: Simulation of Turbulent Flow Inside and Above Wind Farms: Model Validation and Layout Effects, Bound.-Lay. Meteorol., 146, 181–205,, 2013.  a

Wu, Y. T. and Porté-Agel, F.: Modeling turbine wakes and power losses within a wind farm using LES: An application to the Horns Rev offshore wind farm, Renew. Energ., 75, 945–955,, 2015. a

Yılmaz, A. E. and Meyers, J.: Optimal dynamic induction control of a pair of inline wind turbines, Phys. Fluids, 30, 085106,, 2018. a

Short summary
In this work, we conduct a validation study to compare a numerical solver against measurements obtained from the offshore Lillgrund wind farm. By reusing a previously developed inflow turbulent dataset, the atmospheric conditions at the wind farm were recreated, and the general performance trends of the turbines were captured well. The work increases the reliability of numerical wind farm solvers while highlighting the challenges of accurately representing large wind farms using such solvers.
Final-revised paper