Efficient Mann turbulence generation for offshore wind farms with applications in fatigue load surrogate modelling

Surrogate models are vital for offshore wind farm optimisation, or digital twin for rapid execution. Reliable surrogates are difficult to make, especially for large wind farms and for fatigue loads. One limitation is the high resolution of turbulence needed for fatigue calculations and the vast size of a wind farm. We present code optimisations for implementing the Mann model for farm-sized flows, which are infeasible otherwise. Additionally, we present a case study of a full load surrogate for the Lillgrund offshore wind farm using large turbulence box generation in aeroelastic wind farm simulations using HAWC2Farm. Various mappings between the wind field in front of a turbine and the corresponding turbine structural fatigue loads are presented. The best-performing surrogate model, using a proper orthogonal decomposition of the input space with an artificial neural network, is able to significantly reduce the error in fatigue load estimates by as much as a factor of 5 compared to conventional methods. The presented recommendations are crucial to generating reliable wind farm load surrogates for wind farm optimisation purposes.


Introduction
The scale and rate of wind farm installations continue to grow globally, and with it, the need for efficient wind farm designs.To achieve this, optimisations are commonly performed to decide on the wind turbine model, wind farm layout, wind farm control strategy, and many more aspects to maximise the wind energy power output while minimising costs.At the heart of these optimisations are low-fidelity wind farm simulators which can rapidly compute the response of a wind farm design for a wide range of wind conditions.Low-fidelity wind farm simulators are typically steady-state, and therefore unable to capture the dynamic nature of the wind farm.This becomes an issue when fatigue load estimates, which require dynamic time series data, are desired, either to constrain the optimisation or for multi-objective wind farm control [1].It is possible to calculate wind farm fatigue loads using aeroelastic wind farm simulators such as HAWC2Farm [2] or FAST.Farm [3], which can simulate the dynamic structural response of turbines in a wind farm.Surrogate models can then be formulated from these mid-fidelity simulations to quickly estimate difficult-to-calculate quantities within an optimisation loop.The accuracy of a surrogate model relies heavily on that of the underlying data set.Such fatigue load data sets require numerous aeroelastic wind farm simulations using extremely large and high-resolution turbulence boxes, which are difficult and slow to produce using current methods.In this study, we present a methodology for performing aeroelastic wind farm simulations using large, high-resolution turbulence boxes in conjunction with HAWC2Farm, followed by surrogate modelling techniques for predicting turbine fatigue loads.
The study is structured as follows.First, we describe the wind farm simulation setup of the Lillgrund wind farm.Second, we describe the turbulence generation methodology.1000 separate hour-long HAWC2Farm simulations at 100 Hz and a high spatial resolution are performed to generate a comprehensive fatigue load data set.The Mann turbulence model is presented in a way which minimises computational redundancies and exploits parallelisation, allowing for rapid generation of farm-sized turbulence boxes which otherwise would be infeasible to generate.Third, the wind field immediately in front of each turbine is parameterised by including states for the mean wind speed, vertical wind shear, and wakes.The wind states are identified via numerical optimisation and were found to better describe fatigue loads than other parameterisations found in literature.Fourth, a series of surrogate models are trained between the wind states and the fatigue loads.Among the plethora of surrogate models, we have considered quadratic least-squares (QLS) for their simplicity and artificial neural networks (ANN) for their accuracy.Finally, the results and conclusions are presented.The combination of aeroelastic simulation design and surrogate model design demonstrated in this study is important in the field of wind farm optimisation to ensure optimised solutions are grounded in accurate and unbiased wind farm simulation data.

Aeroelastic wind farm simulation database with HAWC2Farm
We present a database of fatigue load results from aeroelastic wind farm simulations of the Lillgrund wind farm.The Lillgrund wind farm consists of 48 wind turbines, each with a nameplate capacity of 2.3MW.The wind farm, with a layout illustrated in Fig. 1, is simulated using HAWC2Farm [2].HAWC2Farm can model the dynamic and aeroelastic response of each wind turbine in a wind farm using the aeroelastic turbine simulator, HAWC2 [4] in conjunction with the dynamic wake meandering model [5].1000 HAWC2Farm simulations with a sampling frequency of 100Hz are performed with a quasirandom sequence of input parameters sampled using the Halton sequence [6].The Halton sequence is chosen to produce a more evenly distributed sample space compared to a pseudorandom sequence.The ambient wind speed, U ∞ , turbulence intensity, TI, wind direction, θ, and shear exponent, α are sampled within the bounds described in Table 1.  ).The dimensions of this 72 km long turbulence box will cause most turbulence generators to run out of memory or take several hours to compute.To overcome this, the four main steps in generating a Mann turbulence box (Fig. 2) were implemented with several innovations to be able to generate these immense boxes as listed below.Using these modifications, a single turbulence box of the aforementioned dimensions can be computed in approximately 5 minutes on the DTU Sophia HPC cluster [10].
1 Coefficient calculation The first step (step 1 in Fig. 2) is to generate the coefficients, the turbulence box of a given size.For small wave vectors (approximately |k| ≤ 3 L as suggested by Mann [7]), it is recommended to approximate an intermediate value, A k by solving: where Φ k is the Mann spectral tensor for wave number, k, and H indicates the Hermitian of the matrix.Eq. ( 1) must be solved numerically, requiring numerous evaluations of the Mann spectral tensor to accurately estimate the integral, which can be computationally costly.For large turbulence boxes, the active region of integration varies greatly depending on the wave vector, making a fixed grid integration inaccurate for this application.We, therefore, use adaptive quadrature to subdivide the integration as required [11].
C k in Fig. 2 can then be determined by decomposing A k from Eq. (1) using, for example, a Cholesky decomposition [12].For large wave numbers, an approximation of C k can be used, requiring only a single decomposition of the Mann spectral tensor, Φ 1/2 k which has an analytical solution.C k can therefore be computed as: where s is a constant term (See Mann (1998) for details [7]).A number of code optimisations are used in this study to reduce the time and memory requirement for this step in high-resolution boxes: • Reuse of spectral tensor stencil: The spectral tensor decompositions for fixed box discretisation, size, Γ and L, can be reused for different random seeds and turbulence intensities.Typically, the box dimensions are known in advance, making it possible to reuse the spectral tensor calculations to quickly generate new turbulence fields.The spectral tensor decompositions can be stored in a 5-dimensional array of size N x × N y × (N z /2 + 1) × 3 × 3 in row-major array storage, ensuring each tensor is stored contiguously for faster memory access.In cases where memory is scarce, the storage of the stencil can be bypassed, and instead can be calculated directly as needed in step 3 .
• Parallelisation: Each coefficient, C k , can be calculated independently of one another, therefore presenting the opportunity for parallel computing.This is particularly the case when using high performance computing clusters, which can provide 32 or more process threads.
2 and 3 Monte Carlo simulation The heart of the turbulence generation algorithm is the Monte Carlo simulation of the correlated turbulence frequency components.The uncorrelated random phases ( 2 ) in r k become correlated in accordance with the Mann model after the linear mapping The evaluation of this step is easily parallelised.
4 Frequency to space domain transformation The final step is to convert the turbulence boxes from the frequency domain to the space domain using an inverse 3D Fourier transform.A number of memory and computational code optimisations can be performed by splitting the task into three rounds of 1D inverse Fast Fourier Transforms (FFT).
• Parallelisation: The 3D inverse FFT algorithm can be parallelised by considering it as three stages of 1D FFTs, also known as a volumetric decomposition [13].In each stage, the 1D signals are independent, allowing the FFT to be performed on each signal in parallel, providing a noticeable speed-up for high-dimensional turbulence boxes.While other specialised parallel 3D FFTs exist which use pencil or slab decompositions [14], the volumetric decomposition was found to be general enough to not consider other optimisations for the presented problem.• Real inverse Fourier transform: The first two rounds of inverse Fourier transforms map from complex to complex domains, whereas the final transform maps to the real numbers.As a result, the final FFT should make use of a real inverse FFT to best exploit the available symmetries.While the algorithm is mathematically equivalent, a real inverse FFT reduces both the computational and memory requirements by half due to symmetry.This additionally speeds up step 4 as only half the coefficients must be computed and stored, hence the dimension of size (N z /2 + 1) (for even-length signals).• Low-prime box dimensions: Past implementations of Mann turbulence generators required dimensions to be powers of 2 [4].This limitation has been relaxed with modern Fourier transform implementations using improved planners [15].Using such methods, any signal length is possible, where dimensions with small prime factors run the fastest.It was found that using box dimensions with the largest prime factor being less than 13 ran in a comparable time to powers-of-2 dimensions.This provides more box discretisation options near the desired threshold, which can noticeably reduce the memory and computational requirements in the turbulence generation procedure.
An implementation of these recommendations is provided as an open-source turbulence generator code [16].

Wind field parameterisation
The wind field parameterisation must be carefully chosen when designing the surrogate model to ensure the fatigue loads can be adequately described.It is common to see wind speed and turbulence intensity as wind field parameters for surrogate models, such as in Riva et al. (2020) [17].While this parameterisation is often adequate for predicting power output, more spatial information is required to predict fatigue loads.For example, when a turbine is partially waked, the fatigue loads of some components can be either increased or decreased depending on the wake centre position.For this reason, other literature has considered both lateral and vertical shear as input parameters [18], as well as several spatial moments across the rotor [19].
In this study, we present two wind field parameterisation methods.The first uses Proper Orthogonal Decomposition (POD), and the second is a skew-normal representation.The first method uses POD on the wind field database to determine a reduced-order basis which can be used to decompose the wind field.POD is commonly used in fluid dynamics to decompose a flow field into orthogonal modes [20].We indicate the ith input variable with x i ∈ R nx .All inputs are collected in the array X = [x 1 , ...x np ] ∈ R nx×np , where index i denotes the point number and j the feature.By decomposing X as where, The first most significant vectors in U can be used as an approximate basis for decomposing wind field observations.The POD operation orders the vectors in order of variance explained, for which only 6 basis vectors are required set to recreate 99.99% of the variations in this data set (Fig. 4).The basis vectors in Fig. 3 capture the mean velocity, vertical shear, partial wakes, full wakes, as well as higher-order dynamics.We use 10 basis vectors as the input for the surrogate model as well as the standard deviation of the extracted wind field over the 10-minute observation is also provided as a proxy measure of turbulence intensity.Finding the reduced-order representation, x, of a single wind field requires a simple matrix multiplication: x = Ũ T x, where Ũ is the truncation of U .The second parameterisation method aims to capture the presence of a vertical shear profile and a wake, described using skew-normal functions [21].The use of a skew-normal distribution in the lateral direction allows for superimposed or stretched wakes to be modelled under the assumption that most compounded wake effects originate from the same row of turbines.With the vertical wake position and skew parameter fixed at zero to simplify the fit, the following expression is used to describe the wind field: where ϕ is a scaled normal distribution function and Φ is the scaled cumulative distribution function in order to generate a skew-normal shape in the y−direction: where erf(x) is the error function.The key parameters are the mean hub wind speed, U hub , shear exponent, α, and a single wake parameterised by a skew-normal distribution with scaling amplitude, A, lateral width, w y , vertical width, w z , lateral location, y 0 , skew parameter, β.Like the POD representation, we append the wind field standard deviation to represent turbulence intensity.The wind field parameterisation used in Eq. (4) is established using numerical minimisation.

Surrogate modelling
For this work, we have considered two types of surrogate models linking the wind states to the DEL: the quadratic least-squares and feed-forward fully-connected neural networks.The quadratic least-squares are extremely fast to train, and provide valuable physical insight into the dependence between input and output, but cannot model highly nonlinear functions.On the other hand, neural networks are complex architectures that are difficult to train, and require orders of magnitude more parameters, while providing a much lower prediction error.Since input and outputs can have different units, and order magnitudes, we have applied a min-max scaler.Each input feature, x, is scaled to an arbitrary range of [−1, +1].
Similarly, the output, y, is scaled to the range [y min , y max ], here set to [−1, +1].For ease of notation we will write x = ψ x ( x) and ȳ = ψ y (y).

Quadratic least-squares
The quadratic least-squares assume parabolic DEL with respect to input features, except for the wake centre y 0 .The DWM model predicts a DEL-y 0 dependence with two peaks, which is symmetric or asymmetric based on the load channel [22].To address this, y 2 0 and y 3 0 were added to the input database and removed y 0 to prevent rank deficiency.The quadratic least-square algorithm produces a global approximation of the high-fidelity model.For a single data point and output channel j, the prediction is given by where c j , b j and A j are the constant, linear and quadratic coefficients, for output channel j.
The training and prediction of this model have been done with the OpenTURNS library [23].

2.4.2.
Artificial Neural Network (ANN) ANN-based surrogates use feed-forward fullyconnected neural networks (FFNN) to capture highly non-linear trends in the structural response of a turbine [24,25].The FFNN has two hidden layers and is described by Eq. 8, where activation functions h 1 and h 2 are used.Grid search is performed over the number of neurons and activation functions to optimise the model.
The FFNN algorithm is implemented using TensorFlow [26] and trains five separate surrogates (one single output, four multi-output) to capture turbine response in power and load channels.The most optimal model architecture for each target channel is presented in Table 2.

Turbulence generation benchmark
First, we compare the computational requirements of the presented turbulence generator, Mann.rs, against other available turbulence generators.Among these are PyConTurb [12], which uses the Sandia method for generating turbulent wind fields [27], and the Mann turbulence generator provided with HAWC2 [4].Turbulence over a square domain with 200m sides is simulated at increasing resolutions on all generators on a single core of a node with 256 GB of RAM in the Sophia HPC [10].The computational time of PyConTurb increases rapidly with the number of grid points and runs out of memory at 128 3 grid points.This is due to the use of the Sandia method, where the number of multiplications scales proportional to N x N 3 y N 3 z , compared to the Mann model, which scales proportionally to N x N y N z .Mann.rs outperforms the HAWC2 turbulence generator by over 1000 times for low-resolution boxes, and approximately 15 times faster for high-resolution boxes.

Wind field parameterisation
We compare the POD parameterisation with 10 basis vectors and the skew-normal wake fitting with a reference case of using the rotor effective wind speed (REWS) to reconstruct the 30 × 30 wind fields.Given that the first basis vector of the POD decomposition is a uniform wind field, the REWS method can be interpreted as the POD parameterisation with one basis vector.The root-mean-square percentage error (RMSPE) is used to quantify the performance in Table 3, showing that the POD and skew-normal parameterisations outperform REWS at reconstructing the wind fields by over 11%, with POD performing the best with an RMSPE of 1.35%.for varying resolutions using PyConTurb [12], the HAWC2 turbulence generator [4], and the presented turbulence generator, Mann.rs [16].Table 3. RMSPE of the wind field reconstructions using various parameterisations.

Surrogate model performance
Surrogate models trained on 800 samples from the database in Section 2.1 were tested on a final 200 samples to evaluate their ability to predict fatigue loads for various turbine components given wind field parameters (Fig. 8).Both skew-normal and POD parameterisations of the wind field outperformed using only REWS and TI by 31% and 50%, respectively, for blade root edgewise DEL (MxBr) using QLS.The ANN further improved predictions by 123% compared to QLS, better approximating the relationship between wind field and fatigue loads.Similar trends were observed for other turbine components, except for tower bottom fore-aft DEL (MyTB), where REWS slightly outperformed skew-normal and POD parameterisations by 1-2%, but still lagged behind the ANN with POD parameterisation.Training and hyperparameter search for the ANN took several hours, while QLS was trained in milliseconds.

Conclusions
We present a procedure for generating load surrogate models for offshore wind farms using aeroelastic simulations of the Lillgrund wind farm.To enable high-resolution turbulence boxes for simulations, we propose a computationally efficient method for generating Mann turbulence.The aeroelastic simulations serve as input data for structural fatigue load surrogate models, with the wind field in front of a turbine as input.The best-performing surrogate model can distinguish spatial features in the wind field, such as the position of incoming wakes, by using proper orthogonal decomposition of the wind in conjunction with an artificial neural network.Both the generated HAWC2Farm data set and the turbulence generator are available online.The presented methods have applications in offshore wind farm layout and control optimisations.

Figure 2 .
Figure 2. Flow diagram of the Mann turbulence generation process, where the block sizes approximately represent the memory requirement of each step.The illustrated steps are: 1 coefficient generation, 2 random vector generation, 3 Monte Carlo simulation, and 4 transformation from frequency to spatial domain.

Figure 3 .
Figure 3. First eight (out of 10) basis vectors of the POD used to parameterise the wind field.

Figure 4 .
Figure 4. Cumulative explained variance versus number of basis vectors used to parameterise the wind field.

Figure 5 .
Figure 5. Examples of 10-minute averaged wind fields of partial wakes (A and B) and a full wake (C), where the sampled fields (upper) are reconstructed (lower) using a fitting of Eq. (4).

Figure 6 .
Figure 6.Flow diagram of the mapping process from rotor wind field to fatigue loads.

Table 2 .
Best model architecture for FFNN surrogates per output channel