Optimal design of passive-adaptive pendulum tuned mass damper for the global vibration control of offshore wind turbines

This paper presents an optimal design procedure for a pendulum tuned mass damper (PTMD) to mitigate the global structural vibrations of offshore wind turbines (OWTs) in the fore – aft and side – side directions. The procedure is tested to the design of a PTMD to be applied to the 5-MW benchmark baseline monopile wind turbine proposed by the National Renewable Energy Lab (NREL). The computation of wind and wave spectra, as well as the evaluation of the hydrodynamic and aerodynamic loads, is conducted by using an in-house built MATLAB® routine working together with an ANSYS® 3-D finite element (FE) global model for evaluating the resultant peak displacement response at the OWT hub by a power spectral density (PSD) analysis. In order to validate the OWT FEM model, a result comparison is made with the NREL OpenFAST, finding good matches between the two codes. An in-house built genetic algorithm (GA) toolbox, coded in MATLAB®, is then used to optimally design the parameters of a PTMD with a simplified 2-degrees-of-freedom (2DOF) model. The chosen GA fitness function targets the minimization of the peak response of the primary structure as evaluated by the 2DOF model. The design parameters of the PTMD are the flexural rigidity and damping, the mass ratio and pendulum length. After the 3-D FE model of the OWT without any control device has been validated, and the PTMD has been optimized by the simplified 2DOF model, the performances of the PTMD are examined on a 3-D global FE OWT + PTMD model in ANSYS®.


| INTRODUCTION
In the last decade, offshore wind turbines (OWTs) became one of the most implemented technologies for renewable-power generation.High structural vibration levels occur in an OWT due to its height, slenderness and the heavy loads acting both on its above-water and its submerged parts.In addition to this, in order to improve the productivity of modern OWTs, support structures continue to increase in slenderness and size, and structural designers continuously work to put forward the size limits of this kind of structures by providing more efficient design solutions regarding the structural dynamics in the environment, the structural/material strength and the operating or vibration control systems.
][7][8][9][10][11][12] Passive control consists of adding one or more non-active devices to the main structure to absorb or transfer part of kinetic energy. 5A passive control system does not require an external power source 13 ; it is easy to implement and is widely used in practical civil engineering structures, 1 such as in wind turbine systems. 14me applications of passive control devices include the tuned mass damper (TMD), which is able to mitigate the structural vibration by transferring the energy between the main structure and an auxiliary/secondary mass.The TMD is composed of a mass-spring-dashpot attached to the primary structure 15 ; usually, it is tuned to the fundamental structure natural frequency, and it is able to reduce substantially the response associated to the first mode of vibration.The conventional design of the TMDs based on the natural frequencies may not be suitable for OWTs, where the complex dynamics of the coupled substructure-rotor assembly may imply non-negligible harmonic contributions to the response at frequencies that are higher than the controlled one (mode shape the TMD is tuned with).In addition, optimal values of the design parameters of the control device may change depending on the wind velocity in operational conditions. 6,16particular type of TMD is the so-called pendulum TMD (PTMD) in which the oscillating secondary mass is a pendulum, whose natural frequency is influenced by the length of the pendulum itself.PTMDs have been investigated and widely applied in high-rise slender buildings due to the simplicity of installation and maintenance. 17,18In comparison to the TMD, the pendulum length is an extra dimension used for tuning, enabling, for the same device, several optimal PTMD configurations.Guimarães et al 19 implemented an inverted PTMD to reduce the vibration of a simplified OWT model.Roffel et al 20 evaluate the performance of a non-linear three-dimensional PTMD (3D-PTMD) reducing the responses of a flexible tower by using finite element (FE) representation.Recently, Sun and Jahangiri 21 developed an analytical model of the National Renewable Energy Lab (NREL) monopile 5-MW OWT 22 coupled with a 3D-PTMD; they showed that the focus on unidirectional vibration attenuation is inadequate for real applications.
As mentioned above, the pendulum length can be considered an additional design/tuning parameter of the PTMD with respect to the classical TMD.This is very important because if coupled with additional systems (e.g., monitoring system or weather forecasting systems), the length of the PTMD can be simply passively adapted, before the specific upcoming weather event reaches the structure, to tune the PTMD to the expected dynamic or the expected weather condition in order to maximize its efficiency.This constitutes a valid passive-adaptive (i.e., adaptive but independent from electric power availability) alternative to the active or semi-active technologies which, as is known, do not operate properly under lack of electrical power.In addition, some modern technologies allow the passive adaptation of the TMD inertial mass, 23,24 opening to possibilities of moving towards complete passive-adaptive PTMDs devices.
In this context, it is very important to develop reliable models and design procedures for the installation of PTMDs in wind turbines, with the goal of exploring the global 3-D behaviour of the PTMD-provided turbine under realistic loads in order to correctly quantify the performances of such kind of vibration control technology.This paper uses as basis the OWT environment provided in Petrini et al, 25 applying it for both parked and operating conditions.The main goal of this paper is to provide a design tool for mitigating OWT vibrations using an optimal 3D-PTMD.The NREL monopile 5-MW baseline OWT 22 coupled to the PTMD is modelled using an integrated MATLAB®/ANSYS® routine.A comparison of the OWT without the PTMD is performed between the present model and the NREL OpenFAST software.An in-house built genetic algorithm (GA) toolbox [26][27][28] finds the optimal PTMD parameters for a simplified 2-degrees-of-freedom (2DOF) model minimizing the hub peak displacements, with a combination of wave and wind (rotationally sampled) spectra.

| NUMERICAL GLOBAL MODELLING OF OWTS AND OF THE DEMANDING ENVIRONMENT FOR STRUCTURAL ANALYSIS PURPOSES
OWTs are located in a complex and high-demanding environment 26 with some non-linear interactions and high variability of loading conditions. 29The schematization of such a complex environment is shown in Figure 1, as well as the relevant macro-geometric parameters of the problem: the mean water depth (h) under the still water level (SWL), the hub height (H hub ) above the SWL, the blade length or rotor radius (R) and the bedrock depth (d). 30 general, the main focus of the structural analysis for PTMD design purposes is the peak or the RMS of the structural response.In order to obtain them, the response is modelled as composed of a mean part, induced by the mean wind and the sea current (eventually including the vortex shedding effects) and by a fluctuating/stochastic part, induced by the wind turbulence and by the sea waves.The first one is obtained by a static analysis, while the second one can be evaluated, in frequency domain, by a power spectral density (PSD) analysis.The actions due to the mean wind and the sea current are then modelled by equivalent static forces, while the actions due to turbulent wind and waves are modelled by their power spectrum.
Beside the wind and wave actions, stochastic models are also capable of describing the wind sampling effects of OWTs due to rotating blades in operating conditions.This is modelled by a rotor torque plus a sampling effect on turbulent wind forces.This section briefly summarizes the F I G U R E 1 Environmental actions over an offshore wind turbine models used for obtaining the peak dynamic response of the turbine under wind-current and wave loads for both parked and operating conditions by the static + PSD analyses.For the complete analytical background of the equations reported in this section, the reader is reminded to the Appendix A of the paper.

| Wind field modelling
The wind instantaneous velocity U(M,t) in a generic point M in space, at height z and time t, is separated into two contributions: (i) the mean wind velocity V m over 10 min in the along-wind direction and (ii) an atmospheric turbulence U 0 = [u,v,w]-characterized by high frequency fluctuations.These components are given in a three dimensional Cartesian reference system as the one (x,y,z) shown in Figure 1 by where î, ĵ and k are the versors of x, y and z; V m (z) is the mean wind velocity at height z in the along-wind direction î, which can be considered a stochastic variable in the mid/long run (daily to yearly cyclic variation); u,v and w are, respectively, the longitudinal, lateral and upward components of the turbulent wind, varying in space and time and that, for ordinary winds, are commonly modelled as three (uncorrelated each other) Gaussian, ergodic, stationary, zero-mean, multivariate/monodimensional stochastic processes.Details of the models used for the above-mentioned components are provided in the Appendix A of this paper.

| Aerodynamic loads of blades during operation
The load components due to the mean wind on the rotating blade (deterministic in the short term) and the wind speed fluctuations (stochastic) can be modelled in different ways.By following the blade-element/momentum (BEM) theory, 31 with the math summarized in Appendix A, the aerodynamic characteristics of the rotor can be evaluated.
For what concerns the first of the above-mentioned contributions (mean wind), the aerodynamic forces on the blade can be applied at different radii once all the coefficients of BEM theory have been estimated.The forces per unit length perpendicular to the plane of rotation (dF X ) and in the direction of blade motion (dF Y ) on an element, known as the out-of-plane and in-plane forces, respectively, are computed as Concerning the second contribution (turbulent) load component, the power spectrum of the single blade root bending moment induced by the turbulent wind velocity component in x direction and in rotating blades conditions is estimated by where Δr is the blade element length and S o u and κ o u r 1 ,r 2 , τ ð Þare the rotationally sampled wind velocity spectrum and the cross-correlation for a pair of points on the blade at radii r 1 and r 2 , respectively, 31 both evaluated as detailed in Appendix A: C l is the lift aerodynamic coefficient of the blade, and α is the wind to blade angle of attack.

| Vortex shedding effect on the tower
At certain critical ranges for the flow mean velocities, the frequency of vortex shedding originating around the tower coincides with the first natural (across-wind direction) frequency of the lateral motion of the tower resulting in lock-in vibrations.The lock-in effect can be considered as the applied maximum across-wind displacement r VS across À Á max given by 32 where D is the diameter of the tubular section under wind action and S t and S c are the Strouhal and Scruton numbers, respectively. 332 | Wave and sea current dynamics and induced loads

| Sea current modelling
The sea currents in shallow water can be characterized by a velocity field in which intensity decreases with depth. 34DNV-OS-J101 35 defines the variation in current velocity V curr (z 0 ) with the depth z 0 (see Figure 1) as follows: where v tide (z 0 ) and v wind (z 0 ) are the current velocities generated at the water depth z 0 by the tide and wind, respectively; v tide0 and v wind0 are the tidal and wind-generated current at SWL; h 0 is the reference depth for the wind-generated current (typically 50 m).Unless data indicate otherwise, v wind0 = kV 10 − 1h , with 0.015 ≤ k ≤ 0.03 and V 10 − 1h the 1-h averaged wind speed at 10-m height.

| Waves dynamics
According to Chandrasekaran, 36 the particles' velocities, accelerations and dynamic pressure as function of the surface elevation of the waves can be computed by a number of wave theories such as the linear of first-order or Airy theory, Stokes second-and fifth-order theory, Solitary wave theory, Cnoidal theory, Dean's stream function theory and the numerical theory by Chappelear.The most appropriate theory can be found based on the relations between the wave height H w , wave period T w and water depth h. 37V-OS-J101 35 defines the range of validity for different wave theories in function of the wave steepness parameter, S w = 2πH w = gT 2 w , the shallow water parameter, μ w = 2πh= gT 2 w , and the Ursell parameter, U r = S w = 4π 2 μ 3 w À Á (both defined in Figure 2), as shown in Table 1.
The horizontal velocity and acceleration of water particles can be written as a function of the surface elevation. 38I G U R E 2 Regular travelling properties (adapted from Brebbia and Walker 34 ) where H w (ω, z 0 ) assumes different shapes depending on the adopted wave theory (see Table A1 in Appendix A).For linear Airy theory,

| Morison's equation for the evaluation of hydrodynamic loads
There are two main hydrodynamic forces acting on an underwater supporting slender structural member of an offshore structure acting along the wave propagation direction: the drag and inertia forces.The drag force is proportional to the square of the water particle velocity and is associated to the effects of viscosity in the fluid.The inertia force is independent of the viscosity and is composed of the virtual mass of the submerged member in motion, which can be an added inertial water mass or a force in opposition to the motion of the member and the inertia force due to the accelerating fluid.
In addition to the drag and the inertial forces, the motion of fluid on slender members can give rise to vortex shedding that generates a lift force on the member (orthogonal to the waves propagation direction), which is sinusoidal-like, at a frequency corresponding to the eddy-shedding frequency.

Wave and current-induced inertial and drag force
By applying the Airy wave theory, the drag and inertial forces per unit length in the x direction at depth z 0 can be computed as a function of the horizontal total particle velocity u t (x,z 0 ,t) = V curr (z 0 )+u w (x,z 0 ,t) by the linearized Morison equation as follows: in which c m , c d and c i = c m +1 are, respectively, the added mass, drag and inertia coefficients of the section with diameter D, and ρ w is the water density (taken as 1025 kg/m 3 ).By substituting Equation 8 in the linearized Morison's Equation 9and by assuming that V curr is constant in time (short-term analysis), the total drag + inertial force can be written as follows.

The use of Morison
If η(x,t) is modelled as a stationary stochastic process, the PSD of the wave force F I+D (x,z 0 ,t) can be obtained with the assumption of uncorrelated nature of the stationary process and its derivative 38 as where S η is the power spectrum of the sea-surface elevation.For the proposed model, S η = S JS (see Appendix A, Equation A18, for JONSWAP formulation).

Lift force due to the underwater vortex shedding effect
For steady flow past a cylinder of diameter D, the Strouhal number S t (z 0 ) = Df vs /V curr (z 0 ) represents the non-dimensional frequency of lift forces at a frequency of vortex shedding f vs .Considering an accelerating flow, the presence of dynamic lift forces is expected for long waves and small cylinder.For the Keulegan-Carpenter number (K C ) K C ≥ 15, Keulegan and Carpenter 40 estimate the occurrence of lift forces.Chakrabarti 41 define the magnitude of lift force as where ω 0 is the angular frequency of the incident wave, ψ n the phase angle of the nth harmonic force, and c n l the lift coefficient for the nth harmonic (and a function of K C ).

| Dynamic structural response
As said before, structural analysis in the frequency domain can be performed to evaluate the dynamic response of OWTs subjected to both the environmental and rotating blade actions.Following Davenport, 42 the peak response r p of the structure is expressed in terms of its mean value (r m ) generated by the mean wind (including the VS effect on the tower) and by sea current and the response standard deviation (σ r ), which is related to the effects of the turbulent wind and of the waves, by making use of the response peak factor g r , as reported in Equation 13.The soilstructure interaction is not considered in this paper (the structure is assumed to be fully restrained at the seabed).
where the peak factor is computed as where ν is the cycling rate of the structural response (taken equal to the first eigen frequency of the system) and T wind is the time interval over which the maximum value is evaluated.
Operatively, the r m is computed by a static analysis with the application of the maximum values of the forces or of the displacements evaluated in previous equations for the mean wind (Equation 2), for the vortex shedding effects (Equations 4, which directly computes the response, and 12 by assuming herein cos [−] = 1) and for the sea current (Equation 10), while σ r can be evaluated by a PSD analysis by applying the power spectra of waves' force (Equation 11) and of turbulent wind components as sampled by rotating blades (Equation 3), as the area underpinned by the response PSD for the response parameter r (Manenti and Pretini 29 ), which in the following sections will be the displacement at the hub of the OWT.

| Analytical background
The method proposed to select PTMD parameters follows the optimal PTMD design applied to high towers using GA, proposed by Colherinhas et al. 28 For optimal design purposes, the OWT dynamics is reduced from multiple DOFs (MDOFs) to a single DOF (SDOF) 15,[43][44][45][46] model.Then, the PTMD is attached to the SDOF falling into a 2DOF discrete system described by Figure 3.The PTMD design parameters are the flexural rigidity (K p ), damping (C p ), the mass ratio between the pendulum and system (μ = M p /M s ) and pendulum length (L p ).
The SDOF dynamic system equivalent to the structure without any control device is governed by the Equation 15 below.
The equivalent mass M s , damping C s and stiffness K s of the SDOF are expressed by where EI is the bending stiffness, L is the monopile tower total lenght ('H hub +h' in Figure 1), c is the viscous damping coefficient, m the mass per unit length of the tower beam and M RNA is the mass of the Rotor-Nacelle Assembly.
The PTMD, which is designed to be implemented in the turbine inside the tower and attached to the nacelle (Figure 4), in the analysis for its optimal design, is attached to the SDOF system resulting in a 2DOF discrete system described by Figure 3.The equations of motion of the 2DOF system are given by where θ(t) and y(t) are the two relevant DOFs of the system.Considering F s (t) = e iωt , y(t) = H y (ω)e iωt and θ(t) = H θ (ω)e iωt , the frequency transfer functions H y (ω) and H θ (ω) of the equivalent-2DOF system are, respectively, computed by where In the optimal design process, the 2DOF system is implemented for optimally designing the PTMD with respect to the fore-aft (FA) displacements of the OWT, and the resulting optimal torsional stiffness and damping of the PTMD are considered as optimal also for the side-to-side (SS) direction, to represent the equivalent 3-D PTMD in successive analyses.
For what concerns the external force in frequency domain F s (ω), it is obtained by a linear combination of the environment spectra and the rotating blade spectra divided by the distance (L aero ) between the hub and the aerodynamic centre of the blade (i.e., S FI + D and S M /L aero ), and it is applied in the FA direction.The linear relationship between the input and the output in the frequency domain can be established as y(ω) = H y (ω) F S (ω), leading to the relationship between the PSD of the response and the stochastic excitation in the form. 38 For a stationary stochastic response analysis, it is reasonable to consider a harmonic excitation with an amplitude of . The use of this excitation gives a physical interpretation to the response of the system that can be computed using a pseudo-excitation method 38,47 as

| GA optimization
The optimization is carried out by hierarchically subdividing the parameters of the PTMD in primary and secondary design variables (DVs). 23Then couples of values for K p and C p (secondary DVs) are fixed, and optimal values of μ and L p as primary DVs are found.An in-house built GA toolbox computes the PTMD optimal parameters of the 2DOF model, following the procedures developed by Colherinhas et al. 28 The goal of this genetic optimization is to minimize the frequency response peak displacements of the tower described in the analytical 2DOF model, considering the influence of the environment and rotating blade actions as specified in the previous sections.The fitness function (only in FA direction), maximizing its inverse, as shown Equation 21.
where i corresponds to each chromosome evaluated in the population of N ind individuals.
The ranges of variation from 0.10 ≤ Lp ≤ 10.00 and 0.01 ≤ μ p ≤ 0.20 are chosen for the primary DVs.The case-study OWT is modelled using the FE method following the detailed specifications of the Phase I (monopile type) of the NREL 5-MW OWT by Jonkman and Musial. 22An integrated MATLAB®/ANSYS® APDL routine creates samples of wave and wind rotational spectra, as well as the hydrodynamic and aerodynamic loads over the structure, following the flowchart specifications presented in Appendix B (Figure A2).The gross properties of the NREL baseline OWT modelled are shown in Table 2, while a schematic representation of the FE model developed in ANSYS is shown in Figure 5.
BEAM188 element type is used for the tower and the submerged monopile, and a CTUBE tapered section connects the bottom and the tip of the tower considering the thickness variation along the height.The rotor-nacelle assembly (RNA) is modelled by two MASS21 elements, one for the hub (connected to the tower top with a rigid CERIG command) and other for the nacelle (with a COMBIN14 connection to a tower top).The PTMD was modelled using a rigid connector (BEAM4) of length L p and a MASS21 connected with a COMBIN14 torsional spring at the tower top.
To analyse the stochastic aerodynamic loads of blades in the frequency domain, the rotationally sampled spectrum 31 (see Appendix A and previous sections) is applied at the hub node.The wind turbulent model along the tower follows the Kaimal spectrum together with its respective coherence and was applied along the 10 nodes of the tower, while the wind-induced vortex shedding effect r VS across (Equation 4) was added to the mean value of the peak response in the across-wind (SS) direction.r VS across was imposed to the node located at the hub height, when the wind velocity V m herein (V hub ) falls within reduced velocity range of 0.8/S t < V R < 1.6/S t , being V R the reduced velocity defined as (DNVGL 39 ) where f n is the first natural frequency of the tower in the SS direction.The wave spectral density S η is represented by the JONSWAP spectrum (Appendix A).The wave force PSD (Equation 11) is applied at the two nodes located at and below the SWL.

| Model validation
The modal frequencies in the FA and SS directions of the presented FEM model agree with those found in literature for the case study structure 6,22 (Table 3).
In order to validate the proposed FEM model of the case study when not provided with any control device, a PSD analysis is performed, and the resultant tower-top peak displacements are compared with those obtained by the NREL OpenFAST model (Figure 6).analysis in the frequency domain, using the same processor specifications (i.e., the expended simulation time is more than 60 times).
Figure 6 presents a response PSD comparison of the FEM, SDOF (Equation 15) and OpenFAST models for V hub = 12m/s, H w = 6 m and T w = 10 s for the FA hub displacements.
The FA and SS tower-top peak displacements of the FEM model are, respectively, 1.04 and 0.34 m versus 1.08 and 0.12 m (OpenFAST).The peak stress at the tower base is S b = 18.16N/mm 2 as evaluated in FEM analysis; the last one will increase when the PTMD is installed due to the additional weight provided by the PTMD mass.

| Optimal PTMD design by the 2DOF model
The following parameters are used in the GA-based PTMD optimal design by the 2DOF model introduced in previous sections: tower with a mass of 532.60 ton; RNA with a mass of m = 350.00ton and flexural stiffness of EI ffi 1,036.24GNm 2 (value obtained from the FEM model); generalized tower mass and stiffness, respectively, M s ffi 882.60 ton and K s ffi 2,532.07kNm 16; the tower damping coefficient was set to 1%. 48For a preliminary investigation, the torsional stiffness and the friction damping of the PTMD are considered, respectively, K p = 12.5 Á 10 5 N/m and C p = 9.0 Á 10 3 Nms (secondary DVs set No. 1), 50 enabling the parametric analysis giving the results shown Figure 7.
The response surface shown in Figure 7 has not a simple curvature, but on the contrary, it presents multiple geometric valleys, something that is due to the rotor dynamics.The dark regions in Figure 7 (and successive Figures 8-11) represent the valleys of this complex surface, that is, the lower response peak displacements, presenting optimal values of the GA optimization.
To find the fastest convergence of the optimization procedure by the implemented GA, the following parameters were established during the adopted GA optimization: • N gen = 100, number of generations; • N ind = 100, number of individuals in the population; • p c = 60%, crossover probability; • p m = 2%, mutation probability; • p elit = 2%, elitism probability; • p dec = 20%, decimation rate; • N dec = 20, step of generation for the occurrence of decimation.Using the same strategy of Colherinhas et al, 28 the secondary DVs (stiffness K p and damping C p ) can be changed to obtain additional suboptimal configurations; that is, by increasing K p the main valley is expected to shift towards the right side on the L p × μ plan, while by increasing C p the response amplitude is expected to decrease, but the surface shape and the valley position are not supposed to be affected.With this in mind, a new set of secondary DVs (set No. 2) is selected with K p = 5 Á 10 5 N/m and C p = 1.5 Á 10 4 Nms, leading to the parametric response surface shown in Figure 9.For this new set of secondary DVs, a thousand optimizations are performed, expending a total time of 2771.51 s (ffi2.8s per optimization).
From this new optimization carried out with the DVs set No. 2, four additional DCs (DC5-DC8) are also selected with the criteria of keeping the same mass ratio of the previous DCs (μ DC5 = 0.10, μ DC6 = 0.05, μ DC7 = 0.03 and μ DC8 = 0.02) as shown in Figure 10.The figure shows that the regression curve shifts to the left when stiffness decreases as expected by. 37The power regression for the set No. 2 takes the form of μ = aL b p , with a = 300.8and b = − 5.629, as also shown in Figure 10.The eight DCs selected above are compared by using the 3-D global FEM model described in previous sections.Due to the particular configuration of the system (the PTMD is located inside the tower), as a prerequisite for a DC to be acceptable, in order to avoid collision of the PTMD mass with the inner surface of the tower, the peak lateral displacements (sway) of the PTMD mass must remain lower than the inner radius of the tubular section of the tower at the location of the PTMD mass.The maximum peak lateral sway experimented by the PTMD mass in the selected DCs is 0.91 m, which is acceptable for avoiding such an event.Figure 11  As expected, due to the variety of conditions explored for performance assessment, there is not a unique optimal DC.For example, in FA direction for V hub = 25 m/s, the DC1 is the best among all explored DCs (−2.53 in reduction).On the contrary, for the wind velocity V hub = 12 m/s, DC2 presents peak displacements reductions higher than DC1, with reductions of −1.61 and −0.93 dB in FA and SS direction, respectively.On the basis of feasibility considerations, it can be stated that, since the mass of the DC1 is huge, it's actual implementation is not suggested (the same applies for the DC5), potentially leading to collateral structural problems in the life-cycle time (e.g., addition of such a mass to the system could increase the vulnerability of the tower to fatigue); the DC2 is possibly the best optimized design in FA direction among those obtained for secondary DVs set No. 1 (K p = 12.5 Á 10 5 N/m and C p = 9.0 Á 10 3 Nms).In fact, compared to DC3 and DC4, the DC2 presents a reduction of −2.11 versus −1.56 and −0.79 dB, respectively, in FA direction, at V hub = 25 m/s, while DC2 performs always lower than DC3 and DC4 in SS direction.
The DCs obtained by the secondary DVs set No. 2 (K p = 5.0 Á 10 5 N/m and C p = 15.0Á 10 3 Nms), in FA direction, present, in general, higher peak reductions compared to set No. 1, as expected by. 28For V hub = 25 m/s, there are reductions, in FA direction, in the order of −2.43 dB (DC6) versus −2.11 dB (DC2) for μ = 0.05 and −2.00 dB (DC7) versus −1.56 dB (DC3) for μ = 0.03.DC6 and DC7 are selected as the best DCs considering the variation of the wind speed and the feasibility of the optimal pendulum parameters.The hard task to select these optimum DCs proves the following Ghassempour's et al. 6 assumptions: The conventional design of TMDs based on natural frequencies may not be suitable for OWTs due to the fact that the optimal values of the DVs strongly depend on the wind velocity at operational conditions.
Figure 14 shows a comparison between the tower-top PSD displacements in FA and SS direction for DC6 and DC7, for wind velocities of V hub = 24 m/s (operational condition).
To further point out on the difficulties in choosing between the two selected DCs, Figure 15 shows the peak response displacement reduction of DC6 and DC7 for both FA and SS direction, in function of V hub (in SS, values of V hub < 11.4 m/s are neglected for scale purposes).It is clear that (especially in the FA direction) the prevalence of one DC on the other one strongly depends on the wind velocity.Higher reductions happen in operating conditions for wind velocity values close to the cut-off speed V hub = 25 m/s.
Finally, regarding the possibility of implementing passive-adaptive control strategies, let us assume that for some reason, the PTMD has been set to the DC8, which (from Figure 13) is shown to work well for both the FA and the SS directions in operating conditions at low wind velocities (e.g., V hub = 12 m/s); in the case that strong winds (V hub ≥ 25 m/s) are forecast in a certain day with FA direction expected to become critical for ultimate strength for such strong winds, then L p and μ parameters (e.g., the last one can be changed, as already said, by implementing an inerter 24,25 ) can be switched in advance to temporary fall in the DC5, which, for the same K p and C p , is more efficient than the DC8 in suppressing the FA vibrations at high wind velocities.This passively adaptive skill, allowed by the presence of L p as additional tuning parameter, can be very useful for optimal life-cycle management purposes of the OWT.

| CONCLUSIONS
Design optimization and a performance assessment procedure are presented for PTMDs in OWT systems for the mitigation of the global vibrations of the tower by installing the PTMD inside the tubular tower section and attached to the hub.A simplified 2DOF model is used for optimizing the OWT + PTMD system for different intensities of the actions by a GA, which aims to minimize the peak displacement making use of the stochastic dynamic theory and the frequency domain (PSD) analysis and leading to the identification of a number of alternative design configurations (DCs).The assessment of the performances and the choice of the best DC among the selected ones is conducted by a global 3-D FE model of the turbine still by conducting a PSD frequency domain analysis and by implementing state of art models for the spectra of the wind (as sampled by the rotor in operating conditions) and the waves actions.
On the basis of the obtained results, the following general conclusions can be provided: • The PTMD is able to mitigate wind-and sea-induced global OWT vibrations both in FA and SS directions.
• It is outlined how the complex shape of the response surface in the primary DVs space leads to a set of suboptimal solution DCs of the optimization problem, which, in turn, are more or less suitable depending on the operating condition of the system.
• The use of the 2DOF model during the optimal design phase allows for significant savings in the computational efforts required during the optimization phase.
• The 3-D PSD FE analysis allows for the fast investigation of the performances at different wind velocities and for the individuation of the more suitable DCs for practical realization and for passively adaptive control strategy purposes.
• Before choosing the more suitable design configuration, the admissibility of the peak sway of the PTMD mass and the increasing of the peak stress at the base of the tower due to the installation of the PTMD additional mass must be checked.
• Moreover, the possibility of passively adapting the pendulum length, if coupled with other passively adaptive skills (rearrangement of K p , C p or μ), allows for the temporary switching from an identified optimal (DC) to different ones, which are more appropriate at different operating conditions (e.g., operating versus parked rotor at different wind intensities).
From the application of the outlined procedure to a case study, a reduction of over 20% of the response peaks for high wind velocities in the FA direction was noticed, which is consistent with the response reduction provided by other semi-active PTMD vibration control systems found in literature, 51 but in the device proposed here, there is the advantage of not requiring electric power in operation.
The paper contributes to provide advanced tools for the optimal design analysis of the PTMD with the purpose of its installation in OWTs as an initial step towards the implementation of passively adaptive vibration control strategies as a convenient alternative to the active and semi-active controls, which are complex to calibrate and manage.
The following research topics are currently under development by the same authors of this paper as further research steps: • development of advanced action models for including non-linear or transitory effects of the environment (hurricane or thunderstorm winds and non-linear waves); • inclusion of the uncertainties in the analysis and probabilistic evaluation of the performances; • consideration of parametric multi-objective optimization using GAs considering the energy produced and the fatigue strength of the system.

PEER REVIEW
The peer review history for this article is available at https://publons.com/publon/10.1002/we.2590.

ORCID
where V hub is the mean wind speed at hub height H hub of the wind turbine and the power law exponent, α 0 , shall be assumed to be 0.2 for normal wind conditions onshore or 0.14 for normal wind conditions offshore locations. 52,53e turbulence components along the three axes are modelled as three independent, multivariate stationary Gaussian ergodic stochastic processes.
To consider the effects of the spatial variation of turbulence in the lateral and vertical direction due to the variation of the wind inside the vertical plan spanned by the moving blades, the spectral description of turbulence must be extended to include information about the crosscorrelations between turbulent fluctuations at points separated laterally and vertically.These correlations decrease as the distance separating two points increases and can be described by 'coherence' functions, which describe the correlation as a function of frequency and separation. 31The coherence C jk (Δr, n) is defined by where n is the frequency, S jk is the cross-spectrum of variation at two points separated by the distance Δr and S jj and S kk are the spectra of variations at each of the points.
Both IEC 61400-1 52 and DNV-OS-J101 35 propose the use of the Kaimal spectrum (modified Karman's PSD 54 ), unless data indicate otherwise, defined by where the integral scale parameter L K and the standard deviation σ K values are given in Table A2.whereσ 1 = I ref (0.75U hub +5.6),I ref is the reference value of the turbulence intensity and Λ 1 the longitudinal turbulence scale parameter. 52A B L E A 1 Velocities and accelerations of first-and second-order wave theories 39

Parameters
Airy wave theory Stokes second-order wave theory General water depth Deep water Horizontal particle velocity, u w πHw T An exponential coherence model may be used in conjunction with the Kaimal autospectrum, taking into account the spatial correlation in the longitudinal velocity component, defined by 52 where C u (Δr, n) is the coherence function defined by the complex magnitude of the cross-spectral density of the longitudinal wind velocity components at two spatially separated points divided by the autospectrum function; Δr is the projection of the separation vector between the two points on to a plane normal to the average wind direction.IEC 61400-1 does not specify the coherence of the other two components to be used in conjunction with the Kaimal model.Burton et al. 31 show that the following expression is often used for those cases A.2 | Blade element momentum theory Blade element momentum (BEM) theory assumes that the forces on a blade element of radius r and length δr (Figure A1) can be calculated by means of two-dimensional aerofoil characteristics using an estimated angle of attack α.
Figure A1A shows that the velocity components at a radial position on the blade can be expressed in terms of the axial and tangential flow induction factors (respectively, a and a 0 ), the wind speed U(1 − a) and the tangential rotational speed of the blade elements (Ωr).Having information about how the aerofoil characteristic coefficients C l and C d vary with the angle α, the forces on the blades can be computed in function of a and a 0 .The resultant relative velocity at the blade is where U ∞ is the mean free stream velocity.
From Figure A1B, the angle of attack is given by α = ϕ − β, where β is the inclination of local blade chord to rotor plane (i.e., blade twist plus pitch angle) and the inflow angle ϕ can be write as follows: Hansen 55 shows that the induction factors can be computed interactively, using the following equations: where σ r ð Þ = c r ð ÞB 2πr is the solidity coefficient evaluated as fraction of the annular area by the rotor disc area, B is the number of blades, c(r) the chord of each blade element at distance r from the rotor horizontal axis, C n = C l cos ϕ+C d sin ϕ and C t (r) = C l (r) sin ϕ − C d (r) cos ϕ are normalized coefficients of the lift and drag projection, and the correction of the assumption of an infinite number of blades is the Prandtl's tip loss factor F, defined as The BEM theory is, then, fully defined by applying the Glauert correction for high values of a (see Hansen 55 ) and can also be extended adding root losses, 31 estimating properly the aerodynamic coefficients, induction factors and angle of attack for each blade element.According to Burton et al., 31 Wilson and Lissaman 56 argue that the drag coefficient should not be included in the calculus of the normalized coefficients.
If drag is excluded from the determination of the flow induction factors a and a 0 , the total torque Q developed by the rotor is defined as where ρ is the air density (taken as 1225 kg/m 3 ), λ tip = RΩ/V m is the tip speed ratio and μ = r/R is the ratio between the blade element radius r and the total radius of the blade R (Figure A1), c(r) is the chord of the blade, C d is the blade drag coefficient and Ω is the rotational velocity of the rotor and W = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi The mechanical power developed by the rotor is P = QΩ and the power coefficient is computed by the following:

A.3 | Aerodynamics during operation
The standard deviation of the blade root bending moment can be computed as where κ u (r 1 , r 2 , 0) is the cross-correlation function κ u (r 1 , r 2 , τ) between the wind fluctuations at radii r 1 and r 2 with the time lag τ set to zero, that is, To find the energy content of the incident wind fluctuations, it is necessary to examine some points on the rotating blade at the blade natural frequencies, provided by a rotationally sampled spectrum.The rotationally sampled spectrum method can be obtained by deriving the power spectrum of the wind seen by a point on a rotating blade, based on the Fourier transform pair where S u (n) is the single sided spectrum of wind speed fluctuations in terms of frequency n (in Hz) and κ u the autocorrelation function for the along wind turbulent fluctuations at a fixed point in space from the corresponding spectrum, for a homogeneous and isotropic turbulence, and an incompressible flow.
After the autocorrelation function, κ u (τ), is found, it is computed the for a point on the rotating blade at a radius r, κ o u r, τ ð Þ, by a derivation (the superscript o denotes a point on a rotating blade).This function is transformed to yield the rotationally sampled spectrum.Burton et al. 31 explain BERTOLLUCCI COLHERINHAS ET AL.
this process by the following three steps (derivation of the autocorrelation function at a fixed point, derivation of the autocorrelation function at a point on the rotating blade and derivation of the power spectrum seen by a point on the rotating bade).
The last step consists in using a DFT to find the integral of the Fourier transform pair of a rotating blade.The transform becomes where N is the number of points in the time series of κ Ão u r, pT=N ð Þ(the asterisk denotes that κ o u r, τ ð Þ is 'reflected' for T > T/2) and the PSD is calculated at the frequencies n k = k/T for k = 0,1,2,…,N − 1.The sum can be evaluated using a fast Fourier transform (FFT) providing for N, a value equals to a power of 2.
After these three steps are applied, the cross-spectrum can be computed for a pair of points at radii r 1 and r 2 on a rotating blade, and the cross-correlation function becomes where s = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U 2 τ 2 + r 2 1 + r 2 2 −2r 1 r 2 cosΩτ q , K ν (x) is a modified Bessel function of the second kind and order ν and Γ(x) is the gamma function.

A.4 | Wave theories
For linear waves, the phase velocity only depends on wave length λ, and it is independent of the wave amplitude as follows: The DNV recommended practice for environmental conditions and loads 39 presents the values of linear wave particle velocity and accelerations in.
where H s is the significant wave height, n p = 1/T p is the peak frequency, T p is the peak period; if H s and T p are expressed in m and s, respectively, then C(γ) = 1 − 0.287 Á lnγ is the normalizing factor, and γ is the peak-shape parameter, defined as for which W (resultant relative wind velocity), a, a 0 (induction factors) and F (tip loss factor) are defined in Appendix A (Equations A6-A9); c(r) is the blade chord (see Appendix A and FigureA1); C n = C l cos ϕ+C d sin ϕ and C t (r) = C l (r)sin ϕ − C d (r) cos ϕ with C l (r) and C d (r) being the lift and drag coefficients of the blade element; ρ is the air density; and N b is the number of blades.By integrating the above equation along the blade, the total mean forces acting on each blade and transmitted to the rotor can be evaluated.
's equation in the estimation of the hydrodynamic loads on slender member should take into account the variation of c d and c m as functions of the Reynolds number (Re = VD/ν, where ν is the kinematic viscosity of the fluid and V is a characteristic velocity of the flow), the Keulegan-Carpenter number (K C = U m T/D, where U m = max[u w ]) and the non-dimensional roughness (Δ = k b /D, where k b is the bottom sea surface roughness height), that is, c d = c d (Re, K C , Δ) and c m = c m (Re, K C , Δ).For an infinite length circular cylinder in an infinite length medium, c m = 1.
Structure with a linear pendulum attached (2DOF) excited by a force F s (t) F I G U R E 4 A schematic description of a PTMD (adapted from Sun and Jahangiri 21 )

T A B L E 2 5
Gross properties of the NREL 5-MW baseline wind turbine 48 Schematic description of the FE model of the OWT developed in ANSYS.Nodes and loads (left); type of FE used (right) T A B L E 3 Natural frequency comparison (values in Hz)

F I G U R E 6
PSD for tower-top displacements of the OWT model.Full frequency range (up); first-peak zoom (bottom) F I G U R E 7 Parametric optimization analysis for the secondary DVs set No. 1 (K p = 12.5 Á 10 5 N/m and C p = 9.0 Á 10 3 Nms) F I G U R E 8 Optimal results of L p and μ obtained with the secondary DVs set No. 1 (colours are associated with the same magnitude scale of previous Figure 7) F I G U R E 9 Parametric analysis of the secondary DVs set No. 2 (K p = 5.0 Á 10 5 N/m and C p = 15.0Á 10 3 Nms) F I G U R E 1 0 Optimal results of L p and μ obtained with the secondary DVs set No. 2(colours are associated with the same magnitude scale of previous Figure 9) The OpenFAST response PSD is computed as the mean of 20 PSDs produced using the MATLAB Toolbox for OpenFAST, 49 from the discrete Fourier transform (DFT) of 20 different time histories obtained by OpenFAST simulations conducted at the same wind intensities.The turbulent stochastic wind velocity fields are generated using the NREL TurbSim software, which works as inputs of the AeroDyn OpenFAST module.The generation of the mean OpenFAST PSD takes a total of 85.52 min with the available computational capacity versus 79.32 s of the ANSYS FEM This way, 400 iteration analyses are carried out by the GA optimization.The results are gathered on the μ versus L p plane, where the power regression line of the optimal solutions found with the secondary DVs set 1 takes the form of μ = aL b p + c, with a = 1.185Á 10 4 , b = − 7.354 and c = 9.249 Á 10 −3 , as shown in Figure 8.Although optimal values fluctuate around the power regression curve of Figure 8, this regression curve (located in a main valley of the surface) gives a good clue for a predesign of PTMDs.Four design cases (DC) are selected (DC1-DC4) from the one obtained considering the secondary DVs set No. 1, with μ DC1 = 0.10 (M p = 85.11t), μ DC2 = 0.05 (M p = 44.88t), μ DC3 = 0.03 (M p = 30.64t) and μ DC4 = 0.02 (M p = 14.55 t).The DC1-DC4 are also indicated in Figure 8.

F I G U R E 1 1
Stresses at the base S b .Secondary DVs set No. 1 (left); secondary DVs set No. 2 (right)

F I G U R E 1 2
Controlled and uncontrolled OWT tower-top peak response displacements for optimal DCs obtained with the secondary DVs set No. 1 (DC1-DC4), in FA and SS directions F I G U R E 1 3 Controlled and uncontrolled OWT tower-top peak response displacements for optimal DCs obtained with the secondary DVs set 2 (DC5-DC8), in FA and SS directions T A B L E 4 Resume of selected optimal DCs (results in dB) Design parameters Secondary DVs set No. 1 (K p = 12.5 Á 10 5 N/m and C p = 9.0 Á 10 3 Nms) Secondary DVs set No. 2 (K p = 5.0 Á 10 5 N/m and C p = 15.0Á 10 3 Nms) shows the peak stresses S b at the base for the different DCs at different V hub values.The local increasing of the response values obtained for each DC at relatively low velocities is due to the wind induced vortex shedding effect occurring at the critical velocity range specified above.A cut-off wind speed V out = 25 m/s has been set for V hub (i.e., for V hub ≥ V out , the rotor is parked), resulting in a local drop of the response from V out on.The decreasing trend of S b by increasing the DC from 1 to 4 and from 5 to 8 is partially due to the contingent decreasing of the pendulum mass.DC1 and DC5 present the higher values of stresses at the base, close to 22 N/mm 2 at moderate wind speeds and for V hub < 8 m/s, due to vortex induced vibrations.In Figures 12 (secondaryDVs set No. 1) and 13 (secondary DVs set No. 2), the peak response displacements at the hub d p for FA and SS directions are compared for all the considered DCs and shown as function of V hub ; at H w = 6 m and T w = 10 s, the response of the uncontrolled configuration is also shown for comparison purposes.For V hub ≥ V out = 25 m/s, that is, cut-out speed, the rotational spectrum is set off (i.e., it is replaced by an ordinary Kaimal wind spectra acting of parked blades), and there is a drop of the peak displacement curves.It also can be noticed an increase of the peak displacements at low wind speed in SS direction, due to vortex induced vibrations.F I G U R E 1 4 PSD of tower-top displacements of DC6 and DC7 for V hub = 24 m/s F I G U R E 1 5 Peak response reduction between the control and uncontrolled OWT In general, it is graphically shown that the DCs corresponding to the secondary DVs set No. 2 present higher peak response displacement reductions than those of the DVs set No. 1.A detailed study for V hub = 12 m/s (operating conditions and out from the vortex-shedding velocity range for the uncontrolled case and for all DCs) and 25 m/s (parked conditions) is performed, and a comparison summary of the peak displacements reduction and the corresponding stress at the base (SB) is presented in Table 4 by considering the OWT with and without PTMD for both FA and SS directions.The comparison of results is presented in dB with the purpose of facilitating the understanding of the performances, that is, for the response peak displacements: r dB p = 10 Á log 10 r PTMD p =r p ; for the stress at the base: σ dB b = 10 Á log 10 σ PTMD b =σ b À Á (where the superscript PTMD denotes the OWT with the PTMD).By adopting these units, better PTMD performances (larger reduction of the response) correspond to lower (negative) values.
A) Description of a blade element; (B) blade element velocities (adapted from Burton et al 31 )

A. 5 |
PSD of wind wavesThe two empirically derived PSDs commonly used are the Pierson Moskowitz (S PM ) and Joint North Sea Wave Project (JONSWAP-S JS ) 0.07 for n ≤ n p and σ = 0.09 for n > n p .