Load Distribution Optimization Method for Fatigue Test of Full-scale Structure of Biaxial Resonant Wind Turbine Blade

: Wind turbine blade is the key component of wind turbine to realize wind energy capture, so it is necessary to verify the rationality of blade structure design through full-scale structural fatigue test. In order to improve the test accuracy and efficiency, this paper proposes to establish the dynamic analysis model of the multi-degree of freedom system of the tested blade by using the finite beam element to realize the calculation method of the load amplitude of the biaxial fatigue test and integrate the dynamic analysis model and the target load calculation method into the particle swarm optimization algorithm to realize the optimization of the load distribution of the biaxial resonant full-size structure fatigue test. Through the design of a biaxial resonant loading scheme of 2.5MW-52.5m wind turbine blades, the results show that this method can quickly and accurately adjust the optimization parameters such as the position of the exciter and the motion quality of the exciter in the flapwise and edgewise on the premise that the fatigue cumulative damage caused by the test load is not less than the cumulative damage of the target load.


Introduction
Wind turbine blade is a key component of wind turbine to convert wind energy into mechanical energy. Its performance directly affects the operational reliability and power generation efficiency of wind turbines [1]. Due to the harsh working environment, blades subjected to random wind load and inertia load coupling effect are very easy to fatigue failure, so the blade fracture caused by fatigue has become the main form of wind turbine blade failure [2]. In order to ensure the structural reliability of blades, the newly developed or process improved wind power blades must be tested for full-scale structural fatigue according to relevant standards before being put into the market, so as to check whether the blades meet the fatigue life specified in the design. At present, wind power equipment manufacturing enterprises and third-party certification institutions mainly use the uniaxial loading mode (as shown in Figure 1) to carry out full-scale structural fatigue tests on blades. The test system excites the blade by the inertia force generated by the rotation of the eccentric mass driven by the motor, so that the blade enters the resonance state and vibrates back and forth in the flap and edgewise, so as to produce the test load. However, with the continuous development of blades to large-scale, uniaxial loading fatigue testing technology has been unable to meet the testing requirements of blade quality and reliability [3]. To solve this problem, it is urgent to carry out the research of blade biaxial loading fatigue test. Before the fatigue test of biaxial resonance full-scale structure, the test loading scheme needs to be designed to ensure that the fatigue cumulative damage generated by the test load in the key area is not lower than that of the target load [4]. Considering the particularity of biaxial resonant loading mode, the design of test loading scheme is essential to achieve the consistency between the test load in the key area of the blade and the target load distribution by adjusting the position of the exciter on the flap and edgewise, the moving quality of the exciter, the position and quality of the fixed counterweight of the blade, which belongs to the category of engineering optimization problem [5]. However, the design of the biaxial resonance loading test scheme changes the mean value of the target load while adjusting the above parameters, which requires recalculation of the target load amplitude according to the current loading scheme and judgment of the consistency between the test load distribution and the target load distribution. The whole optimization process will become relatively complicated, and there is still a lack of effective test load distribution optimization method for biaxial resonance full-scale structural fatigue test.
In view of this, taking the calculation of test load amplitude under biaxial resonant loading as the starting point, the tested blade is simplified into a cantilever beam model, and the multi-degree of freedom system dynamics analysis of the beam model is carried out by using the finite element method to establish the calculation method of test load amplitude, On this basis, the test load amplitude calculation method and the target load calculation method proposed in the previous research work are integrated into the particle swarm optimization algorithm to optimize the fatigue test load distribution of biaxial resonance full-size structure. Finally, the effectiveness of the optimization algorithm is verified by designing the biaxial resonance loading scheme of 2.5MW wind turbine blades.

Calculation method of test load amplitude based on dynamic analysis of multi degree of freedom system
In fact, under the condition of biaxial resonance loading, the mean value of the test load is determined by the blade self-weight, the mass and position of the exciter, and the mass and position of the counterweight. The test load amplitude is determined by the inertial load generated by the reciprocating motion of the above masses. The calculation of the test load amplitude can be regarded as the steady-state response under the blade resonance state. Through the steady-state response analysis results, the variation law of blade displacement with time is established, and the derivative is obtained to obtain the acceleration time relationship, so as to obtain the inertial force of the blade at any time and realize the calculation of the test load amplitude.
Considering that the blade under biaxial resonant loading is a continuously distributed mass system with infinite degrees of freedom, the motion law of the system needs to be described by a function containing both time and space coordinates. It is difficult to directly analyze the motion equation from the second-order ordinary differential equations of the finite degrees of freedom system to partial differential equations. Therefore, in this section, under the assumption that the vibration in the flap and edgewise are independent of each other, the blade model is simplified by using the finite element, and the test load amplitude is calculated by solving the steady-state response of the multi-degree of freedom damped vibration system. The process is divided into six steps: Step 1. Simplify the blade model under biaxial resonant loading The blade is divided into n segments, and each segment represents a finite element. For any element ( ) i , its corresponding beam segment is a beam with equal section. The length, bending stiffness and mass per unit length of the beam segment are ( )  , respectively. The adjacent elements are connected by nodes, with a total of n+1 nodes. Assuming that the cross-section of the cantilever beam before and after bending is always plane and perpendicular to the main axis, and the vibration in the swing and shimmy directions are independent of each other, each node of the finite element has two degrees of freedom, namely transverse displacement i u and plane rotation angle 1, 2, , 1) . At this time, the local coordinate system of the element coincides with the global coordinate system, The above process realizes the transformation from infinite degree of freedom continuous distributed quality system to multi degree of freedom discrete quality system, as shown in Figure 2.
Any element represents the force generated by the unit displacement of the nth degree of freedom on the mth degree of freedom, which can be calculated by Equation (7): The stiffness matrix ( ) i k of the finite element ( ) i is obtained as: Using the same method, any element i represents the force generated by the unit acceleration of the kth degree of freedom on the jth degree of freedom, which can be calculated by Equation (9): The mass matrix ( ) i m of the finite element ( ) i is obtained as: Step 3. Overall stiffness matrix and mass matrix of container After the calculation of each element stiffness matrix and mass matrix is realized in the second step, it is necessary to assemble each element stiffness matrix and mass matrix in a certain order to obtain the overall stiffness matrix and mass matrix. Since the blade cantilever model is divided into n segments, there are n+1 nodes in the multi degree of freedom system, and each node has two degrees of freedom, then the overall stiffness * k and mass * m are (2 2) (2 2) n n + × + matrix. Considering that the elements of each element stiffness matrix and mass matrix correspond to the degrees of freedom of element nodes, after all nodes are uniformly numbered, put the element stiffness matrix and mass matrix into the corresponding positions of the overall stiffness matrix and mass matrix, and add them to obtain the overall stiffness matrix * k and mass matrix * m . Taking finite elements ( ) i and ( 1) i + as examples, the position of the stiffness matrix ( ) k is directly placed in the corresponding position in . The position of the stiffness matrix ( 1) i + k of the element ( 1) i + corresponding to the overall stiffness matrix * k is (2 1) (2 4) i i + × + . Put ( 1) i + k directly into the corresponding position in * k and sum the part coincident with ( ) i k the whole process is shown in Figure 3.
Step 4. The undamped free vibration differential equation of multi degree of freedom system is established, and the natural frequency, main vibration mode, regularization factor and regular vibration mode matrix of the system are calculated Before establishing the differential equation of motion, the support conditions need to be imposed on the system. Considering the degrees of freedom 1 0 n u + = and 1 0 n θ + = at the fixed end of the cantilever beam model in Figure 1, only the degrees of freedom of the nodes in front of the fixed end need to be considered in the dynamic analysis. Therefore, remove the 2n+1 row, 2n+2 row, 2n+1 column and 2n+2 column in the overall stiffness matrix * k and mass matrix * m of the system to obtain the stiffness matrix k and mass matrix m of the multi degree of freedom system, both of which are 2 2 n n × matrices. Then the differential equation of undamped free vibration of multi degree of freedom system is: Since Equation (11) has a non-zero solution, the characteristic polynomial is equal to zero The natural frequency ( 1, 2, , 2 ) i i n ω =  of the multi degree of freedom system is calculated by Equation (12), and the corresponding main vibration mode i u is calculated by substituting i ω into Equation (13) ( ) , the corresponding principal mass matrix M and principal stiffness matrix K are calculated by Equations (14) and (15) respectively: For the regular mode shape matrix Figure 3: Schematic diagram of overall stiffness matrix container loading process.
Step 5. The differential equation of damped forced vibration of multi degree of freedom system is established and solved During the fatigue test of blade full-scale structure, the excitation frequency is close to the natural frequency of the system, and the influence of structural damping on the system vibration will be very significant. However, considering that the damping mechanism of blade structure is very complex, in order to reduce the difficulty of analysis, Rayleigh damping hypothesis [6] is used to describe the blade structural damping: Where, α and β are proportional coefficients, which can be calculated by Equations (20) Where, ζ is the damping ratio, subscripts f and e correspond to the waving direction and swing direction respectively, f ω , e ω , f ζ and e ζ can be measured through the free attenuation vibration experiment of the system [7]. Structural damping is a linear combination of system mass matrix and stiffness matrix, therefore ( ) The differential equation of damped forced vibration of multi degree of freedom system is: is the external force on the multi degree of freedom system: Where, 0 F is the amplitude of the excitation force, ω is the excitation frequency, and the position of 0 F in the vector 0 F corresponds to the node position of the exciter. The damped forced vibration differential equation of multi degree of freedom system is decoupled by using the regular vibration mode matrix ′ u , so that ( ) ( ) t t ′ = q uη is substituted into Equation (23) and multiplied by ( ) Where,  Step 6. Calculate the test load amplitude According to the steady-state response analysis results, the lateral displacement i u of the ith node of the multi degree of freedom system (excluding the fixed end node, 1, 2, , The corresponding node acceleration i a obtained by deriving the lateral displacement i u is: If the mass of each section of the cantilever beam is concentrated to the node, the mass i m′ of each node is: According to Equations (31) and (32), the inertia force of each node is i i m a ′ ⋅ , then the test load amplitude of the section corresponding to any node j ( 2,3, , ) j n = ⋅⋅⋅ is equal to all node inertia forces before node j, and the moment is taken for point J, as shown in Equation (33): Where, ji z is the longitudinal distance from node i to node j.

Test load distribution optimization method based on particle swarm optimization
In this section, an optimization method of biaxial resonance full-scale structure fatigue test load distribution based on particle swarm optimization algorithm is proposed. The particle swarm optimization algorithm is used to control the parameters such as the position and moving mass of the exciter, the quantity, position and mass of the counterweight, realize the optimization of the test load distribution and determine the test loading scheme on the premise that the average value of the target load changes with the optimization parameters.
Since the stiffness of the blade in the flapwise is less than that in the edgewise, and the bending deformation amplitude in the flapwise is usually greater than that in the edgewise, when the blade is installed and fixed, make the edgewise of the blade parallel to the horizontal plane and the flapwise parallel to the vertical plane. At this time, the stress ratio corresponding to the bending deformation in the edgewise of the blade is still -1, while the stress ratio in the flapwise will no longer be equal to -1 due to the influence of factors such as blade weight. On this premise, the optimization process of test load distribution is divided into seven steps: Step 1. Import the data of blade stiffness distribution ( ) i EI , mass distribution ( ) i m and length ( ) i l , and use the calculation method of test load amplitude proposed in this chapter to establish the blade finite element model (the blade finite element model is divided into n segments, with n+1 nodes, node 1 is located at the blade tip and node n+1 is located at the fixed end of blade root, 1, 2, , i n =  ), and determine the position of the section where each node is located.
The fatigue load spectrum of the blade is introduced, and the fatigue cumulative damage of the section where each node of the blade finite element model is calculated by using the fatigue test target load calculation method proposed in Reference [8], and the target load with the stress ratio 1 r = − in the flap and edgewise of the section is generated respectively. Since the stress ratio Step 5. The fitness function obj is defined as shown in Equation (39), which describes the error between the test load amplitude and the target load amplitude generated by the corresponding loading scheme of each particle:  The fitness value of each particle in the current iteration process is compared with the minimum fitness value experienced by the particle. If the former is small, the particle is stored in the pbest (i.e. the individual extreme value of the particle) for updating, otherwise the pbest remains unchanged. Then, the individual extremum stored in the pbest is compared with the fitness value of the particle gbest (i.e. the global extremum of the population) corresponding to the minimum fitness of the current iteration, and the gbest is updated with the individual extremum with the minimum fitness.
Step 6. Update the velocity and position of the qth particle, as shown in Equations (40) and (41): Where, 1, 2, , t T =  ( T is the total number of iterations), u is inertia weight, 1 c , 2 c distribution is individual learning factor and group learning factor, 1 r , 2 r is a random number between 0 and 1, ( ) q pbest is the best position traveled by the qth particle.
Step 7. Return to step 4 and continue the search until the error between the test load amplitude and the target load amplitude reaches the expected requirements or the number of iterations reaches the upper limit. It should be noted that the above test load distribution optimization process does not consider the arrangement of counterweight. In fact, the arrangement of vibration exciter and counterweight on the blade section will inevitably affect the stress-strain distribution near the installation position. The number of counterweights needs to be strictly limited when designing the loading scheme. Therefore, in the actual loading scheme design, the test load distribution without counterweight can be optimized according to the above process, and whether the counterweight needs to be added can be judged according to the results. If the distribution optimization error is large, the twodimensional optimization parameters of counterweight position p z and mass p m can be added to particle position , q t P accordingly. In addition, considering that there is a certain difference between the target load amplitude distribution and the test load amplitude distribution of the blade, if the test load amplitude of all sections is consistent with the target load amplitude, it will bring great difficulty to the optimization problem. Therefore, when calculating the particle fitness, the error between the test load amplitude and the target load amplitude in the key area of the blade can be defined as the fitness function, to reduce the difficulty of optimization problems.

Example verification
This section takes 2.5MW-52.5m blade as an example to verify the effectiveness of the proposed method. According to the requirements of the industry standard DNVGL-ST-0376 [9], the load distribution of biaxial resonance fatigue test in the area of 0~40% length from blade root to blade tip is optimized. The blade is made of glass fiber / epoxy resin composite material. The M-N curve parameter m=10 is determined according to the data provided by the designer.     The blade is divided into 23 sections according to the mass distribution and stiffness distribution of the blade. The blade mode is calculated by using the multi degree of freedom system dynamic analysis model of finite beam element proposed in this section. The results are shown in Table 1. Among them, the measured value of blade natural frequency is obtained by free attenuation test and provided by blade designer. According to the analysis results, the modal prediction results of blade on the flap and edgewise are close to the measured results, and the maximum error is not more than 5%, which verifies the effectiveness of the model. The method proposed in this chapter is used to optimize the load distribution of biaxial resonance loading fatigue test for 2.5MW-52.5m blades. The optimization parameters are defined as the position f z and motion mass f d m of the vibration exciter in the flapwise, the position e z and motion mass e d m of the vibration exciter in the edging direction respectively. The optimization objective is that the target load amplitude distribution is consistent with the test load distribution within the range of 0~40% of the blade spanwise, and the relative error of bending moment is no more than 20%. Table 2 gives the value range and optimization results of optimization parameters. Due to the introduction of external mass, the natural frequencies of flap and edgewise are reduced to 0.77Hz and 1.22Hz respectively, during biaxial resonance loading fatigue test, the vibration exciters in two directions will excite the blade near these two frequencies. Figure 8 shows the distribution of the test bending moment amplitude and the target bending moment amplitude of the blade. From the distribution results of the edgewise, it can be found that the loading scheme does not change the target load amplitude in this direction, and the test load amplitude within 0~40% of the spanwise of the optimized blade is greater than the target load amplitude. At the same time, due to the arrangement of the flapwise exciter at the blade spanwise 22.17m, the introduced static mass of the exciter makes the load amplitude curve of edging test appear a break point at this position. From the results of flapwise distribution, it can be found that the target load amplitude changes due to the introduction of the external mass of the loading scheme. The introduction of external mass increases the average value of the target load of each section within the range of 0~40% of the blade spanwise and makes the section stress ratio 1 r − ＜ . However, considering that the total fatigue cumulative damage of the section remains unchanged, the amplitude of the target load decreases and the amplitude of the decrease from the blade root to the blade tip gradually weakens.  Figure 9 shows the relative error between the spanwise test bending moment amplitude and the target bending moment amplitude of the optimized blade, which is calculated by Equation (42).
From the relative error between the test bending moment amplitude and the target bending moment amplitude, it can be found that the relative error in the edging direction is controlled within 20%, and the relative error is small at both ends and large in the middle within the optimization range of 0~40% of the blade spanwise, while the relative error in the flapwise fluctuates greatly, and the test load amplitude of some sections is lower than the target load amplitude. The reason for this phenomenon is that the target load amplitude distribution is determined by the fatigue cumulative damage of the blade section under 20-year working condition. For the edging direction, the main load leading to the fatigue cumulative damage is the bending moment caused by the blade selfweight, so the target load amplitude distribution is close to the smooth curve. The amplitude distribution of the test load in the edging direction is determined by the bending moment caused by the self-weight of the blade and the inertia force of the reciprocating motion of the external mass. Although it is different from the target load amplitude distribution, it is generally close to the smooth curve. Therefore, the relative error in the edging direction is small at both ends and large in the middle. The fatigue cumulative damage under the condition of blade is mainly caused by random wind load. There are certain differences in the damage size of each section, and the target load amplitude distribution is no longer continuous. However, the test load amplitude distribution in flapwise is still determined by the inertia force of blade self-weight and reciprocating motion of external mass. Therefore, the relative error in flapwise presents a certain fluctuation.

Conclusion
In this paper, the optimization method of fatigue test load distribution of biaxial resonant fullscale structure is studied. Taking the multi-degree of freedom system dynamic analysis model based on finite beam element as the starting point, the calculation method of test load amplitude is established. On this basis, the calculation method of test load amplitude and the target load calculation method are integrated into particle swarm optimization algorithm, the load distribution optimization of biaxial resonance full-scale structure fatigue test is realized. Finally, the effectiveness of this method is verified by designing the biaxial resonance loading scheme of 2.5MW-52.5m wind turbine blades. The results show that the method proposed in this chapter can quickly and accurately adjust the optimization parameters such as the position of the exciter in the flap and edgewise and the motion quality of the exciter on the premise that the fatigue cumulative damage generated by the test load in the specified area of the blade is not less than the cumulative damage of the target load, it can provide a theoretical reference for the distribution optimization of fatigue test load of biaxial resonance full-scale structure.