Sensitivity Study of the Influence of Blade Sectional Stiffness Parameters on the Aeroelastic Response of Wind Turbines

With the development of wind turbines as a result of large-scale and offshore trends, the wind turbine size is becoming increasingly larger. The passive control technique is used to alleviate the increasing loads on the blade for the sake of improving the durability of the wind turbine. The ply design of shells considering the coupling effect of bending and torsion is one of the passive control techniques. The bending torsion coupling stiffness is one of the parameters of the blade section stiffness matrix. In order to fully understand the influence of each blade stiffness parameter on the aeroelastic responses of wind turbines and to consider the influence of structural characteristics on the aeroelastic responses in blade design, the influences and sensitivity of each stiffness parameter in the 6 × 6 stiffness matrix of the blade sections on the aeroelastic responses of the wind turbines are systematically studied under steady wind condition. The aerodynamic forces in the aeroelastic model are calculated by an AeroDyn module based on blade element momentum theory, and the structural dynamic responses of the blade are calculated using generalized Timoshenko beam theory and geometric exact beam theory. The NREL baseline 5 MW wind turbine and blade properties are used in this study, where the diagonal stiffness parameters and non-diagonal stiffness parameters of the matrixes of each blade section are scaled according to certain principles. The results show that the axial stiffness, the flap-wise stiffness, and the torsional stiffness in the diagonal are sensitive to the root loads and tip displacement of the blade. The flap-wise bending torsion coupling stiffness, the flap-wise shear-torsion coupling stiffness, and the edge-wise shear-torsion coupling stiffness in the non-diagonal are also sensitive to the aeroelastic responses. For completeness, the effects of other stiffness parameters on the aeroelastic responses are also analyzed and discussed.


INTRODUCTION
In order to reduce the cost of energy (COE) in the wind industry, the wind turbine tends to be designed to be larger and larger. Long blades mean high slenderness ratios, heavy weight, and high extreme load and fatigue load. In order to meet the needs of high stiffness and long fatigue life of the blade, the design method and structure optimization of the continuously upsized blade are both challenging. The most direct way to achieve this goal is to reduce the load on the wind turbine blade in an aeroelastic response. The popular load control methods can be divided into two categories: active control technology and passive control technology. Because the passive control technology does not need driving equipment and has the advantages of fast response and low cost, it is a supplement to the active control method. The ply design of shells considering the coupling effect of bending and torsion is one of the passive control techniques. When the blade is bent, the coupling effect of bending and torsion is used to change the angle of attack to reduce the load.
In order to analyze the influence of blade flexibility and bending torsion coupling effect, and to meet the design requirements of the flexible blade, advanced modeling technology is needed. The CFD/CSD-coupled wind turbine aeroelastic response analysis method may be the development trend of the future, but at this stage, it is not widely used in engineering due to the limitation of computing resources. The aeroelastic analysis model that combines blade element momentum theory with beam theory is the most widely used in engineering. Based on the assumption of small deformation, linear beam theory is gradually replaced by nonlinear deformation beam theory in the analysis of large deformation flexible blades owing to the calculation accuracy (Kim et al., 2013;Gebhardt and Roccia, 2014). Based on the variational asymptotic method (Berdichevskii, 1979) and the conservation of strain energy, 3D analysis of slender beam structure can be decomposed into nonlinear beam theory and section properties analysis theory, which can provide a complete 6 × 6 stiffness matrix. Larsen et al. (2004) studied the influence of large deformation of the blade on the power of wind turbine. His research results show that the main influence caused by large deformation is to change the effective area of the wind turbine, which can reduce the power under low wind speed conditions, and also reduce the pitch angle under high wind speed conditions, and it can improve the bending moment and reduce the thrust of the wind turbine. Ahlström (2006) studied the influence of different slenderness ratios of a 2 MW blade on the wind turbine according to the same idea. The method of changing the slenderness ratio is to scale the stiffness and mass properties of the blade section equivalently. His research shows that when the tip deformation in the direction of blade waving reaches and exceeds 10% of the blade length, the power of the fan will decrease significantly. Then, Ponta et al. (2016) used the NREL 5 MW to study the influence of wind turbine deformation. The results show that even in the case of small deformation, in order to maintain stable power in the aeroelastic response, a smaller pitch angle is needed when the wind speed is higher than the rated wind speed.
The passive load control technology is realized by the bending torsion coupling effect of the composite blade. The bending torsion coupling effect can be realized by material bending torsion coupling or geometric bending torsion coupling of composite blades (Ashwill et al., 2010;Larwood et al., 2014). The coupling effect of bending and torsion can be reflected in the stiffness matrix of the blade section. Lobitz and Veers (1998) of SNL are the first to study adaptive blade passive load shedding.
Their research shows that the bending torsion coupling has an obvious effect on the fatigue load reduction of the wind turbine. Hayat and Ha (2015) further studied the relationship between the bending torsion coupling degree and the reduction of fatigue load. The bending torsion coupling degree is realized by changing the ply angle and the asymmetry of the ply. Gozcu et al. (2014) also show that the bending torsion coupling effect of the blade can reduce the fatigue load of the wind turbine. Vesel and McNamara (2014) studied the optimal design of blade based on aeroelastic response and bending torsion coupling effect with minimum COE as the optimization objective.
It can be found that there are abundant research results on the influences of large deformation and bending torsion coupling effect of the blade, but it is noted that the flap-wise stiffness and bending torsion coupling stiffness of flexible blade only account for a part of the complete composite beam section stiffness matrix, as shown in Figure 1 (coupling deformation is illustrated by the plate in the figure). However, the effects of other coupling terms on the aeroelastic response are rarely studied.
In order to systematically analyze the influence of each stiffness coefficient in the section stiffness matrix of the composite blade on the load and structural response in the aeroelastic response analysis of wind turbine, this paper takes the NREL 5 MW wind turbine as the benchmark model to adjust the parameters of the blade's stiffness matrix. Based on the BEM method and Aerodyn module, the aeroelastic response of the wind turbine and the load of each component are obtained. The influence of each stiffness coefficient on the aeroelastic response of the wind turbine is obtained by using the Sensitivity Analysis Factor (SAF). The results of this paper reveal the sensitivity of blade stiffness to wind turbine aeroelastic response and load and improve the understanding of wind turbine coupling mechanisms. In the future aeroelastic simulation and blade design, the influence of coupling stiffness should be considered actively to make more effective use of the influence of favorable factors and avoid the influence of unfavorable factors so as to provide theoretical support for the formation of more advanced wind turbine blade design method.

MODEL DESCRIPTION
The aeroelastic simulation analysis model of wind turbine combines blade element momentum (BEM) theory and nonlinear beam theory based on geometrically exact beam theory (GEBT). BEM theory has high calculation efficiency and stable results, and the accuracy can meet the needs of engineering design. The GEBT, also known as the generalized Timoshenko beam theory, can simultaneously consider the large flexible deformation of the blade and the complete 6 × 6 stiffness matrix of the blade sections.

Aerodynamic Model
After years of development and revision, BEM theory is relatively mature. There are many studies about BEM theory; for example, the complete derivation process is discussed in many works related to wind turbine design (Manwell et al., 2009;Tongguang Wang and Qian., 2019). BEM theory consists of two parts: momentum theory and blade element theory. In momentum theory, it is assumed that the wind turbine is a penetrating disk and the flow is stable. The wind turbine absorbs the kinetic energy in the wind and converts it into mechanical energy. Blade element theory, also known as strip theory, simplifies wind turbine blades into a finite number of blade segments superimposed along the radial direction. It is assumed that the flow between each blade element does not interfere with each other, and the aerodynamic force acting on each blade element is only determined by its airfoil and local inflow velocity. The blade element momentum method is to solve the momentum theory and blade element theory simultaneously to obtain the aerodynamic load of each blade element and then obtain the aerodynamic performance of the whole blade.
In order to ensure the accuracy of BEM theory, many scholars have put forward a variety of correction methods for the original theory after years of research. For example, tip loss correction (Glauert, 1935), hub loss correction (Glauert, 1935), correction of large thrust coefficient (Glauert, 1926;Buhl, 2005), skewed wake correction (Pitt and Peters, 1980), etc. When the wind turbine operates in an unsteady state, the angle of attack of the wind turbine section will be unsteady. If the unsteady flow separation occurs, the obvious unsteady aerodynamic load will be caused. The B-L dynamic stall model is adopted, which is based on Leishman and Beddoes (1989). In fact, the rotation of the blade will produce Coriolis force and centrifugal force and then cause stall delay phenomenon, resulting in the so-called threedimensional rotation effect. Based on different theories, scholars have developed a series of three-dimensional rotation effect correction models (Corrigan and Schillings, 1994;Du and Selig, 1998;Chaviaropoulos and Hansen, 2000). Most of these models adopt the following forms to modify the aerodynamic force: where C l,2D and C d,2D are the two-dimensional lift and drag coefficient of the airfoil. C l,3D and C d,3D are the corrected three-dimensional lift and drag coefficient. The second term at the right end of the formula is the quantity to be determined by the modified model. Both dynamic stall and three-dimensional rotation effects lead to stall delay, which is not a linear superposition relationship but a certain coupling relationship (Guntur et al., 2016;Elgammi and Sant, 2017;Zhu et al., 2019).

Structural Dynamic Model
The nonlinear beam model based on geometrically exact beam theory is used in structural dynamics. According to GEBT, the position vector of any point on the beam with initial bending and torsion can be described as follows: where R is the position vector of the point on the reference line of the deformation state; C Bb is the rotation tensor between the deformed state and the undeformed state; ξ is the position vector of any point on the beam section in the undeformed state; w is the warping vector of the beam section. Based on the definition of the position vector, the kinetic energy K and the strain energy U of the beam can be obtained by using the Hamilton principle as follows. The detailed derivation process can be found in Hodges' literature (Hodges, 1990).
where δ is Lagrangian variational operators; δW is the virtual work of external force per unit length of the beam; δA represents the corresponding boundary conditions. Using the definition of generalized strain and velocity, K, U, and δW can be described as where V B and Ω B represent the linear and angular velocities of the beam reference line; γ and κ represent the generalized strain; δq B is virtual displacement; δψ B is virtual torsion angle; f B and m B are distributed forces and moments per unit length of beam; I and S are mass matrix and stiffness matrix of beam section; and these can be complete 6 × 6 matrices. The nonlinear governing equations are solved by the Newton-Raphson method, and the linearized equations described by incremental expressions are obtained. The linear equations are discretized by the finite element method, and the discretized linear governing equations are obtained as follows: where M, G , K, and F represent generalized mass matrix, generalized damping matrix, generalized stiffness matrix, and generalized force matrix respectively. The ordinary differential equations can be solved by the method like generalized-α method.

Stiffness Adjustment Strategy
The reference model of the NREL 5 MW wind turbine is taken as the baseline model. Detailed parameters of the baseline model can be found in the technical report (Jonkman et al., 2009). The GEBT can consider a complete 6 × 6 stiffness matrix. By setting different stiffness coefficients in the stiffness matrix and then calculating the aeroelastic response of the wind turbine, the corresponding influence of each stiffness coefficient on the aeroelastic response of the wind turbine can be obtained. In order to intuitively analyze the sensitivity of stiffness to aeroelastic response and avoid the influence of other factors, the parameters other than the stiffness matrix are fixed. The stiffness coefficient is adjusted in the following two ways.
For the main diagonal of the stiffness matrix, the stiffness factor is set as follows: where K ii ′ represents the stiffness coefficient on the diagonal of the stiffness matrix of the comparison model; K ii represents the stiffness factor on the diagonal in the baseline model. The subscript i represents the position of the coefficient in the stiffness matrix, as shown in Figure 1; α represents the scaling factor. Three groups of comparative analyses are conducted for each stiffness coefficient, and the scaling factor settings of each model are shown in Table 1. Since the influence of flexible blades on wind turbines is mainly investigated, the scaling factor is less than 1.
For the non-main diagonal of the stiffness matrix, the stiffness factor is set as follows: where K 55 represents the flap-wise stiffness of each section in the baseline model; β represents the scaling factor. Three groups of comparative analyses are conducted for each stiffness coefficient, and the scaling factor settings of each model are shown in Table 2. Generally speaking, the coupling stiffness is less than the flap-wise stiffness, so the scaling factor is less than 1. Except for the model corresponding to K 12 ′ , the scaling factors of other models are the same.

SIMULATION ANALYSIS RESULTS OF STEADY WIND CONDITION
To validate the accuracy of the aeroelastic analysis method combined with Aerodyn and GEBT, the aeroelastic response calculation results of the baseline model under steady-state wind conditions are compared with those of the commercial software SIMPACK. The initial rotation speed of the wind turbine is 5 RPM. The wind speed is kept at 12 m/s. The simulation time is 200s. The comparison results of simulation analysis are shown in Figure 2. Figures 2A,B show the comparison of generator power and rotation speed under 200 s steady-state wind conditions. The simulation results of SIMPACK and this method are basically consistent. Figures 2C-F show the comparison between the calculated results of tip deformation and root load. During the start-up period of 0-50 s, the tip deformation fluctuates irregularly, and the wind turbine entered a stable state after 50 s. The amplitude of the tip displacement in the out-ofplane direction of this method is obviously larger than that of SIMPACK, and the calculated results of the mean value in this direction are close to each other. The results of the two methods  are close to each other in the in-plane displacement amplitude, and the mean value of this method is slightly larger than that of SIMPACK. The amplitude and maximum of the blade root bending moment in the out-of-plane direction of the wind turbine calculated by this method are greater than those calculated by SIMPACK, and the blade root bending moment in the in-plane direction calculated by the two methods is almost the same. It can be seen from the results of a steady-state aeroelastic simulation that this method and SIMPACK aeroelastic simulation system can achieve the expected aeroelastic simulation effect by adopting the same unit setting and the same control strategy. The comparison results show that the calculation results of this method are reliable.
The same steady-state wind condition and initial condition are used for sensitivity analysis. In the design of wind turbine components, the load is one of the main input parameters, so the sensitivity of blade stiffness to blade root moment is taken as the main research object.
Taking the effect of flap-wise stiffness on the aeroelastic response of wind turbine as an example, Figure 3 and Figure 4 show the influence of flap-wise stiffness on power and pitch angle, and the influence of flap-wise stiffness on blade tip deformation, respectively. The '5_5_0.3' in the legend means the EIflap in the stiffness matrix of Figure 1 is 0.3 times of the baseline model. It can be seen from Figure 3 that when the flap-wise stiffness is smaller, the output power of the wind turbine is more difficult to maintain stability. Especially, when the flapwise stiffness decreases to 30% of the baseline model, the blade tip deformation is too large to capture enough wind energy, even when the pitch angle is adjusted to 0, it still can not provide enough lift and torque, and the output power can not reach the rated power. The results are consistent with the previous studies (Larsen et al., 2004). SAF (Sensitivity Analysis Factor) is used to quantitatively analyze the influence of each stiffness coefficient on the physical quantity of root moment. SAF is defined as follows: where ΔA/A is the change rate of the dependent variable, such as ΔMx/Mx, ΔUx/Ux etc.; ΔF/F is the change rate of the independent variable, such as ΔK 11 /K 11 , ΔK 23 /K 23 etc. SAF > 0 indicates that the independent variable and the dependent variable change in the same direction, and SAF > 0 means the opposite direction of change. The larger the absolute value of SAF, the more sensitive the dependent variable is to the independent variable, On the contrary, it is not sensitive. Here, we define |SAF|≤0.1 as insensitive; 0.1<|SAF|≤1 represents moderate sensitivity; |SAF|>1 represents high sensitivity. Figure 5 shows the sensitivity of the main diagonal stiffness to the root moment. The number of the legend indicates the position of the stiffness parameter in the stiffness matrix, for instance, '4_3' of the legend is the stiffness coefficient of the fourth row and the third column in the stiffness matrix. It can be seen from Figure 5A that K 33 , K 55 , and K 66 are moderately sensitive to the mean value of root moment Mx, and K 33 is highly sensitive to the range of Mx, as can be seen from Figure 5D.
It can be seen from Figure 5B that all main diagonal coefficients have low sensitivity to the mean value of blade root moment My. For the range of My, K 33 , K 44 , K 55 , and K 66 have moderate sensitivity in Figure 5E.
As can be seen from Figure 5C,F that the average of Mz is very sensitive to K 55 and K 66 . For range of Mz, K 55 has high sensitivity, K 33 and K 44 have medium sensitivity. Figure 6 shows the sensitivity coefficient of non-main diagonal stiffness to blade root moment. It can be seen from Figure 6A,D that K 16 and K 56 have high sensitivity to the mean value of blade root moment Mx. And K 26 , K 36 and K 46 have moderate sensitivity to the mean value of blade root moment Mx. all the non-main diagonal coefficients have low sensitivity to the range value of blade root moment Mx.
As can be seen from Figure 6B,E, K 16 and K 56 have high sensitivity to the My, including mean and range values.
For the mean value of root Mz, except K 13 , K 15, and K 35 , other coefficients are highly sensitive. For the range of Mz, K 15 , K 24 , K 26 , K 45 , and K 46 are moderately sensitive, while other coefficients are insensitive.

SIMULATION ANALYSIS RESULTS OF TURBULENT WIND CONDITION
In order to analyze the simulation results under the condition of turbulent wind, Mann (Mann, 1994) uniform shear model is used. The representative value of the turbulence standard deviation σ 1 is given by the 90% quantile for the given hub height wind speed. The  turbulence wind in this section is set as follows: hub center height is 90 m, average wind speed is 12 m/s, longitudinal turbulence degree is 14%, horizontal turbulence degree is 9.8%, and vertical turbulence degree is 7%. The calculated wind speed in each direction of the turbulent wind is shown in Figure 7.
The aeroelastic responses of the wind turbine are analyzed by using the simulation analysis method and stiffness adjustment strategy in Model description. Taking the flap-wise stiffness as an example, the comparison results of pitch angle, blade root moment My, and blade tip displacement Ux are shown in Figure 8A shows the adjustment of blade pitch angle to ensure stable power output. It can be seen that the rated power cannot be reached even if the pitch angle is adjusted to 0°in the simulation time of 90 s and 170-180 s. But it does not affect the sensitivity research. It also can be seen from Figure 8 that the flap-wise stiffness has an obvious effect on the bending moment and deformation of the blade. In order to quantify the degree of impact, the SAF factor is also introduced to analyze the sensitivity under turbulent wind conditions.
In the case of turbulent wind conditions, the extreme values of blade root moment and tip deformation are taken as the objects for sensitivity analysis. The sensitivity factor of the stiffness coefficient of the main diagonal and non-main diagonal are shown in Figure 9 and Figure 10. The results show that K 33 , K 44 , K 55 , and K 66 have moderate and high sensitivity to the moment and deformation of blade aeroelastic response.
It can be seen from Figure 10 that K 16 and K 56 is highly sensitive to MX. K 16 , K 46 , and K 56 have moderate and high sensitivity to the My. For the Mz of blade root, except K 13 , K 15 , and K 35 , other coefficients are highly sensitive. Coupling stiffness K 56 have high sensitivity to the tip deflection Ux in Figure 10D. Coupling stiffness K 14 , K 16 , K 25 , K 45 , K 46 , and K 56 have high sensitivity to the tip deflection Uy in Figure 10E. K 56 have high sensitivity to the tip deflection Uz in Figure 10F.

CONCLUSION
The influence of blade beam section stiffness coefficients on wind turbine aeroelastic responses is studied by the aeroelastic analysis method based on the aeroelastic analysis module and nonlinear beam theory. The aeroelastic analysis module is Aerodyn based on the BEM theory, and the nonlinear beam theory is GEBT. The section stiffness of the blade is a complete stiffness matrix based on the generalized Timoshenko beam theory. As the flap-wise stiffness in the matrix is one of the most important stiffness in the blade design, the flap-wise stiffness of the baseline model is used as the reference to adjust the stiffness coefficients. By adjusting the single stiffness coefficient while other conditions remain unchanged, the aeroelastic response analysis is carried out. Finally, the sensitivity analysis coefficient is used for quantitative comparison and the aeroelastic simulation results are compared. The conclusions are as follows: a. The sensitivity analysis results of turbulent and steady wind conditions are basically consistent, but the sensitivity of turbulent wind conditions is generally lower than that of steady wind conditions. b. The deformation of the blade has an obvious influence on the power output of the wind turbine. The main reason is that the swept area of the wind turbine is reduced due to the FIGURE 10 | Influence of non-main diagonal stiffness coefficient on medium and high sensitivity.
Frontiers in Energy Research | www.frontiersin.org August 2021 | Volume 9 | Article 707082 deformation of the blade, which is not enough to capture enough wind energy. Even if the pitch angle is adjusted, the power output is still unstable. c. The main diagonal coefficients, which are more sensitive to the root moment Mx (medium and high sensitivity), are K 33 , K 55 , and K 66 . For the mean value of Mx, the change of stiffness is opposite to the change of mean value. For the range value of Mx, the change of K 33 is in the same direction as that of the range value, that is, the decrease of K 33 , will reduce the range value of Mx. d. All main diagonal coefficients are positively sensitive to the mean value of My and are not sensitive. For the range value of blade root moment My, excluding the case of power generation not reaching the rated power, the changes of K 33 , K 44 , and K 66 are in reverse with the range value, and the changes of K 55 are in the same direction with the range value of My. e. For blade root moment Mz, K 55 has high sensitivity and reverse change, while K 33 and K 44 have medium sensitivity and change in the same direction. f. The non-diagonal coefficients K 16 and K 56 are highly sensitive to blade root moment. K 16 and K 56 had reverse medium sensitivity to the mean value of Mx. K 16 and K 56 had reverse high sensitivity to My (including mean and range values). Here, K 56 is the key parameter of the so-called bending torsion coupling effect, and many scholars have published a large number of research results on it. Except K 13 , K 15 , and K 35 , the other non-diagonal coefficients were highly sensitive to the mean value of root Mz.
The load of wind turbines has a direct impact on the structural design of wind turbines. The work of this paper gives the sensitivity analysis method and results of the blade stiffness to the aeroelastic response of wind turbine under certain conditions, which is helpful to improve the understanding of the coupling effect between structure and load and to provide some guidance for the direction of structural adjustment in blade optimization design.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CC performed the calculations and tests and wrote the manuscript. LW checked the calculation results, and TW helped analyze the results. All authors have read and agreed to the published version of the manuscript.