1 Introduction

The span of a large scale wind turbine blade is commonly huge, but the blade chord is quite short. Therefore, a beam model is widely used to calculate the structural dynamic behavior of blades^{[ 1- 2]}. A blade is generally composed of complex thin-walled composite structures. It is difficult to meet the requirements of modern wind turbine in the structural response by the simple beam model^{[ 3]}. The finite element method (FEM) with composite materials has been gradually used to model wind turbine blades^{[ 4]}. In the aspect of wind turbine aerodynamics, the blade element momentum (BEM) method has been widely used to predict overall aerodynamic performance of wind turbines. However, the interactions of vortexes in wake can hardly be considered by the method. In this way, it is hard to obtain unsteady aerodynamic forces using the BEM method with high precision. On the other hand, the three-dimensional (3D) rotational flow can be computed by directly solving the Navier-Stokes/Euler equations of the computational fluid dynamics (CFD) method, furthermore, the unsteady aerodynamic performance can be obtained. However, different turbulent models should be chosen carefully in flow simulations according to different flow states, and a huge amount of computation is required. Therefore, the CFD method is not applicable to daily calculations, especially for a great number of operational conditions as the case of wind turbine. By comparison, the prescribed wake method can afford accurate unsteady aerodynamic loads while the amount of calculation is not huge^{[ 5]}. It is suitable for wind turbine aeroelastic simulations involving the coupling calculation of aerodynamics and structural dynamics of wind turbine. The fast coupled calculation method based on the BEM model and beam model has been widely used in the wind turbine daily aeroelastic simulations due to the simple models and a small number of calculations^{[ 6- 8]}. However, the simplifications and assumptions of the BEM model and beam model both reduce the accuracies of aerodynamic and structural performances computational structural dynamics (CSD) . The aeroelastic method based on the CFD method in conjunction with the CSD method can predict the unsteady aerodynamic performance and structural dynamic responses with a higher accuracy and can be used for static aeroelastic analysis^{[ 4]}. However, a large number of grids in the flow region and a large number of elements in the solid region will be needed to describe the flow state and the blade structures due to the structural complexity^{[ 10]}. The amount of computation in aeroelastic simulation is hard to bear, especially in long-time simulation of dynamic aeroelastic simulations.

In this work, the composite blade structure is modeled with 3D shell elements using the FEM, and the aerodynamic loads are predicted using a prescribed wake method. A nonlinear static aeroelastic method of large scale wind turbines is thus established. Furthermore, the system dynamic equations are linearized at the averaged static aeroelastic status, and the corresponding modes are analyzed considering geometrical nonlinearity of the deformed blade. Besides that, dynamic responses of blade are analyzed through deducting the deformations and loads of static aeroelasticity, and a nonlinear dynamic aeroelastic method is established.

2 Prescribed wake method^{[ 5]}

Wake geometry prescription functions are used to describe the wake shape, and a series of vortexes are arranged at vortex filaments in the wind turbine wake to calculate induced velocities. The local velocity of blade airfoil is calculated by adding inflow velocity and induced velocity of wake vortexes. Then, unsteady aerodynamic is computed by the Beddoes-Leishman theory. The unsteady prescribed wake method is described in detail in Ref. ^[5], and its brief introduction is given here.

2.1 Blade model

In order to provide a detailed assessment of variation in aerodynamic forces along the span, each blade is divided intoN span-wise elements that are considered to be aerodynamic independent, as shown in Fig. 1. A control point which is used to calculate wake-induced velocities, and aerodynamic forces are placed at the quarter chord of the midspan of each element. Boundary points which are used to determine the blade loading and wake shape are placed at each side of the quarter chord of blade element and vortex filaments shed at boundary points.

2.2 Wake geometry prescription

Wind energy transfers from wind to wind turbine when air passes through a turbine, which induces a deceleration of air in the streamwise direction. From the consideration of continuity, the wake radius will increase with the decrease in the axial velocity. The deceleration of air begins ahead of turbine and continues until the air reaches a new equilibrium condition in the far field. Consequently, the turbine wake geometry is divided into the near wake region and the far wake region. In the prescribed method, the wake geometry changes in the near wake region and keeps unchanged in the far wake region. Wake prescription functions are adopted to describe the axial development and the radius development of turbine wake. For discretization, a revolution of the rotor is divided intoN steps equally, vortex filaments with the strength $\Gamma $ trail form each blade element boundary and convect downstream to infinity, as shown in Fig. 2. Wake induced velocities of each element are calculated by the Biot-Savart law, and aerodynamic forces are predicted.

2.3 Unsteady aerodynamic model and dynamic stall model

Wind turbines operate in an unsteady aerodynamic environment. Factors, such as the turbulent flow, yaw angle, and wind shear conditions, have significant effects on the effective inflow angle of blade airfoil. Furthermore, the effective inflow angle can also be changed by the blade vibration speed due to blade elasticity^{[ 11]}. The central premise of blade dynamic responses simulation is to predict the unsteady aerodynamic performance accurately. A semi-empirical model of Beddoes-Leishman is adopted to represent the unsteady aerodynamic forces, such as the lift, the drag, and the pitch moment of blade undergoing the dynamic stall. In the model, unsteady effects are modeled by the attached flow, trailing edge, leading edge separation, and dynamic stall. Attached flow behaviors are simulated by superposition of indicial response, circulator response, and impulsive response. In the separation flow, the associated loss of circulation with development of trailing edge separation may lead to nonlinearity behaviors of lift, drag, and pitch moment, and the trailing edge separation is a major factor to be considered. Also, there is a lag in the normal force coefficient with respect to the change in the angle of attack for dynamic conditions. A first-order model is thus adopted to take account the lag effects. Static separation positions of blade airfoil are solved using the Kirchhoff flow model. The dynamic stall is identified by predicting the onset of leading edge seperation, in which a vortex near leading edge of airfoil separates from the upper surface and is transported downstream over the chord generally.

3 Structural analysis 3.1 FEM and geometrical nonlinearity

The large flexible wind turbine blade deforms in the flap-wise, edge-wise and twist-wise directions under tremendous dynamic wind loads. It is assumed that composite materials of blade are linear materials, but the small transmogrification assumption is no longer established. The geometrically nonlinear strain displacement relation under the finite deformation should be considered^{[ 12]}. According to the theory of elastic, the Green nonlinear strain displacement relation is

(1)

where ${{u}_{i, j}}$ is the partial derivative of $u_{i}$ to $x_{j}$. The strain is rewritten as a linear part and a nonlinear part, i.e., $\varepsilon _{ij}= (e_{i, j}+\eta _{j, i}) $, where the linear part is $ e_{i, j}=u_{i, j}+ u_{j, i}$, and the nonlinear part is $\eta _{j, i}=u_{k, i} u_{k, j}$. According to the Kirchhoff theory, the stress-strain relation is $\sigma _{kl}=D_{klij }\varepsilon _{ij}$, then, the stress is rewritten as an updated form of $\sigma _{kl}^{t+\Delta t}=\sigma _{kl}^{t}+\sigma _{kl}^{\Delta t}$, where $\sigma .{t}$ and ${{\sigma }^{t+\Delta t}}$ are stresses at the timet and the time $t+\Delta t$, respectively, and $\sigma .{\Delta t}_{kl}$ is the increment of stress. The system equations are solved in the updated formula by the updated Lagrange theory, and the equation of the integral form at the time $t+\Delta t$ is

(2)

where ${{Q}^{t+\Delta t}}$ is the external loads.

The system dynamic equation in the finite element scheme is

(3)

where $F=\int_{v}{\sigma _{ij}^{\Delta t}\delta e_{ij}^{t}\text{d}v}$, $ K_\mathrm L $ is the linear part of the stiffness matrix, and $ K_\mathrm {NL} $ is the nonlinear part of the stiffness matrix caused by the geometrical nonlinearity.

3.2 Linearization of system dynamic equations

The composite structure of large-scale wind turbine is quite complex, and therefore the amount of computation will be huge if system dynamic equations are solved directly using the FEM model. In order to reduce the computational cost, the system dynamic equations can be simplified with the modal analysis method, and the structure properties of blade are represented as modal shapes and corresponding frequencies. However, if the modal analysis is carried out at the undeformed position, the truncation error of modal method will be too large to represent large displacements and mechanical properties of blade accurately. System dynamic equations are linearized at the averaged aeroelastic status through the linearization technique^{[ 12]}, and the equations are divided into two parts, i.e., an averaged static part and a dynamic part around the static position. Large deformations will be represented in the static part using the FEM, whose computation cost is endurable because this process will be carried out only once. Moreover, dynamic vibrations of blade will be represented in the dynamic part efficiently using the modal method.

The time-dependent external load $ Q$ is divided into the averaged part $ Q_\mathrm s $ and the dynamic part $ Q_\mathrm d $. The nonlinear large deformations of blade $ u$ are divided into the averaged static part $ u_\mathrm s $ and the dynamic part $ u_\mathrm d $ around the static position,

(4)

(5)

Substitute Eqs. (4) and (5) into Eq. (3) . Then, the vibration equations around the static position due to the dynamic part of the external load can be obtained as

(6)

where $ M_\mathrm T$ and $ K_\mathrm T$ are the tangent mass matrix and the tangent stiffness matrix, respectively. The modal analysis is carried out based on the previous equations, and modal shapes and frequencies considering the geometrical nonlinearity are obtained.

4 Aeroelastic analysis 4.1 Static aeroelastic analysis

Bending deformations and torsional deformations of blade will occur when the blade suffers from external wind loads. The effective inflow angle of blade airfoil changes due to large torsional deformations, which will also change wind loads. Thus, the deformations and wind loads of blade affect each other interactively, and aeroelastic effects are formed. The information transfering between blade structure and aerodynamics is performed through the interpolation method. Static aerodynamic forces are calculated using the steady prescribed wake method, and then averaged static forces are obtained by averaging wind loads over the azimuth angle. The multi-point constrain (MPC) method is used to apply the external forces to the blade structure, the structure analysis is performed, and the distribution of twist angle along the span can be obtained with the NASTRAN software. Extract the incremental inflow angle of blade due to elastic deformations and feed back the incremental angle to the process of aerodynamic calculations, then, aerodynamic forces of deformed blade can be obtained. Calculate aerodynamic forces and structure deformations interchangeably until the convergence is achieved. Thus, a nonlinear static aeroelastic computation method is established.

4.2 Dynamic aeroelastic analysis

After the static aeroelastic analysis is carried out, the nonlinear system dynamic equations are linearized at the static position. Modal shapes and frequencies of blade are calculated by solving eigenvectors and eigenvalues of Eq. (6) . Deduct the averaged aerodynamic forces and averaged static deformations. Then, calculate aeroelastic dynamic responses by the loose coupling method^{[ 13]}. The calculating process is shown in Fig. 3. The steps are as follows: calculate the aerodynamic performance and apply the wind loads to the blade structure with the unsteady prescribed wake method first at the time stept; obtain structural dynamic responses with the superposition method of vibration mode at the time step $t+\Delta t$; feed back the vibration speed and the twist angle to the inflow speed and the inflow angle of blade; predict the aerodynamic performance at the time step of $t+\Delta t$. Cycle the previous calculation processes until the end of simulation.

5 Calculation results and discussion

The blade aeroelastic behavior of wind turbine NH1500^{[ 14]} is studied in this paper. The rotor diameter is 83 m, and the blade span is 40.5 m, the blade weight is 5 600 kg. 3D blade structures are modeled using shell elements, and the hub is simplified as the beam model to connect three blades. The total element number of the rotor is 110 000, as shown in Fig. 4.

Aerodynamic loads of each blade section are predicted and then applied to the control point of blade element in the prescribed wake model, as shown in Fig. 5. In this paper, the wind turbine operates at the steady wind with a speed of 11.0 m/s, and the blade set angle is 0°.

In order to verify the method of this paper, distributions of normal force coefficient $F_\mathrm n $ along the blade span are calculated using the presented method and the BEM method, respectively. And the BEM method is carried out by GH Bladed software. The results are compared as shown in Fig. 6. It is shown that the result of this paper is similar to the result of GH Bladed.

The static aeroelastic analysis is performed using the previous method described in Subsection 4.1. Distributions of blade deformations with the geometrical nonlinearity in the flap-wise direction along the blade span are compared with the one without the geometrical nonlinearity, as shown in Fig. 7. It is shown that the deformation of blade tip with the geometrical nonlinearity is as large as 2.3 m, which is 79.1 percent of the one without the geometrical nonlinearity. It implies that the blade geometrical nonlinearity cannot be neglected when a large deformation of large scale wind turbine blade occurs.

Dynamic aeroelastic behaviors of turbine blade with yaw angles equal to 0°, 10°, and 20° are calculated using the previous method described in Subsection 4.2. Structural responses of blade tip with the geometrical nonlinearity are compared with the one without the geometrical nonlinearity, as shown in Figs. 8 and Figure 9. It is shown that the structural vibration amplitude of blade tip tends to be zero when the yaw angle is equal to 0°, which implies that the aerodynamic performance is steady. The vibration amplitude increases as the yaw angle changes from 0° to 20°. The vibration amplitude of blade tip with the geometrical nonlinearity is smaller than the one without the geometrical nonlinearity in yaw conditions. Vibration amplitudes with yaw angles equal to 10° and 20°, respectively, are 61%$\sim$65% of the one without the geometrical nonlinearity.

6 Conclusions

An integrative aeroelastic method containing the static aeroelastic method and the dynamic static aeroelastic method of large scale wind turbine has been developed in the paper. Effects of geometrical nonlinearity on blade static deformations and dynamic responses of the NH1500 wind turbine are investigated. The results show that the geometrical nonlinearity greatly reduces averaged static deformations and vibration amplitudes of blade dynamic responses, and the geometrical nonlinearity should not be neglected in the aeroelastic analysis of large scale wind turbine.