Periodic stability analysis of wind turbines operating in turbulent wind conditions

In this work, a method for the stability analysis of wind turbines is described. A system identification technique, formulated for handling stochastic disturbances, is used to identify a periodic reduced order model from suitable recorded time histories of the system. Afterwards, such reduced model is analyzed according to Floquet theory. The formulation is model-independent, in the sense that it does not require knowledge of the equations of motion of the periodic system being analyzed, and it is applicable to an arbitrary number of blades and to any configuration of the machine. 5 In addition, as wind turbulence can be viewed as a stochastic disturbance, the method is also applicable to real wind turbines operating in the field. The characteristics of the new method are verified first with a simplified analytical model, and then using a high-fidelity multibody model of a multi-MW wind turbine. Results are compared with those obtained by the well known operational modal analysis approach. 10

of the Coleman transformation is small in general, nor that this approach amounts to some consistent and bounded approximation of a rigorous Floquet analysis. Given the widespread use of the Coleman transformation, and its general excellent behavior, such a proof remains a goal very worth pursuing but, as of today, yet unattained.
-Second, the Coleman transformation unfortunately exists only for a number of blades greater or equal than three. Although this is the most common wind turbine configuration nowadays, a revival of the two-bladed concept is possible.

5
-Third, codes implementing the Coleman transformation require access to the linearized equations of motion of the system. As a consequence, any addition to a simulation code has an impact on the associated stability analysis tool, resulting in extra software maintenance work.
Other possible approaches to the stability analysis of rotors have been formulated in the frequency domain. For example, the estimation of power spectra along with modal frequencies and damping ratios of an operating wind turbine has been addressed 10 by Avendaño-Valencia and Fassois (2014). That paper considered several parametric and non-parametric methods and their application to experimental data, including the periodic autoregressive (PAR) model. In addition, periodic autoregressive moving average (PARMA) models have been considered by Avendaño-Valencia and Fassois (2013). Two subspace algorithms for periodic systems have been presented by Skjoldan and Bauchau (2011) and Mevel et al. (2014), one being used for numerically generated time series, and the other for experimentally measured ones. 15 The operational modal analysis (OMA) has been extended to the periodic case (Allen et al., 2011b), by using the concept of harmonic transfer function (HTF). In that paper, the simple peak-picking method was used for extracting relevant properties from the spectra, while more specialized fitting algorithms were proposed by Allen et al. (2011a). Subsequent applications and developments can be found in Allen (2014, 2012). Although the method is general, the estimation of the quantities of interest for a stability analysis from noisy spectra remains a somewhat delicate operation, as it will be shown later on in the 20 following pages.
In the authors' opinion, there are two desirable goals in the stability analysis of wind turbines that still need further investigation in order to be fully attained: -First, one would like to account with complete rigor for the periodicity of such systems, without introducing approximations of unknown effects.

25
-Second, one would like to formulate the analysis so that it is system-independent. System independence is here intended to mean that a method can be applied to wind turbine models of arbitrary complexity and topology (e.g. any number of blades, horizontal or vertical axis, etc.), and also to real wind turbines operating in the field.
To answer these needs, Bottasso and Cacciola (2015) proposed a periodic stability analysis formulated in terms of inputoutput discrete-time responses. Such time histories could come from "virtual" experiments performed on a given model, from 30 simplified ones to the more advanced contemporary comprehensive multibody-based aero-hydro-servo-elastic models. Using this approach, a reduced periodic auto-regressive with exogenous input (PARX) model is first identified from a recorded response of the system, and then used for conducting a stability analysis according to Floquet theory. On the practical side, this implies that the analysis respects the periodic nature of the problem. Furthermore, one can easily replace the model with a different one, without having to modify or adjust in any way the stability analysis procedure.
Although this approach attains the two goals outlined above, one of its limits is that it can not be used with measurements obtained on a real wind turbine operating in the field, since the effects of wind turbulence are not considered within the PARX model structure. To address this issue, the same approach was extended to account for the presence of turbulence (Bottasso 5 et al., 2014). Using this new technique, one first identifies a periodic autoregressive moving average with exogenous input (PARMAX) model, whose stability is then analyzed according to Floquet. That paper showed only one example related to the first blade edge-wise mode of a wind turbine rotor. Goal of the present paper is to expand and formulate in detail the PARMAXbased method originally proposed by Bottasso et al. (2014). A second goal of this paper is to compare the PARMAX method with the periodic operational modal analysis (POMA) (see Allen et al., 2011a), which is taken here to represent the accepted 10 state-of-the-art for the stability analysis of wind turbines operating in turbulent wind conditions. The article is organized according to the following plan. The problem of the identification of PARMAX models is addressed in Sect. 2. Here, a newly developed algorithm that has its basis on the prediction error method (PEM) is formulated, with particular emphasis on the guaranteed stability of the PARMAX predictor. Section 3 is devoted to POMA theory. After reviewing the concept of HTFs, the treatment proceeds by discussing the method and its use for conducting periodic stability analyses. 15 As the authors are not aware of a reference collecting together all useful background information on Floquet theory and the signal analysis tools needed for POMA, this material is synthetically reviewed in Appendix A, to ease reading. The accuracy of the PARX and POMA identification techniques is then compared against an exact reference in Sect. 4. To this purpose, first a nonlinear wind turbine analytical model is developed. Then, the stability of its linearized version is studied according to Floquet theory, providing a reference ground truth used for comparing PARX and POMA. The equations of such analytical In accordance with Bottasso and Cacciola (2015), the deterministic behavior of a wind turbine measured output z can be modeled with a PARX sequence as A(q; k)z(k) = B(q; k)u t (k), where k is the time index and q the back-shift operator, such that z(k)q −i = z(k − i). The autoregressive and exogenous parts are defined respectively by polynomials A(q; k) and B(q; k) as both being characterized by periodic coefficients a i (k) = a i (k + K) and b j (k) = b j (k + K), where N a and N b indicate the order of the AR-and X-part, respectively, while K is the period of the system. Finally, u t is the input, assumed here to be the turbulent wind.

10
The stochastic nature of the turbulent wind field violates the assumption of a deterministic and fully measurable input u t . To account for this, the actual wind is viewed as a sum of two distinct contributions: a mean wind u(k) and a turbulence-induced perturbation δu t (k). As the spectrum of the atmospheric turbulence is far from being constant, δu t (k) is modeled by means of a shape filter F(q; k) such that 15 where e(k) is a zero-mean, white and Gaussian noise, with periodic variance σ(k) 2 .
Inserting Eq. (3) into (1), the following is derived A(q; k)z(k) = B(q; k)u(k) + G(q; k)e(k), where C(q; k) = B(q; k)F(q; k). Equation (4) is a PARMAX model whose MA-part is represented by polynomial G(q; k), defined as where g w (k) = g w (k + K) are the MA periodic coefficients and N g the MA order. The overall order of the system is defined as n = max(N a , N b , N g ). The resulting PARMAX sequence is then It should be remarked that the present approach does not consider the effects of nonlinearities nor of rotor speed variations 25 induced by turbulence. The former potential problem can be checked a posteriori by looking at the matching between predicted and measured quantities. The latter can be partially solved by averaging the rotor speed over the analyzed time window.
Typically, because of the large inertia of wind turbine rotors, angular speed variations are not expected to be highly significant, especially within the short time windows required by the proposed approach.

State space representation of PARMAX sequences
In order to perform a stability analysis according to Floquet theory (cf. Bottasso and Cacciola (2015) and the review reported 5 in Appendix A), it is necessary to realize the PARMAX sequence (6) into an equivalent state space representation. To this end, consider a linear discrete-time system with time-varying coefficients in observable canonical form where x(k) = (x 1 (k), . . . , x n (k)) T , while the system matrices are given by Including the presence of the MA-part, the input-output sequence of system (7) can be derived as Comparing Eq. (6) with Eq. (9), the following equivalence relations are obtained 20 which readily give the state space system matrices. Once these are known, stability is assessed according to Floquet theory as described in Appendix A.

Identification through the prediction error method
In the present context, a single-input single-output (SISO) PARMAX model must be identified from a sequence of N measurements. Among the plethora of existing estimation methods, which may range from time to frequency domain or from optimization-based to subspace algorithms, the PEM (Bittanti et al., 1994) is chosen here. This method has been frequently used for rotating systems, such as rotorcraft vehicles and wind turbines. For example, the periodic equation-error method was 5 used for identifying a reduced order model of a helicopter rotor by Bertogalli et al. (1999), whereas Bottasso and Cacciola (2015) proposed a periodic output-error method for the identification of reduced wind turbine models.
The estimation problem, formalized according to the PEM, is the one of finding the periodic coefficients a i (k), b j (k) and g w (k) that minimize cost function J defined as the mean value of the square of the prediction error, i.e.

10
Here ε(k) = z(k) −ẑ(k|k − 1) is the prediction error at time instant k, beingẑ(k|k − 1) (hereafter more concisely written aŝ z(k)) the optimal one-step-ahead prediction of z(k) based on knowledge of all data until time step k − 1. According to Bittanti and De Nicolao (1993) and Ljung (1999), the optimal one-step-ahead predictor of process (6) iŝ As previously argued, the presence of the MA part in the PARMAX model allows for a more adequate characterization of 15 the process noise term, at the cost of a more complex estimation procedure. In fact, the optimal predictor of the PARMAX process expressed by Eq. (12) is nonlinear in the parameters, as anyẑ(k) is a function of its previous valuesẑ(k − w), which in turn depend on the parameters. Consequently, the minimization of cost function (11) involves an iterative optimization. If the MA part in Eqs. (6,12) is neglected, a PARX sequence is obtained and the estimation problem reduces to the so-called equation error approach (Bottasso and Cacciola, 2015;Bottasso et al., 2014). 20 Moreover, it is easy to verify that predictor (12) is by itself a PARX dynamic system, in which the autoregressive part is described by coefficients −g w (k), whereas coefficients a j (k) + g j (k) and b j (k) define two X parts with inputs z(k) and u(k), respectively. This fact is not surprising, since it often happens that the poles of the predictor coincide with the zeros of the system to be predicted. As a consequence, it may happen that, during the optimization, coefficients g w define an unstable predictor, jeopardizing the entire identification process (see Bittanti and De Nicolao, 1993).

25
In the literature there are basically two methods to enforce the stability of the MA part. The first is a heuristic approach in which the coefficients g w (k) are perturbed (for example halved) repeatedly until the achievement of a stable predictor.
This method actually corresponds to a re-initialization of the parameters with unpredictable effects on the convergence of the estimation. The second approach is based on the computation of a new predictor, with different coefficients g w but the same autocorrelation of the unstable one. For the time-invariant case, this new canonical model can be obtained using Bauer's 30 algorithm (Sayed and Kailath, 2001), whereas for the periodic case by solving a suitable periodic Riccati equation (Bittanti and De Nicolao, 1993) or through the multivariate Rissanen factorization (Bittanti et al., 1991;Rissanen, 1973).
In this work, an alternative and original method is proposed. The stability of the predictor is enforced by a nonlinear constraint within the estimation process, and the resulting constrained optimization is performed by an interior-point algorithm (cf. Byrd et al., 2000Byrd et al., , 1999Waltz et al., 2006). The estimation problem is then reformulated as where p is the vector of the unknown coefficients and P(p) are the characteristic multipliers of the PARMAX predictor.
The characteristic multipliers that constrain the estimation problem can be computed from the autoregressive part of Eq. (12), i.e.ŷ(k) = w −g w (k)ŷ(k−w), which can be realized into state space form according to Eqs. (7-10), leading to the following dynamic matrix 10 The periodic coefficients a i (k), b j (k) and g w (k) are approximated by using truncated Fourier expansions, i.e.
b c jm cos(mψ(k))+b s jm sin(mψ(k)) , g c wr cos (rψ(k))+g s wr sin (rψ(k)) , where ψ(k) is the rotor azimuth. The unknown amplitudes of such expansions are collected in the vector of parameters p Due to the nonlinear behavior of the predictor, the possible presence of multiple local minima has to be taken into account.
A suitable starting point for the nonlinear problem can be selected by fitting the recorded data with simpler models such 20 as ARMAX or PARX (Bittanti et al., 1991), or by using a recursive extended least-squares algorithm (Avendaño-Valencia and Fassois, 2013; Spiridonakos and Fassois, 2009). In the present work, convergence to the global minimum is ensured by performing several optimization trials from a randomly chosen set of initial conditions.

Theory of periodic operational modal analysis
The OMA is an output-only system identification technique, which has been widely used to conduct modal analyses of different mechanical systems. Recently, a special attention has been devoted in the literature to the application of OMA in the field of wind energy (Carne and James, 2010), and to the related underlying hypotheses (Chauhan et al., 2009;Tcherniak et al., 2010).
An output-only technique specifically tailored for time periodic systems was developed by Allen et al. (2011b). This technique, 5 named periodic OMA (POMA), exploits the particular behavior of an LTP system in the frequency domain, as described by the HTF (see Sect. A2 for details). In the present paper, POMA will be briefly reviewed and then compared to the PARMAX-based stability analysis proposed here.
Consider a strictly proper periodic system and the exponentially modulated periodic (EMP) expansions of its input and output, noted respectively U and Y, as described in Sect. A2. The input-output behavior of the system can be analyzed through 10 the HTF G as with s ∈ C, and G(s) defined according to Eq. (A44). Projecting (17) onto the imaginary axis, each element of the EMP expansion of Y and U can be computed as the Fourier transform of frequency shifted copies of y(t) and u(t) as As reported in Wereley (1991) and briefly reviewed in Allen et al. (2011b), the input-output behavior in the frequency domain can be expressed as Accordingly, the harmonic frequency response function (HFRF) G(ω) is given by where C j,w and C j,w are defined in Eq. (A45) and Eq. (A46) of Sect. A2.

25
The power spectrum of the output, noted S Y Y (ω), can be written in terms of the HFRF G(ω) and the power spectrum of the input S U U (ω) as where (·) H denotes the complex-conjugate transpose. Inserting now Eq. (21) into Eq. (22), the following expression is derived where W j,w,r,t = B j,r S U U B H w,t . Equation (23) can be simplified first by considering a flat expanded input power spectrum W j,r (ω) = B j,r S U U B H j,r , at least in the band of interest of a specific mode, and secondly by assuming that all modes of the system are "suitably separated".
The first requirement was analyzed extensively for wind turbine problems in Tcherniak et al. (2010). There the authors 10 pointed out that the extended input spectrum could be significantly colored, a problem that requires particular care with simplified output-only methods. The second requirement deserves special attention as well. In fact, not only is the separation of the principal harmonics of two modes required, but it is also necessary that all super-harmonics with significant participation are well-separated. For rotary wing systems, this requirement has to be considered carefully especially when looking at the whirling modes, as the principal harmonics of backward and forward modes are typically separated by about 2Ω. This typically 15 creates a crisscrossing of modes in the frequency-rotor speed plane, leading to frequent frequency encounters.
If such conditions are verified, the extended input spectrum W loses its dependency on ω, and the contribution of mode η p + ıqΩ on mode η j + ıwΩ can be neglected when p = j and q = w. Hence, Eq. (23) is simplified to From Eq. (24) one can notice that the peak related to any super-hamonic of a given mode can be viewed as the peak of a 20 linear time-invariant mode. Accordingly, one is allowed to use a standard LTI frequency domain identification technique (e.g. peak-picking, curve-fitting) to compute frequencies, damping factors and modal shapes from the measured spectra.
Moreover, neglecting again the contribution of overlapping modes, one can also estimate the participation by evaluating the power spectra at the peak frequency, since

25
Expressing the product C j,w C H j,w one gets being (·) * the complex conjugate. From Eq. (26), one could envision several criteria for extracting the participation factors for each harmonic belonging to the j-th mode. The simplest one is to compute the central column of the HTF and to pick the amplitudes of the spectra at the frequency of interest. The participation factors are then extracted according to Eq. (A27), 5 reported in Sect. A1, as One can also perform multiple estimations of the participation factors by looking again at the central column of S Y Y . In fact, from Eq. (26), it appears that the amplitudes picked from the th column at frequency ω j + wΩ are equivalent to those picked from the central column at ω j + (w + )Ω. This means also that computing the central column could be sufficient for having an 10 estimation of frequencies, damping and participation factors, as already noticed in Shifei and Allen (2012).
The POMA technique can be then summarized as follows: -Compute the Fourier transforms of the frequency shifted copies of the recorded output y(t), y k (ω) = FFT y(t)e −ıkΩt , and collect them in vector Y (ω) = (. . . , y k (ω), . . .) T .
-Compute the autospectrum S Y Y (ω) using a standard frequency domain analysis method; in the present paper the method 15 of Welch was employed for this purpose.
-Extract from each peak present in S Y Y (ω) the related natural frequency and damping factors using any standard LTI frequency domain estimation tool (Allen and Ginsberg, 2006). In this paper the straightforward peak-picking method was used, as also done by Allen et al. (2011b).
-Reconstruct the Fourier coefficients c jn , and in turn the participation factors, by evaluating the spectrum in correspon-20 dence to each peak.
It is possible to restrict the analysis to the right-half plane just by noticing that Equation (28) is particularly useful for identifying the Fourier coefficients from the peaks of the "reflected super-harmonics", since according to Eq. (28) one can demonstrate that

Application of periodic operational modal analysis to the Mathieu oscillator
As the actual use of POMA and the correct interpretation of all peaks is not a straightforward exercise in general, a simple 5 Mathieu oscillator is analyzed here in preparation to the application of this method to the wind turbine problems studied later on. The dynamics of a Mathieu oscillator is governed by the following equations The parameters in Eq. (30) were set, following Allen et al. (2011b), as ω 2 0 = 1, ω 2 1 = 0.4, ξ = 0.04 and Ω = 0.8. The system 10 was numerically integrated from x(0) = (1000, 0) T , and studied by means of POMA. The results were then compared with those obtained by the full-Floquet theory described in Sect. A1. Figure 1 shows the power spectra of the central column of S Y Y , y k (ω)y H 0 (ω) for k = −4, . . . , 4. The fundamental peak (i.e. the highest one) is found on the 0-shift curve at 0.16 Hz and corresponds to the amplitude c j 0 c j H 0 . At such frequency, all curves show a prominent peak, from which one may easily compute also the damping factors using for example the standard 15 half power bandwidth method. The participation factors are then extracted by looking at the amplitudes of the power spectra using Eq. (27).
Starting from this peak and moving to the right, the subsequent higher peaks are found on the negative-shift curves, first in the −1-shift one at 0.28 Hz, then in the −2-shift one at 0.41 Hz, etc. The opposite happens when moving to the left. Peaks located at negative frequencies appear as reflected in the positive frequency range, but with opposite shifts. This is clear if one 20 looks at the peak located at -0.10 Hz, which has the −2-shift curve as the one with the highest amplitude, whereas the reflected peak at 0.10 Hz is associated with the 2-shift curve. This complex behavior is easily explained by means of Eq. (28), which also states that the information in the negative frequency range can be reconstructed by looking at the curve with the opposite shift in the positive frequency plane.
Frequencies and damping factors computed from such spectra using the peak-picking method are reported in Table 1. The 25 same table also displays the results obtained from the full-Floquet analysis of the system. The comparison shows good accuracy, especially for frequencies and damping factors of the first highest super-harmonics.
The output-specific participation factors are displayed in Table 2. Multiple estimates have been computed from each spectrum peak in the positive frequency plane. The last column shows also the analytical results. As expected, in general superharmonics with lower participation factors are associated with higher estimation errors. 30 Next, a simplified wind turbine model is used for comparing the results obtained with the PARX and POMA approaches. This is useful because it gives a way of comparing the basic performance of the two methods with respect to a known exact ground truth in the ideal case of null disturbances. Later on in this work, the two methods will be compared for the case of a higher fidelity wind turbine model operating in turbulent wind conditions. As no exact solution is known in that case, the preliminary 5 investigation of this section serves the purpose of clarifying whether significant differences exists between the two approaches even at this more fundamental level. Indeed, it will be shown here that some of the underlying hypothesis of POMA are not always fulfilled, and this leads occasionally to some imprecisions in the estimates of modal quantities of interest.
The analytical model is derived in detail in Appendix B, which also gives a schematic sketch of the system in Fig. 12. The model considers the coupled motion of tower and blades subjected to aerodynamic and gravitational forces. The fore-aft and 10 side-side flexibility of the tower is rendered by two equivalent linear springs, whereas each blade is represented as a rigid body connected to the hub through two coincident linear torsional springs, allowing respectively the blade flap-and edge-wise rotations. The characteristics of each element in the model are chosen so as to match the first tower fore-aft and side-side modes and the first blade flap-wise and edge-wise modes in vacuo of a reference 6 MW wind turbine, as computed using a high-fidelity multibody model. The aerodynamic formulation is inspired by the treatment of Eggleston and Stoddard (1987), 15 in which the aerodynamic forces and moments at the blade hinges are computed assuming linear aerodynamics, small flap and lag angles, uniform inflow over the rotor disk and constant rotor speed. The aerodynamic forces induced by tower motion, not present in the treatment of Eggleston and Stoddard (1987), are additionally considered in this paper. The model represents the complete lower spectrum of a wind turbine, including the first side-side and fore-aft tower modes, the first in-plane and out-of-plane blade modes as well as their related whirling modes. 20 After having collected all degrees of freedom in vector ξ = (β 1 , . . . , β B , ζ 1 , . . . , ζ B , y H , z H ) T , being B the number of blades, and the inputs in vector ν = (θ p1 , . . . , θ p B ) T , being θ k the pitch angle of the kth blade, the resulting nonlinear second order implicit system writes System (31) can be integrated in time using any suitable numerical scheme, starting from a consistent set of initial condi-25 tion. This was done for generating the time histories used for PARX and POMA, paying attention not to excite the system nonlinearities, as the reference solution is based on the Floquet analysis of the linearized problem.
Since any mechanical system is linear inξ, one may compute the mass matrix M (ξ,ξ, t) and rewrite the system as M (ξ,ξ, t)ξ = g(ξ,ξ, ν, t). System (31), if asymptotically stable, converges to a periodic trajectoryξ(t) when subjected to a periodic inputν(t). In such a regime, the linearized periodic equations of motion write where the new stateξ(t) and inputν(t) are defined aŝ and the periodic mass, damping, stiffness and input matrices are defined Notice that M (t) is equal to M (ξ,ξ, t) evaluated on the periodic trajectoryξ. These linearized equations of motion about a periodic orbit were then used for developing the analysis according to Floquet, yielding the ground truth solution.
4.1 Stability analysis of a wind turbine analytical model 5 The parameters of the wind turbine analytical model were defined according to Table 3, which loosely represent a conceptual 6 MW wind turbine. The stability of the model is studied in a uniform axial wind of 9 m/s for a collective pitch angle of −0.54°, corresponding to operation towards the end of the partial load region. The linearized periodic system was first studied using Floquet theory (see Appendix A) in order to get the exact natural frequencies, damping and output-specific participation factors. Next, the model was used for generating all outputs needed 10 for performing the PAR(MA)X and POMA analyses by integrating the system forward in time starting from suitable initial non-zero conditions, chosen in order to excite the modes of interest. In this exercise, the wind was considered as stationary, so that the PARMAX identification reduces to the simpler PARX one as the MA-part is not necessary.
Both PARX and POMA estimates were compared with the full-Floquet results in term of relative errors for frequencies and damping factors, and absolute errors for participation factors. Relative errors are defined as v E /v R − 1, while absolute errors 15 as where v is a specific modal parameter and the subscripts E and R refer, respectively, to an estimated and a real (exact) quantity.

Identification of the blade edge-wise mode
The blade edge-wise mode was excited by imposing the initial edge-wise angles of all blades equal to a unique non-zero value, whilst all other states were set to zero at the initial time. This way the blade in-plane mode was excited, while avoiding the 20 onset of the whirling modes.
Considering first the POMA approach, the harmonic power spectrum for the second blade edge-wise angle, ζ 2 , was computed with frequency shifts from −2Ω to +2Ω. The results obtained this way are reported in Fig. 2.
Clearly, the 0-shift PSD shows a prominent peak at ω E = 0.86 Hz related to the blade in-plane mode, from which one may easily extract the frequency and damping factor of the principal harmonic. The peak-picking method could in principle be 25 applied to any of the peaks displayed in the figure; however, one may observe that most of the peaks are of a low amplitude and often barely noticeable from the side-band of the principal harmonic. For example, the super-harmonic at 0.67 Hz, even if visible within the 0-shift curve, has not enough energy to allow one to estimate its modal quantities to any reasonable accuracy.
Therefore, it was preferred to compute frequency and damping factors only by looking at the highest peaks: the frequency and damping factor of the super-harmonic at ω E + Ω were extracted from the peak at 1.05 Hz of the -1-shift curve, while those of the super-harmonic at ω E + 2Ω from the peak at 1.24 Hz of -2-shift curve, and similarly for the other super-harmonics. For the same reason, participation factors were obtained only by looking at the PSD amplitude at ω E . In fact, at this frequency all curves show peaks that are prominent and distinct enough to compute the participation factors according to Eq. (27).
Next, the PARX analysis was considered. As long as only the blade in-plane mode is significantly excited, as indicated from 5 the 0-shift curve in Fig. 2, the order of the AR-part may be set as N a = 2. A first order X-part (N b = 1) was considered as the inputs (wind speed and pitch angle) are constant in this case. Finally, the number of harmonics for the Fourier series expansion of both the AR-an X-parts, N Fa and N F b , were both set equal to 1. The matching between predicted and simulated output, not reported here for the sake of brevity, showed excellent correlation, proof of the fact that the identified model captures very well the dynamics of interest.
10 Table 4 reports the Floquet modal parameters, assumed as ground truth, as well as the errors obtained by the two methods considered here.
Looking at the results, it appears that both the PARX and POMA methods are able to capture the relevant dynamics related to the principal harmonics, as frequencies, damping and participation factors are of good quality. In particular, damping and participation factors are slightly better estimated by PARX. 15 The estimation of the super-harmonic modal parameters deserves a special mention. The PARX method is able to provide a good matching for all modal parameters of all harmonics: frequencies and participation factors have negligible errors, whereas damping factors show an error lower than 1%. On the other hand, the error of the POMA super-harmonic estimates is typically quite large especially for the damping factors, even though the principal harmonic is well captured. This fact has mainly two possible explanations. First, the hypothesis of well separated modes is here not fully satisfied, as the 20 side band of the tower principal harmonic affects all super-harmonic peaks. The lower the rotor speed, the more this effect is pronounced, as the frequency separations among super-harmonics coincide with multiples of the rotor frequency. Second, but more importantly, according to the dynamics of a periodic system all harmonics belonging to a specific mode descend from a sole characteristic multiplier. Therefore their frequencies and damping factors are strictly connected to each other. This relation is totally ignored by POMA (cf. Allen et al., 2011b), as it considers each peak in the frequency response as a stand-alone mode.

Identification of other low-damped modes
The tower side-side and blade in-plane whirling modes were excited by imposing different initial conditions for each blade edge-wise angle and a suitable lateral displacement of the tower. Figure 3 shows the harmonic power spectral density (HPSD) for the tower side-side displacement y H , with frequency shifts from −2Ω to +2Ω. Here again, the zero-shift curve shows three distinct peaks: the tower side-side mode, the backward and 30 forward in-plane whirling modes, respectively at 0.34 Hz, 0.68 Hz and 1.1 Hz. Accordingly, the PARX complexity was set as N a = 6, N b = 1, N Fa = 1 and N F b = 1. As for the previous case, the matching between predicted and simulated output, not reported here, is excellent.
Comparisons among the exact and identified modal parameters are displayed in Tables 5 through 7.
Figure 3 clearly shows that a good mode separation is here not fully achieved, as whirling super-harmonics interact with each other. This is not due to the specific wind turbine or condition considered here, as in fact any rotating blade system will always have the principal harmonics of its whirling modes separated by about 2Ω. In addition, it also appears that the second  According to the PARMAX-based stability analysis, the system should be perturbed so as to induce a significant response of one or more modes of interest. Among the many possible ways of exciting a specific wind turbine mode, as for example the use of pitch and torque actuators  or of eccentrical masses (Thomsen et al., 2000), impulsive forces were used in this work. Such forces could be realized in practice by pyrotechnic exciters. The rotor angular speed is averaged over 30 the length of the recorded history and used to compute the system period. Afterwards, the signal is re-sampled in order to have an integer number of steps within a period.
The selection of the model complexity deserves special care. As the order of the AR-part N a is strictly related to the number of system modes, it can be estimated by looking at the number of principal-harmonic peaks present in the output PSD. This heuristic approach for the problem at hand turned out to be simple and effective and was preferred to more sophisticated criteria (Skjoldan and Bauchau, 2011;Avendaño-Valencia and Fassois, 2014). As described in Sect. 2.1, the input wind speed was considered as the sum of two contributions, a constant deterministic part and a turbulence-induced one. As long as the 5 deterministic input is considered to be constant, one is allowed only to estimate an X-part with order N b = 1. The MA-part order (noted N g ) as well as the number of harmonics used to model the periodicity of the coefficients (noted N Fa , N F b and N Fg ) were set with a trial an error approach, until the achievement of satisfactory results.
After having performed the estimation for different wind conditions and therefore at different rotor speeds, the results of the analyses in term of frequency, damping and participation factors were fitted using low-order polynomials, computed by means 10 of the robust bi-square algorithm (Kutner et al., 2005). The fitting process was applied only to the frequency and damping of the principal harmonic, indicated with the subscript (·) 0 . The corresponding characteristic exponent was then computed as The super-harmonics were finally obtained by means of Eq. (A17). On the other hand, the participation factors of all superharmonics were fitted with the same bi-square algorithm. This setting allows for the modeling of three periodic modes.
The result of an identification executed at the rated rotor speed is shown in Fig. 4. The excellent superposition of the curves 20 indicates a reduced order PARMAX model of very good quality, capable of modeling the signal behavior despite the small nonlinearities and rotor speed variations characterizing the system that generated the data.
To draw the Campbell diagram, eight different identifications were made in order to cover the entire range of angular speeds of the machine. The results are shown in Fig. 5, where red dots indicate each specific identification, whereas lines refer to their quadratic fits. The gray bands are the 2σ non-simultaneous functional prediction bounds, and measure the confidence level of 25 the fitting curves. From the gray bands one can infer that each frequency and damping factor identified at a specific rotor speed is coherent with the others, as all the estimates define a clear trend. On the other hand, a significant but acceptable uncertainty still characterizes the participation factors.
Similar analyses were conducted by Bottasso et al. (2014), where a different turbulence intensity (IEC level "B", instead of a uniform 5% over the whole wind speed range) was used, caeteris paribus. As the Campbell diagram is similar in both works, 30 one may conclude that the PARMAX-based analysis does not appear to be significantly influenced by turbulence level.
Much longer portions of the time histories analyzed with PARMAX were then processed with the POMA method. In Fig. 6, the HPSD obtained for a wind field with a 3 m/s average speed is shown (notice the similarities with Fig. 2). For this case the turbulence intensity was quite low, and the HPSD lines present well defined peaks. However it was found that, for increasing wind speed, while the n = 0 lines remain well defined, the quality of the peaks associated with the super-harmonics progressively degrades, making the estimation of damping (and, in some cases, also of frequency) increasingly more difficult.
The Campbell diagram obtained from POMA is displayed in Fig. 7. Comparing this figure with the PARMAX plot shows that frequencies are well identified, but the high dispersion of damping factors masks the expected trend. Several differences 5 can also be noticed between the plots with respect to the participation factors. While both approaches indicate that the principal harmonic is the most important in the response, they however detect a markedly different behavior as a function of rotor speed.
In addition, POMA overestimates the participation factors of the ±2 super-harmonics.

Tower side-side mode
The tower side-side mode was excited with a chirp-shaped force applied at the tower top. The frequency band of such signal 10 was set in order to excite only that single mode. The tower base side-side moment was then recorded and used as output. As only the tower side-side peak is visible in the PSD of the response, then N a was set equal to 2. The other coefficients were set as N Fa = 1, N b = 1, N F b = 1, N g = 2 and N Fg = 1.
The agreement between the output predicted with the PARMAX reduced model and the measure, not shown here for the sake of brevity, is very good. The left plot of Fig. 8 shows the Campbell diagram obtained with the PARMAX approach. In 15 this diagram the results of the identifications are approximated with straight lines. Looking at this plot, it appears that at 0.8Ω r the principal harmonic intersects the 2×Rev. For the PARMAX identification this is not particularly problematic, and in fact only the participation factor has been slightly underestimated. On the other hand, this poses a major problem for POMA. In fact, when the signal is frequency-shifted by +2Ω, its average value is transported over the principal peak, making it difficult to estimate the mode shape and the damping of the tower side-side mode. 20 The Campbell diagram obtained from POMA identifications is shown in the right plot of Fig. 8. The plot clearly shows that the damping of the principal harmonic estimated with the half-power bandwidth is double the one estimated by PARMAX.

Backward and forward whirling in-plane modes
The backward and forward whirling in-plane modes were excited with a tower top side-side doublet, whose amplitude and duration were selected such that the input force spectrum is almost flat in the frequency range of interest. The three blade root 25 edge-wise bending moments M 1 , M 2 and M 3 were recorded and the multi-blade coordinate transformation was used to yield the direct and quadrature moments, noted respectively M d and M q . The spectra of M d , displayed in Fig. 9, show well defined peaks. The PARMAX reduced model was set with the following choice of parameters: N a = 8, N Fa = 1, N b = 1, N F b = 1, N g = 2 and N Fg = 1. Both the backward and forward whirling in-plane modes, as well as the side-side tower mode, were nicely visible in the frequency plot of the perturbed time histories. Thus, for each wind speed, only one reduced model capable of representing the behavior of all these three modes was identified. Figures 10 and 11 show at left the periodic Campbell diagram obtained using the PARMAX approach, and at right the one 5 computed with POMA, respectively for the backward and the forward whirling in-plane modes. It should be noticed that both approaches provide the same results in terms of frequencies. The overall trend of the principal harmonic damping factors as functions of the rotor speed is similarly captured. In particular, the PARMAX results are characterized by a lower uncertainty for the backward mode and a higher uncertainty for the forward one. The rising of the damping factors with the angular speed, for these two modes, has been observed also in Skjoldan and Hansen (2013), although for an isotropic condition.

10
Once again, the damping of the super-harmonics obtained with the POMA technique are not well estimated, as already noticed in Sect. 4. Moreover the participation factors of the ±2 super-harmonics are typically too high: for example, in the right plot of Fig. 11 one may notice that the participation of super-harmonic +2 of the forward whirling mode is higher than that of the principal one. This strongly overestimated participation is due to the nearly 2Ω spacing of the whirling modes, which causes their super-harmonics to nearly overlap.

Conclusions
In this paper we have considered a model-independent periodic stability analysis capable of handling turbulent disturbances.
The approach is based on the identification of a PARMAX reduced model from a transient response of the machine. The full-Floquet theory is then applied to the reduced model, yielding all modal quantities of interest. As only time series of measurements are necessary, the method appears to be suitable for the application to real wind turbines operating in the field. 20 In order to assess the validity of the proposed method, the well known POMA was implemented and used for comparison.
Tests were performed first with the help of a wind turbine analytical model, whose exact solution can be obtained by the theory of Floquet, and then with a high-fidelity wind turbine multibody model operating in turbulent wind conditions.
Based on the results obtained in this study, one may draw the following considerations: -Both methods are able to characterize the relevant behavior of the wind turbine in turbulent wind conditions. However, 25 the results provided by the proposed PARMAX analysis are in general more accurate than those given by the POMA technique, especially if one looks not only at the principal harmonics but also at the super-harmonics.
-Often the underlying hypothesis of POMA are not exactly fulfilled, and this leads to inaccuracies especially in terms of damping and participation factors. These effects are more visible for the whirling modes, as they are separated by about 2Ω, which means that there will always be a perfect overlap between the super-harmonics of these two modes at some 30 angular velocity. The PARMAX analysis is less prone to such problems. conditions, where the rotor speed is hardly constant (which, on the other hand, is a fundamental hypothesis of both methods).
The development of the present SISO-PARMAX approach suggests a number of extensions, which are currently under investigation:

5
-The use of multiple outputs in a multiple-input multiple-output (MIMO) PARMAX framework could improve the quality of the results.
-Due to the stochastic nature of turbulence, a multi-history PARMAX applied to different realizations of the same experiment could provide more robust modal results, along with the associated variances.
-The peak-picking method is rather simple and it is unable to exploit all the informational content available in the HPSD, 10 especially in the presence of noisy peaks. Fitting algorithms have been preliminary explored (see Allen et al., 2011a), but their application to the multiple output case has not yet been attempted.

A1 Floquet theory in continuous time
A generic SISO LTP system in continuous time can be written in state-space form as where t is time, x, u and y the state, input and output vectors, respectively, while A(t), B(t), C(t) and D(t) are periodic system matrices such that for any t. The smallest T satisfying Eq. (A2) is defined as the system period. Vector u contains the wind turbine control inputs (i.e. blade pitch angles, electrical torque, possibly the yaw angle) as well as exogenous inputs related to the wind states (e.g. wind speed, vertical or lateral shears, cross-flow, etc.).
To study the stability of Eq. (A1a), its autonomous version is considered together with the associated initial conditions: The state transition matrix Φ(t, τ ) maps the state at time τ , x(τ ), into the state at time t, x(t): and it obeys a similar equation with its associated initial conditionṡ where I is the identity matrix. It can be shown that in the continuous-time case the transition matrix is always invertible (Bittanti and Colaneri, 2009).
An important role in the stability analysis of periodic systems is played by the state transition matrix over one period 5 Ψ(τ ) = Φ(τ +T, τ ), termed monodromy matrix. By definition, the monodromy matrix relates two states separated by a period; consequently, a generic state that is sampled at every period, notedx τ (k) = x(τ + kT ), obeys the following linear-invariant discrete-time equatioñ The system is asymptotically stable if all the eigenvalues of the monodromy matrix, called characteristic multipliers and noted 10 θ j , belong to the open unit disk in the complex plane. It can be shown that the eigenvalues of the monodromy matrix and their multiplicity are time-invariant even if the monodromy matrix is periodic (Bittanti and Colaneri, 2009). For this reason, one can ignore the time lag τ when referring to the characteristic multipliers. The eigenvalues θ j and associated eigenvectors s j are obtained by the spectral factorization of the monodromy matrix, i.e.
In order to determine the frequency content of a periodic system, it is necessary to introduce the so called Floquet-Lyapunov transformation. The Floquet-Lyapunov problem is the one of finding a bounded, periodic and invertible state-space transformation z(t) = Q(t)x(t) such that the resulting governing equatioṅ is time-invariant, i.e. the Floquet factor matrix R is constant. Since R = Q(t)A(t)Q −1 (t) +Q(t)Q −1 (t), the periodic transformation Q(t) must obey the following matrix differential equatioṅ whose solution is Exploiting the periodicity condition Q(τ + T ) = Q(τ ), one gets the relationship between monodromy matrix and Floquet factor, which writes The eigenvalues of the Floquet factor, called characteristic exponents and noted η j , are computed by the spectral factorization of R: with V = [. . . , v j , . . . ]. Inserting Eqs. (A7) and (A12) into (A11), the following result is derived 5 which shows that V = Q(τ )S and, more importantly, that characteristic multipliers and characteristic exponents are related as Notice that there is an infinite number of Floquet factors, and therefore an infinite number of Floquet-Lyapunov transformations. In fact, one can choose any invertible initial condition Q(τ ). In addition, computing characteristic exponents from 10 multipliers by inverting Eq. (A14) leads to a multiplicity of solutions, as in fact where ∈ Z is an arbitrary integer. This indeterminacy, however, does not affect the real frequency content of the response, since the transition matrix is uniquely defined. This aspect of the problem will be further analyzed later on in these notes.
Given Q(τ ) and R, the transition matrix is readily obtained from Eq. (A10) as where the periodic matrix P (t) = Q(t) −1 is termed periodic eigenvector.
Consider now, for each mode, one of the infinite solutions of Eq. (A15), for example the one with = 0, notedη j . Introducing Ω = 2π/T , any other characteristic exponent η j could be computed fromη j as η j =η j + ı nΩ, n ∈ Z.
(A17) 20 Inserting Eq. (A12) into (A16), one can express the state transition matrix as the following modal sum where Z j (t, τ ) = P (t)V I jj V −1 Q(τ ), while I jj is a matrix with the sole element (j, j) equal to 1 and all others equal to 0.
Because of the particular definition of I jj , matrix Z j (t, τ ) is of unitary rank ∀(t, τ ) and it is also equal to ψ j (t)L j (τ ) T , where 25 ψ j (t) = col j (Ξ(t)), with Ξ(t) = P (t)V . Equation (A18) can be now reformulated as Exploiting the periodicity of ψ j (t), Eq. (A20) becomes where ψ j n are the amplitudes of the harmonics of the Fourier expansion of ψ j (t). This expression of the state transition matrix 5 can also be found in Skjoldan and Hansen (2009);Allen et al. (2011b); Bottasso and Cacciola (2015).
From Eq. (A21) it appears that, for each mode, an infinite number of exponents (playing the role of eigenvalues of the LTI system) participates in the response of the system. Furthermore, a single frequency is not sufficient for completely characterizing that mode. All exponents have imaginary parts that differ by integer multiples of Ω and have the same real part; thus, all exponents of a given mode are either stable or unstable. This fact is not surprising, as the stability of the system is just 10 determined by the characteristic multipliers, which are uniquely defined.
For the LTP system, the exponentsη j + ı nΩ play the role of the eigenvalues of the LTI case, as they yield the frequencies ω j n = |η j + ı nΩ| and damping factors ξ j n = −Re(η j )/ω j n of each mode. To describe this situation, this infinite multiplicity of frequencies is termed a fan of modes (cf. Bottasso and Cacciola, 2015). Each harmonic in a fan contributes to the overall response according to its associated "modal shape" ψ j n . The relative contribution of the nth harmonic to the jth mode is 15 measured through its participation factor, defined as φ j n = ψ j n n ψ j n .

(A22)
The triads {ω j n , ξ j n , φ j n } describe completely the behavior of a periodic mode. The participation factors can be defined also as functions of the Frobenius norm of the harmonics of Z j (t, τ ), Z jn = ψ j n L j (τ ) T , as shown in Bottasso and Cacciola (2015): The two definitions are exactly equivalent as, in this specific case, Z j n F = ψ j n L j (τ ) and L j (τ ) stays the same for all Often, although not always, the harmonic with the highest participation is very similar in terms of frequency and damping to the one that would results from the invariant analysis of periodic systems based on the Coleman transformation (Coleman and Feingold, 1958;Hansen, 2004). As suggested by Bottasso and Cacciola (2015), such harmonic may be called the principal one, while the others may be termed super-harmonics. Furthermore, any one of these harmonics could resonate with external excitations.

5
In order to understand how each harmonic appears in a specific output of the system, the output-specific participation factor can be defined. To this end, consider an output of the autonomous system (A3), Inserting Eq. (A20) into (A24) the following is derived Exploiting now the periodicity of the product C(t)ψ j (t), Eq. (A25) can be rearranged as where c j n are the harmonics of the Fourier expansion of C(t)ψ j (t). The output-specific participation factor can finally be defined as φ y j n = |c j n | n |c j n | .
(A27) 15 A2 The harmonic transfer function and the harmonic frequency response function The forced response of system (A1), named y F (t), can be computed as where is the impulse response. From Eq. (A28), it appears that the periodicity of C(t), B(t) and Φ(t, τ ) causes the input-output behavior of a LTP system to be far from being describable as a LTI-like one. In particular, it can be shown that a LTP system subjected to an input at a given frequency may respond at an infinite number of frequencies, which in addition to the input frequency itself include also the integer multiples of the system frequency (Bittanti and Colaneri, 2009;Wereley, 1991). This is also the reason why any output of a wind turbine subjected to a constant-in-time wind (i.e. at the zero frequency) is characterized In the frequency domain, the input-output relation can be expressed by means of the HTF (cf. Bittanti and Colaneri, 2009;Wereley, 1991), which can be interpreted as the extension to periodic systems of the standard time-invariant transfer function.
To this end, the so-called exponentially modulated periodic (EMP) signal is defined as where s ∈ C. According to definition (A30), any v k can be also viewed as the Laplace transformation of v(t) evaluated at It can be shown that a periodic system subjected to an EMP admits an EMP regime (Bittanti and Colaneri, 2009), and that in such a regime its states are EMP signals. In order to exploit this property, one has first to define two doubly infinite-dimensional vectors containing respectively the EMP harmonics u(t) and y(t), as 10 Y(s) = · · · y −1 (s) y 0 (s) y 1 (s) · · · T , (A32a) Next, the doubly-infinite Toeplitz matrices A, B, C and D, containing the Fourier expansions A k , B k , C k and D k of the corresponding system matrices, are defined as 15 and similarly for the B, C and D matrices. Finally, by inserting the EMP expansions of y and u and the Fourier expansions of the system matrices into Eq. (A1), summing up all terms at the same frequency, the input-output relationship is derived as where the HTF is defined as 20 with N = blkdiag{ıkΩI, k ∈ Z}, being I and I identity matrices of suitable dimensions.
The HTF can also be represented by means of the impulse response of the system (Bittanti and Colaneri, 2009). From Eq. (A29), it is easily verified that function h(t, t − r) for a fixed time lag r is periodic and, consequently, that it can be expanded in a Fourier series as The output equation can be then written according to the following convolution which leads to the input-output relation in the Laplace domain where Y (s), U (s) and H k (s) are respectively the Laplace transforms of y, u and h k . Equation (A38) can be evaluated for each element of the EMP output signal Y by substituting the complex number s with the exponentially modulated periodic one s + ıkΩ with k ∈ Z, leading to the following relationship Consequently, since Y (s + ıkΩ) = y k (s) and U (s + ıkΩ) = u k (s) because of Eq. (A31), the HTF can be written as Inserting (A21) into Eq. (A29), one can derive the following expression where the product L T j (τ )B(τ ) and D(τ ) have been expanded in Fourier series, being l j m and d k the related amplitudes. After 15 some manipulations (see also Wereley, 1991;Wereley and Hall, 1990), the Laplace transformation of h k (t − τ )e −ınΩ(t−τ ) can be finally written as Consider now the row-index ∈ Z and the column-index r ∈ Z of the HTF, defined such that the element with = r = 0 (noted G 0,0 ) corresponds to the median element H 0 (s), and the element with = r = −1 (noted G −1,−1 ) to H 0 (s − iΩ). 20 Hence, according to such definitions and thanks to Eq. (A42), the following holds Consequently, the HTF can be computed as and D = D.
From a practical standpoint, the use of the harmonic input-output relation expressed by the HTF implies that one has to 10 consider a truncated finite dimensional approximation of G(s), which corresponds to the use of truncated versions of the EMP input and output signals. The convergence of truncated HTFs has been discussed in Sandberg et al. (2005).

A3 The discrete-time case
In this section the stability analysis of periodic discrete-time systems is briefly reviewed. For a more comprehensive treatment, the reader is referred to Bittanti and Colaneri (2009) and Bottasso and Cacciola (2015). 15 The autonomous dynamic equation of a generic LTP system in discrete time and its initial conditions are where k is a generic time instant and A(k) is a periodic matrix of period K such that A(k + K) = A(k), ∀k. Similarly, the transition matrix obeys the following equation with its initial conditions In this work we consider only reversible systems, i.e. those for which det(Φ(k, κ)) = 0, ∀(k, κ).
For reversible discrete-time systems, the state transition matrix Φ(k, κ) can be decomposed in periodic and contractive parts as where P (k) is periodic and R is constant. Here again, the system is stable if the characteristic multipliers θ j , i.e. the eigenvalues of the monodromy matrix Ψ(κ) = P (κ)R (K) P (κ) −1 , belong to the open unit disk in the complex plane. The relationship between characteristic multipliers and characteristic exponents is In the discrete-time case, the apparent multiplicity of the characteristic exponents manifests itself as a phase indetermination where = 0, . . . , K −1 is an arbitrary integer. As in the continuous-time case, this does not in reality generate any inconsistency as frequencies, damping and participation factors of the various harmonics are unaffected by this apparent arbitrariness.
Following the same approach of the continuous-time case, the transition matrix can be rewritten as where Ξ(k) = P (k)V and After having expanded ψ j (k) in Fourier series, one gets where now ψ j n are the amplitudes of the harmonics of the Fourier expansion of ψ j (k). Coherently, the multiplication of Eq. (A54) with C(k) leads to being c j n the harmonics of the Fourier expansion of C(k)ψ j (k). This shows that the jth mode is characterized by K exponents 20 with the same modulus and different phases. Each exponent can be transformed into the continuous one using the following expression (cf. Franklin and Powell, 1980) where ∆t is the sampling time and subscripts (·) c and (·) d refer, respectively, to the continuous and discrete-time cases. Once the continuous-time exponents are computed, frequencies, damping and participation factors can be readily obtained as in the The simplified upwind horizontal-axis wind turbine model used in this work, depicted in Fig. 12, considers the coupled motion of tower and blades. The tower fore-aft and side-side flexibility are rendered by two equivalent linear springs and dampers.
Each blade is modeled as two rigid bodies connected to each other by means of two equivalent revolute joints, which allow respectively the blade flap and edge-wise rotations. The inner part of the blade is rigidly connected to the hub. Each joint is 5 associated to a rotational spring and a rotational damper. The inertial and structural characteristics of each element are chosen so as to match the first tower fore-aft and side-side mode and the first blade flap-wise and edge-wise modes in vacuo, computed using a high-fidelity multibody model of the wind turbine.
The reference frame used for the derivation of the equations of motions has its origin located at the hub, the x axis directed downward, the z axis directed from the tower to the rotor, and the y axis selected so as to form a right handed triad. To simplify 10 the notation, in the following subscript k, denoting the blade number, will be dropped together with the time dependence whenever possible.
The contribution of the two blade parts to the total energy can be developed separately. Thus, let r U and r D indicate respectively the dimensional abscissa along the inner and the movable parts of the blade, respectively. The position of a generic blade point is given by when the point belongs to the inner part of the blade, and by e cos ψ + r D cos β cos(ψ + ζ) y H + e sin ψ + r D cos β sin(ψ + ζ) when it belongs to the movable part. The kinetic energy of the whole rotor is obtained by summing up the kinetic energy of the hub, T H , and of both the inner and the movable parts of the kth blade, respectively noted T D k and T U k , resulting in where and being ρ(r) the blade mass per unit span.
All springs and gravity contribute to the potential energy of the system as where the potential energy of the side-side and fore-aft springs are defined respectively as V y H = 1/2K y y 2 H and V z H = 1/2K z z 2 H , while that of the flap-wise and edge-wise springs as V β k = 1/2K β β 2 k and V ζ k = 1/2K ζ ζ 2 k . Finally the contribution of gravity can be expressed as The damping function D follows a rather similar procedure, where The aerodynamic model is based on a linearized BEM approach with constant aerodynamic properties along the blade, mostly taken from Eggleston and Stoddard (1987), with the addition of the hub velocity (ẏ H ,ż H ) to the inflow and cross-flow 15 terms, but neglecting the yaw rate. Table 8 gives the meaning of some symbols used in the following equations.
The hub shear force in the fore-aft direction is The hinge out-of-plane moment is 20 The hub shear force in the direction parallel to the chord of the blade, and pointing towards the leading edge, is The hinge moment in the edge-wise direction is This aerodynamic model assumes that the wind velocity varies linearly over the rotor disc, and therefore it is not suited to 5 simulate turbulent wind fields.
The virtual work of the aerodynamic forces and moments results to be The generalized forces follow directly from the previous expression.