Nonlinear modelling and adaptive control of smart rotor wind turbines

This paper develops a nonlinear mid-ﬁdelity aeroservoelastic model for smart rotor wind turbines and studies the turbulent load alleviation of the wind turbines with trailing edge ﬂaps (TEFs) actuated by a novel proportional-derivative model-free adaptive control (PD-MFAC) algorithm. This nonlinear model contains a structural model for the tower and blades represented by geometrically non-linear composite beams and an aerodynamic model for the rotor using a vortex panel method coupled with a stall delay process. The capability of the new aerodynamic model to deal with ﬂow separations and analyze the detailed ﬂow ﬁeld enables it to simulate the dynamic response of the wind turbine blade sections with TEFs with arbitrary size and deﬂection angles. It is shown that the TEF alters the aerodynamic coeﬃcients in a complex manner which could be explained by the evolution of the detailed vortical ﬁeld. Furthermore, three independent PD-MFAC TEF controllers are designed to alleviate the turbulent load acting on the wind turbines. The eﬀectiveness of the controller in terms of turbulent load alleviation is evaluated by the root mean square value of the blade root-bending moment (RBM). A traditional Gain-scheduled PI (GS-PI) controller is also designed as a comparison to the PD-MFAC controller. Simulation results show that to reduce the RBM and blade tip deﬂection (BTD) caused by external disturbances, the PD-MFAC ﬂap controller shows more eﬀective performance than the GS-PI ﬂap controllers.


Introduction
To reduce greenhouse gas (GHG) emissions and achieve energy sustainability, renewable energies are becoming one of the most popular energy sources in the world. Among different renewable energy types, wind power is standing out as one of the fastest-growing and competitively commercialized energy sources [1]. It is expected to increase its global demand by up to 60% by the year 2050. To 5 increase the average power capacity and decrease the Levelized Cost of Energy (LCOE), continuous efforts have been made in the design of larger wind turbines with an increasing rotor diameter, and hub height [2]. To this end, it is inevitable to have a more flexible and slender design for longer and heavier blades. As is well established, due to the inflow turbulence, wind shear, tower shadow, wake effects, gravity, mass and aerodynamic imbalances, and other sources of disturbance, 10 wind turbines are subject to significant fluctuating loads in different frequency domains. Some of these unsteady loads may cause blade structural resonance and fatigue damage, some of those may contribute to structural failure and lead to a reduction of life expectancy [3]. The increasing size of wind turbine blades is accompanied by an increase in the susceptibility to gravitational effects and in wind speed variations across the rotor disk, which add difficulties in mitigating the fluctuating 15 loads acting on the flexible blades. New challenges have been brought about in both modeling and control design of wind turbines. The development of higher fidelity models and the optimization of various flow-control devices to reduce the fluctuating loads, prolong the service life of the wind turbines, will be of great significance for the following years [4].
In terms of numerical models of horizontal axis wind turbines (HAWTs), a structural model 20 and an aerodynamic model are combined to build a full aeroservoelastic model. There are mainly two categories of structural models according to the discretization methods used to present the flexible bodies, i.e. three-dimensional (3D) finite-element method (FEM) and one-dimensional (1D) composite beam models. 3D FEM is an extraordinary method to obtain detailed stress distributions within a structure. However, it is too computationally expensive to be used in testing the control 25 strategies. There are linear composite beam models and nonlinear composite beam models. The most commonly used linear beam formulations are the classical Euler-Bernoulli beam theories and the Timoshenko beam theories for straight and uniform beams, containing the assumption of small displacements [5]. The latter contains shear deformation effects, while the former does not. The wind turbine blades normally have a thin and slender structure where the shear deformations effects 30 are negligible, thus the Euler-Bernoulli beam theory is more often used for simplicity. Furthermore, linear beam models are no longer sufficient to determine the large deformations and highly nonlinear structural response on a smart rotor wind turbine structure. Thus geometrically-exact composite beam theory, developed from the original Euler-Bernoulli beam theory by enabling non-uniform deforming of the cross-section within the linear regime, was developed [6]. There are a wide variety 35 of approaches to modeling HAWT aerodynamics, ranging from the commonly used blade element momentum (BEM) method [7], the actuator type models [8], the free vortex wake (FVW) models [9,10], and the unsteady vortex lattice methods (UVLM) [11,12,13] to higher fidelity methods, such as the CFD methods [14]. CFD models are the most accurate while the most computationally intensive models among others. BEM is the least computational cost while with the lowest accuracy. 40 To provide a good balance between the computational cost and the accuracy, this work will propose a nonlinear aeroservoelastic model which includes a nonlinear Euler-Bernoulli beam theory and a discrete vortex method.
In terms of the control strategies in the load alleviation of wind turbines, there are two main categories, i.e. passive and active load alleviation methods. Examples of passive load alleviation 45 strategies include morphing structures on wind turbine blades [15]. Such passive load alleviation systems are normally not easily maintainable [16]. Among the active control methods, individual blade pitch control (IPC) which alters the aerodynamic properties of the blades is the most widely studied in the past [16,17,18]. Rezaeiha et al. [19] studied the total fatigue loads on the wind turbines caused by different sources and showed that wind turbulence results in over 65% of flapwise 50 fatigue, while gravity contributes to over 80% of edgewise fatigue. The low-frequency deterministic loads, caused by wind shear, tower shadow or gravity, could be well reduced by IPC, while for mitigating the high-frequency non-deterministic loads, e.g. turbulent load, on wind turbine blades, IPC is less effective due to its reaction speed restricted by the actuator limitations. Moreover, using IPC to mitigate the blades' fatigue load is compensated by an increase in pitch activity, 55 thus the design of IPC is a trade-off between the fatigue damage of the blades and blade pitch actuators [20]. Thus, the concept of 'smart rotor control' has been proposed for a more precise and responsive turbulence load alleviation of wind turbine blades as a compliment to IPC [21,22].
The main advantage of the smart rotor is dealing with high-frequency loads caused by turbulence rather than the low-frequency deterministic load. Smart rotors refer to the distributed control 60 surfaces equipped on turbine blades that can be actuated rapidly. Among different types of 'smart rotors', the trailing-edge flaps (TEFs) demonstrated excellence in mitigating loads on wind turbines.   in a PID control scheme depend on the specific operating conditions, therefore the control efficiency will be inevitably reduced in off-design conditions.
In our previous paper [29], a novel proportional-derivative model-free adaptive control (PD-MFAC) scheme has been developed and tested on a preliminary low-fidelity model of a wind turbine configuration. In the PD-MFAC, in order to build a linearized equivalent dynamical model along

Numerical method
The aeroservoelastic model of a wind turbine [30] with smart rotors will be adapted from the low-fidelity model developed in our previous work [29]. Although the proposed aeroelastic model and C M 0 ) could be found in our previous paper [29].
In the numerical simulations, the wind turbine is operating in a turbulent wind field with a mean inflow speed V ∞ and a specific turbulent intensity. The turbulence inflow data w δ , in accordance with the International Electrotechnical Commission (IEC) standard 61400-1 [31], will be generated by a well-known open-source stochastic inflow turbulence tool (TurbSim) developed by National

145
Renewable Energy Laboratory (NREL). As the considered model wind turbine has a rated wind speed of 11.4 m/s, mean wind speeds of 11m/s and 14m/s are chosen to cover the typical operating speed range. As the turbulence intensity at the mean wind speed of 11 -14 m/s ranges from 5% to 20% [32], three typical turbulent intensities from this range are chosen, they are 6%, 10% and 17.5%, respectively. The von Karman turbulence model and the Normal Turbulence Model (NTM) 150 with different turbulence seeds are chosen in TurbSim. Yaw angles are not considered since it is out of the scope of this work.

Structural Model
The same structural model as in the previous work [29] is adopted here. It is a full 6-degree of freedom (DOF) geometrically-nonlinear Euler-Bernoulli beam model where the tower and the blades are defined as a collection of geometrically nonlinear composite beams respectively. The wind turbine configuration and its discrete multi-body composite beam representation are shown in Fig. 1. The intrinsic beam formulations are closed as follows according to Hamilton's principle, , ,  [29]. In this work, a bigger flap size is considered for better control performance than the small flaps. The TEF configuration considered here takes up 20% of the wingspan and 15% of the local chord as shown in Fig. 2. We will model the dynamics of the TEF 160 by modifying the shape of the wing sections (the airfoils) continuously, instead of modifying the local aerodynamic coefficients of the lifting surface. In Fig. 2, the local angle of attack is denoted as α, while the flap deflection angle is denoted as β.

Aerodynamic model
The aerodynamic model is a vortex panel method with a stall delay algorithm, where the sepa- In subsection 2.2.1, we will recall the standard VPM, and in subsection 2.2.2, we will introduce the dynamic stall delay model to determine the dynamic separation point. Since the state-space description of the aerodynamic model is the same as in our previous paper [29], we omitted the 175 detailed formulas here.

Vortex panel model (VPM)
Given the specific geometry of the airfoil deformation and the dynamic separation position, the problem described in Fig. 3 becomes a standard VPM. In this model, the N +3 unknown parameters, Here t and t − 1 refer to the current time step and the last time step respectively, and ∆l SEP , ∆l T E are the length of the separating vortex sheet and the wake sheet respectively. Note that once the vorticity is shed to the flow field as point vortices, their strength remains constant in time until they are merged by coalescence. A coalescence criterion in the far wake has been used to restrict the total number of point vortices. In each time step, the strength of the vortex distributions are resolved, the induced velocity field and the velocity potential are then obtained. The reader could refer to Katz (1981) [33] for more details about solving the standard VPM. The local aerodynamic coefficients of the lifting surface C L , C D , and C M are obtained by integrating the pressure distribution over the airfoil given by the unsteady Bernoulli Equation This VPM has been validated with classical cases such as the starting flow around a thin airfoil at an incidence of 45 o .

A dynamic stall model to determine the dynamic separation point 180
The steady-state separation point for a 2D airfoil with TEFs as a function of the AoA (α) and the TEF deflection angle (β) could be obtained in several ways, such as XFOIL tool, the extended Kirchhoff's Law, or experimental measurements. For simplicity, here we will use an extended Kirchhoff's Law to determine the separation point of the airfoil with TEF. Here the steady-state lift coefficient is given by With the steady-state separation points given, two time lags, i.e. the time-lag due to the LE pressure (T p ) and the time-lag due to the dynamics of the boundary layer (T f ) are introduced here to calculate the dynamic (stall-delayed) separation points. The dynamic separation point after the first delay is determined by where C p L (t) is the unsteady lift coefficient for attached flow, and T p = 1.7 is the time constant for the separated flow [35]. C p L (t) is the lift coefficient with a time-lag to the attached flow lift coefficient C p L (t). Likewise, the dynamic separation point after the second time lag is determined by where T f = 3 is the time constant for the separated flow which can be determined from unsteady aerofoil data [35].       edge vortex, which reduces the lift coefficient through an upwash effect on the lower surface of the TE. When the angle attack is at the lowest value (point '2' in Fig. 8 (b)) or approaching the mean 215 value in the up-stroke (point '3' in in Fig. 8 (b) Similarly, the comparison of the drag coefficient hysteresis curves with different flap dynamics 225 β = 10 o sin(ωt + φ) (φ = 0, π/2, and π) are shown in Fig. 9. The comparison results for the moment coefficients are shown in Fig. 10. The drag coefficients and the moment coefficients could be analyzed similarly to the lift coefficients. For the relative small angle of attack case (α = 8 o ± 10 o sin(ωt)) in Fig. 9 (a) and Fig. 10 (a), we can see that the drag coefficient hysteresis curves rotates clockwise and the the moment coefficient hysteresis curves rotates counterclockwise

High-pass filter
Gain-Scheduled PI controllers are three flap pitch angles. A closed-loop block diagram of the fluid-structure coupling and control system is depicted in Fig. 11.
The derivation and the stability analysis of the PD-MFAC algorithm could be found in our previous work [29]. Here we only give the final expression of this algorithm and the parameters.
For a SISO nonlinear system expressed in the time discrete form y (k + 1) = f y (k) , ..., y k − n y , u (k) , ..., u (k − n u ) , where u (k) is the input signal and y (k) is the output of the nonlinear SISO system at discrete time k, n y and n u are the known orders, which are positive integers, a virtual control input (the input of the equivalent system shown in Fig. 11) is defined as where K 1 and K 2 are constant control parameters, ∆y (k) = y (k) − y (k − 1) is the time difference of the system output, and ∆t is the time step for the discrete system. By introducing the PPD vector φ * L , the system could be linearized along the operation point as follows The PD-MFAC algorithm can be derived as time difference of historical virtual control input vector U * L (k). The virtual control input vector is obtained by U * L (k) = u * (k) , ..., u * (k − L + 1) . For our considered case, L = 3. The step factors ρ i ∈ (0, 1] and η ∈ (0, 2], and the penalty factors satisfy λ > 0, µ > 0. The parameters in the PD-MFAC are tuned by using a particle swarm optimization (PSO) algorithm. In the PSO algorithm, the number of swarms is set to be 20 and the number of maximum 250 iterations is set to be 50. The searching process will continue until the residual converges to 1e-12.
Here is the results tuned for this PD-MFAC scheme: k 1 = 9.8366 × 10 −08 , K 2 = 3.91 × 10 −09 , λ = 8 × 10 8 , ρ = 6.684 × 10 −3 , η = 9 × 10 −4 , µ = 9 × 10 −3 , φ As a comparison to the proposed PD-MFAC scheme, a traditional GS-PI controller is also designed. In this controller, the proportional and integral terms are defined as where y is the flapwise RBM and y ref is the reference flapwise RBM when there is no inflow turbulence. The three parameters k p , k i , and k g in the GS-PI controller are also tuned by a PSO 260 algorithm. The parameters obtained by the PSO algorithm are k p = 5.36e − 02, k i = 8.79e − 08, and k g = 9.8e08.

Numerical Simulations
The open-loop and closed-loop numerical simulation of the NREL 5-MW reference wind turbine with TEFs will be carried out in this section. Using the aeroservoelastic modeling method described purely casued by turbulence, respectively. We can see from these figures that the density bars for the open-loop simulations are more centralized than those for the GS-PI controller simulations.
And the density bars for the GS-PI controller simulations are more centralized than those for PD-MFAC simulations. Thus we can come to a conclusion that, an obvious reduction in the blade root        reduction of the full range of frequency load. It is shown in Table 2 that the proposed PD-MFAC has an evident advantage over the traditional model-free adaptive GS-PI controller concerning the 340 reduction of flapwise RBM caused by turbulence. Additionally, as a bonus, the PD-MFAC controller also outperforms the GS-PI controller in the alleviation of BTD caused by turbulence.

Discussion
An aeroservoelastic model for a wind turbine with TEF configurations has been developed and a model-free PD-MFAC smart rotor controller has been designed for the load alleviation of this

Conclusions
This work has studied the modeling and load alleviation control method of smart rotor wind In the future work, (a) the optimal size of the TEF which offers a better control performance could be studied; (b) the wake effect and the three-dimensional effect could be included; (c) the load 395 alleviation of a larger wind turbine (e.g. 20MW), could be studied using the proposed model and control strategy.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.