Fatigue design load calculations of the offshore NREL 5MW benchmark turbine using quadrature rule techniques

A novel approach is proposed to reduce, compared to the conventional binning approach, the large number of aeroelastic code evaluations that are necessary to obtain equivalent loads acting on wind turbines. These loads describe the effect of long-term environmental variability on the fatigue loads of a horizontal-axis wind turbine. In particular Design Load Case 1.2, as standardized by IEC, is considered. The approach is based on numerical integration techniques and, more specifically, quadrature rules. The quadrature rule used in this work is a recently proposed"implicit"quadrature rule, which has the main advantage that it can be constructed directly using measurements of the environment. It is demonstrated that the proposed approach yields accurate estimations of the equivalent loads using a significantly reduced number of aeroelastic model evaluations (compared to binning). Moreover the error introduced by the seeds (introduced by averaging over random wind fields and sea states) is incorporated in the quadrature framework, yielding an even further reduction in the number of aeroelastic code evaluations. The reduction in computational time is demonstrated by assessing the fatigue loads on the NREL 5MW reference offshore wind turbine in conjunction with measurement data obtained at the North Sea, both for a simplified and a full load case.


Introduction
For the certification of the design of a wind turbine, various load cases have been formulated in the IEC 61400 standard for both onshore [19] and offshore [20] wind turbines. These load cases vary from regular power production in commonly observed environmental conditions, to extreme conditions simulating storm and failure. All cases have one property in common: the aeroelastic simulation code that models the wind turbine must be evaluated a large number of times, which is computationally costly. In this work arguably one of the most costly load cases is considered: Design Load Case (DLC) 1.2, describing fatigue loading for a power producing offshore wind turbine under regular environmental conditions. As the time span of this load case is the full lifetime of the turbine, a large number of model runs is necessary to assess the statistics of the fatigue loading.
The conventional approach to assess the effect of environmental variability on the turbine lifetime is to firstly split the domain of the variables describing the environmental conditions in bins, secondly run the aeroelastic model several times in each bin (the so-called seeds), and finally determine the quantity of interest (e.g. the weighted equivalent load) incorporating the probability of occurrence of each bin. This approach is suggested in the aforementioned standard [19,20] and has been successfully applied in previous research [1,14,31].
If applicable, binning is a versatile and effective tool to assess the effect of parameter variability. However, a major disadvantage of binning is that the number of necessary evaluations of the aeroelastic model is often prohibitively large, which limits its applicability. Several approaches have been suggested over the years to alleviate the high computational cost, e.g. lumping as discussed in Kühn [22,Section 8.3]. This often requires non-trivial preprocessing steps and restrictive assumptions on the aeroelastic model.
Other approaches with similar flexibility exist. Examples include for example quasi Monte Carlo methods [9,29], surrogate methods [2,28,37], or numerical integration techniques by means of quadrature approaches [4,7,12,16,23]. These approaches have in common that regularity in the model (i.e. continuous responses to variability in the input parameters) is leveraged such that fast convergence can be obtained for sufficiently smooth functions. They have been successfully used to assess uncertainties in models in wind energy [6,13,18,26,27,32], but a rigorous approach to model standardized load cases that directly competes with binning is still lacking. Furthermore, since the number of bins is fixed, it is difficult to assess the accuracy of the obtained estimation of the quantity of interest (the loads in the cases studied in this work). Moreover, explicit distributions of the parameters must be provided a priori. These are often unavailable, as there are various distributions that could potentially fit the measurements of the uncertain environment [25].
The goal of this work is to use quadrature rule approaches to assess fatigue loads described by DLC 1.2, leading to a competitive and attractive alternative approach to binning. Quadrature rules have the advantage that they are tailored to calculating integrals, which is the key goal of the load case. Moreover, they are versatile and robust by providing fast convergence for smooth functions, while still yielding accurate results for non-smooth functions. The quadrature rule used in this work is the implicit quadrature rule [5], which has the additional advantage over conventional quadrature rules [12,16,23] that no explicit knowledge about the distribution of the environmental parameters is required: only measurements of the environmental conditions are necessary. Furthermore, it provides an error estimate of the estimated load.
This article is structured as follows. In Section 2 fatigue load calculation is introduced, consisting of describing DLC 1.2, the on-site conditions considered in this article, and the wind turbine under consideration. The propagation of parametric uncertainty through an aeroelastic model is discussed in Section 3. Here, the benchmark binning procedure is briefly explained and the quadrature rule algorithm and its most relevant mathematical properties are discussed in detail. The results obtained by binning and quadrature rules are compared in Section 4, demonstrating the applicability and prospect of quadrature rules. The article is concluded in Section 5.

Fatigue load calculation
The DLCs as specified in the design requirements for offshore wind turbines by IEC [19,20] can be globally put into two categories: ultimate and fatigue load cases. The ultimate load cases describe the simulation of failure of a component due to rare events. Often these cases require nontrivial extrapolation procedures to assess the statistics. On the other hand, the fatigue load cases describe the simulation of regular environmental conditions over a long period of time to assess the effect of wear of the turbine. The last group is the focus of this article, as this type of cases requires a large number of evaluations of an aeroelastic code. In particular, the focus is on DLC 1.2 in combination with DLC 6.4. DLC 1.2 accounts for environmental conditions that result into power production and DLC 6.4 describes similar environmental conditions, with the sole exception that the turbine is idling due to too small or too large wind speeds 1 . There are various other fatigue load cases describing specific scenarios, e.g. start up (DLC 3.1), shutdown (DLC 4.1), and fault occurrence (DLC 2.4), but these are not the primary focus of this work since assessing these cases is significantly less computationally costly.
In this section, the specific details of the load case considered in this article are outlined. In Section 2.1 the details of the DLC are discussed, as standardized by IEC [19,20]. In Section 2.2 and 2.3 the environmental conditions and the wind turbine under consideration are discussed. To keep the simulation as realistic and as reproducible as possible, freely available measurement data [36] and the NREL 5MW reference turbine [21] are used to describe the environment at the offshore site and model a wind turbine that should perform well at the location that is considered. In this work, the aeroelastic code BLADED is employed for the load calculations, which is briefly discussed in Section 2.4.

Design Load Case 1.2: fatigue analysis
The goal of DLC 1.2 and DLC 6.4 is to assess the effect of wear on a wind turbine over a long period of time under normal design situations. All conditions are equal for both load cases, except for the wind condition that describes an idling turbine in DLC 6.4. Basically DLC 6.4 covers the cases where the wind speed at hub height averaged over 10 minutes, denoted as V hub , is smaller than the cut-in speed or larger than the cut-out speed and DLC 1.2 covers all other cases. It is convenient to bundle both cases in one single scenario and make no assumption about V hub .
Given the mean wind speed, the turbulence intensity is defined by means of the normal turbulence model (NTM), which implies that the turbulence standard deviation σ 1 is given by the 90% quantile for the given wind speed at hub height, i.e.
The turbulence intensity, which should be used to generate an inflow wind field using the methods from Veers [34] or Mann [24], then equals σ 1 /V hub . The lateral component of the turbulence standard deviation, denoted as σ 2 , and the upward component, denoted as σ 3 , are defined as σ 2 = 0.8 σ 1 and σ 3 = 0.5 σ 1 respectively (here, it is assumed that the Kaimal model is used). I ref denotes the expected value of the turbulence intensity at 15 m/s, with I ref = 0.16 for the case described in this article. The wind direction θ wind with respect to the turbine is assumed to be uniformly distributed between −12°and 12°, as the fatigue load case under consideration does not model extreme yaw misalignment. Notice that wind direction 0°is equivalent to no yaw misalignment. The wind direction and wind speed are independently distributed random variables.
The sea state is incorporated using a joint probability distribution of the significant wave height H s , the peak spectral period T p , and the wind at hub height V hub . A normal sea state is considered, so no anomalous events are incorporated and currents are ignored. The wind and wave directionality (and their dependence) have to be incorporated as co-directional or multi-directional (i.e. they are either equal or not). In this work, the latter option is chosen. For this purpose, the wave direction θ wave (with respect to the turbine) is used to determine the wind/wave misalignment: The distribution of the wind/wave misalignment can be determined by subtracting the stochastic wave direction (which should be inferred from the site conditions) from the uniformly distributed wind direction. Combining the parameters describing the wind and sea conditions yields a five dimensional parameter space, called Ω in this work, consisting of all possible combinations of those parameters, called x ∈ Ω in this work. It consists of the wind speed at hub height denoted by V hub , the wind direction denoted by θ wind , the significant wave height denoted by H s , the peak spectral period denoted by T p , and the wind/wave misalignment denoted by θ wind/wave . Hence x = (V hub , θ wind , H s , T p , θ wind/wave ) T . The concrete goal of DLC 1.2 and DLC 6.4 is to assess the effect of regular, long-term variability of these parameters on the forces acting on the wind turbine.

Meteorological mast IJmuiden measurements
The environmental conditions used for the load case calculations are based on the measurements done using meteorological mast IJmuiden [36]. This meteorological mast is situated approximately 85 km from the Dutch coast (at coordinates 52 • 50 53 N 3 • 26 8 E). At the measurement site the water depth is approximately 28 meters, see Figure 1 for an illustration. The measurement campaign started in November 2011 and lasted until March 2016. Up to some interruptions, the wind conditions at various heights, wave conditions, and ocean current conditions at various depths have been obtained using a lidar and a buoy throughout that period of time.
For the purpose of this work, we are mainly interested in the statistical behavior of the wind in relation with the sea state. To this end, a straightforward preprocessing step is performed to combine the hourly measured sea state with the 10-minute averages of the wind conditions at that specific moment. After removal of all invalid measurements (due to measurement errors or interruptions), 24 650 samples of the wind speed and wind direction at hub height (which is 90 m, see the next section) and significant wave height, wave direction, and peak spectral period are obtained. These samples of the combined wind and wave conditions constitute 88% of all available measured wave directions (measured hourly) and 10% of all available measured wind directions (measured every 10 minutes).
The obtained measurements are visualized in Figure 3. Notice that Figure 3a clearly shows that there are no measurements with V hub = 0 m/s, i.e. the measurement devices have a small tolerance resulting in a measured wind speed that is unconditionally positive. Figure 3b demonstrates that the dominant wind and wave direction is south west. There is also a significant number of waves from north, possibly due to the open connection with the North Sea at the offshore site of the meteorological mast. Nonetheless, waves from north rarely occur in combination with wind from south west, as the wind/wave misalignment has approximately zero mean.

NREL 5MW reference offshore wind turbine
The NREL 5MW reference wind turbine for offshore system development is a fictional offshore wind turbine designed to support concept studies aimed at assessing offshore wind technology [21]. It is a horizontal-axis wind turbine rated at 5 MW with a rotor diameter of 126 m. Its hub height is 90 m, which allows for the straightforward usage of the measurement data discussed in the previous section, as measurements at 90 m are available. The cut-in, rated, and cut-out wind speed of the turbine equal 3 m/s, 11.4 m/s, and 25 m/s respectively. As the mean wind speed of the meteorological mast at 90 m is 10.13 m/s, simulation of the NREL 5MW wind turbine at the location of the IJmuiden mast describes a realistic scenario of a power producing wind turbine. The geometry of the turbine is depicted in Figure 2.

Aeroelastic code BLADED
The simulation of the turbine is performed using the aeroelastic wind turbine modeling tool BLADED (version 4.6). One run of BLADED to obtain the forces acting on the turbine over a period of 10 minutes takes approximately 50 minutes on the hardware used for the numerical experiments, which significantly dominates the total computation time. All other algorithms discussed in this article (such as determining the bins, the quadrature rules, and all post processing procedures) are negligible in this regard.
Given specific values of the environmental parameters, BLADED is used firstly for simulating the wind turbine to obtain the forces and secondly for post processing the obtained results. Obtaining forces consists of two steps. Firstly, a turbulent wind field using the average wind speed at hub height and the turbulence intensity is generated. Then secondly, using this turbulent wind field, the random sea state, and the design of the NREL 5MW wind turbine, a time history of the forces acting on various components of the turbine is calculated. By repeatedly doing such calculations for various environmental conditions, and incorporating the frequency of occurrence of each specific set of environmental parameter values, weighted equivalent loads are calculated for various representative slopes of the S-N curve. An equivalent load is a representative load of a specific load time history resulting in the same damage. A weighted equivalent load is an equivalent load where the frequency of occurrence of several time histories is taken into account. Determining these loads is the goal of this work, which is a standardized way of reporting and assessing wind turbine performance.
An equivalent load is calculated using the time history of the forces acting on the turbine. It is out of the scope of this article to fully discuss all details, we refer interested readers to Burton et al. [8] and the references therein for a general introduction or the Theory Manual of BLADED [11] for the specific implementation details within BLADED.
There exists variability in the output of BLADED due to the random nature of the turbulent wind field and sea state. To model these, it is common to generate several wind fields and sea states and average the time histories of the forces over these realizations. These so-called seeds are also incorporated in the framework in this article.
It is important to emphasize that the procedure of determining equivalent loads and calculating the seeds is similar for both binning, which is the well-known conventional approach, and numerical integration by means of quadrature methods, which is the novel approach requiring significantly less BLADED evaluations suggested in this article. No modifications to the code or input files of BLADED are necessary, except for the frequency of occurrence of each individual simulation. This is a major advantage of the approach suggested in this article, as all existing well-developed frameworks can be used without significant changes (i.e. the approach is non-intrusive).

Numerical integration for load calculations
The procedure to determine all equivalent loads that are necessary to obtain weighted equivalent loads is the most costly part of DLC 1.2, as it requires the calculation of a large number of time histories of loads. In this section, it is demonstrated why this happens in case binning is used and a numerical integration framework is proposed as efficient and accurate alternative.
To introduce equivalent loads in our framework, let u : Ω → R be the map that yields an equivalent load (the quantity of interest) for specific values of the environmental conditions x ∈ Ω (see Section 2.1). To incorporate the uncertain nature of the environment, the equivalent loads are combined into a single weighted equivalent load, which is ideally calculated as follows: Here, ρ(x) is the probability density function modeling the uncertainty of the environment. The integer m is the slope of the S-N curve, where it is implicitly assumed that the S-N curve can be represented by a monomial with power m. The notation L eq is used to represent this (ideal) weighted equivalent load.
In the case studied in this article, it is not assumed that there is a known probability density function of the environmental parameters, but instead measurements are directly used (see Section 2.2). In this case, L eq can be introduced as follows: where y k ∈ Ω are all measurement points available (as depicted in Figure 3) and K = 24 650 depicts the number of measurement points. Each y k is a 5-dimensional vector describing the measured parameters (see Section 2.2). Both (3.1) and (3.2) can be used in the framework proposed here, i.e. either measurement data or a distribution can be used straightforwardly.
Usually it is not viable to evaluate (3.1) exactly, as u is a map that involves evaluating a complex aeroelastic model. Also (3.2) cannot be used, as it is intractable to evaluate the aeroelastic model for each measurement point. Therefore, the integral is replaced by a weighted sum, obtaining the following approximation: where x 0 , . . . , x N are N + 1 specific representative values of the environmental parameters, called nodes in this article, and w 0 , . . . , w N the weights. The weights model the frequency of occurrence of the nodes. The notation L eq N is used to represent a weighted equivalent load using N runs of BLADED. Ideally, we would have that L eq The standardized approach to determine the nodes and weights is by means of binning. In that case, x k are equidistantly spaced and w k equal the frequency of occurrence of each environmental condition. This approach is discussed briefly in Section 3.1.
One major disadvantage of binning is the large number of nodes (and concomitant model runs) that are obtained for the five-dimensional parameter space considered in DLC 1.2. An alternative is to determine the nodes and weights such that they form an interpolatory quadrature rule. This yields an accurate estimation with relatively small number of nodes. The challenge is to construct the rule that incorporates the measurement data without reconstructing a (possibly inaccurate) probability density function and without requiring model evaluations at all samples y k . The method used in this work is discussed Section 3.2.
In (3.1), (3.2), and (3.3) it is assumed that the function u can be evaluated exactly. This is generally not the case, as the obtained equivalent loads still depend on the generated wind field and random sea. This is often solved by using seeds, which consists of constructing various wind and sea states with similar environmental parameters and averaging the obtained equivalent loads. Since the obtained equivalent loads are used in a weighted sum, it is a natural idea to use more seeds for nodes with high weight and vice versa. With this idea, seeds can be incorporated in the quadrature rule framework. Details are discussed in Section 3.3.

Binning
Binning is the conventional approach to assess variability in the wind and waves. The idea is to evaluate the aeroelastic code at equidistantly placed locations and post process the results using the probability density function of the variables under consideration. The bins can be interpreted as quadrature rule nodes, yielding a rule with equidistant x k . The weights w k describe the frequency of occurrence of x k .
If a distribution ρ(x) is known, the weight of the k-th bin b k is typically determined as follows [19, Annex G]: If measurements are considered, the weights are the fraction of measurements in the bin of x k , i.e.
where K k is the number of measurements in the bin of x k and K is the total number of measurements. An example of a quadrature rule obtained by applying binning is illustrated in Figure 4a, where the same measurement points as in Figure 3a have been used to determine the bins. The weights are illustrated by means of colors and nodes with zero weight are not depicted. The mean wind speed at hub height is binned with bin sizes of 2 m/s and the significant wave height is binned with bin sizes of 0.5 m. Notice that many weights equal zero, requiring no concomitant model runs. Nonetheless, the majority of the nodes (117 out of 210) has possibly small, but positive weight.
The advantage of this approach is that it is versatile and straightforward to implement. Moreover it is supported by many software packages for wind turbine simulation. However, if multiple environmental parameters are considered the number of bins increases exponentially. To see this, notice that if d parameters are considered with B bins in each direction, the total number of bins is proportional to B d , even if bins with zero weight are neglected. This severely limits the applicability of binning if a large number of parameters is considered.

Interpolatory quadrature rules
As an alternative to binning, we consider interpolatory quadrature rules. The nomenclature and mathematical framework of these rules is introduced in Section 3.2.1. In Section 3.2.2 the quadrature rule used in this work is explained.

Nomenclature
The key property of an interpolatory quadrature rule is that it is constructed such that it integrates all polynomials up to a certain degree exactly. Hence if x 0 , . . . , x N and w 0 , . . . , w N are the nodes and weights of a quadrature rule of degree D, then for all polynomials ϕ of degree D or less. The degree of a multivariate polynomial is determined as the sum of its powers, e.g. ϕ(x) = V 3 hub θ 7 wind has degree 10. The probability density function ρ does not appear in the definition of the quadrature rule directly, but is incorporated in the nodes and the weights. If nodes x 0 , . . . , x N are given, an interpolatory quadrature rule of degree N can be obtained by solving the following linear system for the weights w 0 , . . . , w N : Here, ϕ j are monomials, i.e. polynomials of only one term (e.g. the V 3 hub θ 7 wind example above). In the univariate case (d = 1), it holds that ϕ j (x) = x j . Multivariate polynomials are sorted graded lexicographically [37] to accommodate quadrature rules of any number of nodes N (we will call such a rule, with a little abuse of nomenclature, a rule of degree N ). The values µ j are the raw moments of the distribution, defined as follows: which uses (3.1). If measurements are used, the raw moments are defined as follows: where y k constitute all available measurements. The latter definition, which is employed in this work, uses (3.2).
A quadrature rule of high degree is not necessarily accurate if applied to functions that are not a polynomial. On the other hand, if the linear system from (3.4) is satisfied (i.e. the quadrature rule has degree N ) and all weights are non-negative (i.e. w k ≥ 0 for all k), then the following inequality holds [7]: This inequality is commonly known as the Lebesgue inequality. The right hand side of the inequality, i.e. min ϕ − u ∞ , is the minimal distance between the aeroelastic model u and any polynomial ϕ. It is notoriously difficult to determine this distance [35] (even if the model is known analytically). This is not a major issue, since the bound tells that the error made by integration using an interpolatory quadrature rule is at most twice as bad as the smallest distance between the model and any polynomial, even if that distance is unknown.
In other words, if the computational model can be approximated by means of a polynomial, a quadrature rule is a viable numerical integration methodology. A model can be approximated by a means of a polynomial if it is continuous and if its domain is bounded, which is the case for BLADED and the parameter space Ω considered in this article. Many more convergence results considering quadrature rules and polynomial approximation exist. The interested reader is referred to Brass and Petras [7] and Watson [35] and the references therein for more information.
There exist various approaches to construct interpolatory quadrature rules with positive weights. For univariate cases, the Gaussian quadrature rules [17] and Clenshaw-Curtis quadrature rules [10] are arguably the best-known rules. In higher dimensional spaces, constructing quadrature rules has similar limitations as when using binning: the number of nodes grows exponentially as the number of parameters grows. To alleviate this, the Smolyak sparse grid can be used [30,33], although this approach requires that all parameters are independently distributed (which is clearly not the case, as described in Section 2.1).
In this work we use the implicit quadrature rule to determine the nodes and the weights [5]. This quadrature rule is constructed as a subset of a provided sample set, so it is very suitable for the measurement data that is considered in this work. It has the advantage that its number of nodes does not grow exponentially fast if the number of parameters is increased, which is the case with binning. Moreover, it has positive weights and therefore yields an accurate estimation of E[u], whose error can be efficiently estimated. It is briefly explained in Section 3.2.2.
An example of the implicit quadrature rule, determined using the measurement points from Figure 3a, is depicted in Figure 4b. The rule is chosen to consist of 117 nodes (the number of bins with non-zero weight). It is clearly visible that the nodes have positive weights and to a certain extent represent the distribution of the measurement points.

Implicit quadrature rule
Consider K sample points y 0 , . . . , y K−1 to be given. The goal is to construct, for a desired number of nodes N (with N K), a quadrature rule with nodes x 0 , . . . , x N (with x k ∈ {y k }) and positive weights w 0 , . . . , w N such that Notice that the right hand side of this equality is formed by the raw moments from (3.6), which are simply the averages of all ϕ j determined over the samples. As ϕ j are polynomials, this estimate can be determined for very large values of K without any significant computational cost.
The key idea of the implicit quadrature rule is to interpret the right hand side of (3.8) as a quadrature rule with K nodes x k = y k and weights w k = 1/K for all k. Then the Vandermonde-matrix of (3.4) is not a square matrix, but is the following (N + 1) × K-matrix: T . Such a null vector can be exploited to construct a smaller quadrature rule [4], by noticing that where α ∈ R is an arbitrary constant (any multiple of a null vector is a null vector). A smaller quadrature rule can be obtained by selecting α such that w k − αc k = 0 for one k. The goal is to keep all weights positive, i.e. w k − αc k ≥ 0, which is achieved by choosing α as follows: This value is firstly such that w k − αc k ≥ 0 for all k, secondly such that there exists a k = k 0 with w k0 − αc k0 = 0, and finally such that {w k − αc k } can be used as weights of the quadrature rule, since c is a null vector. Therefore, the k 0 -th node and weight can be removed from the quadrature rule by updating the weights using w k ← w k − αc k (as w k0 − αc k0 = 0). Repeatedly applying this procedure to the quadrature rule with K nodes leads eventually to a quadrature rule of degree N with N + 1 nodes, which was the original goal. Note that there are always two possible nodes that can be removed, as both c and −c are null vectors and (3.10) only uses the positive values of the null vector. Eventually the rule consists of N + 1 nodes and the matrix from (3.9) becomes a square matrix. In that case no more nodes can be removed from the rule without also removing a function ϕ j that will be integrated accurately by the rule. This yields the nodes x k and weights w k that can be used in (3.3) to accurately estimate equivalent loads.
The procedure to remove nodes can also be used to approximate the error made by a quadrature rule. The trick is to, given a quadrature rule, iteratively remove a function ϕ j (the last row of (3.9), this decreases the degree of the rule) and subsequently remove a node from the quadrature rule. Then a sequence of quadrature rules with positive weights of decreasing degree and decreasing number of nodes is obtained, such that the following expression can be used to estimate the accuracy of the rules: Here, n denotes the number of nodes of the rule upon removal of nodes (namely N −n nodes). Determining e n for several increasing n (with n < N ) gives insight in the decay of the error made by the quadrature rule, as e n is bounded by the errors made using A n and A N : (3.11) The left hand side of the inequality describes that e n is larger than the difference of the errors made by A n and A N . The right hand side of the inequality describes that e n is smaller than the sum of both errors. If the quadrature rule is applied to a function u that can be numerically integrated using a quadrature rule (i.e. (3.7) is a good error estimate), a small e n implies that the error of the quadrature rule is small. An even better estimation of the error can be obtained by repeatedly determining e n for different sequences of quadrature rules, and taking the average of the obtained values. This is possible as multiple sequences exist given the freedom in choosing c or −c.

Seed balancing
Both binning and quadrature rule approaches yield nodes that describe values of the environmental parameters that should be evaluated using the aeroelastic code. These evaluations have been denoted by the function u(x). However, these evaluations still contain a degree of randomness: for each set of parameters, a certain number of seeds must be evaluated. This involves three steps: firstly various turbulent wind fields and sea states are constructed, secondly the loads are determined for each wind field and sea state, and finally the obtained loads are averaged. Hence, if S k seeds are used to determine u(x k ), this can be denoted as follows: where u (s) (x k ) is an evaluation of the aeroelastic code using the s-th wind field and s-th sea state. Ideally, we would like to accurately approximate E[u ∞ ], with The approximation used in this work is A N [u S ] for finite N and S. Here, S = k S k denotes the total number of BLADED evaluations used in this quadrature rule estimate. It yields the following error estimate: The focus in this section is on the seed error, since the quadrature error has been discussed extensively in Section 3.2. It is common to use a fixed number of seeds for each node (i.e. all S k are equal). However, an error of similar order of magnitude determined with a smaller total number of seeds can be obtained by selecting the number of seeds per node based on the frequency of occurrence of the node. In this way, the accuracy in the overall approximation of the integral is preserved with a reduced cost. In this section this intuition is derived mathematically. First, notice that the seed error behaves as follows [9]: Throughout this section we write A ∼ B if the growth or decay of A is asymptotically the same as B. The constant of proportionality is the variance of the integrand, which typically depends on x k and m (the inverse S-N slope, recall that the integrand in (3.3) is u m ). It is omitted from this estimation, but can straightforwardly be incorporated in the procedure discussed in this section, provided that its dependence on x k and m is explicitly known. The computational time of BLADED does not depend on the specific wind field or sea state under consideration, so the total computational time scales linearly with the number of seeds. Hence given an accuracy goal E, the goal is to determine ε S0 , . . . , ε S N that firstly minimize the total computational time S, denoted as the total number of seeds used for all calculations, and secondly achieve the accuracy goal, denoted by the seed error: Notice that for N → ∞, the asymptotic error of this expression decays to 0. Even though the error estimate is asymptotic, the relation S k = ε −2 S k will be used throughout this work. The problem sketched in (3.14) is a constrained minimization problem, that can be solved using the method of Lagrange multipliers [3], which yields the solution where A is a scaling constant such that N k=0 ε S k w k = E. Hence the number of seeds for the k-th node should ideally equal The number of seeds S k is larger if the associated weight w k is larger and no seeds (hence no BLADED runs) are necessary if w k is zero. The latter is straightforward to observe: if the weight of a node is zero, it does not contribute to the quadrature rule approximation. Notice that the relation between S k and w k is non-linear due to the non-linear relation between ε S k and S k . Only in rare cases it holds that all values of S k as defined by (3.15) are integers. It is difficult to solve the optimization problem stated here with the restriction that all S k are integer, so in this article the number of seeds is simply rounded up such that the accuracy goal is reached in all cases.
The effect of the seeds and their dependence on the weights is illustrated in Figure 5 and summarized in Table 1, where the quadrature rules from Figure 4 are reconsidered. Here it is assumed that by default five seeds are used per bin or node, which is equivalent to an accuracy goal of E = 1/ √ 5. If binning is used, the total number of seeds reduces from 585 (five per bin with non-zero weight) to 323, which is only 56% of the original cost. If the proposed quadrature rule from Section 3.2.2 is used, the total number of seeds reduces to 341, which is 58% of the original cost.

Numerical results
In order to demonstrate the benefits of the proposed quadrature rule, the five-dimensional DLC 1.2 is considered and the approach is compared to binning. However, it is infeasible to assess the complete load case using binning, as it results in an excessive number of runs of BLADED. Therefore, first two simplified test cases are considered for the comparison. The first test case, which is discussed in Section 4.1, consists of assessing the accuracy of the quadrature rule with respect to binning if the "model" u is a known test function. The second test case, which is discussed in Section 4.2, consists of the DLC 1.2 load case, where only the wind conditions are considered uncertain (and the sea state is fixed at nominal conditions). The full five-dimensional load case is then assessed purely by means of a quadrature rule and contains the five aforementioned parameters, modeling both the wind and the sea state. The results of this case (which follows the IEC standard [19,20] to a large extent) are discussed in Section 4.3.

Test functions
The five-dimensional Genz test functions [15] are considered, which are six functions u 1 , . . . , u 6 constructed specifically to assess the accuracy of integration methods (see Table 2). The goal is to calculate E[u] defined by measurement data (see (3.2)), i.e. the goal is to use binning and an implicit quadrature rule to estimate The data under consideration are the measurements from the IJmuiden meteorological mast (see Section 2.2), scaled to the domain of definition of the Genz test functions. Therefore,ŷ k denotes y k after scaling of all data to Ω = [0, 1] d (with d = 5). As the functions u f are known exactly, no seeds are employed in this test case.
Binning (see Section 3.1) is performed by fixing a number B = 1, . . . , 7, constructing B bins in each direction (obtaining B 5 bins), and removing all bins with zero weight. Then the implicit quadrature rule (see Section 3.2.2) is constructed such that its number of nodes is equal to the number of obtained bins. The vectors a and b, that define the Genz test function under consideration, are selected randomly in the unit hypercube and a is scaled such that a 2 = 5/2. The experiment is repeated 100 times (i.e. with Table 2: The test functions from Genz [15]. All d-variate functions depend on the d-element vectors a and b. The vector b is an offset parameter to shift the function. The vector a describes the degree to which the family attribute is present.

Integrand Family Attribute
Corner Peak Here, u f is the f -th Genz test function constructed with the k-th pair of random vectors a and b and N denotes the number of bins with non-zero weight or, equivalently, the number of nodes of the quadrature rule. The obtained results are depicted in Figure 6.
First, notice that the quadrature rule consistently outperforms binning. The convergence behavior of binning is approximately equal for all functions, as binning does not leverage smoothness of the function due to its local approximation. Moreover a large number of bins is necessary to obtain an accurate estimation. On the other hand, the convergence behavior of the quadrature rule clearly reflects the attributes of the test functions (see Table 2). The functions u 1 , u 2 , u 3 , and u 4 are smooth and can be approximated accurately using a polynomial of low degree. The integration error of these functions decays fastest. The fifth and sixth Genz test functions are not smooth and it is not straightforward to approximate these functions by means of a polynomial. This can also be observed in the convergence behavior of the quadrature rule, as it clearly converges slower than the other test functions. 3 6.00 · 10 5 N 5.93 · 10 5 N 5.87 · 10 5 N 5 5.63 · 10 5 N 5.57 · 10 5 N 5.49 · 10 5 N 10 6.41 · 10 5 N 6.37 · 10 5 N 6.22 · 10 5 N 12 6.76 · 10 5 N 6.70 · 10 5 N 6.53 · 10 5 N

Two-dimensional verification load case
The two-dimensional case under consideration is a case in which the wind speed and wind direction vary, but the sea state is considered to be fixed. The distribution of the wind speed and wind direction are as described in Section 2.1 and 2.2, so the wind speed follows the data and the wind direction is uniformly distributed between −12°and 12°. The parameters describing the sea state are fixed at their mean values obtained from the data, which are H s = 1.46 m, T p = 6.76 s, and θ wind/wave = −2.11°.
In this case, the conventional approach is to bin the wind speeds using bins of 2 m/s and to bin the wind direction into three bins of 8°. The wind speed varies between (approximately) 0 m/s and 30 m/s, yielding 15 bins in this dimension (these are also depicted in Figure 4a), resulting in a total number of 45 bins. The frequency of occurrence of each bin is determined using the number of measurements in the bin, which yields a non-zero weight for any bin. For a fair comparison, the implicit quadrature rule is constructed such that it contains 45 nodes based on the same two parameters.
For both cases, five seeds per bin (or node) are used and by means of rainflow cycle counting the weighted equivalent loads are determined. For the quadrature rule, the seeds are balanced using the procedure from Section 3.3. Seed balancing is not applied to binning, as we want to keep binning as close to the conventional approach as possible. The number of required BLADED runs of the balanced seeds in conjunction with a quadrature rule decreases from 225 to 173, which is (approximately) 77 % of the original cost. Some obtained loads for several components calculated using these three approaches (i.e. binning and quadrature rule with five seeds, and a quadrature rule with balanced seeds) are tabulated in Table 3.
The difference between the loads determined with a quadrature rule (both determined with five seeds or balanced seeds) and those determined with binning is less than 5%. Strictly speaking, it is not evident whether the quadrature rule or binning yields a more accurate answer in this case, as no exact values of the expected values of the loads are known. However, based on the results of Section 4.1, we have a strong indication that the quadrature rules yield more accurate estimates. The largest variation is in the loads on the yaw bearing, which is likely due to its large sensitivity to the varying environmental conditions. The contrary is true for the blade root, which is significantly less sensitive to environmental variations. This will also be observed in the next section, where the five-dimensional case is considered.
One main advantage of the quadrature rule approach is that its accuracy can be assessed without any additional costly aeroelastic calculations using the procedure discussed in Section 3.2.2. To this end, equivalent loads have been determined using N = 1, . . . , 44 nodes, and these are compared with the value from Table 3. As discussed in Section 3. Here, L eq N,k refers to calculating the equivalent load using N nodes of the k-th quadrature rule sequence. The error is scaled with the loads from Table 3, such that loads of different order of magnitude and with different characteristics can be compared.
The obtained quadrature errors are depicted in Figure 7. Notice that all reported errors decay approximately three orders of magnitude, to a point where the error is smaller than 1 % of the equivalent load. The errors decay rapidly for small N and the rate of decay reduces for increasing N . This is likely due to the fact that the error decays algebraically (as u 5 did in Section 4.1), i.e. its decay is approximately a straight line if depicted on a log-log plot. However, significantly more BLADED runs are necessary to confirm this claim numerically (e.g. see the number of nodes considered in the construction of Figure 6). The left figure depicts the convergence for the equivalent loads of the rotating hub for various slopes of the S-N curve. It demonstrates that the behavior of the error is similar for different values of the S-N slope, as all lines decay rapidly. The right figure demonstrates that the behavior of the error is more or less independent from the specific component under consideration.

Five-dimensional Design Load Case
The five-dimensional load case consists of DLC 1.2 as described in the IEC standard [19,20]. The random parameters are in this case the five aforementioned parameters, consisting of wind speed, wind direction, wind/wave misalignment, significant wave height, and peak spectral period.
Applying binning to this case is infeasible, as the number of bins with non-zero frequency of occurrence equals 3 410, which would result in approximately 15 000 evaluations of the aeroelastic simulation code. Even after optimizing the number of seeds using seed balancing from Section 3.3 the number of evaluations is still 8 961, which remains intractable.
Instead, an implicit quadrature rule is determined consisting of 100 nodes. This number of nodes is larger than the number of nodes from the previous section, to account for the higher dimensionality of the problem. Furthermore the number of nodes is sufficient to obtain a small error (which can be assessed during the simulation using the methods discussed in Section 3.2.2). The seeds are optimized using the method from Section 3.3, resulting in 419 individual aeroelastic simulation code evaluations. This is a major reduction compared to the number of bins, since the number of BLADED evaluations is less than 5% of that of binning. The equivalent loads acting on various components are summarized in Table 4.
The loads of the rotating hub and blade root are for the wind turbine and offshore site under consideration relatively close to the values that were calculated in the two-dimensional load case. These loads are apparently less sensitive to environmental variations which could be due to their high position on the turbine. The loads acting on the yaw bearing are significantly larger if a random sea is considered.  As the yaw bearing is the component that directly links the tower (which is connected to the sea) and the rotating hub (which is connected to the blades), its loads are sensitive to variations in both the wind and the sea. Contrary to the two-dimensional case, the sea state is random, so these variations are much larger.
To assess the accuracy of the obtained loads, the convergence is assessed in a similar way as done in the previous section: by removal of nodes from the quadrature rule with 100 nodes, sequences of quadrature rules are constructed that are used to assess the accuracy. There is a minor subtlety in this case: the number of seeds is optimized, so the nodes of the smaller quadrature rules are evaluated using an incorrect number of seeds. We do not correct for this, so the estimated error might be slightly larger. The results are gathered in Figure 8.
Compared to the two-dimensional test case, the errors of the quadrature rules are slightly larger, which is either due to the larger dimensionality of the problem or due to the unbalanced seeds used to construct these figures. Nonetheless, the errors clearly decay to a relative error smaller than 1%. For larger slopes of the S-N curve, the error is larger, which is likely due to the higher power used in the expression that is integrated (this increases the constant of proportionality discussed in Section 3.3). Again, the convergence behavior is possibly algebraic, though significantly more simulations are necessary to fully quantify the error. This is out of the scope of this work.
Concluding, this test case demonstrates that the proposed quadrature rule is capable of accurately approximating weighted equivalent loads in standardized form using less than 500 BLADED evaluations, which is a significant reduction compared to the approximately 15 000 BLADED evaluations that binning requires.

Conclusion
A novel approach based on quadrature rules has been proposed to calculate wind turbine fatigue loads as described by Design Load Case 1.2, which is standardized by IEC [19,20]. The quadrature rule under consideration is the implicit quadrature rule, with the key properties that it has positive weights and can be constructed using measurement data. It is based on polynomial approximation, which leverages smoothness in the model to achieve high accuracy. To demonstrate the efficiency of the new approach, it has been compared to binning, the conventional approach. In both approaches, the number of seeds per bin or node can be balanced to maximize accuracy in the available computational time. The environmental parameters are based on real offshore measurements in the North Sea and the wind turbine under consideration is the NREL 5MW reference wind turbine.
Both quadrature and binning approaches have been applied to a simplified two-dimensional load case, for which it has been demonstrated that the accuracy of the quadrature rule is comparable to that of binning. Moreover, the error of the quadrature rule converges algebraically upon addition of nodes.
The main advantages of the quadrature rule, i.e. a significant reduction of computational time with similar accuracy, have been demonstrated by considering the full load case, which is governed by five random parameters. In this case, binning is infeasible. The accuracy of the quadrature rule has again been assessed numerically, confirming that the error of the rule decays rapidly.
Throughout this work, all results are generated based on five seeds per node. A possible improvement is to vary the number of seeds per node, for which a further study of the significance of the seed error compared to the quadrature error needs to be performed.
Overall, the results demonstrate that for fatigue load cases our proposed quadrature rule forms a highly promising alternative to binning with significantly lower computational cost and similar accuracy. Therefore, a topic for future research is to extend the numerical integration framework to other load scenarios, e.g. fatigue load cases with more uncertain parameters (where the benefit of using a quadrature rule is even larger) or ultimate load cases (where no rainflow cycle counting is applied).