Effect of uncertainty sources on the reliability level of wind turbine rotor blades

1. Abstract Wind turbine rotor blades are large composite structures performing most of their design life under random cycle loading patterns. Concurrently material properties of the constituent plies exhibit inherent variability. In order to ensure a safe and cost effective design, uncertainty related to the basic variables (material properties, loads etc.) should be quantified and taken into account in the design calculations. Recently the blade design standard DNVGL-ST-0376 was released, showing the trend of using dedicated probabilistic analysis for wind turbine blade design. A critical evaluation of the new edition of the standard is performed herein, particularly in terms of the ability of the suggested safety factors to satisfy the target failure probability level of 1E-04. To this end, probabilistic analysis methodology is employed, starting with the measurement uncertainty for the static and fatigue properties of the composite material and going all the way up to the blade failure probability on a layer by layer basis. The application is performed on the INNWIND.EU reference 10MW rotor blade of 90m length. Furthermore, the current reliability level of the specific blade design, following the new standard DNVGL-ST-0376, is estimated considering the various failure modes i.e. fibre failure (short term and fatigue strength), buckling and inter fibre failure of the composite laminates, while taking into account sources of variability that contribute to the physical, statistical, measurement and model uncertainty. Results indicate that while for the static (extreme) analysis deterministic results are conservative, the opposite is observed for fatigue analysis.


Introduction
Designing the 20MW wind turbine of the future with about 250m rotor diameter is a challenging task. The task becomes more complicated considering that wind turbine structures should continuously operate in a fully stochastic environment for their entire design life. To provide a safe and cost effective design, uncertainty quantification and probabilistic analysis become mandatory.
Reliability analysis of wind turbine rotor blades under extreme and variable amplitude loading has already been performed in the past [1]- [5]. The main effort, however, was put on the uncertainty quantification related to the environmental conditions e.g. wind speed and turbulence as well as the statistical treatment of the sectional stress resultants, with less emphasis placed on the uncertainty on the material properties. When mechanical properties were considered as random variables only part of the uncertainty was quantified; namely the physical and statistical uncertainty related to the limited sample size of experimental tests. Measurement uncertainty was totally ignored. The available standards for composite materials do not provide any indication of the measurement uncertainty obtained through application of the experimental method and therefore, this is difficult to consider in combination with design guidelines and standards for the blade design. This, in turn, leads to the introduction of partial safety factors, which are based on engineering judgement rather than experimental evidence. This is clearly evident in the new edition of DNVGL-ST-0376 design standard for wind turbine blades.
In the present work, measurement uncertainty for both static and fatigue properties of the composite material is quantified following principles of metrology and based upon appropriate experimental data [6]. This information is of high importance if comparisons are to be made between candidate materials or the effect of a parameter, whether environmental or other, on the material properties is to be sought. Moreover, to identify the role of measurement uncertainty on the reliability level of the blade, the probability of failure under static and variable amplitude loading is estimated based on probabilistic analysis tools developed for application on wind turbine blades in [7]- [9]. The application is performed on the 90m INNWIND.EU reference 10MW wind turbine blade [10]. The probability of failure against static and variable amplitude fatigue loading is estimated at the material ply level on several locations on each section of the blade taking into account the various sources of uncertainty i.e. physical, statistical, measurement and model uncertainty.
Having the methodology and the tools to reasonably quantify several types of uncertainty and perform reliability analysis at the blade structural level, namely the building ply of each laminate, the reliability level of the specific blade design is quantified following the new standard DNVGL-ST-0376. Since the new edition encompasses the difference between deterministic and probabilistic design it provides the perfect platform to locate the disparities between the two methodologies. In parallel, it allows also a critical evaluation relevant to the ability of the suggested safety factors by the new standard to provide a blade design with a target annual failure probability 1E-04.

Rotor blade
The case study is based on the INNWIND.EU 10MW reference (WT) [10], a horizontal, variable speed, collective-pitch controlled wind turbine with a three-bladed upwind rotor. Each blade expands to a radius equal to 89.17m and has a nominal mass of 41,788kg.
The structural design of the blade defined in [10] is based on a load carrying box girder with two shear webs and a cylindrical root. A third shear web is present close to the trailing edge from 21.8m to the tip. The material used is E-Glass/Epoxy composite and balsa wood sandwich panels. One unidirectional and two different multidirectional laminates were used to form the plies of numerous lamination schemes, namely a uniaxial with 0° fibre orientation (along the blade length), a biaxial with ±45° fibre orientation and a triaxial with combined 0° and ±45° fibre orientations. Details on the blade geometry and materials are given in [10]. For the purposes of this work, small modifications on the blade design were performed considering each ply of the multidirectional fabrics separately. This was done to allow application of the Puck failure criterion, which is directed for uniaxial laminas (layers). The assumption is supported by the fact that the layers used in the design of the blade are non-crimp stitched and not woven.

Mechanical properties of the composite ply
Lamina properties that need to be determined for further use in the structural design are: the elastic modulus in the fibre direction and perpendicular to the fibre direction, E1 and E2, respectively, the major Poisson's ratio, ν12, and the in-plane shear modulus G12, as well as the strength in tension and compression along the fibre direction and perpendicular to it, XT and XC, YT and YC, respectively, and the in-plane shear strength S. The present analysis is focused on the property determination based on the ISO standard series. For the fatigue properties, constant amplitude (CA) fatigue tests are performed to estimate the parameters of the SN curve. According to wind turbine blade design standard, DNVGL-ST-0376 [11], the power law model (SN curve) is fitted to the experimental data. Its formulation is given herein by: Where Nf the number of cycles to failure, Smax the maximum (or absolute minimum) stress level, ε the error of the model following normal distribution, N(0,se 2 ), with se the standard deviation of the error, a=e C , and ε1=e ε . The fatigue experimental data are derived following ISO 13003 [12], with the exception of the formulation of the SN curve.
Although the standard practice for model fitting involves the least square fitting method, for the purposes of this work, the maximum likelihood estimators (MLE) are used to determine parameters of the power law model C, b and se. It should be noted that MLE estimators are equivalent to the least square ones when the normality assumption for the error ε holds. The normality assumption was investigated by implementing the Kolmogorov-Smirnov test on the calculated residuals of the fitting procedure.
For both static and fatigue properties implementing the ISO standard series, measurements are performed for the quantities shown in Table 1. Further, more complex measurands i.e. quantities that are not measured directly but are rather determined through the above mentioned measurements are shown in Table 2.

Statistical uncertainty
In Table 2 the subscripts 1 and 2 refer to the principal axes of the composite lamina, i.e. along the direction of the fibres and perpendicular to them, respectively, while subscripts x and y in shear stress/strain and shear modulus equations refer to the loading direction of the ±45 o specimens and transverse to it. Furthermore, the subscript F stands for the value of the quantity in Table 2 at the failure of the specimen and the exponents (1) , (2) correspond to different measurements. More specifically in Eq. (6) σi (1) is the stress in the i direction at axial strain value of εi (1) =0.05% while σi (2) is the axial stress at a strain value of εi (2) =0.25%. Moreover, in Eq. (8), the shear stress τ12 (1) corresponds to a shear strain value of γ12 (1) =0.1% while τ12 (2) to a shear strain value of γ12 (2) =0.5%.
Shear stress MPa Shear strain % 12 12 = − Modulus of elasticity GPa Poisson's ratio --- Shear modulus GPa The mean values of the material properties calculated during the respective test series are reported, along with the corresponding standard deviation.
Depending on the starting scale at which uncertainties are built up in the composite materials, three main approaches are encountered [13]. Uncertainties can be modelled from a micro-scale level where the macroscopic material properties of the lamina are obtained using various micromechanical models. Uncertainties can also be quantified from the meso-scale or the ply level and it is related with the uncertainties introduced when mechanical characterization of the lamina is performed. Finally, the macro scale level involves uncertainties determined by experimental tests similar to that required for the ply characterization but for more general layup schemes for the laminates or for subcomponents tests. Herein, uncertainties were modelled at the meso-scale level based on the experimental data derived by coupon tests.
To quantify uncertainty in the material properties, except from the experimental data, detailed information about the performance of the tests is needed. Due to the limited information related to the INNWIND.EU material properties, the OPTIDAT database [14] is used instead. To be more specific, testing results with full information about the testing procedure followed for their derivation, for which the author had complete access, concerned the measurements from one laboratory reported in [15]. Tests were performed on unidirectional specimens with stacking sequence [0]4 (or [90]7 for the properties in the transverse direction) from one batch of material. For the determination of the shear properties, axial testing on [±45]S coupons were performed. Therefore, the outcome of the present statistical analysis can be viewed only as a favourable scenario missing some part of the introduced uncertainty, as e.g. plate to plate variability (coupons cut from one plate) or laboratory to laboratory variability.
For the fatigue properties of Glass/Epoxy material, a database to perform the various statistical analyses was provided by Fraunhofer-IWES from a previous industrial project [6]. More specifically, data concerned measurements of the number of cycles to failure on several stress levels for three different ratios of minimum to maximum load (R-ratio) namely R=0.1, -1 and 10 for constant amplitude fatigue tests. Tests were performed on unidirectional specimens with stacking sequence [0]4. More information about the constituents of the material as well as the derivation of the experimental data can be found in [6].
It is assumed that during the performance of the tests every effort was made to identify and correct possible systematic effects for the creation of both databases.
Summarising, for the probabilistic calculations implemented in this work, the OPTIDAT database is used to quantify uncertainty i.e. to estimate the coefficient of variation (CoV) of each static material property, while the UD properties of the INNWIND.EU design [10] are assumed and assigned as the mean values of the material properties. Regarding the fatigue properties, i.e. parameters of the power law model, both mean and CoV values are derived by statistically analysing the database in [6].

Uncertainty quantification methodology
The uncertainty analysis is performed based on the principles of the Guide for Uncertainty Measurement (GUM) [16]. The combined standard uncertainty up i.e. the uncertainty of the measurands, is determined by the uncertainty propagation law given by: where p the quantity (measurand) determined by combining more than one of the measured (input) quantities xj, uxj the standard uncertainty of xj i.e. the uncertainty of the result of a measurement expressed as a standard deviation, f the mathematical model between the measurand p and input xj and r the correlation coefficient between xj and xk. The assumption behind Eq. 9 concerns the approximation of the function f using Taylor series expansion while keeping only the linear terms i.e. linearization of the function f.
The sources of uncertainty that were taken into account in the uncertainty quantification analysis are listed in Table 1. That is, uncertainties due to test instruments (resolution of calliper), specimen (dimensions, tolerances), test procedure (standard's recommendations) and statistical uncertainty due to the limited number of measurements e.g. thickness measurements. It should be highlighted that GUM [16] suggests a uniform treatment for all uncertainty components i.e. uncertainties due to systematic and random effects. That is, the result of a measurement after correction for recognized systematic effects is still only an estimate of the value of the measurand. Thus, GUM distinguishes only between methods used to evaluate uncertainty components i.e. type A and type B uncertainties and not the components themselves. The former are uncertainties that can be determined by analysing a series of independent repeated measurements while the latter are uncertainties that cannot be determined by repeated observations, but are rather based on scientific judgement.
Thus, for the first three sources of uncertainty (type B) of Table 1, a rectangular distribution was assumed, while a normal probability distribution for the (type A) statistical uncertainty. That is, regarding to the thickness/width measurement of the test coupons, one of the sources of uncertainty presented in Table 1 is related to the precision of the measuring instrument. For a calliper with accuracy ±0.02 mm, a rectangular distribution can be assumed implying equal probability for any value in the interval of precision. Those uncertainties propagate to the more complex measurands presented in Table 2. Because, in practise, series of tests is performed and the mean values of the measurand of Table 2 along with their standard deviations are reported, the associated probability distributions to the mean values (and the variances) of the measurand of Table 2 were assumed to be a normal distribution based on the Central Limit Theorem approximation. The analytical formulation of the combined uncertainty up of the static and fatigue material properties of Table 2 was developed and partially presented in [9], [17]- [19]. A brief discussion is performed herein, selectively for topics of interest.

Uncertainty on the mean value of a series of measurements
Following [20], the combined standard uncertainty on the mean value of a series of n measurements for a material property p can be derived by Eq. (9). Its formulation is given by: where sp is the sample standard deviation and up is the combined standard uncertainty related to the static material properties and is estimated as explicitly described in [17] while the various sources of uncertainty considered are presented in Table 1. Because for a given property p, e.g. XT every single measurement pi, where i=1,…n corresponds to the number of the consecutive measurements for the specific property, results in slight different up value, in Eq. (10), the average value of the combined standard uncertainty of the measurand is used.

Uncertainty on the variance of a series of measurements [9]
In a similar approach with the derivation of the combined standard uncertainty uμp, Eq. (10), the combined standard uncertainty for the variance of a series of n measurement for a material property p can be derived by Eq. (9). Its formulation is then given by: where μp, sp are the sample mean and standard deviation respectively and up is the combined standard uncertainty related to the static material properties i.e. measurands and is estimated as described in [17]. Again, the average value of the combined standard uncertainty of the measurand is used as for Eq. (10).

Uncertainty on measurement of applied load (fatigue) [18]
The force, F, is measured through a load cell on the testing machine, which should be verified and calibrated for example according to ISO 7500-1 [21], to achieve the requested accuracy specified within the ISO material testing standards for static tests. For the CA fatigue tests, ISO 13003 [12] does not make reference to any particular calibration standard, requiring only that the load reading device should be accurate to within ±2% of the full scale value for fatigue cases. Assuming a rectangular distribution, the type B uncertainty on the load is derived by: where Ffs is the full scale value of the load cell. According to the study in [18], it was found that the requirement of the standard for ±2% accuracy of the full scale value of the load reading device covers several uncertainties among them the ones that are related to the mean value and the amplitude on the load command.

Uncertainty on fatigue cycles [18]
In case the fatigue cycles are continuously counted and displayed up to the failure of the specimen, depending on whether the testing machine counts the cycle at the beginning or end of the fatigue cycle and whether the specimen fails at the beginning or end of the fatigue cycle, the uncertainty can be assumed in the range 1 cycle. Thus, assuming a rectangular distribution, the uncertainty in cycle counting is: Uncertainty introduced in the fatigue cycles due to test parameters such as frequency, heating of the coupon during test etc. might be neglected assuming that ISO 13003 [12] is strictly followed.
The uncertainty on the fatigue cycles as estimated by Eq. (13) exhibits a quite low value, especially when cycles to failure in the order of magnitude 1E+04-1E+06 are considered and therefore it can be ignored.

Uncertainty on fatigue parameters [18]
The slope and intercept parameters of Eq. (1) are estimated by means of MLE, as already discussed. Assuming the normality for the residuals, the analytical equations of the estimated parameters with respect to the input data i.e. pairs of maximum/minimum stress level (Smax) and the number of cycles to failure (Nf) can be found in several statistical books and reports as e.g. [20]. However, due to the introduction of complex mathematical expressions as well as due to the difficulty of directly applying Eq. (9) on these expressions, the combined standard uncertainty for the two parameters of the power law model (C, b) was estimated following numerical procedures and more specifically, the finite difference method. According to this method [20], the term ( ) of Eq. (9) can be approximated by quotients of finite differences and thus the product ( ) of Eq. (9) can by approximated by the values Δi given by: Where p corresponds to the slope (or intercept) parameter of the power law model, input quantities xi stand for the maximum/minimum stress level (Smax) that the tests were performed and u(xi) is the combined uncertainty for xi in which the uncertainties of Table 1 are taken into account and propagate up to the stress level. Eq. (9) can be rewritten in terms of Δi and is given by: Unless input quantities are correlated, the second sum on the right-hand side of Eq. (15) can be omitted. As already mentioned, the uncertainty on the fatigue cycles estimated by Eq. (13) were ignored due to the quite low values. Table 3 forms the basis for the development of probabilistic models for the material properties necessary for the reliability analyses performed. Under the heading Data the mean values, μp, as derived by the INNWIND.EU blade design report [10] and the coefficient of variation of each property, i.e. the ratio of the standard deviation over the mean value, CoV(p), of the material properties, as determined through the OPTIDAT database are presented. Applying the methodology of section 4.1 and more specifically, Eq. (10)-(11) on the OPTIDAT experimental data, the uncertainty of the mean value, expressed as a coefficient of variation, CoV(μp), and the uncertainty of the standard deviation CoV(sp) of each property is determined. Results are reported in Table 3 under the heading Propagation Law.

Elastic and static strength properties
In the reliability analysis, the physical uncertainty of the basic variables is taken into account by considering only the statistics of the experimental data, μp and CoV(p) of Table 3. By considering the asymptotic properties of the maximum likelihood estimators (ap-MLE) for the mean and standard deviation of the experimental data the physical and statistical uncertainty can be taken into account. Results are shown in Table 3 under the heading ap-MLE and they were derived by evaluating the inverse of the Fisher information. The correlation coefficient between the two parameters (μp, sp) indicated a quite low value for all the material properties and thus it was further ignored. To consider physical, statistical and measurement uncertainty, statistics under the heading Propagation Law of Table 3 are considered into the reliability calculations. Due to the limited data, the correlation coefficient between the parameters (μp, sp) was not possible to be determined implementing the approach of section 4.1. For this reason, the correlation coefficient was assumed equal to the one derived by ap-MLE method and thus for this case was ignored.
In line with JCSS [22], all material properties might be assumed to follow the lognormal distribution. Correlation between the material properties was taken into account in the reliability analysis by implementing the Nataf transformation and adopting the formulation developed in [23]. Part of the correlation matrix between the material properties was derived directly from the available OPTIDAT data. Only five from the thirty six required correlation coefficients were possible to be evaluated (ρE1-v12=0.21, ρE1-XT=0.55, ρE2-YT=-0.04, ρv12-XT=0.15, ρG12-S=0. 24). The remaining ones were assigned a low correlation value, equal to 0.1, ensuring a positive semi definite matrix. An example for the empirical cumulative distribution function (CDF) considering the various types of uncertainty of Table 3 can be seen in Figure 1a and Figure 1b for the tensile strength in the fibre direction (XT) and the shear strength (S), respectively. Both property distributions are depicted with the continuous black line. To estimate the CDF of each property introducing the quantified measurement uncertainty (or statistical uncertainty), a simulation procedure was implemented. That is, assuming that μp, sp are normally distributed independent random variables with mean values as indicated in Table 3 under the heading Data and CoV values given either under the heading Propagation Law (or under the heading ap-MLE depending on whether measurement uncertainty is considered or not), NB sample points are generated for the mean value and standard deviation of each property. Further, lognormal parameters are estimated for every set of μp, sp following basic statistical theory. Eventually, samples of size Nf for the property are generated for each set of the lognormal parameter values while the correlation matrix of the material properties is also taken into account. The derived sample (NB*Nf) is considered to be representative for the physical, statistical and/or measurement uncertainty and can be further fed into simulation techniques such as Monte Carlo to perform reliability analysis of wind turbine rotor blades. The various types of uncertainty namely physical, statistical and measurement uncertainty mainly affect the tails of the distribution.
A comparison of the distribution of the material property estimated following the uncertainty quantification analysis presented in section 4 as well as the standard, DNVGL-ST-0376 [11], is carried out. Tolerance bound theory was used along with Monte Carlo simulation as described in [24] to obtain the i th % fractile with confidence level 95% assuming normal distribution with unknown variance. Results are added in Figure 1 (the dashed line curve). Additionally, the distribution of the material property following GL2010 [25] corresponding to the lower confidence limit of the sample mean in order to estimate the i th % fractile with confidence level 95% assuming normal distribution with known variance was also superimposed (grey colour curve). The characteristic value of the material property i.e. 5% fractile, according to both DNVGL-ST-0376 and GL2010 is given by the intersection of the red straight line with the respective distribution curve. Depending on the quality of the experiments and the measurement uncertainty obtained, it is clearly seen that in some cases, the DNVGL-ST-0376 approach might be very conservative, as e.g. for the tension along the fibre direction (XT). In other cases, however, as for the shear strength, which has an increased uncertainty, above 1% for the mean value (CoV(μp)), the DNVGL-ST-0376 approach is almost coinciding with that obtained when calculating the measurement uncertainty. The effect of this difference, nevertheless, should be more systematically analysed with respect to the final result on the blade design.

Fatigue properties
Similar to the properties derived from the database of the static experiments, the results for the fatigue properties are presented in Table 4. The intercept C and slope b were estimated from the experimental data of [6]. Results are reported under heading Data, together with the variance se 2 , of the error ε of Eq. (1). To remind that tests were performed for three different R ratio (0.1, -1 and 10) i.e. the ratio of the minimum to maximum load during testing as recommended by DNVGL-ST-0376. More specifically, the standard suggests the use of the Shifted Goodman and piecewise linear constant life diagrams (CLD). In order to construct those diagrams, fatigue tests are necessary to be performed for R=-1 (tension-compression region) for the first type of CLD while additional tests at least for R=0.1 (tension-tension region) and R=10 (compressioncompression region) for the second type of graph according to the standard. Considering the asymptotic properties of the maximum likelihood estimators (ap-MLE) for the two parameters C, b and the variance of the error se 2 , the physical and statistical uncertainty due to the limited number of experimental tests can be taken into account. Results are shown in Table 4 under the heading ap-MLE. The procedure of section 4 was further applied on the experimental data of [6] and the coefficient of variation of the intercept, CoV(C), and the slope CoV(b) is determined. Results are reported in Table 4 under the heading Propagation Law. It is pointed out that the only source of uncertainty leading to the CoV values of Table 4 is the uncertainty introduced by the stresses as a result of uncertainty in measurements both of the load and the coupon area (width and thickness). The CoV(se 2 ) of the variance of the error as well as the correlation coefficients among the parameters of the power law model were assumed to be equal to the ones estimated by the ap-MLE procedure. Table 4 comprises the base to develop probabilistic models for the fatigue parameters as was the case of static properties. That is, when considering physical uncertainty of the fatigue properties only the error ε of Eq. (1) is taken into account as a stochastic variable while slope and intercept parameters are assumed constants and equal to the MLE estimations presented in Table 4. Furthermore, in a reliability analysis considering both the physical and the statistical uncertainty, the information under heading ap-MLE is used. Finally, assuming that the measurands (C, b) are characterized by a normal probability distribution with mean values the MLE estimations and standard deviation the combined standard uncertainty indicated under the heading Propagation law of Table 4, physical, statistical and measurement uncertainty are considered in reliability calculations. Samples representative of the physical, statistical and/or measurement uncertainty are generated and fed into a Monte Carlo simulation in a similar way to that expounded for the case of the static properties.  From Table 4, it is obvious that when introducing measurement uncertainty, the CoV values for the parameters of the power law model are increased in comparison to the respective ones when only statistical uncertainty is taken into account. The effect is more pronounced for the uncertainty budget of the slope parameter in R=10, that is in case of CA fatigue experiments in the compression-compression region. This reflects the general difficulty in performing fatigue tests in the specific sector and highlights also the need for the construction of a robust database concerning the specific experimental tests.
The effect of considering statistical and measurement uncertainty directly on the S-N curve formulation can be seen in Figure 2. The average S-N curve, as well as the characteristic S-N curve, i.e. 5 th % fractile with 90% confidence level according to DNVGL-ST-0376 are plotted in Figure 2 with a black continuous and a black dashed line respectively. The results were reduced by one standard deviation in the intercept C, Figure 2(a), and the slope parameter b, Figure 2(b). S-N curves considering only the statistical uncertainty are plotted in blue, while those by considering both the statistical and measurement uncertainty are plotted in red. Note that the 5% fractile for a normal distribution is obtained by reducing the mean value by 1.96 times the standard distribution. To facilitate clarity in Figure 2 the average value was only reduced by 1.0 standard deviation.
The consideration of the measurement uncertainty affects significantly the derived S-N curves.

Loads on the blade
The strength of the blade under extreme loading and elastic stability is examined for the power production state of the wind turbine, DLC 1.1, realising the Normal Turbulence Model for all wind speed bins in the interval between cut-in and cut-out wind speeds. More specifically, one 10 minute time series of stress resultants at various cross-sections along the blade length for twenty different wind-speed bins was provided by POLI-Wind group, Politecnico di Milano [26].
Load extrapolation was performed on the sectional stress resultants time series namely the flap and edgewise bending moments, the respective shear forces and the torsional moment following IEC 61400-1 ed.3 [27]. Axial force is mainly dictated by gravitational and centrifugal forces and thus it was considered as a deterministic variable. Local maxima were extracted from the 10 min simulations adopting the peak over threshold method. The threshold value was set according to IEC i.e. equal to the mean value plus 1.4 standard deviation of the time series for each wind speed bin. A time separation of 10 seconds was also imposed to ensure statistical independence between successive maxima. The 3p-Weibull was selected to fit local maxima for all the wind speed bins. The long-term exceedance probability of the extreme sectional stress resultants for the 10 min and 1 year reference period (T) was estimated.
Correlation between extreme stress resultants was taken into account in the reliability calculations by implementing the Nataf transformation and adopting the formulation developed in [23]. Correlation was estimated based on the load time series. More specifically, one correlation matrix was estimated for each of the wind speed bin considered in DLC 1.1. As an approximation, the final correlation matrix used in the reliability calculations was estimated as the average of the derived correlation matrices and is presented in Table 5. It is observed that Further, flap and edge moment stress resultants are neither uncorrelated nor fully correlated. Similar, behaviour was observed in other research works e.g. in [28]. Due to the limited number of available time series, statistical uncertainty was taken into account. The Bootstrap technique was used to quantify the uncertainty as described in [8]. Moreover, model uncertainties related to the aero-elastic load calculations, considered in [5], were also included in the probabilistic model and is given by: where Xdyn stands for the uncertainty related to modeling of the dynamic response for the wind turbine including uncertainty in damping ratios and eigenfrequencies, Xexp for the uncertainty related to the modeling of the exposure such as the terrain roughness and the landscape topography, Xst for the statistical uncertainty related to the limited amount of wind data, Xaero is related to the uncertainty in assessment of lift and drag coefficients, Xstress for the uncertainty related to the computation of stresses from the wind load and Fext for the physical and statistical uncertainty. For the extreme edge moment stress resultants, the annual long-term exceedance probabilities can be seen in Figure 3 considering all types of uncertainty. The load values have been normalized with respect to the maximum load value from all provided simulations. The characteristic value with a recurrence period of 50 yrs. as IEC standard suggests is also depicted in Figure 3. It is very interesting to notice the great influence of the model uncertainties of [5] to the annual distribution of the extremes. The CoV values for the annual loads considering all uncertainty types were estimated between 26-30% when the respective ones considering only physical uncertainty are 4-7%.
Concerning the analysis under variable amplitude loading, additionally to DLC 1.2 that is the normal operation of the wind turbine throughout its lifetime, the parked condition of the wind turbine, DLC 6.4, is also considered while load time series were provided in the range of [26 m/s, 40 m/s].
It should be highlighted that for the design of a wind turbine rotor blade all the DLCs specified in the standard [11] should be considered. Additionally, due to the multi-axial loading of the rotor blade and the different load carrying capacity of an airfoil on the various load directions, the entire load envelop of the bending moments (in the flap and edge direction) should be examined to identify the critical load directions. More specifically, according to [11], the load envelope should be discretised in at least 12 equally distributed bending moment directions while verification analyses should be performed for each one of them. Thus, the limitations of the present work concerns the examination of very specific DLCs. On the other hand, because of the unidimensional load extrapolation performed herein and in combination with the correlation matrix of Table 5, only one load directions was mainly investigated. This means that there may be a load combination of the flap and edge moment in other directions more critical than the one examined in this example. The multiaxial load modelling problem for wind turbine rotor blades has been dealt both in deterministic and probabilistic terms in e.g. [29] and [30] respectively.

Safety factors
In a classic deterministic design, the material properties are reduced while the loads are increased through safety factors to account for parameters which may deviate from the design assumptions during operation of the wind turbine blade. In a probabilistic design, safety factors are avoided by directly considering the variability of the basic variables. However, there are certain factors, the effect of which is not included in the experimental databases selected for the statistical characterisation of material properties.
The general material partial safety factor, γm, suggested by DNVGl-ST-0376 takes into account the dependence on the type of material, the processing, component geometry and influence of the manufacturing process on the strength. It should be noted that a similar formulation is also given in IEC 61400-1 [27], where the minimum recommended partial safety factor takes into account possible unfavourable deviations/uncertainties of the strength of material from the characteristic value, uncertainties in geometrical parameters etc.
In the present work the uncertainty of the component geometry (assumed invariable), or the influence of the manufacturing processing on the strength (assumed equal from plate to blade) were not considered. Therefore, to estimate as accurate as possible the reliability level of the blades it was decided to take into account specific partial safety factors for the material properties as shown in Table 6. The purpose is not to double count uncertainties but to complete the missing uncertainties of this analysis by considering the appropriate safety factors. In this work, the estimated reliability level following this approach is defined as the 'current' reliability level of the blade.
Safety factors suggested by DNVGL-ST-0376 are the base factor γm0, the partial reduction factor for criticality of failure mode γmC, the factor for irreversible long-term degradation γm1, the factor for temperature effects γm2, the factor for manufacturing effects γm3, the factor for the accuracy of the analysis method γm4 and the factor for the accuracy of the load assumptions γm5. The partial reduction factors used for each analysis type are shown in Table 6. It was assumed that γm3, γm4 and γm5 are equal to unity for all analyses performed. It is highlighted that accuracy of the structural analyses i.e. γm4, is considered by taking into account appropriate model uncertainties explained in section 8. According to DNVGL-ST-0376 the product of γm0 and γmC factors results in a design with an annual failure probability of 1E-04. Normally, both factors should be ignored in a structural reliability analysis. However, because in the multiplication γm0*γmC several uncertainties are hidden that partially were not accounted for in the coupon experiments e.g. uncertainties in geometrical parameters, while their effect cannot be identified and separated, as a compromise, the γmC factor was set equal to unity. The last column of Table  6 gives the reduction factor γm i.e. the product of all the assumed partial reduction factors that will be applied to the material properties sample values. Regarding the loads, the aero-elastic simulation results already include safety factors following GL2010 [25], which are not different than those in the new standard DNVGL-ST-0376 [11]. The partial safety factor values used to multiply loads in the aero-elastic calculations are presented in Table 7. It should be noted that for different limit states i.e. different types of analysis (short term fibre failure, fatigue etc.) and different techniques for estimating the maximum loads, different safety factors should be applied. More specifically, for serviceability and fatigue limit states a safety factor equal to 1.0 should be applied while for ultimate limit states, the partial safety factor considering DLC 1.1 should be equal to 1.25 if load extrapolation is performed. To compensate for this and be aligned with the new standard, the load sample values generated in the Monte Carlo simulation were reduced approximately 35% for the inter fibre failure and fatigue analysis and 10% for fibre failure and buckling analysis.
For the estimation of the reliability level of the blade, the model uncertainties for loads as described in [5] and implemented in the developed probabilistic tool, Eq. (16), were not considered to avoid double counting of uncertainties.

Probabilistic tool
The probabilistic analysis was performed using PROBUST [8], a dedicated tool for the reliability and structural analysis of composite rotor blade sections under ultimate and fatigue loading. The tool is operating in MATLAB environment by modifying FORTRAN subroutines of THIN [7] and CUFSM [31] into the Matlab Executable as well as developing new source code.
PROBUST consists of four distinct processes: The first process concerns the stochastic representation of the basic variables i.e. the in-plane thermo-mechanical properties of the unidirectional composite ply, the ply thickness, the developed sectional stress resultants as well as the gradient of operating and curing temperature and the moisture content. Concerning fatigue properties, the uncertainty analysis presented in the previous section was fully implemented in the current code.
The second process concerns the structural analysis of the rotor blade. The program can be used for stress and instability analysis in multi-cellular orthotropic beam sections using thin wall theory and finite strip method (FSM) respectively. Additionally, fatigue analysis is performed by selecting a constant life diagram and applying the rainflow counting method. Specifically, two CLDs were implemented in the code namely, the Shifted Goodman diagram and the Piecewise Linear CLD as suggested by DNVGL-ST-0376 [11]. Thus, the cross section into PROBUST is discretized using 2-node elements. The effective thermo-elastic properties for every different laminate in the blade section are evaluated by means of Classical Lamination Theory (CLT) and assigned to every element. Cross-section properties are then calculated providing the necessary structural input for an aero-elastic beam element analysis. Sectional stress resultants are computed performing aero-elastic simulations. The derived loads are further used to calculate the normal and shear stress resultants of each node in the laminated shell elements. The solution of the abovementioned calculations is based on the Euler-Bernoulli kinematic assumptions. The stress-strain fields are evaluated for each layer for all the nodes in the blade section using CLT.
Regarding buckling analyses, only the axial normal stresses are taken into account in the calculations. The specific type of analysis is performed iteratively for a user-defined range of half-wave lengths to determine the critical load factor and the corresponding mode shape.
For the fatigue analysis, two types of constant life diagrams were implemented in the code as recommended by DNVGL-ST-0376 namely, the Shifted Goodman diagram and the Piecewise Linear CLD. The number of allowable load cycles Nf adopting the shifted Goodman diagram (without safety factors) is given by: where Sa, Sm the amplitude and mean stress that the allowable number of cycles Nf need to be specified, a, b, ε1 the parameters of Eq. (1) for R=-1 and XT, XC the static tensile and compressive strengths in the fiber direction respectively. For the Linear Piecewise CLD the allowable number of cycles is given by: and 1 ( , +1 , , , +1 , 1 , 1 +1 ) + +1 − 2 ( , , , , 1 ) +1 − 3 ( , +1 , , +1 , 1 +1 ) +1 = 0 (19) where a1TT, b1TT, ε1TT the parameters of the power law model for the first known R-ratio R1TT in the tensile region, a1CC, b1CC, ε1CC the parameters of Eq. (1) for the first known ratio R1CC in the compression region, ai, bi, ε1i and ai+1, bi+1, ε1i+1 the parameters of Eq. (1) for ratios Ri and Ri+1 (when moving counter clockwise in the CLD graph) respectively for which experimental data were provided and r= (1+R)/(1-R). It is noted that Eq. (19) is a nonlinear equation with respect to Nf and it has to be solved numerically. When parameters of the power law model are considered random variables in Eq. (17)- (19), the probabilistic expression of the allowable number of cycles to failure is obtained. Assuming constants the intercept and slope parameters and ε1 equal to unit, Eq. (17)- (19) represent the classical deterministic expressions for Nf. The Palmgren-Miner damage index (D) is finally calculated at the ply level of the composite structure. It should be mentioned that only the axial stresses along the fibre direction are taken into account in the fatigue calculations presented herein.
For the study case in the present work, the section was discretized with 136 2-node elements. The material properties used to perform the stress calculations are given in section 4.2 and 4.3 for static and fatigue analysis respectively. It is reminded that for the probabilistic analysis, the property values for the UD ply of the INNWIND.EU blade design [10] were used as the mean values of each static material property making the assumption of separated UD plies of the BIAX and TRIAX fabrics. Furthermore, OPTIDAT database was used for uncertainty quantification i.e. the estimation of the CoV value of each property. For the fatigue properties, the IWES database was used to estimate both mean and CoV values. Property values were presented in Table 3 and Table 4 for static and fatigue properties respectively.
In the third process, the limit state function is determined for every analysis type (static, buckling, or fatigue). For the strength analysis under ultimate loading, the limit state function is given by: where X the vector of the basic variables, f1(X) a function expressing either the strength ratio i.e. the ratio of allowable stress over the applied one of each failure criterion or the stretch factor fs for Puck criterion, XR a random variable quantifying the model uncertainty related to the adopted failure criterion and XD a random variable expressing the model uncertainty related to the structural model used to evaluate deflections.
For the elastic stability, the limit state function is described by: where λ1 the critical eigenvalue and Xλ a random variable quantifying the model uncertainty related to the evaluation of the load factor. For the fatigue analysis the limit state function is given by where = ∑ , with ni the applied load cycles for the i-th bin of the Markov matrix derived by the rainflow counting method and Nfi the allowable number of cycles for each bin calculated by Eq. (17)- (19). In Eq. (22) Xd is a random variable quantifying the model uncertainty related to the calculation of the linear damage index D. For all the limit state functions when g(X)>0, the composite ply (or the blade section for the buckling analysis) is in the safe mode.
The fourth process in PROBUST contains probabilistic methods for reliability estimation. In the current exercise, the crude Monte Carlo (MC) method was used for strength analysis under variable amplitude loading. Sample values for every random variable were generated and repetitive simulations were performed. The failure probability (Pf) for every ply is estimated as the ratio of the number of failures to the total number of simulations. A major issue of concern in implementing the MC method is the random number generators. For correlated non-normal variables, the formulation developed in [23] was adopted.

Model uncertainty
The model uncertainties discussed in the formulation of the limit state functions (Eq. (20)- (22)) are based in the work performed within the INNWIND.EU project [26]. Six organizations participated in the benchmark on structural analysis tools for wind turbine blades. As all partners were provided with the same information regarding the blade external structure, the internal material lay-out as well as loading that should be imposed on the blade it was possible to compare the output and assess the variability of the results. The results were derived by each partner using their own tools (fully in-house developed or combined with commercially available ones), with its specific structural analysis approach (thin wall theory and finite element (FE) models using beam, shell or solid elements) and their preferable analysis type (linear or geometrical non-linear).
For the section of interest i.e. the section with the maximum chord, coefficient of variations of the respective model uncertainties XD, Xλ and Xd were estimated in [26] and are shown in Table 8. Work on the quantification of the model uncertainty related to failure criteria of composite materials, XR, has already been performed in [32] based on experimental data. Results were adopted herein and presented in Table 8.
According to JCSS [22], model uncertainty is characterized by the choice of a suitable probability distribution offering mathematical convenience to the analyst. For this reason, due to the large CoV value e.g. in case of XR, negative values may be introduced and thus quantities that should always be positive become negative ones e.g. the strength ratio f1 of Eq. (20). To circumvent the difficulty a lognormal distribution is selected for these cases.

Strength analysis under extreme loading
Having developed appropriate probabilistic models for the static and fatigue material properties, reliability analysis of the blade is performed. According to [26], the critical section of the blade was found to be close to the maximum chord section. The Puck criterion was considered to identify failure for each layer and perform the verification analysis that DNVGL-ST-0376 requires for the composite laminates i.e. fibre and inter-fibre failure. The annual failure probability was estimated by implementing Monte Carlo with 2E+06 iterations. It should be mentioned that the β index is related to the failure probability by Pf=Φ(-β), with Φ() the standardized normal distribution N(0,1). The analysed section consists of 136 2-node laminated elements with 9 plies per element.
The section is considered to be a series system of elements and in turn each element is assumed to be a series system of layers while each layer fails due to several failure modes. Assuming first ply failure for each laminated element and positively correlated behaviour of the layers in the stacking sequence, the annual β value of the element is given by the minimum annual β value for every lamination scheme, i.e. the worst case of all plies in each element, and is depicted in Figure 4. Additionally, the estimated annual β index value of the layer corresponds to the worst case among the five failure modes of the Puck criterion (namely tensile fibre failure, compressive fibre failure, mode A, B and C for the inter fibre failure mode). Because of the different loading conditions at every point in the multi-cellular blade section, the entire β index distribution, element by element, is depicted in Figure 4. Due to the sample size of 2E+06, a lower bound in the estimated probability of Pf=5E-07 (β=4.89) is implied. Also, note that the target probability of failure referenced in DNVGL-ST-0376 is 1E-04, corresponding to a β index value of 3.72. The green curve, with x marks, refers to the reliability level of the blade for which physical, statistical and measurement uncertainties for the material properties are considered. Introducing the model uncertainties of Table 8, XD and XR the blue (with triangular marks) and red (with dashes) curve are obtained respectively. It is reminded that XD corresponds to the model uncertainty related to the structural model used to evaluate deflections while XR is related to the failure criteria for composite material. Finally, the 'current' reliability level of the blade design, as explained in section 6, is depicted with the black (with circular marks) curve when appropriate safety factors given in Table 6 are considered in the analysis. The results on the complete section of the blade in Figure 4, fluctuate from low reliability levels (even failure) to expected reliability levels (i.e. β>3.72). It should be stated that the effect considering the measurement uncertainty of the material properties on the reliability level of the blade is minor for this particular example, were the material database used was quite large (each property defined by at least 25 experiments versus the traditional required 7). When model uncertainties are introduced, the β index distribution is considerably affected, especially, for the model uncertainties related to the formulation of the failure criterion of composite materials (red with dashed curve). Taking into account the safety factors as discussed earlier, the reliability level of the blade is even more reduced. Clearly, the target probability level of β ≥ 3.72 is not achieved for a large number of elements. The low β index values of Figure 4 correspond mainly to the A and B inter-fibre failure modes of Puck criterion indicating the initialization of matrix cracks in the [±45] layers. It is important to say that the failure probabilities are zero for all the elements of the blade section considering fibre failure or the catastrophic failure mode C. The outcome of this work does not imply that the blades spinning around the word should start fail (with a significant failure rate) but there is a high probability that off-axis layers in the laminates will develop matrix cracks. The development of these matrix cracks does not necessarily lead to blade collapse. However, it will certainly affect operational life under cyclic loading. That is, these matrix cracks in combination with the cycle loading of the blade can lead to more substantial failure modes (e.g. extensive matrix cracks, delamination etc.) in the laminates of the blade.
In Figure 5, probabilistic and deterministic results for the same case are set side by side. Specifically, the distribution on the blade section of the 'current' reliability level, as defined in section 6, is compared against the stretch factor fs distribution of the Puck failure criterion (worst case among the failure modes) as calculated following the DNVGL-ST-0376. fs distribution is depicted by the blue continuous curve while its values are plotted against the right hand side axis of the graph. In the same graph the critical stretch factor value of 1.0 is shown. fs values less than unity indicate failure for the specific element. It is clearly seen that the stretch factor values for some elements e.g. the element #33 (fs=0.93) are only slightly below the limit and thus the specific design could be considered as an adequate one according to DNVGL-DT-0376. This, however, means that the specific design should result in an annual failure probability approximately equal to 1E-04. The probabilistic approach of this work, however, indicates a much larger annual failure probability. The result can be attributed mainly to the considerable effect of the failure criterion model uncertainties in the analysis on top of the safety factors required by the standard, as seen in Figure 4.

Buckling analysis
Similarly to the analysis presented for the strength of the blade under extreme loading in this section a buckling analysis of the specific blade section is discussed. Results are depicted in Table 9. Specifically, the first row of Table 9 indicates the annual β index value of the blade section when measurement uncertainty for the stiffness properties is taken into account. The β values were estimated implementing the response surface method as described in [8] combined with the Monte Carlo. The second row corresponds to the reliability level when the model uncertainty of Table 8 Xλ is also considered. The third row of Table 9 corresponds to the 'current' reliability level, as defined in section 6, of the blade section considering the appropriate safety factors of Table 6. As for the static case a drastic change on the estimated β values is observed when model uncertainties and appropriate safety factors are considered. Still, note that the acceptable β index value for an annual failure probability of 1E-04 is 3.72. This result was, nevertheless, anticipated because the deterministic design following DNVGL-ST-0376 indicated a load factor equal to 0.52 as shown in the last row of Table 9. A load factor less than unity indicates buckling of the section. It is clear that the blade should be reinforced at this area for example by inserting thicker sandwich panels.

Strength analysis under variable amplitude loading
Regarding strength under variable amplitude loading, DNVGL-ST-0376 defines two approaches to take into account the mean stress effect; namely the shifted Goodman diagram and the linear piecewise constant life diagram. The β index distribution for the blade section at the maximum chord was estimated for both cases and is presented in Figure 6. In order to make this comparison and have failure probability values other than 1E-06, i.e. the limit of the calculation procedure, the number of operational years of the wind turbine was increased more than the specified design life of 20-30 years. It is obvious that the β index distribution is significantly affected by the CLD implementation. For the shifted Goodman CLD, the critical parts of the section are the spar cup and the leading edge panel in the suction and pressure side of the blade. For Linear Piecewise CLD, critical parts of the blade section are the trailing panel and the cup in the pressure side but not in the suction side of the blade. Both approaches indicate almost the same reliability level but with totally different β-index distribution along the section.
Finally, the effect of measurement uncertainty for the fatigue properties directly on the reliability level of the blade section is presented in Figure 7 using the shifted Goodman CLD and taking also into account the appropriate safety factors of Table 6 and the model uncertainty  on the linear damage index value of Table 8. The estimated β index values correspond to the most critical element of the section, that is, the element #69 at the spar cup in the suction side of the blade.
In Figure 7 the β values are plotted against the years of operation of the wind turbine The labels "Mat(P)", "Mat(P+S)", "Mat(P+S+M)" correspond to different formulations for the material fatigue properties, namely: physical only, physical and statistical and physical, statistical and measurement uncertainty, respectively. The effect of considering measurement uncertainty of the static and fatigue material properties on the blade reliability is quite obvious. When measurement uncertainty is ignored, the failure probability of the blade section in terms of fatigue analysis is zero as indicated by the curve with circles (black) or the line with crosses (green) in Figure 7. The 'current' reliability level for T=25 yr, as defined in section 6, is given by the bold marker and it has a β index value equal to 3.6. The annual β index value of the blade section under variable amplitude loading is much higher than the one required from the standard i.e. β=3.72.
In the same graph the linear damage index as calculated following DNVGL-ST-0376 is also presented. The D index values are given on the right hand side axis of the graph. It is noted that, for an acceptable design in deterministic terms the D index should be less than unity. Quite small index values were observed. For the presented design assumptions, the fatigue was not a limiting failure mode for the blade. Yet, as seen for the estimations on the 25 years of operation, the deterministic approach clearly overestimates the fatigue capacity of the blade, leading to non-conservative design values, when compared to the probabilistic one.

CONCLUSIONS
The results presented above are not general for all wind turbine blades. They are based on a specific structural design, as well as connected with the available material databases and load conditions. For this specific design, the highest annual probability of failure was observed for the elastic stability case (buckling), while the annual probability of failure against fatigue was found to be negligible. This being said, some trends can be clearly identified.
Uncertainty analysis was performed according to the principles of metrology using the OPTIDAT database and the database in [6]. Regarding static properties, variability on the sample mean was found to fluctuate between 0.8-10%. For the standard deviation, the CoV values were much higher at 19-22%. Measurement uncertainty affects mainly the tails of the distribution of the material property while, depending on the CoV value of the sample mean, it could have a more pronounced effect on these. The new standard DNVGL-ST-0376 captures cases with higher variability, yet is quite conservative for cases with low fluctuations.
For the fatigue cases, variability on the intercept parameter of the power law model, Eq. (1), was found to fluctuate between 1.29-2.36% depending on the R-ratio that CA fatigue tests were performed. For the slope parameter, the CoV value was much higher at 12.5-66.2%. The greater CoV values i.e. CoV>20%, are mainly related to the data derived by CA fatigue tests in the compression-compression region. This is not appropriately captured by the blade design standard, as shown in Figure 2.
For the blade, when considering all sectional stress resultants and the multi-axial state of stress of the laminate and the layers, the reliability of the blade section under extreme loading is not driven by the unidirectional layers on the lamination sequence, but rather from the offaxis layers. In contrast to that, for the fatigue case, where only the axial stresses (the stresses along the fibre direction) are taken into account, the reliability is driven by the caps of the blade (and the unidirectional layers). The uncertainty relevant to the failure criteria models for composite materials is still quite high, driving the results for the extreme (static) load cases. It is recommended that this uncertainty is reduced through verification of appropriate failure criteria by experimental evidence. This significant uncertainty relevant to the failure criteria, is already known and commented in the literature within the World Failure Exercise, among others. This could also be the background for the conservative approach of the new blade design standard in the static case. Nevertheless, a compensation for advanced probabilistic methods incorporating these types of uncertainty should be foreseen in the standard. Currently this is not the case.
The fatigue formulation for multi-axial state of stress for composite materials is in its infancy relevant to variable amplitude fatigue loading. This is more pronounced in the wind energy sector, where also the sectional stress resultants coming from various load cases show a strong variability making even more difficult concurrent stress counting on the constituent layer of the blade section. Work is recommended in this direction, although this falls outside the scope of the development of probabilistic methods for structural reliability estimation.
Regarding the effect of the CLD implementation, each formulation resulted in a different β index distribution indicating also differences in the locations of the critical parts of the blade section. Uncertainty related to CLD selection affects greatly the reliability estimations.
Nevertheless, opposite to the static case, the deterministic approach of the design standard leads to non-conservative results for the fatigue case. Measurement uncertainty in experiments performed in fatigue was proven quite important, while up to now it has been neglected. This might be a reason for the large deviation and uncertainty of the fatigue failure criteria employed in composite structural design, such as wind turbine blades.
In light of the above findings, in combination with the structural model uncertainties obtained in [26], it is also recommended to revisit the partial safety factors used in the blade design through a dedicated safety factor re-calibration.