0
Research Papers

Parametric Instability and Localization of Vibrations in Three-Blade Wind Turbines PUBLIC ACCESS

[+] Author and Article Information
Takashi Ikeda

Department of Mechanical Systems Engineering,
Hiroshima University,
1-4-1, Kagamiyama,
Higashi-Hiroshima, Hiroshima 739-8527, Japan
e-mail: tikeda@hiroshima-u.ac.jp

Yuji Harata

Department of Mechanical Engineering,
Aichi Institute of Technology,
1247 Yachigusa, Yakusa-cho,
Toyota, Aichi 470-0392, Japan
e-mail: y-harata@aitech.ac.jp

Yukio Ishida

Institute of International Education and Exchange,
Nagoya University,
Fro-cho, Chikusa-ku,
Nagoya Aichi, 464-8601 Japan
e-mail: ishida@nuem.nagoya-u.ac.jp

1Corresponding author.

Contributed by the Design Engineering Division of ASME for publication in the JOURNAL OF COMPUTATIONAL AND NONLINEAR DYNAMICS. Manuscript received July 6, 2017; final manuscript received March 22, 2018; published online May 17, 2018. Assoc. Editor: Brian Feeny.

J. Comput. Nonlinear Dynam 13(7), 071001 (May 17, 2018) (11 pages) Paper No: CND-17-1297; doi: 10.1115/1.4039899 History: Received July 06, 2017; Revised March 22, 2018

Nonlinear vibration characteristics of three-blade wind turbines are theoretically investigated. The wind turbine is modeled as a coupled system, consisting of a flexible tower with two degrees-of-freedom (2DOF), and three blades, each with a single degree of freedom (SDOF). The blades are subjected to steady winds. The wind velocity increases proportionally with height due to vertical wind shear. The natural frequency diagram is calculated with respect to the rotational speed of the wind turbine. The corresponding linear system with parametric excitation terms is analyzed to determine the rotational speeds where unstable vibrations appear and to predict at what rotational speeds the blades may vibrate at high amplitudes in a real wind turbine. The frequency response curves are then obtained by applying the swept-sine test to the equations of motion for the nonlinear system. They exhibit softening behavior due to the nonlinear restoring moments acting on the blades. Stationary time histories and their fast Fourier transform (FFT) results are also calculated. In the numerical simulations, localization phenomena are observed, where the three blades vibrate at different amplitudes. Basins of attraction (BOAs) are also calculated to examine the influence of a disturbance on the appearance of localization phenomena.

Renewable, natural power sources, such as wind and solar, have generated much attention because using them reduces the consumption of fossil fuels and thereby lessens the burden on the environment. Wind power is one of the most promising renewable, natural power sources, and wind turbines are being used in many countries. Wind turbines consist of a tower and blades, which vibrate due to both wind velocity and the rotation of blades. These vibrations can lead to severe accidents, including the collapse of the tower and/or breakdowns [1]. Furthermore, as the size of wind turbines continues to increase, as a means to reduce the cost of power generation, the stiffness becomes comparatively small. This results in a higher risk of vibration problems and their serious consequences.

Many studies on the vibrations in wind turbines have been conducted. Larsen and Nielsen [2,3] investigated parametrically excited vibrations, and their stabilities in a wind turbine blade subjected to sinusoidal nacelle motion. Kondo [4] presented a mathematical model of a wind turbine blade and confirmed its validity by comparing the numerical results with experimental data. Chopra and Dugundji [5] investigated nonlinear vibrations, caused by self-excited flutter oscillations, for a wind turbine blade with three degrees-of-freedom (3DOF). Inoue et al. [6] showed the occurrence of superharmonic resonance in an elastic rotating cantilever. Ichikawa [7] discussed the instability of an asymmetric blade in a horizontal-axis wind turbine using Routh's criterion. Yamane [8] analyzed the stability of a wind turbine composed of a flexible two-blade rotor and a flexible tower based on Floquet theory. Saravia et al. [9] used the finite element method to calculate dynamic instability regions, due to parametric excitation, for rotating, thin-walled composite beams. Ishida et al. [10] investigated unstable vibrations, caused by both wind velocity and parametric excitation, in a rigid blade supported by a rotational spring. Acar et al. [11] analyzed the in-plane blade-hub dynamics in a three-blade, horizontal-axis wind turbine and showed the response curves of superharmonic and subharmonic resonances due to external and parametric excitations created by gravity. Chao et al. [12] showed that the equations of motion are quite similar to those of centrifugal pendulum vibration absorbers. Ikeda et al. [13] investigated unstable vibrations due to the rotation of an asymmetric blade in two-blade wind turbines. To avoid unstable vibrations, teetering two-blade wind turbines have been invented. Janetzke and Kaza [14] performed a stability analysis of whirl flutter, and Yamane et al. [15] and Motoi [16] conducted experiments for teetering two-blade wind turbines. The studies mentioned above focus on analyzing the vibrations of wind turbines, while there are a few studies on controlling such vibrations. Ikeda et al. investigated the vibration control of wind turbines using tuned liquid dampers [17,18] and tuned mass dampers [19]. Although previous studies mentioned earlier [27,9,10] have considered the vibrations of wind turbine blades, the nonlinear vibrations of wind turbines, due to the effect of the tower motion on the blades, have not been thoroughly investigated.

This paper investigates nonlinear vibration characteristics of three-blade wind turbines. The wind turbines are modeled as a coupled system, composed of a flexible tower with two degrees-of-freedom (2DOF), and three blades, each with a single degree-of-freedom (SDOF). These wind turbines form nonlinear systems due to the nonlinear restoring moments acting on the blades. The blades are subjected to winds, and the wind velocity increases proportionally with height due to vertical wind shear. The free vibrations of the wind turbine are theoretically analyzed to obtain the natural frequency diagram. The corresponding linear system with parametric excitation terms is analyzed to determine the rotational speeds at which unstable vibrations appear in the blades. The nonlinear system, considering the nonlinear restoring moments acting on the blades, is then analyzed, and swept-sine tests are applied to the equations of motion for the wind turbine to obtain the frequency response curves. Stationary time histories and their fast Fourier transform (FFT) results are also shown. In the numerical simulations, it is found that “localization phenomena” occur in the blades, which vibrate at different amplitudes. Basins of attraction (BOAs) are also shown to examine the influence of a disturbance on the appearance of localization phenomena.

Equations of Motion.

Figure 1 shows an illustration of a three-blade wind turbine. A nacelle, in which a generator and a gear box are set, is installed on top of a flexible tower. Three blades are attached to the generator shaft. In the stationary Cartesian coordinate system O-x0y0z0, the z0-axis coincides with the shaft axis, and the center of the nacelle coincides with the origin O. Figures 2(a) and 2(b) show the top and side views, respectively, of an analytical, 5DOF model of a wind turbine. The 5DOF model is composed of a 2DOF system, which represents the tower, and three SDOF blades. Each blade is assumed to be a solid bar with out-of-plane, i.e., flapwise vibration characteristics. The 2DOF system consists of mass mT, springs with spring constants kz and kx, and dashpots with viscous damping coefficients cz and cx in the z0- and x0-directions, respectively. The nacelle and the shaft are assumed to be a lumped mass which is included in the system's mass. The motion of mass mT is restricted to the x0z0-plane. Blade “i” (i = 1, 2, 3), with mass mi, length li, width bi, thickness hi, and density ρi, is connected to the shaft by a rotational spring, with stiffness ki, and a damper, with viscous damping coefficient ci. The inclination angle of blade “i” is designated as θi. Blades “2” and “3” are positioned at angles α2 and α3, respectively, from the angular position α1(=0 deg) of blade “1.”

The moving Cartesian coordinate system Ob-xyz is fixed to the nacelle and is shown in Fig. 3, where the z-axis coincides with the center of the shaft, and the y-axis is vertical, through the supporting point of the blades. Let ω be the rotational speed of the blade around the z-axis. The coordinates of the infinitesimal element dsi of blade “i” are given by Display Formula

(1)xi=x0+sicosθicos(ωt+αi)yi=sicosθisin(ωt+αi)zi=z0sisinθi}

where (x0, z0) is position of the 2DOF system's mass, and si is the distance of the element from the origin Ob. The total kinetic energy T of the wind turbine is given by Display Formula

(2)T=12mT(x˙02+z˙02)+i=130li12ρibihi(x˙i2+y˙i2+z˙i2)dsi

Substituting Eq. (1) into Eq. (2), we obtain the following equation: Display Formula

(3)T=12mT(x˙02+z˙02)+i=1312mi(x˙02+z˙02)+i=1312mi[13li2θ˙i2liθ˙ix˙0sinθicos(ωt+αi)liθ˙iz˙0cosθi+13li2ω2cos2θiliωx˙0cosθisin(ωt+αi)]

where mi = ρibihili. The total potential energy U of the wind turbine includes the components generated from the 2DOF system's stiffness, gravity, and the rotational spring supporting the blade. Therefore, U is given by Display Formula

(4)U=12kxx02+12kzz02+i=1312kiθi2+i=130liρibihigsicosθisin(ωt+αi)dsi=12kxx02+12kzz02+i=1312kiθi2+i=1312miglicosθisin(ωt+αi)

where g is the acceleration of gravity.

It is assumed that the wind velocity increases proportionally with height, due to vertical wind shear, as shown in Fig. 4. Let V0 be the wind velocity for y = 0 at origin Ob, and V0±Vi at y=±li, where it is assumed that |V0|>|Vi|. The nonlinear moment Ni acting on blade “i” around origin Ob due to the wind velocity is given by Display Formula

(5)Ni=0li(V0+Vi×yili)2ρAbisicos2θidsi=ρAbili2cos2θi[12V02+23V0Vicosθisin(ωt+αi)+14Vi2cos2θisin2(ωt+αi)]

where ρA is air density. When the nonlinear moment acts on blade “i,” the virtual work δWi generated by the virtual displacement δθi is given as Display Formula

(6)δWi=Niδθi

Therefore, the generalized moment that acts on blade “i” is Ni. Substituting Eqs. (3)(5) into Lagrange's equations, and incorporating the viscous damping terms, we obtain the equations of motion in the x0- and z0-directions of the 2DOF system in Eqs. (7a) and (7b), respectively, and the three blades in Eq. (7c) as follows: Display Formula

(7a)(mT+m1+m2+m3)x¨0+cxx˙0+kxx0i=13mili2[(θ¨isinθi+θ˙i2cosθi+ω2cosθi)cos(ωt+αi)2ωθ˙isinθisin(ωt+αi)]=0
Display Formula
(7b)(mT+m1+m2+m3)z¨0+czz˙0+kzz0i=13mili2(θ¨icosθiθ˙i2sinθi)=0
Display Formula
(7c)mili23θ¨i+ciθ˙i+kiθi+mili2ω23sinθicosθimilig2sinθisin(ωt+αi)mili2{z¨0cosθi+x¨0sinθicos(ωt+αi)}ρAbili2cos2θi[12V02+23V0Vicosθisin(ωt+αi)+14Vi2cos2θisin2(ωt+αi)]=0(i=1,2,3)

The following dimensionless parameters are introduced: Display Formula

(8)x0=x0xst,z0=z0xst,li=lixst,bi=bixst,ci=cxp0M,ci=czp0M,ci=ciMp0xst2,kx=kxp02M,kz=kzp02M,ki=kiMp02xst2,V0=V0xstp0,Vi=Vixstp0,μA=ρAxst3M,μi=miM,μT=mTM,ω=ωp0,t=p0t(i=1,2,3)}

where M=mT+m1+m2+m3, xst=Mg/kz, and p0=kz/M. It should be noted that all primes will be omitted hereafter for simplicity although the quantities are still dimensionless in the theoretical analysis and results. Applying Eq. (8) to Eq. (7) results in the equations of motion for the 2DOF system and blades in dimensionless form as follows: Display Formula

(9a)x¨0+cxx˙0+kxx0i=13μili2[(θ¨isinθi+θ˙i2cosθi+ω2cosθi)cos(ωt+αi)2ωθ˙isinθisin(ωt+αi)]=0
Display Formula
(9b)z¨0+czz˙0+kzz0i=13μili2(θ¨icosθiθ˙i2sinθi)=0
Display Formula
(9c)μili23θ¨i+ciθ˙i+kiθi+μili2ω23sinθicosθiμili2sinθisin(ωt+αi)μili2[z¨0cosθi+x¨0sinθicos(ωt+αi)]μAbili2cos2θi[12V02+23V0Vicosθisin(ωt+αi)+14Vi2cos2θisin2(ωt+αi)]=0(i=1,2,3)

Note that parametric excitation terms and nonlinear terms such as sinθi and cosθi are included in Eq. (9).

Natural Frequencies.

In order to determine the natural frequency of the wind turbine, Eq. (9) is linearized and expressed without considering damping terms, parametric excitation terms, and wind velocity as follows: Display Formula

(10)x¨0+kxx0=0(a)z¨0+kzz0i=13μili2θ¨i=0(b)μili23θ¨i+kiθi+μili2ω23θiμili2z¨0=0(i=1,2,3)(c)}

Note that Eq. (10a) is independent of Eqs. (10b) and (10c). The solutions of the free vibrations for Eqs. (10a)(10c) are assumed as follows: Display Formula

(11)x0=Axcos(pt+β)z0=Azcos(pt+β)θi=Aicos(pt+β)(i=1,2,3)}

where p and β represent the natural frequency and the phase angle, respectively. Substituting Eq. (11) into Eq. (10) and equating the coefficients on both sides of the resulting equations, we obtain the simultaneous equations for unknown variables Aj (j = x, z, 1, 2, 3) in matrix form as follows: Display Formula

(12)B{AxAzA1A2A3}T=0

where B is a 5 × 5 matrix. In order to have non-trivial solutions of Eq. (12), the following equation must hold: Display Formula

(13)det(B)=0

Equation (13) is the frequency equation, and from it, the natural frequencies of the 5DOF system can be numerically obtained for an arbitrary value of the rotational speed. When the values of μi, li, and ki are identical to μ, l, and k, respectively, Eq. (13) is reduced to Display Formula

(14)(p2kx)[μl2(p2ω2)3k]2{μl2(49μ)p44[μl2(kz+ω2)+3k]p2+4kz(μl2ω2+3k)}=0

The solutions of Eq. (14) are defined as Display Formula

(15)px=kx(a)pj=ω2+3k/(μl2)(j=2,3)(b)μl2(49μ)p44[μl2(kz+ω2)+3k]p2+4kz(μl2ω2+3k)=0p1,p4(c)}

where px is the natural frequency of the 2DOF system in the x0-direction and is constant. Here, p1 to p4 are the natural frequencies of the 2DOF system in the z0-direction, coupled with the three blades. Note that pj (j = 2, 3) are double roots.

Unstable Regions.

In this section, the corresponding linear system with parametric excitation terms is analyzed to determine the rotational speeds where unstable vibrations appear due to parametric excitation. The regions where these unstable vibrations appear can be used to predict at what rotational speeds the blades may vibrate at high amplitudes in a real wind turbine, which could cause it to breakdown. Equation (9) is linearized and expressed without considering wind velocity as follows: Display Formula

(16)x¨0+cxx˙0+kxx0=i=13μiliω22cos(ωt+αi)(a)z¨0i=13μili2θ¨i+czz˙0+kzz0=0(b)μili23θ¨iμili2z¨0+ciθ˙i+[ki+μili2ω23μili2sin(ωt+αi)]θi=0(i=1,2,3)(c)}

Equation (16a) is independent of Eqs. (16b) and (16c), does not include any parametric excitation terms, and hence does not affect the instability of the wind turbine. Equations (16b) and (16c) are coupled with each other, and because Eq. (16c) includes parametric excitation terms of the frequency ω, unstable vibrations may appear. Each unstable vibration has a predominant frequency n(ω/2) (n: positive integer), hereafter designated as [n(ω/2)], and the rotational speeds at which it occurs constitute an unstable region. Furthermore, periodic vibrations of frequencies n(ω/2) (n = 2m; m: positive integer) may appear on the boundaries of unstable regions. The solutions of these periodic vibrations are assumed based on Floquet theory [20] as follows: Display Formula

(17)z0=u1cosωtv1sinωt+u2cos2ωtv2sin2ωt+u3cos3ωtv3sin3ωt+R1θ1=u4cosωtv4sinωt+u5cos2ωtv5sin2ωt+u6cos3ωtv6sin3ωt+R2θ2=u7cosωtv7sinωt+u8cos2ωtv8sin2ωt+u9cos3ωtv9sin3ωt+R3θ3=u10cosωtv10sinωt+u11cos2ωtv11sin2ωt+u12cos3ωtv12sin3ωt+R4}

Substituting Eq. (17) into Eqs. (16b) and (16c) and equating the coefficients of the terms of the frequencies (i = 1, 2, 3) and the constant terms on both sides of each resulting equation, we obtain the following equation: Display Formula

(18)A0X=0

where A0 is a 28 × 28 matrix whose elements consist of the system parameters, and Display Formula

(19)X=(u1,v1,,u12,v12,R1,R2,R3,R4)T

The condition det(A0) = 0 must be satisfied so that Eq. (18) has nontrivial solutions for the periodic vibrations. The boundaries of the unstable regions for Eq. (17) can then be determined from det(A0) = 0.

Periodic vibrations of frequencies n(ω/2) (n = 2m − 1; m: positive integer) may also appear on the boundaries of unstable regions, and their solutions are assumed as follows: Display Formula

(20)z0=u1cos(ω/2)tv1sin(ω/2)t+u2cos(3ω/2)tv2sin(3ω/2)t+u3cos(5ω/2)tv3sin(5ω/2)t+R1θ1=u4cos(ω/2)tv4sin(ω/2)t+u5cos(3ω/2)tv5sin(3ω/2)t+u6cos(5ω/2)tv6sin(5ω/2)t+R2θ2=u7cos(ω/2)tv7sin(ω/2)t+u8cos(3ω/2)tv8sin(3ω/2)t+u9cos(5ω/2)tv9sin(5ω/2)t+R3θ3=u10cos(ω/2)tv10sin(ω/2)t+u11cos(3ω/2)tv11sin(3ω/2)t+u12cos(5ω/2)tv12sin(5ω/2)t+R4}

Equation (20) is substituted into Eqs. (16b) and (16c) in the same manner as Eq. (17), which leads to Display Formula

(21)B0X=0

where B0 is a 28 × 28 matrix whose elements consist of the system parameters. The boundaries of the unstable regions for Eq. (20) can then be determined from det(B0) = 0.

Natural Frequencies.

Figure 5 shows the natural frequency diagram of the three-blade wind turbine, calculated from Eq. (15). The horizontal and vertical axes represent the rotational speed ω of the blades and the natural frequency p of the wind turbine, respectively. The values of the wind turbine's parameters are μT = 0.94, μi = 0.02 (i = 1, 2, 3), kx = kz = 1.0, ki = 0.5, li = 40.0, bi = 2.0, α1 = 0 deg, α2 = 120 deg, and α3 = 240 deg. Natural frequencies px, p1, p2, p3, and p4 are indicated by the solid lines and are defined in Eq. (15), where p1 > p2 = p3 > p4. Note that px is kx=1.0 obtained from Eq. (15a), because the equation of motion of the 2DOF system in the x0-direction, Eq. (10a), is independent of Eqs. (10b) and (10c). The natural frequencies, p2 and p3, correspond to the double roots of Eq. (15b). As the rotational speed increases, p2 and p3 also increase, due to centrifugal force, and asymptotically approach p = ω. The natural frequencies, p1 and p4, are obtained from Eq. (15c). The natural frequency p1 is close to px at lower rotational speeds and asymptotically approaches p = ω as the rotational speed increases. By contrast, p4 is close to p = ω at lower rotational speeds, and asymptotically approaches p = 1 as the rotational speed increases. The natural frequencies of the blades correspond to p4 for ω < 1 and p1 for ω > 1. The dashed lines represent p = n(ω/2). Points A, B, C, and D represent the rotational speed at which unstable vibrations [n(ω/2)] for n = 3, 4, 5 and 6 may appear, respectively. Note that p = (1/2)ω is omitted, because it does not intersect p4 for ω < 1, or p1 for ω > 1, and thus unstable vibration [(1/2)ω] does not appear in the blades.

Unstable Regions.

Figure 6 shows the unstable regions caused by parametric excitation. The horizontal and vertical axes represent the rotational speed ω of the blades and the blade length li, respectively. The values of the parameters are same as those in Fig. 5, and in addition, cx = cz = 0.03 and ci = 0.04. Note that as the blade length li increases, blade mass mi also increases, because the values of ρi, bi, and hi are constant. Here, it is assumed that the values of the dimensionless blade stiffnesses ki are constant in order to avoid overly complex results. Unstable vibrations [n(ω/2)] appear in the shaded regions. The black circles represent unstable vibrations; the white circles represent stable vibrations with zero amplitude. They are numerically calculated from Eq. (16). Note that Fig. 6 is calculated for ω = 0.05–0.25 because some discrepancies exist between the theoretical and simulation results near ω = 0, and unstable regions do not appear for ω > 0.25. It is found that the theoretical results are mostly in agreement with the numerical results. Unstable vibrations [ω], [(3/2)ω], [2ω], [(5/2)ω], and [3ω] appear when li > 77.4, li > 39.4, li > 39.9, li > 39.9 and li > 40.6, respectively. Note that the unstable vibration [(1/2)ω] does not appear as previously mentioned in Fig. 5. Unstable vibrations [(3/2)ω], [2ω], and [(5/2)ω] appear when li = 40. The horizontal, dotted line (Fig. 8) show the rotational speeds at which resonant peaks are expected to appear in the nonlinear system. As blade length li is increased, the unstable regions where these vibrations appear widen, and unstable vibration [ω] is expected to appear when li > 77.4.

Figure 7 shows the influence of the stiffness of the rotational spring ki on the unstable regions for li = 40. Unstable vibrations [ω], [(3/2)ω], [2ω], [(5/2)ω], and [3ω] appear when 0 < ki < 0.142, 0 < ki < 0.530, 0 < ki < 0.513, 0 < ki < 0.509 and 0 < ki<0.485, respectively. As the stiffnesses of the rotational springs are decreased, the unstable regions where these vibrations appear widen. Figures 6 and 7 suggest that real wind turbines with longer length blades and/or smaller rotational spring stiffnesses may be at a higher risk for breakdowns.

Frequency Response Curves.

Figures 8(a) and 8(b) show the frequency response curves of the translations x0 and z0 for the 2DOF system, respectively; Figs. 8(c)8(e) for blade inclinations θi (i = 1, 2, 3), respectively. These curves are obtained by applying the swept-sine test to Eq. (9), which includes the nonlinear terms in real wind turbines. The values of the wind turbine's parameters are μΑ = 3.2 × 10−8, μT = 0.94, μι = 0.02, kx = kz = 1.0, ki = 0.5, li = 40.0, cx = cz = 0.03, ci = 0.04, bi = 2.0, α1 = 0 deg, α2 = 120 deg, α3 = 240 deg, V0 = 6.0, Vi = 0.6, and λ = ±3.0 × 10−8. The values of li and ki are same as those in Figs. 6 and 7, respectively. Here, λ represents the acceleration of the sine sweep excitation, and ωt in Eq. (9) is replaced by (1/2)λt2+ω0t in the numerical calculations, where ω0 represents the initial rotational speed for the swept-sine test. The maximum amplitudes of the time histories are calculated, and these values are plotted at each rotational speed in Fig. 8, when the rotational speed is either increased or decreased. The solid (black) line represents the results when the rotational speed is increased from 0.05 to 0.20 (λ= 3.0 × 10−8, ω0 = 0.05), and the dashed (red) line represents the results when the rotational speed is decreased from 0.20 to 0.05 (λ= –3.0 × 10−8, ω0 = 0.20). Resonant peaks appear due to wind and parametric excitation near the rotational speeds indicated by A, B, and C, which correspond to the points of the same names in Figs. 57, where unstable vibrations appear in the linearized system. Note that the response curves for the three blades are identical. When the rotational speed is decreased from ω = 0.2, the blades vibrate at increasing higher amplitudes from point A2 and peak at point A3, which exists out of the frame in Figs. 8(c)8(e). Jump phenomena occur, and the amplitudes of the blades jump from point A3 to point A1. The blades then vibrate at extremely low amplitudes until point B2. The amplitudes of the blades increase again and peak at point B3, also out of the frame in Figs. 8(c)8(e). Jump phenomena occur, and the amplitudes of the blades are almost zero at point B1. When the rotational speed is increased from ω = 0.1, the amplitudes of the blades peak at comparatively high amplitudes at point B3 near rotational speed B, and at lower amplitudes at point A3 near rotational speed A. The frequency response curves for the blades near rotational speeds A and B bend to the left, exhibiting softening behavior, because each blade is supported by a rotational spring, and thus behaves like a pendulum. In real wind turbines, starting and stopping them and/or a change in wind velocity affect the rotational speed of the blades. As demonstrated by Fig. 8, a decrease in the rotational speed may cause the blades to vibrate at higher amplitudes.

Response Near Rotational Speed A.

Figures 9(a)9(d) show the stationary time histories of the 2DOF system translations x0 and z0, and blade inclinations θi (i = 1, 2, 3) for different initial conditions at ω = 0.180 near rotational speed A in Fig. 8. Eight different patterns of blade vibrations are possible and categorized according to their characteristics as patterns I, II, III, or IV, in Table 1. In pattern I, all three blades vibrate at high amplitudes. In pattern II, two of the blades vibrate at high amplitudes, and the remaining blade vibrates at low amplitudes, and thus three patterns exist, II-1, II-2, and II-3, depending on which of the three blades vibrates at low amplitudes. In pattern III, two of the blades vibrate at low amplitudes, and the remaining blade at high amplitudes, and thus three patterns also exist, III-1, III-2, and III-3, depending on which of the blades vibrates at high amplitudes. In pattern IV, all three blades vibrate at low amplitudes. Figure 9(a) shows pattern I, Figs. 9(b) and 9(c) show pattern II-1 and III-1, respectively, and Fig. 9(d) shows pattern IV. Figures 9(b) and 9(c) demonstrate “localization phenomena,” i.e., even though the blades are identical, their vibration amplitudes are different. Localization phenomena were also observed in identical multicantilevers, attached to a disk [2124], and in structures with multiple, identical nonlinear dynamic dampers [25,26]. Note that the time histories for x0 in the 2DOF system are different depending on which pattern of the blades appear, as shown in Figs. 9(a)9(d). This is because x0 is parametrically and nonlinearly coupled with the blade inclinations in Eq. 9(a). The time histories for z0 in the 2DOF system are also different in Figs. 9(a)9(d) and are more complicated than those of x0 because Eq. (9b) includes both parametric and nonlinear terms.

Figure 10 shows the FFT results at ω = 0.180 in Fig. 8, corresponding to pattern II-1 in Fig. 9(b). The vertical axes represent the amplitudes of vibrations in logarithmic scales, and the horizontal arrows represent the constant amplitudes. It is found that the frequency (3/2)ω is predominantly included in z0 and θi because unstable vibration [(3/2)ω] appears in the linearized system, Eq. (16), near rotational speed A in Fig. 5. Note that other frequencies are also included in z0 and θi due to the nonlinear terms of Eq. (9).

Figures 11 and 12 show the enlarged frequency response curves near rotational speed A of Fig. 8, and they are obtained by applying the swept-sine test in the same manner as in Fig. 8 when different initial conditions are given at the initial rotational speed ω0 = 0.180. The values of the parameters are the same as those in Fig. 8. Both Figs. 11 and 12 demonstrate the occurrence of localization phenomena in the vibration patters of the blades. Figure 11 is an example of vibration pattern II-1, where blades 1 and 2 vibrate at high amplitudes in Figs. 11(c) and 11(d), respectively, and blade 3 vibrates at low amplitudes in Fig. 11(e). Pattern III-1 is shown in Fig. 12, where blade 1 vibrates at high amplitudes in Fig. 12(c), and blades 2 and 3 vibrate at low amplitudes in Figs. 12(d) and 12(e), respectively.

Depending on the initial conditions, the blades may exhibit vibration patterns II or III, in which localization phenomena occur. The influence of the initial conditions on the system can be investigated by giving a disturbance to one of the blades and calculating BOAs. Figure 13 shows BOAs when disturbances are given to blade 1 near rotational speed A. A disturbance is defined by the coordinate system (θ¯1,θ¯˙1) which represents a deviation from the steady-state solution of blade 1. Disturbances are given to a grid of 200 × 200 starting points in the (θ¯1,θ¯˙1) plane. The white domain represents the BOA where vibration patterns I or IV appear. The gray domain represents the BOA where vibration patterns II or III appear. Figures 13(a) and 13(b) show BOAs at ω = 0.170 and 0.180, respectively. The origin of the coordinate system corresponds to the steady-state solution on branches A1A2, and the BOAs demonstrate point symmetry about the origin. In Fig. 13(a), the white domain is more dominant, and pattern IV is more likely to occur than pattern III, and thus, localization phenomena seldom occur. As the rotational speed increases, as shown in Fig. 13(b), the total area of the gray domain increases, and thus pattern III, i.e., localization phenomena, is more likely to occur. Figures 13(c) and 13(d) show BOAs at ω = 0.170 and 0.180, respectively, and the origin of the coordinate system corresponds to the steady-state solution on branches A2A3 in Fig. 8. In Figs. 13(c) and 13(d), the gray domain is more dominant than in Figs. 13(a) and 13(b), and is most dominant in Fig. 13(c). Thus, pattern II, i.e., localization phenomena, is most likely to occur.

Response Near Rotational Speed B.

Figures 14(a) and 14(b) show the stationary time histories and FFT results, respectively, at ω = 0.116 near rotational speed B in Fig. 8. Eight different patterns of blade vibrations are also possible and categorized in the same manner as those for the responses near rotational speed A. Figure 14 corresponds to pattern III-1: blade 1 vibrates at high amplitudes, while blades 2 and 3 vibrate at low amplitudes. The FFT results show that the frequency component 2ω is predominantly included in the vibrations of all three blades because unstable vibration [2ω] appears in the linearized system, Eq. (16), near rotational speed B in Fig. 5.

Figure 15 shows BOAs at ω = 0.116 when a disturbance is given to blade 1. The origins of the coordinate systems in Figs. 15(a) and 15(b) correspond to the steady-state solutions on branches B1B2 and B2B3 in Fig. 8, respectively. In Fig. 15(a), the white domain, corresponding to pattern IV, is dominant, and thus pattern III, i.e., localization phenomena, seldom occur. On the other hand, the gray domain, corresponding to pattern II, is dominant in Fig. 15(b), thus localization phenomena often occur, and may cause the turbine to breakdown.

The nonlinear vibration characteristics of a three-blade wind turbine have been investigated. The unstable regions, frequency response curves, stationary time histories, FFT results, and basins of attraction were numerically calculated. The results can be summarized as follows:

  1. (1)In the corresponding linear system of the three-blade wind turbine, unstable vibrations of frequencies (n/2)ω (n = 2, 3, 4, …) appear at specific regions of rotational speeds due to parametric excitation. These unstable regions can be used to predict the rotational speeds where resonant peaks may appear in the nonlinear system.
  2. (2)In the nonlinear system, the nonlinear restoring moments acting on the blades are considered, and resonant peaks appear at constant amplitudes in the frequency response curves near the rotational speeds where unstable vibrations appear in the linearized system.
  3. (3)The frequency response curves exhibit softening behavior and bend to the left due to the nonlinear restoring moments acting on the blades. The vibrations of the blades peak at higher amplitudes when the rotational speed is decreased.
  4. (4)Localization phenomena may occur due to the nonlinear restoring moments acting on the blades, even when identical blades are used.
  5. (5)Basins of attraction demonstrate that the occurrence of localization phenomena depends on the magnitude of the disturbance and the rotational speed.

For future work, experiments will be conducted using a wind turbine prototype, and the theoretical results will be compared with experimental data, in order to confirm the validity of the theoretical analysis presented in this paper.

  • bi =

    width of blade “i” (i = 1, 2, 3)

  • ci =

    viscous damping coefficients of blade “i

  • cx, cz =

    viscous damping coefficients of the 2DOF system in the x0- and z0-directions, respectively

  • g =

    acceleration of gravity

  • hi =

    thickness of blade “i

  • ki =

    spring constant of rotational spring which supports blade “i

  • kx, kz =

    spring constants of the 2DOF system in the x0- and z0-directions, respectively

  • li =

    length of blade “i

  • M =

    total mass of the 2DOF system and three blades (= mT + m1 + m2 + m3)

  • mi =

    mass of blade “i

  • mT =

    mass of the 2DOF system

  • Ni =

    nonlinear moment acting on blade “i” due to the wind velocity

  • Ob-xyz =

    moving Cartesian coordinate system (see Fig. 3)

  • p0 =

    natural frequency of the 2DOF system motion in the z0-direction

  • Vi =

    variation of wind velocity at edge of blade “i

  • V0 =

    wind velocity at center of nacelle

  • t =

    time

  • (x0, z0) =

    stationary Cartesian coordinate system (see Fig. 1)

  • αi =

    angular position of blade “i

  • θi =

    inclination angle of blade “i

  • λ =

    acceleration of sine sweep excitation

  • ρA =

    air density

  • ρi =

    density of blade “i

  • ω =

    rotational speed of blade

Ishihara, T. , Yamaguchi, A. , Takahara, K. , Mekaru, T. , and Matsuura, S. , 2005, “ An Analysis of Damaged Wind Turbines by Typhoon Maemi in 2003,” Sixth Asia-Pacific Conference on Wind Engineering (APCWE VI), Seoul, South Korea, Sept. 12–14, pp. 1413–1428. http://windeng.t.u-tokyo.ac.jp/ishihara/paper/2005-5.pdf
Larsen, J. W. , and Nielsen, S. R. K. , 2006, “ Non-Linear Dynamics of Wind Turbine Wings,” Int. J. Non-Linear Mech., 41(5), pp. 629–643. [CrossRef]
Larsen, J. W. , and Nielsen, S. R. K. , 2007, “ Nonlinear Parametric Instability of Wind Turbine Wings,” J. Sound Vib., 299(1–2), pp. 64–82. [CrossRef]
Kondo, H. , 1988, “ Aeroelastic Response and Stability of Wind Turbine Generators: Behavior of the Turbine Blade,” JSME Int. J. Ser. 3, Vib., Control Eng., Eng. Ind., 31(4), pp. 719–726.
Chopra, I. , and Dugundji, J. , 1979, “ Non-Linear Dynamic Response of a Wind Turbine Blade,” J. Sound Vib., 63(2), pp. 265–286. [CrossRef]
Inoue, T. , Ishida, Y. , and Kiyohara, T. , 2012, “ Nonlinear Vibration Analysis of the Wind Turbine Blade (Occurrence of the Superharmonic Resonance in the Out of Plane Vibration of the Elastic Blade),” ASME J. Vib. Acoust., 134(3), p. 031009. [CrossRef]
Ichikawa, T. , 1984, “ Some Consideration on Aeroelastic Instabilities Caused by Coupling Between Propeller-Type Rotor and Its Supporting Structure,” National Aerospace Laboratory, Tokyo, Japan, Technical Report No. TR-804 (in Japanese).
Yamane, T. , 1986, “ Coupled Rotor/Tower Stability Analysis of Horizontal-Axis Wind Turbine,” Trans. JSME, 52(477), pp. 1514–1519 (in Japanese). [CrossRef]
Saravia, C. M. , Machado, S. P. , and Cirtunez, V. H. , 2011, “ Free Vibration and Dynamic Stability of Rotating Thin-Walled Composite Beams,” Eur. J. Mech. A/Solids, 30(3), pp. 431–441. [CrossRef]
Ishida, Y. , Inoue, T. , and Nakamura, K. , 2009, “ Vibrations of a Wind Turbine Blade (Theoretical Analysis and Experiments Using a Single Rigid Blade Model),” JSME J. Environ. Eng., 4(2), pp. 443–454. [CrossRef]
Acar, G. , Acar, M. A. , and Feeny, B. F. , 2016, “ In-Plane Blade-Hub Dynamics in Horizontal-Axis Wind-Turbines,” ASME Paper No. DETC2016-60344.
Chao, C.-P. , Lee, C.-T. , and Shaw, S. W. , 1997, “ Non-Unison Dynamics of Multiple Centrifugal Pendulum Vibration Absorbers,” J. Sound Vib., 204(5), pp. 769–794. [CrossRef]
Ikeda, T. , Harata, Y. , and Ishida, Y. , 2012, “ Unstable Vibrations of a Wind Turbine Tower With Two Blades,” ASME Paper No. DETC2012-71551.
Janetzke, D. C. , and Kaza, K. R. V. , 1983, “ Whirl Flutter Analysis of a Horizontal-Axis Wind Turbine With a Two-Bladed Teetering Rotor,” Sol. Energy, 31(2), pp. 173–182. [CrossRef]
Yamane, T. , Matsumiya, H. , Kawamura, S. , Muzutani, H. , Nii, Y. , and Gotanda, T. , 1992, “ Vibration Characteristics of an Experimental Wind Turbine With a 15m Teetered Rotor,” JSME Int. J. Ser. 3, Vib., Control Eng., Eng. Ind., 35(2), pp. 268–273.
Motoi, H. , 1993, “ Prevention of Parametric Excitation in the Two-Bladed Wind Turbine With Teetering Mechanism,” Asia-Pacific Vibration Conference'93, Kitakyushu, Japan, Nov. 14–18, pp. 450–455.
Ikeda, T. , Takahashi, H. , Harata, Y. , and Ishida, Y. , 2011, “ Vibration Control of Wind-Turbine Towers Using Tuned Liquid Dampers,” JSME Chugoku-Shikoku Branch 49th Conference, pp. 149–150 (in Japanese).
Ikeda, T. , Takahashi, H. , Harata, Y. , and Ishida, Y. , 2011, “ Vibration Suppression of a Wind Turbine Tower Using a Cylindrical Tuned Liquid Dampers,” Dyn. Des. Conf., 2011, p. 537 (in Japanese).
Ikeda, T. , Harata, Y. , Sumida, J. , and Ishida, Y. , 2012, “ Vibration Control of Wind-Turbine Towers by Dynamic Absorbers,” JSME Chugoku-Shikoku Branch 50th Conference, Paper No. 125-1 (in Japanese).
Nayfeh, A. H. , and Mook, D. T. , 1979, Nonlinear Oscillations, Wiley, New York.
Vakakis, A. F. , and Cetinkaya, C. , 1993, “ Mode Localization in a Class of Multidegree-of-Freedom Nonlinear Systems With Cyclic Symmetry,” SIAM J. Appl. Math., 53(1), pp. 265–282. [CrossRef]
King, M. E. , and Vakakis, A. F. , 1995, “ A Very Complicated Structure of Resonances in a Nonlinear System With Cyclic Symmetry: Nonlinear Forced Localization,” Nonlinear Dyn., 7(1), pp. 85–104.
Dick, A. J. , Balachandran, B. , and Mote, C. D., Jr. , 2008, “ Intrinsic Localized Modes in Microresonator Arrays and Their Relationship to Nonlinear Vibration Modes,” Nonlinear Dyn., 54(12), pp. 13–29. [CrossRef]
Olson, B. J. , Shaw, S. W. , Shi, C. , Pierre, C. , and Parker, R. G. , 2014, “ Circulant Matrices and Their Application to Vibration Analysis,” ASME Appl. Mech. Rev., 66(4), p. 040803. [CrossRef]
Ikeda, T. , 2010, “ Bifurcation Phenomena Caused by Multiple Nonlinear Vibration Absorbers,” ASME J. Comput. Nonlinear Dyn., 5(2), p. 021012. [CrossRef]
Ikeda, T. , 2011, “ Nonlinear Responses of Dual-Pendulum Dynamic Absorbers,” ASME J. Comput. Nonlinear Dyn., 6(1), p. 011012. [CrossRef]
Copyright © 2018 by ASME
View article in PDF format.

References

Ishihara, T. , Yamaguchi, A. , Takahara, K. , Mekaru, T. , and Matsuura, S. , 2005, “ An Analysis of Damaged Wind Turbines by Typhoon Maemi in 2003,” Sixth Asia-Pacific Conference on Wind Engineering (APCWE VI), Seoul, South Korea, Sept. 12–14, pp. 1413–1428. http://windeng.t.u-tokyo.ac.jp/ishihara/paper/2005-5.pdf
Larsen, J. W. , and Nielsen, S. R. K. , 2006, “ Non-Linear Dynamics of Wind Turbine Wings,” Int. J. Non-Linear Mech., 41(5), pp. 629–643. [CrossRef]
Larsen, J. W. , and Nielsen, S. R. K. , 2007, “ Nonlinear Parametric Instability of Wind Turbine Wings,” J. Sound Vib., 299(1–2), pp. 64–82. [CrossRef]
Kondo, H. , 1988, “ Aeroelastic Response and Stability of Wind Turbine Generators: Behavior of the Turbine Blade,” JSME Int. J. Ser. 3, Vib., Control Eng., Eng. Ind., 31(4), pp. 719–726.
Chopra, I. , and Dugundji, J. , 1979, “ Non-Linear Dynamic Response of a Wind Turbine Blade,” J. Sound Vib., 63(2), pp. 265–286. [CrossRef]
Inoue, T. , Ishida, Y. , and Kiyohara, T. , 2012, “ Nonlinear Vibration Analysis of the Wind Turbine Blade (Occurrence of the Superharmonic Resonance in the Out of Plane Vibration of the Elastic Blade),” ASME J. Vib. Acoust., 134(3), p. 031009. [CrossRef]
Ichikawa, T. , 1984, “ Some Consideration on Aeroelastic Instabilities Caused by Coupling Between Propeller-Type Rotor and Its Supporting Structure,” National Aerospace Laboratory, Tokyo, Japan, Technical Report No. TR-804 (in Japanese).
Yamane, T. , 1986, “ Coupled Rotor/Tower Stability Analysis of Horizontal-Axis Wind Turbine,” Trans. JSME, 52(477), pp. 1514–1519 (in Japanese). [CrossRef]
Saravia, C. M. , Machado, S. P. , and Cirtunez, V. H. , 2011, “ Free Vibration and Dynamic Stability of Rotating Thin-Walled Composite Beams,” Eur. J. Mech. A/Solids, 30(3), pp. 431–441. [CrossRef]
Ishida, Y. , Inoue, T. , and Nakamura, K. , 2009, “ Vibrations of a Wind Turbine Blade (Theoretical Analysis and Experiments Using a Single Rigid Blade Model),” JSME J. Environ. Eng., 4(2), pp. 443–454. [CrossRef]
Acar, G. , Acar, M. A. , and Feeny, B. F. , 2016, “ In-Plane Blade-Hub Dynamics in Horizontal-Axis Wind-Turbines,” ASME Paper No. DETC2016-60344.
Chao, C.-P. , Lee, C.-T. , and Shaw, S. W. , 1997, “ Non-Unison Dynamics of Multiple Centrifugal Pendulum Vibration Absorbers,” J. Sound Vib., 204(5), pp. 769–794. [CrossRef]
Ikeda, T. , Harata, Y. , and Ishida, Y. , 2012, “ Unstable Vibrations of a Wind Turbine Tower With Two Blades,” ASME Paper No. DETC2012-71551.
Janetzke, D. C. , and Kaza, K. R. V. , 1983, “ Whirl Flutter Analysis of a Horizontal-Axis Wind Turbine With a Two-Bladed Teetering Rotor,” Sol. Energy, 31(2), pp. 173–182. [CrossRef]
Yamane, T. , Matsumiya, H. , Kawamura, S. , Muzutani, H. , Nii, Y. , and Gotanda, T. , 1992, “ Vibration Characteristics of an Experimental Wind Turbine With a 15m Teetered Rotor,” JSME Int. J. Ser. 3, Vib., Control Eng., Eng. Ind., 35(2), pp. 268–273.
Motoi, H. , 1993, “ Prevention of Parametric Excitation in the Two-Bladed Wind Turbine With Teetering Mechanism,” Asia-Pacific Vibration Conference'93, Kitakyushu, Japan, Nov. 14–18, pp. 450–455.
Ikeda, T. , Takahashi, H. , Harata, Y. , and Ishida, Y. , 2011, “ Vibration Control of Wind-Turbine Towers Using Tuned Liquid Dampers,” JSME Chugoku-Shikoku Branch 49th Conference, pp. 149–150 (in Japanese).
Ikeda, T. , Takahashi, H. , Harata, Y. , and Ishida, Y. , 2011, “ Vibration Suppression of a Wind Turbine Tower Using a Cylindrical Tuned Liquid Dampers,” Dyn. Des. Conf., 2011, p. 537 (in Japanese).
Ikeda, T. , Harata, Y. , Sumida, J. , and Ishida, Y. , 2012, “ Vibration Control of Wind-Turbine Towers by Dynamic Absorbers,” JSME Chugoku-Shikoku Branch 50th Conference, Paper No. 125-1 (in Japanese).
Nayfeh, A. H. , and Mook, D. T. , 1979, Nonlinear Oscillations, Wiley, New York.
Vakakis, A. F. , and Cetinkaya, C. , 1993, “ Mode Localization in a Class of Multidegree-of-Freedom Nonlinear Systems With Cyclic Symmetry,” SIAM J. Appl. Math., 53(1), pp. 265–282. [CrossRef]
King, M. E. , and Vakakis, A. F. , 1995, “ A Very Complicated Structure of Resonances in a Nonlinear System With Cyclic Symmetry: Nonlinear Forced Localization,” Nonlinear Dyn., 7(1), pp. 85–104.
Dick, A. J. , Balachandran, B. , and Mote, C. D., Jr. , 2008, “ Intrinsic Localized Modes in Microresonator Arrays and Their Relationship to Nonlinear Vibration Modes,” Nonlinear Dyn., 54(12), pp. 13–29. [CrossRef]
Olson, B. J. , Shaw, S. W. , Shi, C. , Pierre, C. , and Parker, R. G. , 2014, “ Circulant Matrices and Their Application to Vibration Analysis,” ASME Appl. Mech. Rev., 66(4), p. 040803. [CrossRef]
Ikeda, T. , 2010, “ Bifurcation Phenomena Caused by Multiple Nonlinear Vibration Absorbers,” ASME J. Comput. Nonlinear Dyn., 5(2), p. 021012. [CrossRef]
Ikeda, T. , 2011, “ Nonlinear Responses of Dual-Pendulum Dynamic Absorbers,” ASME J. Comput. Nonlinear Dyn., 6(1), p. 011012. [CrossRef]

Figures

Grahic Jump Location
Fig. 2

Analytical model: (a) top view and (b) side view

Grahic Jump Location
Fig. 3

Coordinate system for blade “i

Grahic Jump Location
Fig. 6

Unstable regions showing the influence of blade length li

Grahic Jump Location
Fig. 7

Unstable regions showing the influence of spring stiffness ki

Grahic Jump Location
Fig. 8

Frequency response curves calculated by the swept-sine test when μΑ = 3.2 × 10−8, μT = 0.94, μi = 0.02 (i = 1,2,3), kx = kz = 1.0, ki = 0.5, li = 40.0, cx = cz = 0.03, ci = 0.04, bi = 2.0, α1 = 0 deg, α2 = 120 deg, α3 = 240 deg, V0 = 6.0, Vi = 0.6, and λ = ±3.0 × 10−8. (a) Translation x0 of the 2DOF system; (b) translation z0 of the 2DOF system; and (c)–(e) blade inclinations θi (i = 1, 2, 3).

Grahic Jump Location
Fig. 14

Stationary time histories and FFT results of pattern III-1 at ω = 0.116 in Fig. 13: (a) time histories and (b) FFT results

Grahic Jump Location
Fig. 15

Basins of attraction when blade 1 is disturbed at ω = 0.116: (a) on branches B1B2, (b) on branches B2B3 in Fig. 8. BOA for patterns III and II in (a) and (b), respectively.

Grahic Jump Location
Fig. 4

Definition of wind velocity

Grahic Jump Location
Fig. 5

Natural frequency diagram for the three-blade wind turbine when μT = 0.94, μi = 0.02 (i = 1,2,3), kx = kz = 1.0, ki = 0.5, li = 40.0, bi = 2.0, α1 = 0 deg, α2 = 120 deg, and α3 = 240 deg

Grahic Jump Location
Fig. 9

Stationary time histories at ω = 0.180 in Fig. 8: (a) pattern I; (b) pattern II-1; (c) pattern III-1; and (d) pattern IV

Grahic Jump Location
Fig. 10

Fast Fourier transform results of pattern II-1 at ω = 0.180 in Fig. 9(b)

Grahic Jump Location
Fig. 11

Frequency response curves demonstrating pattern II-1

Grahic Jump Location
Fig. 12

Frequency response curves demonstrating pattern III-1

Grahic Jump Location
Fig. 13

Basins of attraction when blade 1 is disturbed at (a) ω = 0.170 and (b) ω = 0.180 on branches A1A2, and at (c) ω = 0.170 and (d) ω = 0.180 on branches A2A3 in Fig. 8. BOA for pattern III in (a) and (b), and those for pattern II in (c) and (d).

Tables

Table Grahic Jump Location
Table 1 Vibration patterns of blades

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In