Nonlinear Flap-Wise Vibration Characteristics of Wind Turbine Blades Based on Multi-Scale Analysis Method

Nomenclature

1  Introduction

The blade plays a critical role in the operation of wind energy conversion by wind turbines, and its vibration is of great concern to the related researchers. It is a typical thin-walled composite beam structure, and the low-order vibration modes are mainly characterized by flapping and oscillation during its operation. The energy-rich region of the blade during vibration is mainly concentrated in the first and second orders due to the low rotational speed. Previous research has primarily focused on the impact of different parameters on the vibration characteristics. Some scholars simplify the blade as a rotating rigid body molding blade modal analysis to study the impact of different hub speeds on its vibration characteristics [1]. CFD method to study the vibration performance of the blade at different wind speeds, the centrifugal force produces a nonlinear effect that increases the intrinsic frequency and changes the period and amplitude [2,3]. The blade is simplified as an Euler beam model, which can be solved by using the finite element method, the vibration pattern superposition method, and the direct integration method, but the vibration pattern superposition method for modal analysis of the kinetic equations is only applicable to linear calculations [4]. Hamilton’s principle is used to represent the nonlinear kinematics of the blade; under the internal resonance condition, the equations may be analytically solved by the multiscale technique, and the accuracy of the perturbation method and the numerical solution can be compared [5–7].

A nonlinear Reduced-Order Model (ROM) based on large deformation theory and the Galerkin method is proposed for the vibration characteristics of horizontal axis wind turbine blades considering the coupling of waving and pendulum vibration as well as blade rotational rigid body motion [8]. The Euler-Lagrange control equations for dynamic systems are derived based on Hamilton's principle, and the vibration control equations are solved using the Extended Galerkin Method (EGM) by separating the structural variables in space and time and defining the assumed mode shapes to satisfy the basic boundary conditions [9,10]. The fluid-solid interaction effects are used to linearize the simulation of homogeneous blade in-plane vibrations. Van der Pol oscillators are used to characterize the time-varying characteristics of blade vortices and can be used to simulate fluid-solid coupling processes through their interaction with the blade equations of motion [11–15]. Moreover, transition sets made up of bifurcation and hysteresis sets may be built to examine the nonlinear vibration of the blade by including system factors such as singularity theory and van der Pol damping.

Multiple Scale Analysis (MSM) refers to the decomposition of the motion of a system into ingredients on different time scales, thus revealing the dynamics of the system [16]. Through multi-scale analysis, various dynamic characteristics at different time and spatial scales can be discovered in a two-dimensional nonlinear vibration system, such as periodic vibration, chaotic vibration, and twin phenomena. These characteristics not only help us understand the motion law of the system but also provide important references for the control and optimization of the system [17–19]. When studying the stability of the steady-state solution of blade vibration in 1:1 internal resonance using the multi-scale method, the Routh-Hurwitz criterion can be used [20].

The analysis shown above demonstrates that classical finite element methods are typically used in the study of blade flap-wise direction vibration, whereas analytical and semi-analytical approximation approaches are still seldom used [21,22]. The nonlinear characteristics of the centrifugal force and the nonlinear vibration characteristics caused by excitation and hub speed are therefore investigated by varying the excitation amplitude and frequency based on the periodic characteristics of wind-induced blade action [23]. In this work, the modal function and Galerkin method are used to truncate the blade flap-wise vibration control equation, and the vibration response is solved using multi-scale method analysis and Runge-Kutta to investigate the flap-wise vibration displacement at different hub speeds, excitation amplitude, and excitation frequency fluctuation blade vibration characteristics, revealing the influence of parameters such as hub speed on the blade flap-wise vibration.

2  Structural Equations of Motions

2.1 Derivation of the Vibration Governing Equation

The blade’s characteristic of short chord length and long span-wise length is simplified as being rotated around a hub with a radius of R, at an angular velocity ω rad/s. The density of the blade material is ρ and a cross-sectional moment of inertia is I. Accounting for the nonlinear response attributed to centrifugal forces, the wind loads are transformed into the transverse periodically distributed load Fcos⁡(ωt), where F and ω are the excitation amplitude and circular frequency, respectively. The model of the rotating beam of the blade and the deformation of blade segments are illustrated in Fig. 1.

Figure 1: Blade rotating beam model (left) and blade Micro-segment deformation (right)

Drawing upon the geometric theory of large deformations and the beam model theory of blade segments, and limiting our considerations solely to flap-wise vibration, we can derive the differential equation governing blade vibration. This equation can be expressed as follows:

∂2M∂x2+ρA∂2v∂t2−∂∂x(N∂v∂x)=f(x,t)(1)

In Eq. (1), A, M and N represent the cross-sectional area, bending stiffness, and axial force of the blade, respectively, while v denotes the displacement deformation in the y direction. ∂2M∂x2, ∂N∂x and ∂2v∂x2 signify the second-order and first-order partial derivatives of A, M and v with concerning x. It is expressed by the Kelvin-Voigt equation describing the blade intrinsic structure relationship as:

σ=Eε+ηε˙(2)

The geometric nonlinear relationship is:

ε=12(∂v∂x)2−y∂2v∂x2(3)

The centrifugal force in the section x is represented by N1 and can be written as:

N1=ρΩ2∫xlA(R+x)dx(4)

According to the assumption of the Kelvin-Voigt blade constitutive relationship, the axial force N can be obtained:

N=N1+∫∫Aσdydz=N1+EA2(∂v∂x)2+ηA∂2v∂x∂t∂v∂x=ρΩ2AR(l−x)+ρΩ2A(l2−x2)2+EA2(∂v∂x)2+ηA∂2v∂x∂t∂v∂x(5)

By the substitution of Eq. (5) into Eq. (1), the nonlinear differential equation for blade vibration can be obtained. The formula is as follows:

EIv,xxxx+ρAv,tt−(EAv,xx+ηA(v,xxtv,x+v,xtv,xx)−ρΩ2A(R+x))v,x−(EA2v,x2+ηAv,xtv,x+ρΩ2AR(l−x)+12ρΩ2(l2−x2))v,xx=f(x,t)(6)

In Eq. (6), v,x, v,xx and v,xxxx represent the first, second, and fourth partial derivatives of x, respectively, and v,xxt represents the second-order partial derivatives of x and the first partial derivative of t. Boundary conditions for the fixed support of the leaf root and free boundary conditions at the leaf tips of the blade are expressed as:

v(x,t)|x=0=∂v(x,t)∂x|x=0=0v(x,t)|x=l=∂3v(x,t)∂x3|x=l=0(7)

Assuming the blade cantilever beam model, for a homogeneous equal-section blade free vibration, E, I, A and ρ are constants, f(x,t)=0, then the blade vibration Eq. (5) can be simplified as follows:

EIv,xxxx+ρAv,tt−(EAv,xx+ηA(v,xxtv,x+v,xtv,xx)−ρΩ2A(R+x))v,x−(EA2v,x2+ηAv,xtv,x+ρΩ2AR(l−x)+12ρΩ2(l2−x2))v,xx=f(x,t)(8)

By introducing the dimensionless variable v∗=v/l, x∗=x/l, t∗=ρAEIl2, A1=EIT4ρAl4, B1=ET2ρl2, C1=D1=ρAEIρη, G1=ρAL3Ω2EI, H1=Al32I, J1=ηlAρEI, N1=ρAl4RΩ2EI, Q1=ρAl5Ω22EI. Substituting the dimensionless variable into Eq. (8) and removing (∗), the dimensionless vibration differential equation is obtained as:

v,tt+A1v,xxxx−B1v,xxvx−C1v,xxtvx2−D1v,xtv,xxv,x+G1(R+lx)v,x−H1v,x2v,xx−J1v,xtv,xxv,x−N1R(1−x)v,xx−Q1(1−x2)v,xx=0(9)

2.2 Discretizing Equations Using the Galerkin Method

In engineering practice, the nonlinear factors have less influence on the vibration mode shapes, but more influence on the amplitude and frequency. The original vibration differential equation of the blade is coordinate transformed, and the blade vibration differential equation is transformed into a modal equation expression. With the coordinate transformation, the expression is:

v(x,t)=∑i=1nqi(t)Yi(x)(10)

In Eq. (10), qi(t) is the modal coordinate or main coordinate of the blade, and also the time modulation of the i modal. Employing the modal of a linear system for blade nonlinear vibration characteristics research, the blade vibration differential equations are converted into spatial known, and then the blade flapping vibration characteristics are studied in the time domain. According to the completed finite element model of the composite blade, the Block Lanczos method is used to solve the modal. The full fixed constraint is adopted for the blade root part, and the number of modal extractions is set. The free modal solution of the blade is obtained, and the first three-order flapping vibration modal shapes of the blade are obtained as shown in Fig. 2.

Figure 2: The previous three vibration modes of the blade

Following the assumption of the rotating cantilever beam model of the blade beam the vibration mode function can be selected as [18]:

Yi(x)=cos⁡βixl−cosh⁡βixl+ξ(sin⁡βixl−sinh⁡βixl)(11)

In Eq. (11), ξ=sin⁡(βi)−sinh⁡(βi)cos⁡(βi)+cosh⁡(βi), the transcendental equation 1+cos⁡(βi)cosh⁡(βi)=0 is satisfied by βi, and the first three modal shapes of the blade are obtained by solving the transcendental equation and normalizing their maximum amplitudes, which are shown in Fig. 3.

Figure 3: The previous three vibration modes of the blade

Vibrating system in which the first several modes prevail, the first-order truncation of the blade flapping vibration differential equation is executed utilizing the Galerkin method following the vibration displacement modal hypothesis. Taking selected modal functional and coordinates transformation into the vibration differential equation (Eq. (8)) obtained:

q¨(t)+ω02q(t)+bq2(t)+cq˙(t)q2(t)+dq3(t)=0(12)

In Eq. (12), the coefficients were truncated using the Galerkin method and integrated over the entire length of the blade, resulting in the equation being expressed as:

a=a1,ω0=(a2+a6(R+lx)+a9(1−x)−a10(1−x2))/a,b=(−a3)/a,c=(−a4+a5−a8)/a,d=(−a7)/a,a1=∫0lY2(x)dx,a2=∫0lA1Y(x)Y‴′(x)dx,a3=∫0lB1Y(x)Y″(x)Y′(x)dx,a4=∫0lC1Y″(x)Y′2(x)Y(x)dx,a5=∫0lD1Y′2(x)Y″(x)Y(x)dx,a6=∫0lG1Y′(x)Y(x)dx,a7=∫0lH1Y″(x)Y′2(x)Y(x)dx,a8=∫0lJ1Y″(x)Y′2(x)Y(x)dx,a6=∫0lG1Y′(x)Y(x)dx,a7=∫0lH1Y″(x)Y′2(x)Y(x)dx,a8=∫0lJ1Y″(x)Y′2(x)Y(x)dx,a9=∫0lN1Y″(x)Y(x)dx,a10=∫0lQ1Y″(x)Y(x)dx.(13)

3  Equations and Mathematical Expressions

Multi-scale method analysis conceives representation response features to be expanded as several time scale functions. Setting the results of the differential equations of blade vibration as a measure of multiple independent variables, y is considered as a function of t and ςt. Introducing the independent variable Tn=ςnt,n=0,1,⋯, the derivative concerning time t thus becomes a partial derivative expansion concerning T, i.e., [19].

ddt=dT0dt∂∂T0+dT1dt∂∂T1+⋯=D0+ςD1+⋯d2dt2=D02+2ςD0D1+ς2(D12+2D0D1)+⋯(14)

In Eq. (15), Dk=∂/∂Tk⋯(k=0,1). The differential equation for blade vibration is solved with t substituted, and the form of the equation solution can be written as:

q=q0(T0,T1,T2,⋯)+ςq1(T0,T1,T2,⋯)+ς2q2(T0,T1,T2,⋯)+⋯(15)

Based on perturbation, the vibration of the blade is studied, and it is noted that ω=(1+γ)ω0, γ is the parameter for the frequency mismatch of stimulation. By substituting Tn=ςnt,n=0,1,⋯ and Eq. (15) into Eq. (12) and comparing the coefficients of the same order on both sides of the equation, we obtain a set of approximate partial differential equations for each order, expressed as:

ς0:D02q0+ω02q0=0ς1:D02q1+ω02q1=−2D0D1q0ς2:D02q2+ω02q2=−2D0D1q1−D12q0−2D0D2q0−bq02−cq02D0q0−dq03(16)

By solving Eq. (16) can be obtained:

q0=Z(T1,T2)eiω0T0+Z¯(T1,T2)e−iω0T0(17)

By substituting Eq. (18) into Eq. (16) and eliminating the long-term term, can be obtained:

D1Z=0,q1=eiω0T0+e−iω0T0(18)

Substituting Eqs. (17) and (18) into Eq. (16), the expression is obtained:

D02q2+ω02q2=−2iω0D2Zeiω0T0−D12Zeiω0T0–bZ2e2iω0T0−2bZZ¯−ciω0Z3e3iω0T0−2ciω0Z2Z¯iω0T0−cZ2Z¯iω0eiω0T0−dZ3e3iω0T0−3dZ2Z¯iω0T0+cc(19)

In Eq. (19), cc is the long-term term, eliminating the long-term term can be obtained:

−2iω0D2Z−D12Z−2ciω0Z2Z¯−cZ2Z¯iω0−dZ3e3iω0T0−3dZ2Z¯=0(20)

To solve for Z, let its complex form be expressed as:

Z=12αeγi(21)

Eq. (21) involves functions α and γ concerning time, and the derivative of Z concerning time can be obtained:

Z˙=12α˙cos⁡γ−12αγ˙sin⁡γ+i(12α˙sin⁡γ+12αγ˙sin⁡γ)(22)

Substitute Eq. (22) into Eqs. (21) and (23) for comparison, separate the real part and imaginary part, the blade nonlinear vibration characteristics for stable solutions, then there is a˙=0 and γ˙=0. Hence, the derived blade slowly varies amplitude and variable phase to satisfy the relationship that is the blade bifurcation response equation as follows:

(−38da3−aω0σ)2−a2ω02=(14ω0)2(23)

Nonlinear terms in the vibration control equations consist of multiple modal time-domain functions, and the coefficients are integrals containing multiple modal forms, which require more complicated calculations. The main resonance response bifurcation response equations of the blade nonlinear vibration system are solved with numerical methods, by giving a value that both cos⁡γ and sin⁡γ are satisfied in the interval [−1, 1]. Then, the first-order approximate solution of the system can be expressed as:

q=αcos⁡(ω0t+γ)(24)

In Eq. (24), the amplitude is denoted by α and the phase difference is denoted by γ.

4  Simulation Results and Discussion

The blade vibration equation is solved by the numerical method of Runge-Kutta to study the nonlinear flap-wise vibration characteristics of the NREL 5 MW wind turbine blade. Where the wind turbine hub radius R=1.5 m, the total length of the blade is 61.5 m, the density ρ=1.25 kg/m3, and the modulus of elasticity E=39 GPa. With hub speed as the variable, vibration characteristics of the blade at mini rotating speed of 4.1 rpm, rated rotating speed of 12.1 rpm, and hub rotating speed of 22.1 rpm are studied according to the actual working condition of the blade. Blades vibration dimensionless displacement response with hub to rated speed is shown in Fig. 4.

Figure 4: The vibration displacement response curve of the unexcited hub at a speed of 12.1 rpm

In Fig. 4, the blade vibration displacement response curve is shown at a hub speed of 12.1 rpm, and the absolute value of the dimensionless displacement amplitude varies in the range of 0.127–0.72, with a minimum amplitude of 0.127. The overall vibration dimensionless displacement response exhibits periodic changes, and the nonlinear influence leads to a local small-scale response fluctuation in the response curve. Phase diagrams showing the evolution trajectory of different states over the vibration system of the blade are shown in Fig. 5.

Figure 5: The phase diagram when the unexcited hub speed is 12.1 rpm

In Fig. 5, the motion state is a closed trajectory, indicating that the blade vibration is a stable motion. The two ends of the dx axis are consistent with the vibration displacement response value at the hub speed of 12.1 rpm, the leftmost is −0.72, the rightmost is 0.72, and the middle change value is the displacement response change value change area is −0.72–0.72. dx/dt changes response range is −0.65–0.65, and the points in the phase diagram are the regular points. The existence of nonlinear terms will make the relationship curve between amplitude and frequency show special patterns such as multi-peaks and distortions, and the relationship between amplitude and frequency of different wheel speeds is shown in Fig. 6.

Figure 6: Relationship between blade vibration amplitude and frequency response under different wheel speeds

The amplitude-frequency response curves of blade-free vibration under different wheel speeds are shown in Fig. 6. The multi-peak phenomenon varies with different speeds, and peaks occur simultaneously at frequencies of 0.159 and 0.34 for all speeds. When the frequency is 0.159, the maximum amplitude of the low speed is 16.41, followed by the rated speed with an amplitude of 14.93, and the minimum amplitude is at the high speed, which is 14.46. Points D, E, and F have peaks at a frequency of 0.34, with the first peak at 5.19 for the high speed, the second at 5.16, and point F is the peak for the low-speed amplitude.

FEM simulation software is used to solve the Campbell diagram analysis of blade vibration. The constraint conditions and loads are set to be exactly the same as the numerical simulation, and the finite element simulation results are compared with the numerical settlement results. The Campbell diagram of blade vibration is drawn as shown in Fig. 7.

Figure 7: Campbell diagram of blade vibration

It can be seen in Fig. 7 that the two vertical lines are the cut-in speed and the rated speed of the wind wheel, respectively. The ray starting from the origin represents the rotation frequency of the blade, and the two horizontal lines are the first and second-order frequencies of the blade vibration. The rated speed of the wind wheel is 12.1 rpm, and the blade rotation frequency 1P is 0.202 Hz and 3P is 0.606 Hz. The first-order natural frequency of the blade is 0.887 Hz, which is 77.5% and 32.5% different from the rotation frequency 1P and the rotation frequency 3P of the blade. The first-order and second-order natural frequencies of the blade do not intersect with the blade rotation frequency in the wind wheel cut-in speed and rated speed range, that is, the blade does not resonate in the actual working process. That the curve direction and trend are in general agreement with the numerical results, of the same magnitude, indicates that the results are credible.

When the blades are rotated with the rotor hub, the effect of wind on the blades is shown as periodic changes, which can be simplified as periodic load Fcos⁡(ωt). Taking it the same period as that of hub rotation, the load amplitude varies with wind speed. As a matter of engineering practice, different excitation amplitude F is taken under the hub rotation speed of 12.1 rpm to survey the influence of excitation on blade vibration. The vibration displacement response curve of the blade at the excitation amplitude F = 1.23 N is shown in Fig. 8.

Figure 8: The vibration displacement response curve when the excitation amplitude F = 1.23 and the hub speed is 12.1 rpm

In Fig. 8, with the blade loaded periodic load and hub speed of 12.1 rpm, the maximum magnitude of the dimensionless displacement of flapping vibration is 3.16, with an increased magnitude of 2.44 compared to that without considering wind load, which is a 71% increase. It is hard to visualize the movement space evolution trajectory from the vibration displacement response curve, and the phase trajectory after the introduction of wind load was visualized for the system movement space evolution as shown in Fig. 9.

Figure 9: The phase diagram when the excitation amplitude F = 1.23 and the wheel speed is 12.1 rpm

In Fig. 9, the equilibrium position and amplitude nonlinear calculation results are larger than the results when wind load is not considered. dx axis values at both ends increase from −0.72 to −3.53, the maximum value changes from 0.72 to 3.53, and the middle change value is the displacement response change value. dx/dt changes response range, −3.16 to 3.16. With the rated rotational speed of the hub of 12.1 rpm, the effect of different excitation amplitudes on the nonlinear vibration of the blade is investigated, and the relationship between the amplitude of the frequency of the period is shown in Fig. 10.

Figure 10: The amplitude-frequency response relationship of different F values when the wheel speed is 12.1 rpm

In Fig. 10, as the hub is at rated speed, the amplitude-frequency response is bimodal, and the excitation amplitude has a remarkable influence on amplitude so that when the excitation amplitude increased with 1, the corresponding vibration amplitude response increased 10 times. The excitation amplitude has little effect on the frequency. The first peak frequency in the rated speed in Fig. 9 is consistent with that without external excitation in Fig. 6, which is 0.159 Hz. The second peak has a slight change, which is 0.34 without excitation. The frequency is reduced to 0.279 under external load excitation.

When the external excitation load frequency fluctuated minutely, both the frequency and amplitude for the nonlinear vibration dimensionless displacement response varied. Adopt excitation loads fundamental frequency and hub-rated speed consistent at 12.1 rpm, ω=(ω0+γ)t. When γ=0.01, γ=0.1, and γ=1, the non-linear vibration dimensionless displacement response is as shown in Fig. 11.

Figure 11: Displacement response curves with different values of γ when the wheel speed is 12.1 rpm

In Fig. 11, the vibration displacement response of the blade changes obviously with the increase of γ, and the maximum amplitude decreases significantly. When γ=0.01, the maximum amplitude of the vibration response is 2.27, and the maximum amplitude of γ=0.1 is 1.69, which is 26% lower than that of γ=0.01. When γ=1, it is reduced to 0.61, which is 73% lower than that of γ=0.01, and the slight fluctuation of the external excitation load frequency has little effect on the response period.

To visually observe the spatial state of blade vibration motion changing with the fluctuation of external excitation frequency, the spatial phase diagram at γ=0.1 is compared and analyzed with the phase diagram under the wind load at a hub rotation speed of 12.1 rpm (See Fig. 7). The evolutionary trajectory of the motion space at this time is shown in Fig. 12.

Figure 12: Wheel speed of 12.1 rpm and phase diagram at γ=0.1

In Fig. 12, it can be observed by comparing Figs. 7 and 11 that when there is a fluctuation in the excitation frequency, there is a significant change in the spatial evolution trajectory of blade vibration, and the maximum response amplitude decreases from 3.53 to 2.27, which is reduced by 36%. The maximum value of dx/dt also decreases from 3.16 to 2.04, which is a reduction of 35.4% compared to when there is no fluctuation. Furthermore, significant changes in the evolutionary trajectories of many movements can be observed in the intermediate region of the phase diagram, as shown in Figs. 9 and 10.

In the phase space, the displacement dy and the velocity dy/dt are regarded as the conjugate variables to calculate the local maximum of the steady state process. The characteristics and variations of the system motion are presented by intercepting the Poincare cross-section to explore the trajectory of the evolution of the blade vibratory motion. The Poincare section is visualized by taking the rated speed of 12.1 rpm and no external load excitation, the external excitation amplitude of F = 10.12 but no fluctuation of the external excitation frequency, the external excitation amplitude of F = 10.12 and the external excitation frequency fluctuation of a = 0.1, as shown in Fig. 13.

Figure 13: Poincare section plot

In Fig. 13, it can be determined whether the vibration is chaotic. It can be seen from Fig. 11 that when there is no external excitation on the blade, the maximum value of dy is 0.753 and the maximum value of dy/dt is 2.5 × 10−3. As external load excitation is introduced, there is a 57.67% increase in the maximum value in dy, which is 1.779, and the maximum value of dy/dt is increased to 8.6 × 10−3, with an increased ratio of 70.93%; the frequency fluctuation a = 0.1 of the external load excitation increases the dy/dt value by 0.077 and decreases the dy/dt value by 0.410−3 compared to no frequency fluctuation. All of the discrete points on the cross-section are minority, without obvious fractal structure, and the blade vibration is periodical movement.

5  Conclusions

In this work, the wind turbine blade is simplified as a cantilever beam spinning around the hub. Following the finite element simulation analysis, the suitable vibration mode function is chosen, shortened using the Galerkin method, solved using the Runge-Kutta numerical method, and compared to the finite element simulation findings. The study’s findings indicate that:

(1) There is a significant variable effect of hub speed change in the blade vibration amplitude with time. Various periphery speed blade free vibration amplitude frequency response is separate, although, at the frequency of 0.159 and 0.34, each speed appears at the simultaneous peak; frequency is 0.159, minor speed amplitude is the maximum, then the rated speed amplitude, magnitude at high speed is minimum, respectively, 16.41, 14.93, 14.46.

(2) At rated speed, when harmonic external excitation is introduced, the maximum amplitude value becomes 3.16, with an amplification ratio of 71%. Amplitude frequency response is bimodal, and the excitation amplitude influence on that amplitude is notable, when the excitation amplitude increased by 1, the corresponding vibration amplitude response increased 10 times. The frequency of the first peak is consistent with that without external excitation, which is 0.159 Hz, and the second peak has a slight change, decreasing from 0.34 without excitation to 0.279 with a reduced rate of 18%.

(3) As the external excitation frequency of the load fluctuates, the blade vibration spatial evolution trajectory alters, and the maximum value of the response decreases from 3.53 to 2.27, which is reduced at a ratio is 36%. The change in the motion evolution trajectory can be easily seen from the middle area of the phase diagram.

Acknowledgement: The authors would like to acknowledge the generous help from the School of Mechanical and Electronic Engineering and the College of Automation.

Funding Statement: This work was supported by the National Natural Science Foundation of China (No. 51965034).

Author Contributions: The authors confirm their contribution to the paper as follows: study conception and design: Yuqiao Zheng, Qifa Lang; data collection: Qifa Lang; analysis and interpretation of results: Qifa Lang, Chenglong Shi; draft manuscript preparation: Tiancai Cui, Heyu Zhang. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: All data generated or analyzed during this study is included in this published article.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Sun, R., Li, C., Ding, Q., Chen, W. P. (2018). Modal and vibration analysis for the large wind turbine blade. Journal of Engineering for Thermal Energy and Power, 33(10), 400–407. [Google Scholar]

2. Algolfat, A., Wang, W. Z., Albarbar, A. (2022). Comparison of beam theories for characterisation of a NREL wind turbine blade flap-wise vibration. Proceedings of the Institution of Mechanical Engineers Part A: Journal of Power and Energy, 236(7), 1350–1369. [Google Scholar]

3. Karimi, B., Hamed, M. (2018). Nonlinear kinematics analysis and internal resonance of wind turbine blade with coupled flap-wise and edgewise vibration modes. Journal of Sound and Vibration, 116, 390–408. [Google Scholar]

4. Ju, D. Y., Sun, Q. (2017). Modeling of a wind turbine rotor blade system. Journal of Vibration and Acoustics, 139(5), 1–15. [Google Scholar]

5. Evgueni, S., Sudhanva, K. C. (2019). Determination of natural frequencies and mode shapes of a wind turbine rotor blade using Timoshenko beam elements. Wind Energy Science, 4(1), 57–69. [Google Scholar]

6. Behnam, K., Hamed, M. (2018). Nonlinear kinematics analysis and internal resonance of wind turbine blade with coupled flap-wise and edgewise vibration modes. Journal of Sound and Vibration, 435, 390–408. [Google Scholar]

7. Jokar, H., Mahzoon, M., Vatankhah, R. (2020). Dynamic modeling and free vibration analysis of horizontal axis wind turbine blades in the flap-wise direction. Renewable Energy, 146, 1818–1832. [Google Scholar]

8. Jokar, H., Vatankhah, R., Mahzoon, M. (2021). Nonlinear vibration analysis of horizontal axis wind turbine blades using a modified pseudo arc-length continuation method. Engineering Structures, 247, 1103–1131. [Google Scholar]

9. Farsadi, T. (2022). Variable thickness thin-walled rotating blades made of functionally graded porous materials. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 236(14), 7674–7689. [Google Scholar]

10. Farsadi, T. (2021). Enhancement of static and dynamic performance of composite tapered pretwisted rotating blade with variable stiffness. Journal of Vibration and Acoustics, 143(2), 021009. [Google Scholar]

11. Wang, D., Hao, Z. F., Chen, Y. S., Zhang, Y. X. (2018). Dynamic and resonance response analysis for a turbine blade with varying rotating speed. Journal of Theoretical and Applied Mechanics, 56(1), 31–42. [Google Scholar]

12. Navadeh, N., Goroshko, I., Zhuk, Y., Etminan, M. F., Soleiman, F. A. (2021). Finite element analysis of wind turbine blade vibrations. Vibration, 4(2), 310–322. [Google Scholar]

13. Huang, J. D., Xia, H. J., Li, D. Y., Guo, K. X. (2020). Nonlinear aeroelastic modal analysis of large wind turbine flexible blades. Journal of Mechanical Engineering, 56(14), 174–182. [Google Scholar]

14. Li, Z. G., Song, C. Y., Gao, Z. Y., Zhang, L. R., Wang, J. W. (2017). Development of aeroelastic coupling numerical model for large-scale megawatt wind turbine and its application. Proceedings of the CSEE, 1–13 (In Chinese). [Google Scholar]

15. Yin, Y. J., Liao, M. F., Lv, P. (2017). Analysis of transient respnse of wind turbine blade based timsheko beam. Taiyangneng Xuebao/Acta Energiae Solaris Sinica, 38(9), 2552–2558 (In Chinese). [Google Scholar]

16. Rabiee, F. H., Jafari, A. A. (2022). Nonlinear forced vibration of rectangular plates by modified multiple scale method. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 236(9), 4621–4630. [Google Scholar]

17. Li, Z. G., Cheng, Y. Y., Gao, Z. Y., Zhang, L. R., Wang, J. W. (2022). Numerical model of structural dynamic response for wind turbine flexible blad. Journal of Engineering Thermophysics, 43(9), 2374–2381. [Google Scholar]

18. Luo, R. (2021). Generalized formulation for free vibration of elastic solids with static loads and application to rotating tapered cantilever beam vibration. Journal of Vibration Testing and System Dynamic, 5(2), 149–168. [Google Scholar]

19. Farzaneh, R., Asghar, A. J. (2022). Nonlinear forced vibration of rectangular plates by modified multiple scale method. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 236(9), 4621–4630. [Google Scholar]

20. Malte, R., Ould, M. E. (2020). A numerical method to compute global resonant vibrations of ships at forward speed in oblique waves. Applied Ocean, 14(1), 95–166. [Google Scholar]

21. Waseem, M., Ajmal, M., Kim, T. W. (2015). Development of a new composite drought index for multivariate drought assessment. Journal of Hydrology, 527, 30–37. [Google Scholar]

22. Mohanraj, M., Pandiyarajan, R., Jayakumar, K. S., Antony, P. D., Sabarish, S. et al. (2022). Determination of vibration analysis in single and double cracked cantilever numerical beam. Materials Today: Proceedings, 52, 462–466. [Google Scholar]

23. Liang, B., Zhang, T. T. (2018). Numerical optimization and cyber-physical- social computing for vibrations of the elliptical treadmill based on GSO-BPNN model. IEEE Access, 6, 10886–10895. [Google Scholar]