Project:Sandbox

Mehrle : Template of Abstract for Mechatronics Master’s Proceedings

Wind turbine, CFD, FEM, For-aft tower oscillation, Pitch angle control, Multi-system closed-loop.

= Introduction =

work is focussed on the damping of the tower axial oscillations of a wind turbine for nominal operation conditions. For this case, an active pitch controller is designed in order to reduce the influence of structural oscillations acting on the tower top. The controller adjusts the blade pitch angle so as to modify the aerodynamic thrust, with the objective of actively damping or compensating unwanted for-aft oscillations, particularly at undesirable tower vibrational modes.

At first, the system is defined. On one side, the wind turbine modelled is 5.XM series offshore wind power plant provided by REpower Systems, Senvion. It is designed using the CAD software-package CATIA (with the help of MATLAB) and latter modelled using the numerical Computational Fluid Dynamics approach in ANSYS-FLUENT. In addition, the wind turbine offshore supporting structure has a basic mono-pile type.

Then an appropriate control system is developed. Thus, a simplified non-linear mathematical model of the tower structural dynamics and the turbine performance is made, which is then linearised around the nominal operation conditions and working points. In addition, an excitation force oscillates sinusoidally at the tower first mode of vibration representing undesired vibrational forces during operation. The closed-loop transfer function is then derived and a controller is developed and tested. Finally, the connection and interface between both systems is searched for, which will lead to the desired multi-physics simulation system.

Wind Turbine Aerodynamics
Wind turbines perform at given relations between the rotational and wind speed together with an optimal angle of attack which maximize the power extraction during operation. Each turbine’s performance is represented via the torque, thrust and power coefficient performance curves. The wind turbine modelling comes from the lift and drag forces which are extracted from the aerodynamic characteristics of every single differential blade element: $$\begin{aligned} dL &=& \frac{1}{2}\rho c_{L}(\alpha,\lambda) v_{rel}^2 c(r)\partial{r}\\ dD &=& \frac{1}{2}\rho c_{D}(\alpha,\lambda) v_{rel}^2 c(r)\partial{r} \end{aligned}$$ The factor tip speed ratio is a key parameter which defines the ratio between the rotational speed at the tip of the blades and the undisturbed wind speed far upstream; $$\begin{aligned} \lambda = \frac{\omega R}{U_{\infty}}\end{aligned}$$ Therefore, by introducing the afore parameter and inserting the torque and thrust coefficients: $$\begin{aligned} T &=& \frac{1}{2}\rho AU_{\infty}^{2}c_{T}\\ M &=& \frac{1}{2}\rho ArU_{\infty}^{2}c_{M}\end{aligned}$$ In which the thrust coefficient $$c_{T}(\alpha,\lambda)$$ and the moment coefficient $$c_{M}(\alpha,\lambda)$$ are introduced, which primarily are dependant on the differential lift and drag coefficients of the turbine, however modelled to depend solely on the optimized angle of attack, $$\alpha$$ and on the tip speed ratio, $$\lambda$$, providing a distinct simplification of the turbine model. The power extracted, when no gear box is defined: $$\begin{aligned} P = \frac{1}{2}\rho AU_{\infty}^{3}c_{P}\end{aligned}$$ where the power coefficient $$c_{P}$$ is defined as: $$\begin{aligned} c_{P} = \lambda c_{M}\end{aligned}$$ The forces acting over a differential blade element are depicted in the following Figure [fig:AIRF].

(airfoill) at (0,0) ;

(v1) at (0,0) ; (v2) at (-7,0) ; (v3) at (7,0) ; (v4) at (0,7) ; (v5) at (0,-7) ; (v6) at (-6,-2.6) ; (v7) at (8,3.5) ; (v8) at (-7,-6) ; (v9) at (-7,0) ; (v10) at (-7,0) ; (v2)edge(v3); (v4) edge (v5); (v6) edge (v7);

(v8) edge (v1);

(v11) at (-2,-10) ; (v12) at (-1,-10) ; (v13) at (0,-10) ; (v14) at (1,-10) ; (v15) at (2,-10) ; (v16) at (-2,-8) ; (v17) at (-1,-8) ; (v18) at (0,-8) ; (v19) at (1,-8) ; (v20) at (2,-8) ;

(v11) edge (v16); (v12) edge (v17); (v13) edge (v18); (v14) edge (v19); (v15) edge (v20); (v16) edge (v20);

(-4.5,0) arc (180:201:5.0cm); (-4.1,-1.87) arc (201:215:5.0cm); (-6.5,0) arc (180:215:7.5cm);

at (-4.7,-2.9) $$\alpha_{c_{L,max}}$$ ; at (-5.1,-1.1) $$\beta(r)$$ ; at (-7.3,-2.2) $$\gamma(r)$$ ;

at (-2.8,-3.7) $$v_{rel}$$ ; at (0,-7.5) $$Wind$$ $$flow$$, $$U_\infty$$ ;

(v21) at (6,-0.4); (v22) at (2,-0.4); (v21) edge (v22); at (4,-1) $$rotation$$ $$direction$$ ; at (4,0.3) $$rotation$$ $$plane$$ ; at (6.2,3.1) $$chord$$ $$reference$$ $$line$$ ; at (-0.5,3.5) $$rotation$$ $$axis$$ ;

(v21) at (7,6) ; (v22) at (-6,7) ; (v1) edge (v21); (v1) edge (v22); at (-4.5,6.5) $$L$$ ; at (5,5) $$D$$ ;

Basics of Computational Fluid Dynamics
Computational Fluid Dynamics (CFD) is the science of numerically predicting fluid flow, heat and mass transfer, chemical reactions, and related phenomena. CFD solves primarily equations for conservation of mass, momentum and energy. These governing equations are formulated in integral form and can solely be calculated iteratively, employing the finite volume method (FVM) as the main spatial discretisation approach.

Navier-Stokes Conservation Equations
The conservation equations in the fluid mechanics are derived from the principles of mass, momentum and energy conservation in fluid volumes. Defining a control volume as shown in Figure [fig:conv], the rate of change of every intensive fluid quantity, $$\phi$$, per unit of time (transient term) results in the rate of the quantity which trespasses the surface of the control volume (convective term), the magnitude diffusing rate (diffusive term) and the variable volumetric sink and sources (generation term). The representation of the convective flux in a control volume can be observed in Figure [fig:conv].

,[width = 0.7]

(O) at (0.,0.); (1) at (-3,-0.9); (2) at (2.6,0.6); (3) at (-0.8,1.2); (4) at (0.2,-1.4); (upp) at (-0.5,3); (btt) at (0.5,-2.5);

(1) to[out=-29,in=190](4); (4) to[out=5,in=270](2);

(2) to[out=125,in=-20](3); (3) to[out=175,in=80](1);

(1) to[out=85,in=210](upp); (upp) to[out=20,in=99](2);

(1) to[out=-85,in=200](btt); (btt) to[out=20,in=-85](2);

(r1) at (0.2,2.8); (r2) at (0.75,2.3); (l1) at (-0.7,2.5); (l2) at (-0.15,2.0);

(l1) to[out=20,in=-85](l2); (l1) to[out=20,in=180](r1); (r1) to[out=20,in=-85](r2); (l2) to[out=20,in=195](r2);

(r1u) at (0.2,5.3); (r2u) at (0.75,4.8); (l1u) at (-0.7,5); (l2u) at (-0.15,4.5);

(l1u) to[out=20,in=-85](l2u); (l1u) to[out=20,in=180](r1u); (r1u) to[out=20,in=-85](r2u); (l2u) to[out=20,in=195](r2u);

(r1)edge(r1u); (r2)edge(r2u); (l1)edge(l1u); (l2)edge(l2u);

(v1) at (0.1,2.4); (v2) at (-1.3,3.8);

(v1)edge(v2);

(n) at (-1.5,4.1) $$\vec{n}$$ ; (n) at (1.4,4) $$\vec{v}\partial{t}$$ ; (n) at (2.9,1.3) $$\xi_{c}$$ ; (xi) at (1.2,1.8) $$\partial{\xi}$$ ; (n) at (-0.2,-1) $$V_{c}$$ ;

(xi)edge(v1);

The so-called transport equation for whichever fluid magnitude in a closed control volume results in :

$$\begin{aligned} \underbrace{\frac{\partial{}}{\partial{t}}\int_{V_{c}}\rho\phi dV}_\text{transient term} + \underbrace{\oint_{\xi_{c}}\rho\phi\vec{v}\vec{n}d\xi}_\text{convective term} = \underbrace{\oint_{\xi_{c}}\rho\Gamma\nabla{\phi}d\xi}_\text{diffusive term} + \underbrace{\int_{V_{c}}q_{\phi}dV}_\text{generation term}\end{aligned}$$ Where $$\Gamma$$ is the diffusive constant and $$q_{\phi}$$ the generation constant.

By inserting certain values for the intensive fluid magnitudes $$\phi$$, the different conservation equations are found:

\frac{\partial{}}{\partial{t}} \int_{V_{c}}\rho dV + \oint_{\xi_{c}} \rho \vec{v} \vec {n} d\xi = 0 \end{aligned}$$ For a steady and incompressible flow it results in the divergence of the flow velocity field: $$\begin{aligned} \nabla*\vec{v} = 0 \end{aligned}$$ \frac{\partial{}}{\partial{t}} \int_{V_{c}}\rho \vec{v} dV + \oint_{\xi_{c}} (\rho \vec{v}) \vec{v} \vec {n} d\xi = \underbrace{\oint_{\xi_{c}} \bar{\bar{\sigma}}\vec{n}d\xi}_\text{surface forces} + \underbrace{\int_{V_{c}} \rho \vec{a} d V}_\text{volumetric forces} \end{aligned}$$ In which $$\bar{\bar{\sigma}}$$ denotes the stress tensor, whereas $$\vec{a}$$ refers to the vector of the volumetric forces, such as the gravity, acting over the control volume. \frac{\partial{}}{\partial{t}} \int_{V_{c}}\rho T dV + \oint_{\xi_{c}} (\rho T) \vec{v} \vec {n} d\xi = \underbrace{\oint_{\xi_{c}} (\vec{p} + (\frac{\kappa}{c_{p}}\nabla T))\vec{n}d\xi}_\text{convective and conductive terms} \end{aligned}$$
 * For $$\phi = 1$$ it gives the continuity equation : $$\begin{aligned}
 * For $$\phi = \vec{v} $$ it provides the conservation of momentum equation in vectorial form : $$\begin{aligned}
 * For $$\phi = h$$ - the enthalpy $$h = c_{p}T$$ - an example for the energy equation for heat transfer in solids is obtained : $$\begin{aligned}

Wind Turbine Controller Design State-of-the-Art
The wind turbine controller system is essentially composed of a group of sensors and actuators which are combined, and controlled via a software-based system, to function at a particular desired operating point. The principal goal of a wind turbine control system is to extract the maximum power maintaining certain structural and operative limits and for variable wind and rotational speeds. Furthermore, another distinct objective of the control system is to reduce fatigue, stresses and vibrations due to mechanical loads which act on the wind turbine during operation.

The aerodynamic model of the turbine is required so as to obtain the appropriate values for torque and wind thrust. The simplified closed-loop system model, also defined as feedback control system, in frequency response is represented in Figure [fig:cloop].



In which $$G_{c}(s)$$ represents the controller transfer function, the $$G(s)$$ stays for the linearised plant model and $$H(s)$$ is the feedback transfer function, usually $$1$$. Besides, $$Y(s)$$ is Laplace transform of the measured output, $$R(s)$$ the reference point, $$E(s)$$ the error signal and $$U(s)$$ the input signal (the angle of attack-pitch value to control). The closed-loop $$T(s)$$ transfer function which shows the ratio between the reference input and the output yields: $$\begin{aligned} T(s) = \frac{Y(s)}{R(s)} = \frac{G_{c}(s)G(s)}{(1 + H(s)G_{c}(s)G(s))}\end{aligned}$$ An optimal PID controller is of interest for pitch control systems. An ideal PID controller transfer function is provided by: $$\begin{aligned} G_{c}(s) = K_{p}(1 + \frac{1}{T_{i}s} + T_{d}s)\end{aligned}$$ In practice, the pure derivative action is never used, due to the “derivative kick” produced in the control signal for a step input, and to the undesirable noise amplification. It is usually replaced by a first-order low pass filter. Thus, the Laplace transformation representation of the approximate PID controller is given by : $$\begin{aligned} G_{c}(s) = K_{p}(1 + \frac{1}{T_{i}s} + \frac{T_{d}s}{1 + \tau T_{d} s})\end{aligned}$$ Where $$K_{p}$$ refers to the proportional controller gain, $$T_{i}$$ denotes the integrative compensator and $$T_{d}$$ stays as the derivative term. In addition, $$\tau$$ is the filter term.

Adding further terms to the control loop can improve the behaviour of the feedback control system, such as notch filters or lead-lag compensators. This will certainly improve its operation around undesired scopes as high frequency modes of vibration.

= System Overview =

The coupled system developed is composed of three distinct parts. The CFD plant model of the turbine in ANSYS-FLUENT, which provides the system aerodynamic loads; the MATLAB-SIMULINK model representing the operation behaviour of the tower dynamics, giving information regarding the tower for-aft oscillations and of the rotor angular speed when needed, and the pitch angle control system designed for the linearised rotor-tower model. All together comprise the desired closed-loop system which actively damps the system oscillations. In Figure [fig:sys_rep] a representation of the overall system can be observed.



In the complete coupled system, the CFD Turbine Plant obtains the wind thrust force $$T$$ and the generated turbine torque $$M_{t}$$. This information is employed by the rotor-tower model (or the FEA) in SIMULINK to calculate the resulting tower for-aft oscillation, $$v_{tower}$$ and the actual rotor speed $$\omega_{rotor}$$ (solely for the full rotor-tower system). The controller then obtains the required angle of attack to compensate the tower oscillation (by means of using the upper node tower velocity, $$v_{tower}$$, as the measured output), feeding it back as blade pitch angle to the ANSYS-FLUENT simulation and therefore closing the loop.

The wind speed, $$U_{\infty}$$, starting turbine position, the oscillating/wave forces, $$F_{w}$$, the electrical power extracted, $$P_{electric}$$ and the optimal tip speed ratio, $$\lambda_{opt}$$ are given by the nominal turbine operation conditions. In Figure [fig:for_aft_tower] a representation of the tower for-aft oscillations is shown.

(O) at (0,0) ; (upp) at (0,3.5) ; (btt) at (0,-2.6) ;

(O) to[out=60,in=265](upp); (upp) to[out=255,in=110] (O); (btt) to[out=80,in=285](O); (O) to[out=255,in=110] (btt);

(0,0) ellipse (0.4cm and 0.15cm)[draw=black, fill = gray!15!white];

(a) [rectangle, fill = gray!15!white, shape border rotate=0, shift = (0.1,0) ,draw, minimum height=1.2mm, minimum width=2.2mm] ;

(a) [rectangle, fill = gray!15!white, shape border rotate=0, shift = (0.3,0) ,draw, minimum height=1.7mm, minimum width=3.7mm] ; (tu1) at (0.6,-0.2); (tu2) at (1,-0.2); (tb1) at (0.45,-6.5); (tb2) at (1.15,-6.5);

(tu1)edge(tb1); (tu2)edge(tb2); (tb1)edge(tb2);

(b1) at (-0.55,-6.5); (b2) at (2.15,-6.5); (b1)edge(b2);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(O) at (0,0) ; (upp) at (0,3.5) ; (btt) at (0,-2.6) ;

(O) to[out=60,in=265](upp); (upp) to[out=255,in=110] (O); (btt) to[out=80,in=285](O); (O) to[out=255,in=110] (btt);

(0,0) ellipse (0.4cm and 0.15cm)[draw=black, densely dashed, fill = gray!0!white];

(tu1) at (0.6,-0.2); (tu2) at (1,-0.2);

(a) [rectangle, fill = gray!0!white, shape border rotate=0, densely dashed, rotate = (22), shift = (-0.90,0.102) , draw, minimum height=1.2mm, minimum width=2.2mm] ; (a) [rectangle, fill = gray!0!white, shape border rotate=0, densely dashed, rotate = (22) , shift = (-0.65,0.1) , draw, minimum height=1.7mm, minimum width=3.7mm] ;

(tu1) to[out=-60,in=105](tb1); (tu2) to[out=-60,in=105](tb2);

(O) at (0,0) ; (upp) at (0,3.5) ; (btt) at (0,-2.6) ;

(O) to[out=60,in=265](upp); (upp) to[out=255,in=110] (O); (btt) to[out=80,in=285](O); (O) to[out=255,in=110] (btt);

(0,0) ellipse (0.4cm and 0.15cm)[draw=black, densely dashed, fill = gray!0!white];

(tu1) at (0.6,-0.2); (tu2) at (1,-0.2);

(a) [rectangle, fill = gray!0!white, shape border rotate=0, densely dashed, rotate = (-15), shift = (1.1,0.185) , draw, minimum height=1.2mm, minimum width=2.2mm] ; (a) [rectangle, fill = gray!0!white, shape border rotate=0, densely dashed, rotate = (-15) , shift = (1.35,0.195) , draw, minimum height=1.7mm, minimum width=3.7mm] ;

(tu1) to[out=260,in=85](tb1); (tu2) to[out=260,in=85](tb2);

Wind Turbine Characteristics
The wind power plant type selected is a 5.XM series Senvion offshore three-blade wind turbine with a nominal power extraction of $$5.075\ \mathrm{MW}$$. Its design is made in CATIA. Its main operation characteristics are listed in the following Table [tab:turbine_param] extracted from.

ccc

Property &amp; Value &amp; Unit

Nominal Power &amp; $$5.075$$ &amp; $$\mathrm{MW}$$

Nominal Wind Speed &amp; $$14$$ &amp; $$\mathrm{m/s}$$

Cut-in Wind Speed &amp; $$3.5$$ &amp; $$\mathrm{m/s}$$

Cut-off Wind Speed &amp; $$30$$ &amp; $$\mathrm{m/s}$$

Nominal Rotational Speed &amp; $$12.1$$ &amp; $$\mathrm{rev/min}$$

Nominal Angular Speed &amp; $$1.267$$ &amp; $$\mathrm{rad/s}$$

Nominal Tip Speed Ratio &amp; $$5.7$$ &amp; $$\mathrm{-}$$

[tab:turbine_param]

The wind turbine is fundamentally composed of four principal components, namely the three-blade rotor, the nacelle, the spinner and the supporting - simple mono-pile type - tower structure. No other component is defined so as to simplify the turbine design and consequently the final CFD simulation. The Wortmann FX63-137 airfoil type is selected. Its aerodynamic and geometric properties are shown in Table [tab:aproper].

ccc

Characteristic &amp;Value &amp; Unit

Maximum Thickness/Chord Length &amp; $$0.137$$ &amp; $$\mathrm{-}$$

Position of Maximum Thickness/Chord Length &amp; $$0.297$$ &amp; $$\mathrm{-}$$

Maximum Lift Coefficient &amp; $$2.01$$ &amp; $$\mathrm{-}$$

Angle of attack at maximum lift coefficient &amp; $$10.5$$ &amp; $$\mathrm{deg}$$

Drag Coefficient at maximum lift coefficient &amp; $$0.05$$ &amp; $$\mathrm{-}$$

Aerodynamic Efficiency &amp; $$91.7$$ &amp; $$\mathrm{-}$$

Zero Lift Angle &amp; $$-9$$ &amp; $$\mathrm{deg}$$

[tab:aproper]

Then, the rotor blade is constructed in order to maintain the maximum lift coefficient, $$c_{L,max}$$, by means of keeping $$\alpha_{c_{L,max}}$$ fixed for nominal operation conditions along the whole blade. The Schmidt theorem is employed for the blade angular and chord distribution design, according to. Its properties shown in Table [tab:blade_prop].

ccc

Property &amp; Value &amp; Unit

Net Blade Length &amp; $$58$$ &amp; $$\mathrm{m}$$

Starting Twist Angle &amp; $$33.2726$$ &amp; $$\mathrm{deg}$$

Tip Twist Angle &amp; $$-3.829$$ &amp; $$\mathrm{deg}$$

Starting Chord Length &amp; $$5.7916$$ &amp; $$\mathrm{m}$$

Average Chord Length &amp;$$3.5362$$ &amp; $$\mathrm{m}$$

Tip Chord Length &amp;$$1.1758$$ &amp; $$\mathrm{m}$$

Maximum Chord Length &amp;$$6.394$$ &amp; $$\mathrm{m}$$

Radial Position at Maximum Chord Length &amp;$$8.19$$ &amp; $$\mathrm{m}$$

[tab:blade_prop]

The spinner and the nacelle are also modelled in CATIA as simple geometries. The supporting structure properties are provided by. For this case, the tower is not modelled as a CAD design for the CFD simulation, nonetheless solely its dynamic behaviour.

Rotor &amp; Tower Mathematical Model
Having defined the system stiffness $$k$$ and the tower first mode of vibration $$\omega_{r}$$, the tower is modelled as a one degree of freedom mass-spring-damper system excited by the impact of an oscillating force, which represents the vibrational forces acting on the turbine during operation. It essentially defines the complete tower as a single-degree of freedom linear spring (as the aerodynamic thrust behaves as the damper) and leaves aside the other higher vibration modes of the flexible tower structure. Figure [fig:msd] shows a representation of the mass-spring-damper system.

(a) at (2.5,-1.5); (b) at (2.5,1.5); (c) at (-2.5,-1.5); (d) at (-2.5,1.5);

(a)edge(b); (b)edge(d); (c)edge(d); (c)edge(a);

(s1) at (-5.5,0.5); (s2) at (-2.5,0.5);

(s3) at (-5.0,0.9); (s4) at (-4.75,0.1); (s5) at (-4.5,0.9); (s6) at (-4.25,0.1); (s7) at (-4.0,0.9); (s8) at (-3.75,0.1); (s9) at (-4.0,0.9); (s10) at (-3.75,0.1); (s11) at (-3.5,0.9); (s12) at (-3.25,0.1); (s13) at (-3.0,0.9); (s14) at (-2.75,0.1);

(s1)edge(s3); (s3)edge(s4); (s4)edge(s5); (s5)edge(s6); (s6)edge(s7); (s7)edge(s8); (s8)edge(s9); (s9)edge(s10); (s10)edge(s11); (s11)edge(s12); (s12)edge(s13); (s13)edge(s14); (s14)edge(s2);

(d1) at (-5.5,-0.6); (d2) at (-2.5,-0.6);

(d3) at (-4.5,-0.25); (d4) at (-3.5,-0.25); (d5) at (-4.5,-0.95); (d6) at (-3.5,-0.95); (d7) at (-4.2,-0.25); (d8) at (-4.2,-0.95);

(d9) at (-3.5,-0.60); (d10) at (-4.2,-0.60);

(d3)edge(d4); (d5)edge(d6); (d4)edge(d6); (d7)edge(d8);

(d1)edge(d10); (d2)edge(d9);

(f1) at (2.5,0.2); (f2) at (5.5,0.2);

(f1)edge(f2);

(x1) at (0,2.2); (x2) at (3.5,2.2); (x3) at (0,2.6); (x4) at (0,1.8);

(x1)edge(x2); (x3)edge(x4);

(s) at (-4,1.4) $$k$$ ; (d) at (-4,-1.4) $$c_{d}$$ ; (d) at (0,0.2) $$m$$ ; (f) at (4,0.65) $$F(t)$$ ; (x) at (2,2.75) $$x$$,$$\dot{x}$$,$$\ddot{x}$$ ;

(b1) at (-0.55,-6.5); (b2) at (2.15,-6.5); (b1)edge(b2);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

System which has the following characteristic differential equation: $$\begin{aligned} m\ddot{x} + c_{d}\dot{x} + kx = F(t)\end{aligned}$$  $$\begin{aligned} m\ddot{x} + kx = \frac{1}{2}\rho \pi R^{2}(U_{\infty} - \dot{x})^2c_{T} + F_{W}\end{aligned}$$ The excitation force $$F(t)$$ s composed of $$F_{W}$$, the waves excitation force and $$T$$, the thrust force. The former is defined as an sinusoidal load which oscillates at the resonance system frequency $$\omega_{r}$$ and with amplitude given by a constant plus the system oscillating mass. Therefore: $$\begin{aligned} F(t) = F_{W}(t) = m_{osc}\hat{F} \sin{w_{r}t}\end{aligned}$$ Together with the application of the conservation of angular momentum to the rotor: $$\begin{aligned} I \dot{\omega} = M_{T} - M_{e}\end{aligned}$$ $$\begin{aligned} M_{T} = \frac{1}{2}\rho \pi R^{3} (U_{\infty} - \dot{x})^{2} c_{M}\end{aligned}$$  $$\begin{aligned} M_{e} = \frac{c_{P,opt}*\rho \pi R^{5}}{2 \lambda_{opt}^{3}}\omega^{2}\end{aligned}$$ Finally, the coupled system of non-linear differential equations which models the turbine-tower performance is given by, with the introduction of the constant factor $$q_{cte}$$: $$\begin{aligned} \ddot{x} &=& \frac{1}{m} \left(\frac{1}{2}\rho \pi R^{2}(U_{\infty} - \dot{x})^2c_{T} + F_{w}\right)\\ \dot{\omega} &=& \frac{1}{I} \left(\frac{1}{2}\rho \pi R^{3} (U_{\infty} - \dot{x})^{2} c_{M} - q_{cte}\omega^{2}\right)\end{aligned}$$ Thus, Figure [fig:matmodeltower] represents the forces and torques which are acting on the wind turbine and which are defined in our mathematical model.

(O) at (0,0) ; (upp) at (0,3.5) ; (btt) at (0,-2.6) ;

(O) to[out=60,in=265](upp); (upp) to[out=255,in=110] (O); (btt) to[out=80,in=285](O); (O) to[out=255,in=110] (btt);

(0,0) ellipse (0.4cm and 0.15cm)[draw=black, fill = gray!15!white];

(a) [rectangle, fill = gray!15!white, shape border rotate=0, shift = (0.13,0) ,draw, minimum height=1.11mm, minimum width=2.3mm] ;

(a) [rectangle, fill = gray!15!white, shape border rotate=0, shift = (0.3,0) ,draw, minimum height=1.72mm, minimum width=3.37mm] ;

(massr) at (2.3,2.4) $$m_{osc}$$ ; (massc) at (2.0,2.1);

(massc)edge(0.8,0); (0) at (0.8,0) ; (tu1) at (0.6,-0.2); (tu2) at (1,-0.2); (tb1) at (0.45,-6.5); (tb2) at (1.15,-6.5);

(tui1) at (0.7,-0.2); (tui2) at (0.9,-0.2); (tbi1) at (0.55,-6.5); (tbi2) at (1.05,-6.5);

(tu1)edge(tb1); (tu2)edge(tb2); (tb1)edge(tb2);

(tui1)edge(tbi1); (tui2)edge(tbi2);

(b1) at (-0.55,-6.5); (b2) at (2.15,-6.5); (b1)edge(b2);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(b3) at (-0.55,-6.65); (b4) at (-0.45,-6.5); (b3)edge(b4);

(H1) at (1.55,-6.45); (H2) at (1.55,-0.2); (h1) at (2.55,-6.45); (h2) at (2.55,-5.0);

(H1)edge(H2); (h1)edge(h2);

(H) at (2,-3.25) $$H$$ ; (h) at (2.90,-5.75) $$h$$ ; (T1) at (1.3,.2) ; (T2) at (3.5,0.2) ; (Fw1) at (1.3,-0.2) ; (Fw2) at (3.5,-0.2) ;

(T1)edge(T2); (Fw1)edge(Fw2);

(T1) at (1.3,0.2) ; (T2) at (3,0.2) ; (Fw1) at (1.3,-0.2) ; (Fw2) at (3,-0.2) ;

(k) at (2.4,0.7) $$T$$ ; (r) at (2.4,-0.7) $$F_{W}$$ ;

(Mt1) at (-0.4,.2) ; (Mt2) at (-2.9,0.2) ; (Mt22) at (-3,0.2) ; (Me1) at (-0.4,-0.2) ; (Me2) at (-3,-0.2) ; (Me12) at (-0.5,-0.2) ;

(d) at (-1.9,0.5) $$M_{t}$$ ; (b) at (-1.9,-0.5) $$M_{e}$$ ;

(Mt1)edge(Mt2); (Mt1)edge(Mt22); (Me2)edge(Me1); (Me2)edge(Me12);

(-2.5,-5.0) sin (-2,-4.95) cos (-1.5,-5.0) sin (-1,-5.05) cos (-0.5,-5.0)sin (0,-4.95) cos (0.5,-5.0)sin (1,-5.05) cos (1.5,-5.0)sin (2,-4.95) cos (2.5,-5.0)sin (3,-5.05) cos (3.5,-5.0); (w1) at (-0.5,-5.1) ; (w2) at (-0.5,-5.3) ; (w3) at (-0.5,-5.5) ; (w4) at (-0.5,-5.7) ; (w5) at (-0.5,-5.9) ; (w6) at (-0.5,-6.1) ; (w7) at (-0.5,-6.3) ; (wf1) at (0.37,-5.1) ; (wf2) at (0.3,-5.3) ; (wf3) at (0.23,-5.5) ; (wf4) at (0.16,-5.7) ; (wf5) at (0.07,-5.9) ; (wf6) at (-0.03,-6.1) ; (wf7) at (-0.15,-6.3) ;

(w1)edge(wf1); (w2)edge(wf2); (w3)edge(wf3); (w4)edge(wf4); (w5)edge(wf5); (w6)edge(wf6); (w7)edge(wf7);

(wf8) at (-0.16,-6.38) ; (wf0) at (0.4,-5.05) ; (wf0) to[out=255,in=65](wf8);

(w1) at (-2.5,-5.6) $$F_{water}$$ ;

= Controller Design =

System State-Space Representation
Prior to system linearisation, the non-linear model defined in by Equations (29) and (30), is comprised into a sole state-space model with the following state vector, $$\vec{x}$$, the input vector, $$\vec{u}$$, and the output vector, $$\vec{y}$$:

$$\begin{aligned} \vec{x} = \begin{bmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{4} \end{bmatrix} = \begin{bmatrix} x\\ \theta\\ \dot{x}\\ \dot{\omega} \end{bmatrix}\end{aligned}$$ $$\begin{aligned} \vec{u} = \begin{bmatrix} u_{1}\\ u_{2} \end{bmatrix} = \begin{bmatrix} \alpha\\ F_{W} \end{bmatrix}\end{aligned}$$

$$\begin{aligned} \vec{y} = \begin{bmatrix} y_{1}\\ y_{2} \end{bmatrix} = \begin{bmatrix} \dot{x}\\ \dot{\omega} \end{bmatrix}\end{aligned}$$

where $$x$$ is the for-aft tower displacement, $$\dot{x}$$ represents the for-aft tower velocity, $$\theta$$ the azimuthal radial position of the rotor and $$\omega$$ being the angular speed of the rotor blades. $$\alpha$$ is the system s controllable variable and the $$F_{W}$$ the system s uncontrollable input variable which varies sinusoidally at the first tower eigenfrequency.

Thus, the coupled differential system is rearranged, represented in the standard form of the state derivatives and finally has the following state evolution equations: $$\begin{aligned} \dot{\vec{x}} = \begin{bmatrix} x_{3}\\ x_{4}\\ -\frac{k}{m}x_{1} + \frac{\rho \pi R^{2}}{2m}(U_{\infty} - x_{3})^{2}c_{T} + \frac{F_{W}}{m}\\ \frac{1}{I}(\frac{1}{2}\rho \pi R^{3}(U_{\infty} - {x_{3}})^{2}c_{M} - \frac{c_{P,opt}\rho \pi R^{5}}{2 \lambda_{opt}^{3}}x_{4}^{2}) \end{bmatrix}\end{aligned}$$ Once the state formulation is defined, the next step is to obtain the linearised model.

Linearisation
Prior to linearisation, it is first necessary to set a working point, or linearisation point, for the state, $$\vec{x}$$, input, $$\vec{u}$$, and output vectors, $$\vec{y}$$. Therefore, given the generalised state and output non-linear functions are: $$\begin{aligned} \dot{\vec{x}} = \vec{f}(\vec{x},\vec{u},t)\\ \vec{y} = \vec{h}(\vec{x},\vec{u},t) \end{aligned}$$ Applying Taylor series expansion, the matrices of this linear state-space model are given by the partial derivatives of the system and output equations with respect to the state and inputs, yielding: $$\begin{aligned} \dot{\vec{x}} = \vec{f}(\vec{x}_{0},\vec{u}_{0},t) + \underbrace{\frac{\partial \vec{f}}{\partial \vec{x}}\Big|_{\vec{x}_{0},\vec{u}_{0},t}}_\text{A} \varDelta \vec{x} +\underbrace{\frac{\partial \vec{f}}{\partial \vec{u}}\Big|_{\vec{x}_{0},\vec{u}_{0},t}}_\text{B} \varDelta \vec{u}\\ \dot{\vec{y}} = \vec{h}(\vec{x}_{0},\vec{u}_{0},t) + \underbrace{\frac{\partial \vec{h}}{\partial \vec{x}}\Big|_{\vec{x}_{0},\vec{u}_{0},t}}_\text{C} \varDelta \vec{x} +\underbrace{\frac{\partial \vec{h}}{\partial \vec{u}}\Big|_{\vec{x}_{0},\vec{u}_{0},t}}_\text{D} \varDelta \vec{u} \end{aligned}$$ It yields a Multiple Input Multiple Output system (MIMO). The system comprises two controlled variables, the inputs (one controllable) and two manipulated variables, the outputs, and therefore the system has to show the interaction between each input variable with the correspondent manipulated variable. Thus, the system transfer function $$G(s)$$ results in : $$\begin{aligned} G(s) = \begin{bmatrix} G_{v_{t},\alpha}(s) & G_{v_{t},F_{W}}(s) \\ G_{w_{rotor},\alpha}(s) & G_{w_{rotor},F_{W}}(s) \end{bmatrix}\end{aligned}$$

Close-Loop Model Transfer Function
Adjusting the aerodynamic rotor damping accordingly is achieved by controlling the blades pitch angle $$\varphi$$ from its original position in order to reach the appropriate angle of attack, $$\alpha$$, which provides the required aerodynamic damping force. Therefore, the system is simplified into a Multiple Input Single Output (MISO), as there is no intention to control nor manipulate the rotor angular speed, $$\omega_{rotor}$$. See Figure [fig:cloop_sys].

At this point, the closed-loop transfer function of the system can be defined. First, the output-inputs plant relation results: $$\begin{aligned} Y_{1}(s) = G_{v_{t},\alpha}(s)U_{\alpha}(s) + G_{v_{t},F_{W}}(s)U_{FW}(s)\end{aligned}$$ Where the closed-loop transfer function, $$T(s)$$, giving the relation between the wave forces input and the tower velocity output, yields: $$\begin{aligned} T(s) = \frac{Y_{\alpha}(s)}{U_{FW}(s)} = \frac{G_{v_{t},F_{W}}(s)}{1 + G_{c}(s)G_{v_{t},\alpha}(s)} \end{aligned}$$

PID Controller
The general ideal PID transfer function, Equation (14), is inserted into the closed-loop TF and its characteristic equation is obtained: $$\begin{aligned} \phi(e) = s^{3}(1 - 0.4311K_{p}T_{d}) +s^{2}(1.07 - 0.4311K_{p}-\nonumber\\ -0.2216K_{p}T_{d}) + s(1.838 - 0.2216K_{p} - 0.4311\frac{K_{p}}{T_{i}})\nonumber\\ + (0.8598 - 0.2216\frac{K_{p}}{T_{i}})\end{aligned}$$ The resulting third-order system is then forced to follow a PT2-PT1 behaviour, which : $$\begin{aligned} G_{PT2*PT1} = K(s^{2} + 2\xi \omega s + \omega^{2})(s + \omega)\end{aligned}$$ Introducing values for the damping ratio $$\xi$$ and a desired settling time $$t_{s}$$, knowing that its relation is given by: $$\begin{aligned} t_{s} = \frac{4}{\xi \omega},\end{aligned}$$ then different PID factors can be obtained. With the objective that the system should respond rapidly to the input signal given by the wave forces, relatively high damping ratios and low settling times are desired. For $$t_{s} = 0.4$$ and $$\xi = 0.7$$, the PID compensator terms result in: $$K_{p} = -80.2433$$, $$T_{i} = 0.0694$$, $$T_{d} =  0.0029$$. Thus, the linearised system response and control variable are depicted in Figure [fig:lin_0.4]:



Implementation of the Controller into the System Mathematical Model
The chosen PID controller is therefore introduced into the tower-rotor mathematical model. Figure [fig:nonlinearmodel].



Essentially, the non-linear plant defining the dynamic behaviour of the wind turbine calculates the for-aft displacement, $$x_{tower}$$ and speed $$v_{tower}$$, and the rotor angular position $$\theta$$ and velocity $$\omega_{r}$$ depending on the values of wind speed $$U_{\infty}$$, optimal tip speed ratio $$\lambda_{opt}$$, optimal power coefficient $$c_{P,opt}$$, thrust and torque coefficients, $$c_{T}$$ and $$c_{M}$$ respectively, and finally the excitation wave force $$F_{W}$$. The former three input values are constant as are provided by the turbine nominal conditions, Table [tab:turbine_param]. On the other hand, the rotor coefficients depend on the $$\alpha$$ and on the $$\lambda_{new}$$, which directly depends on the system outputs tower for-aft speed and rotor angular velocity. Its formula is given by: $$\begin{aligned} \lambda_{new} = \frac{\omega_{rotor}R}{v_{real}}\end{aligned}$$ Being $$v_{real}$$, the real or effective wind speed. It is basically given by the difference between the undisturbed wind speed $$U_{\infty}$$ and the tower for-aft speed $$v_{tower}$$: $$\begin{aligned} v_{real} = U_{\infty} - v_{tower}\end{aligned}$$ The controller, on the other hand, reads the information given by the measured output, $$\Delta v_{tower}$$, and provides the according angle of attack, $$\Delta \alpha$$, which compensates the tower oscillations. This $$\Delta \alpha$$ is then sum up with the starting optimal angle of attack, $$\alpha_{c_{L,max}}$$, Table [tab:aproper]. The combination of the new tip speed ratio together with the new angle of attack provide the new according turbine torque and thrust coefficients.

= 3D Wind Turbine CFD Setup =

The computational mesh has been detached into 5 separated grid zones. The stream grid and the disc grid, which gathers the bulk of the turbine geometry (rotor blades plus spinner). On the other hand, in order to introduce the influence of the pitch blade adjustment, three more regions are created within the disc grid. In Figures [fig:grid_setup] and [fig:PlaneYZ], a representation of the system grids defined XY and YZ planes is shown.

(airfoill) at (0,0) ;

(O) at (0.0,0.0); (y) at (0.0, 6); (z) at (-6.1, 0); at (1.1, 6.1) $$Y_{0}$$ ; at (-6.8,-1.3) $$X_{0}$$ ;

(O) – (z) node [near end] (0,0) arc (-150:150:1 and 1);  ; at (-5, 3.5) $$\omega_{rotor}$$ ;

(a) at (0.2,-1.105) ; (a) at (0.2,-1.105) ; (c1) at (0.05,-1.25); (c3)at (0.35,-1.25); (c4) at (0.05,-1); (c2) at (0.35,-1);

at (1.8,-0.5) $$Z_{0}$$ ;

(O) edge (y); (O) edge (z);

(xx1) at (-11,0); (xx2) at (25,0); (xx1) edge (xx2);

(xxx1) at (-8.3,7); (xxx2) at (-3,7); (xxx1) edge (xxx2);

(xxx2) at (-4.5,8.1) $$Flow\ Direction$$ ;

(d1) at (0.9,5.0); (d2) at (0.9,-5.0); (d3) at (-0.9,5.0); (d4) at (-0.9,-5.0);

(d1) edge (d2); (d2) edge (d4); (d3) edge (d4); (d3) edge (d1);

(DISC) at (6,-3.5) $$Disc\ Grid$$ ; (a) at (2.8,-3.5); (b) at (0.2,-2.5); (a) edge (b);

(cv1) at (23.0,11.0); (cv2) at (23.0,-11.0); (cv3) at (-10.0,11.0); (cv4) at (-10.0,-11.0);

(cv1) edge (cv2); (cv2) edge (cv4); (cv3) edge (cv4); (cv3) edge (cv1);

(CV) at (15.9,8.0) $$Stream\ Grid$$ ;

(cota1) at (-11.5,11.0); (cota2) at (-11.5,-11.0);

(cota1) edge (cota2);

(cota3) at (-10.0,12.0); (cota4) at (23,12.0);

(cota3) edge (cota4);

(cota5) at (-1.3,5.0); (cota6) at (-1.3,-5.0);

(cota5) edge (cota6);

(ccc1) at (-13,0.0) $$4D$$ ; (ccc2) at (7.0,13) $$4D$$ ; (ccc3) at (-2.1,2.8) $$D$$ ;

(v1) at (-5,-12.5); (v2) at (5,-12.5); (v1) edge (v2); (ccc1) at (0,-13.8) $$v_{tower}$$ ;

(airfoill) at (0.00,-0.05) ;

(O) at (0.03,-1.105); (y) at (0.03, 5); (z) at (-6.1, -1.105); at (0.6, 5) $$Y_{0}$$ ; at (-5.4,-0.7) $$Z_{0}$$ ; (a) at (0.2,-1.105) ; (c1) at (0.05,-1.25); (c3)at (0.35,-1.25); (c4) at (0.05,-1); (c2) at (0.35,-1); (c1) edge (c2); (c3) edge (c4); at (1,-0.5) $$X_{0}$$ ;

(O) edge (y); (O) edge (z);

(g1di) at (-0.3,-0.85); (g1dd) at (0.4,-0.85); (g1ui) at (-0.3,3.6); (g1ud) at (0.4,3.6);

(g1di) edge (g1ui); (g1dd) edge (g1ud); (g1di) edge (g1dd); (g1ud) edge (g1ui);

(g1di) at (-0.3,-0.85); (g1dd) at (0.4,-0.85); (g1ui) at (-0.3,3.6); (g1ud) at (0.4,3.6);

(g1di) edge (g1ui); (g1dd) edge (g1ud); (g1di) edge (g1dd); (g1ud) edge (g1ui);

(g1di) at (-0.3,-0.85); (g1dd) at (0.4,-0.85); (g1ui) at (-0.3,3.6); (g1ud) at (0.4,3.6);

(g1di) edge (g1ui); (g1dd) edge (g1ud); (g1di) edge (g1dd); (g1ud) edge (g1ui);

(g1di) at (-0.3,-0.85); (g1dd) at (0.4,-0.85); (g1ui) at (-0.3,3.6); (g1ud) at (0.4,3.6); (O) at (0.03,-1.105); (y) at (0.03, 5); (O) – (y) node [near end] (0,0) arc (-150:150:1 and 1);  ; at (-1.5, 5.5) $$Y^{'}_{blade\ grid\ 1}$$ ; (g1di) edge (g1ui); (g1dd) edge (g1ud); (g1di) edge (g1dd); (g1ud) edge (g1ui);

(y) at (-1.8, 2.7) $$\dot{\varphi}$$ ;

(g1di) at (-0.3,-0.85); (g1dd) at (0.4,-0.85); (g1ui) at (-0.3,3.6); (g1ud) at (0.4,3.6); (O) at (0.03,-1.105); (y) at (0.03, 5); (O) – (y) node [near end] (0,0) arc (-150:150:1 and 1);  ; at (0.5, 5.5) $$Y^{'}_{blade\ grid\ 2}$$ ; (g1di) edge (g1ui); (g1dd) edge (g1ud); (g1di) edge (g1dd); (g1ud) edge (g1ui);

(y) at (1.3, 2.5) $$\dot{\varphi}$$ ;

(g1di) at (-0.3,-0.85); (g1dd) at (0.4,-0.85); (g1ui) at (-0.3,3.6); (g1ud) at (0.4,3.6); (O) at (0.03,-1.105); (y) at (0.03, 5); (O) – (y) node [near end] (0,0) arc (-150:150:1 and 1);  ;; at (1.5, 3.8) $$Y^{'}_{blade\ grid\ 3}$$ ; (g1di) edge (g1ui); (g1dd) edge (g1ud); (g1di) edge (g1dd); (g1ud) edge (g1ui);

(y) at (0.52, 4.4) $$\dot{\varphi}$$ ;

(4.3,2.0) arc[x radius=4.74cm, y radius =4.74cm, start angle=35, end angle=85]; at (3.3,3.9) $$\omega_{rotor}$$ ; (0.0,2.0) arc[x radius=3cm, y radius =2cm, start angle=98, end angle=140]; at (-1.0,2.2) $$\theta$$ ;

at (-2.5,5.5) $$Blade \ Grid \ 1 $$ ; (-2,4.7) – (-0.1,3)[]; at (-5.2,-5) $$Blade \  Grid \ 2 $$ ; (-5,-4.8) – (-4,-3.6)[]; at (4.5,-5.5) $$Blade \  Grid \ 3 $$ ; (4.5,-5.2) – (3.8,-3.6)[];

It is important to point out that the grid tetra elements are simplified and converted to polyhedra prior to the simulation start.

Overall, the present simulation is set as a transient with the Spalart-Allmaras turbulent model, with different rotating and moving domains, relatively to each other. The sliding mesh approach is employed defining a mesh motion with matching mesh interfaces. Regarding the time step size setting, $$2 \ \mathrm{deg}$$ are allowed per time step for the fastest rotational speed occurring. It is finally rounded to $$0.025\ \mathrm{s}$$.

SIMULINK Tool Implementation
The interface between both programs ANSYS-FLUENT and MATLAB-SIMULINK has been developed in and adapted for the present. Fundamentally, the CFD plant model for every case is composed of the following characteristics:


 * 4 general inputs directly to the FLUENT case: wind speed and X, Y &amp; Z wind speed vector components.
 * 3 UDF inputs: pitch angular, rotor angular and tower for-aft speed.
 * 2 outputs: rotor thrust and torque.

Furthermore, 5 User Defined Functions (UDF) are employed, using the ones developed by as the basis, in order to allow the different mesh motions to take place during the simulation.

= Results =

Wind Turbine Simulation
Finally, the simulation is carried out with the use of the SIMULINK Tool for nominal conditions without feedback. Figures [fig:velocitycontoursx] and [fig:velocitycontoursz] the velocity contours at $$10\ \mathrm{s}$$ of the simulated wind turbine, for the YZ and XY planes respectively. Furthermore, In Table [tab:sima] the force, torque and power obtained are shown.

ccc

Property &amp; Value &amp; Unit

Aerodynamic Thrust &amp; $$1.225$$ &amp; $$\mathrm{MN}$$

Aerodynamic Thrust without Nacelle nor Spinner &amp; $$1.220$$ &amp; $$\mathrm{MN}$$

Rotor Torque &amp; $$6.335$$ &amp; $$\mathrm{MNm}$$

Power &amp; $$8.027$$ &amp; $$\mathrm{MW}$$

[tab:sima]

System Response without Controller and with Wave Forces
The non-controlled system is tested with the excitation forces and its state responses are given by [fig:states_ncont_fw].

System Response with Controller and Wave Forces
The controlled system is tested with the excitation forces and its states are shown in Figure [fig:states_04], while the controlled variable is given in Figure [fig:dalpha_04].





When the rotor angular speed reaches figures within operational conditions (just before $$1.267\ \mathrm{rad/s}$$), the controller is switched on and starts compensating the system. The adjusting of the angle of attack with figures between $$+ 3.602$$ and $$- 3.057\ \mathrm{deg}$$, following the wave forces signal frequency, modifies the aerodynamic thrust so that the system is virtually capable of compensating its effects. Nevertheless, the system yet oscillates within a certain small range of $$+ 0.19$$ and $$- 0.17\ \mathrm{m/s}$$ in comparison with the excitation of $$\pm 3\ \mathrm{m/s}$$ the uncontrolled system undergoes in Figure [fig:states_ncont_fw].

In addition, when observing the control variable behaviour over time, Figure [fig:dalpha_04], the effect of the controller switch-on, generating an starting peak, can be seen. This performance brings forward the need for the controller variable limitation.

Coupled System with Mathematical Rotor-Tower Model
At this point, the closed-loop complete rotor-tower model is implemented in series with the turbine CFD simulation plant. The simulation time has been defined to $$5\ \mathrm{s}$$. Figure [fig:fluentt] depicts the simulation results for the system states in comparison to the results acquired for the same case without the interference of the CFD model.



It yields a non-welcoming outcome, as the CFD model of the turbine requires a great deal more time to reach the nominal rotor angular speed of $$\omega_{0}$$ when compared to the precise period for the closed-loop mathematical model on its own. In addition, problems have arise when coupling more complex cases, not being able of obtaining appropriate results apart from the simpler situations.

Comparison of Stress Loads at the Root of the Tower for the Controlled and Non-Controlled Cases
Essentially, the maximum bending stress is given by : $$\begin{aligned} \sigma_{b} = \frac{M_{bending}}{I}z\end{aligned}$$ In which $$z$$ is the outer diameter tower cross section at the root and $$I$$ is the moment of inertia of the tower along its revolution axis. When the tower is subjected to a sinusoidal loading over its structure, as such is the present case with the introduction of the wave forces, an alternating stress appears which limits are constraint within the maximum bending stress. Therefore this alternating stress yields : $$\begin{aligned} \sigma_{a} = 2\sigma_{b}\end{aligned}$$ Figure [fig:load3] compares the controlled and uncontrolled system alternating stress widths.



Thus, given an oscillating force with an alternating width of $$141\ \mathrm{kN}$$ at the tower top, the alternating stress width results in $$75.45\ \mathrm{MPa}$$ for the uncontrolled system, while the controlled model yields the low figure of $$21.83\ \mathrm{MPa}$$.

The reduction for this undesired case shows the great improvement the controller provides to the system.

= Conclusion =

A controller algorithm has been successfully developed being capable of damping the strong waves oscillations modelled, also showing compensated alternating bending stresses that the tower root undergoes. The improvement reached results clear-cut even for the worst case simulated.

However, the outcome has not been successfully achieved, as the CFD plant has behaved very sensitively to high rotation relative speeds and for-aft motions of the sliding mesh when implemented the system within the coupled model.

Further research and work has to be carried out in order to obtain a final properly working version of the coupled system. The introduction of far smaller time steps sizes, finer meshes, the implementation of other mesh motion methods to simulate the tower for-aft oscillation are solely different possibilities that have to be attempted to make the model functional and valid.

= Acknowledgment =

I would principally like to thank my supervisor Andreas Mehrle for his great support and guidance during the development of this thesis, contributing with his great knowledge to implement the different steps of the present.

Furthermore, I would like to give thanks to Mario Heis, Manuel Berger and Alexander Fuchs for their support in providing all the possible information and data they had available with respect to their respective theses.

I also want to acknowledge all of those who supported me during the realisation of the thesis. With special attention to Jannik Gade, for his patience and stylistic advises, and to Bernhard Hollaus, for the support and great hints he provided me during the controller design.

Without their aid, this work would have not even been possible to begin with.

Alejandro Secades Rodríguez is currently carrying out the Master Study in Mechatronics &amp; Smart Technologies at MCI Innsbruck/Austria. He is also working part-time for the same University in a distinct variety of projects with special focus on CFD, Aerodynamics and Control Engineering.