Wind turbine blade residual life prediction method based on unscented Kalman filter

Abstract

The invention discloses a wind turbine blade residual life prediction method based on unscented Kalman filtering, and the method comprises the following steps: S1, building a blade damage model, operating a wind turbine, and enabling a blade to crack and gradually expand; s2, converting blade crack propagation into a discrete accumulation process to obtain a state transition equation and an observation equation; s3, determining annual average wind speed, cut-in wind speed, cut-out wind speed and Rayleigh distribution. S4, obtaining observation data of the crack length of the blade; s5, defining a damage state process, and setting a covariance matrix of system noise and crack length measurement errors; s6, performing Bayesian prediction by using an unscented Kalman filtering method, and performing Bayesian updating; and S7, setting a blade crack propagation damage threshold value, judging whether the updated state exceeds the set damage threshold value or not, and calculating the remaining service life. According to the prediction method, the problem that errors and calculation complexity are difficult to coordinate and reduce in the prior art is solved.

Description

Technical Field

[0001] The invention belongs to the technical field of wind power generation, and in particular relates to a method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering. Background Art

[0002] Wind turbines are composed of multiple components, and their blades are affected by many factors such as wind shear, tower shadow effect, turbulent wind, gravity load, centrifugal load, and surface erosion. Therefore, blade failure accounts for about 23% of the total failures of wind turbines, and replacing blades will lead to long downtime and increase maintenance costs. Currently, wind turbines are equipped with a large number of sensors to collect operating data of wind turbines. These data are not only used to monitor wind turbines, but also to evaluate the performance and health of various components of the unit based on prediction and health management systems. Among them, the accurate prediction of the remaining service life of the unit blades is particularly critical, which plays an important role in ensuring the normal operation of the unit and preventing accidents.

[0003] The key to predicting the remaining service life of wind turbine blades is to predict the evolution of fatigue damage and solve the nonlinearity and uncertainty of crack propagation. The commonly used methods are the extended Kalman filter and the particle filter. Although the extended Kalman filter has shown certain effects in estimating nonlinear systems, it also has certain limitations. Specifically, when the nonlinearity of the system is low, the error introduced by the extended Kalman filter is not obvious; however, for strong nonlinear systems, the extended Kalman filter directly ignores high-order terms for linearization, which will lead to large errors. In contrast, the particle filter is also applicable to nonlinear systems and can effectively solve the problem of large errors. However, the particle filter method has a relatively large amount of calculation and the implementation process is relatively complicated. There is no feasible solution to the contradiction between reducing errors and computational complexity. Summary of the invention

[0004] The purpose of the present invention is to provide a method for predicting the remaining life of a wind turbine blade based on an unscented Kalman filter, which solves the problem in the prior art that it is difficult to coordinate the reduction of error and the computational complexity.

[0005] The technical solution adopted by the present invention is a method for predicting the remaining life of a wind turbine blade based on an unscented Kalman filter, comprising the following steps: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; S3, determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6 and continue iterative operation.

[0006] The present invention is also characterized in that: The equation for blade crack growth in S2 is: (1); In the formula, a is the current blade crack length; R is the stress ratio, i.e., the ratio of the root mean square value of the minimum stress to the maximum stress in the stress cycle; A , m All are blade material related parameters; is the stress intensity factor, The calculation of requires blade cycle stress, which is obtained by the rainflow counting method.

[0007] The system noise in S2 is , the crack length measurement error is , and All obey a normal distribution with a mean of 0, that is, , ,in, is the covariance matrix.

[0008] The state transfer equation in S2 is: ; The observation equation is: .

[0009] The specific process of determining the annual average wind speed under the standard wind turbine category in S3 is as follows: According to the latest version of the wind power design standard IEC61400-1Y2019 issued by the International Electrotechnical Commission (IEC), the annual average wind speed under the standard wind turbine category is divided into three levels: 7.5m / s, 8.5m / s and 10m / s.

[0010] The process of determining the cut-in wind speed and cut-out wind speed of the wind turbine in S3 is as follows: According to the regulations of the National Renewable Energy Laboratory (NREL) of the United States, the cut-in wind speed of a 5MW wind turbine is 3m / s and the cut-out wind speed is 25m / s.

[0011] The probability density function of the Rayleigh distribution in S3 is: (2); In the formula, is the average wind speed, unit: m / s.

[0012] The specific process of S4 is: S4.1. Update the turbulent wind speed data that conforms to the Rayleigh distribution every 10 minutes, and calculate the time series load of the blades according to the maximum wind energy tracking and constant power control strategy. After a complete cycle, use the rain flow counting method to perform statistical analysis on these time series load data to obtain the number of cycles and the corresponding time series load amplitude range; S4.2, input the data in S4.1 into the blade crack growth model to calculate the blade crack growth rate; S4.3, multiply the rate in S4.2 by the time of the cycle to obtain the crack extension length of the blade in this cycle; S4.4, adding the extension length in this cycle to the blade crack length in the previous cycle, to obtain the crack length at the current moment; S4.5. Repeat S4.1-S4.4 to finally obtain the observation data of the blade crack length in the entire time period.

[0013] The Bayesian prediction in S6 assumes that the posterior probability density function of a given state is , then the state prior distribution at the current k moment is for: (3); In the formula, are independent observed variables; is an unknown variable of the system and obeys a first-order Markov process; is the state transition probability density function of the system; is the posterior probability density obtained at the previous moment; Bayesian updating is to obtain the observation value at the current k moment through the system observation model , update and adjust the system's prior probability density function , realize the transformation from the prior probability density function to the posterior probability density function at time k Derivation of: (4); Where: is the likelihood function, which describes the possibility of the unknown variable appearing under given observations; is the normalization constant.

[0014] The specific process of Bayesian prediction using the unscented Kalman filter is as follows: Step 1: Obtain the Sigma point set and its corresponding weights according to formula (5); (5); Where: is a scaling function used to reduce the total prediction error; is an estimated value; Step 2: Calculate the one-step prediction of 2n+1 Sigma point sets according to formula (6); (6); Step 3: Calculate the one-step prediction and covariance matrix of the system state quantity according to equations (7) and (8); (7); (8); Where: is the weight; Step 4: Based on the one-step prediction value, use the untraceable transformation again according to formula (9) to generate a new Sigma point set; (9); Step 5: Substitute the predicted Sigma point set in S4 into the observation equation according to formula (10) to obtain the predicted value of the Sigma point set; (10); Step 6: According to equations (11), (12) and (13), the predicted value of the Sigma point set is obtained from S5, and the mean and covariance of the system prediction are obtained by weighted summation; (11); (12); (13); Step 7: Calculate the Kalman gain matrix according to formula (14); (14); Step 8: Calculate the state update and covariance update of the system according to equations (15) and (16); (15); (16).

[0015] The beneficial effects of the present invention are: The wind turbine blade remaining life prediction method based on unscented Kalman filtering provided by the present invention collects blade flapping bending moment data in the entire wind speed range, observes crack propagation under maximum wind energy tracking and constant power control strategies according to the three wind speed levels specified in the IEC standard, and combines various conditions of actual operation of the wind turbine to make the observation of crack length more realistic; unscented Kalman filtering is used to predict the evolution of fatigue damage and solve the nonlinearity and uncertainty problems of crack propagation. Compared with other methods, the method has a small amount of calculation, a small error, and a more accurate prediction of the remaining service life. BRIEF DESCRIPTION OF THE DRAWINGS Figure 1 It is a schematic diagram of the rain flow counting method of the wind turbine blade remaining life prediction method based on unscented Kalman filtering of the present invention; Figure 2 is the wind speed Rayleigh distribution diagram under different annual average wind speeds of the present invention; Figure 3 is a blade crack observation flow chart of the present invention; Figure 4 is a blade crack prediction flow chart of the present invention; Figure 5 is a stress range distribution diagram of the rain flow counting method in Example 6 of the present invention; Figure 6 is a comparison diagram of crack observations in a constant wind speed interval and a variable wind speed interval in Example 6 of the present invention; Figure 7 This is a comparison diagram of crack observations of three dangerous nodes of blades in Example 6 of the present invention; Figure 8 is a blade crack prediction diagram based on unscented Kalman filtering in embodiment 6 of the present invention; Fig. 9 is a blade life distribution diagram based on unscented Kalman filtering in Example 6 of the present invention; Fig.10 It is a blade crack prediction diagram of different initial lengths in Example 6 of the present invention. DETAILED DESCRIPTION

[0016] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0017] Example 1 The remaining life prediction method of a wind turbine blade based on unscented Kalman filtering proposed in this embodiment includes the following steps: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; S3, determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6 and continue iterative operation.

[0018] Example 2 The remaining life prediction method of a wind turbine blade based on unscented Kalman filtering proposed in this embodiment includes the following steps: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; The equation for blade crack growth is: (1); In the formula, a is the current blade crack length; R is the stress ratio, i.e., the ratio of the root mean square value of the minimum stress to the maximum stress in the stress cycle; A , m All are blade material related parameters; is the stress intensity factor, The calculation of requires blade cycle stress, which is obtained by rain flow counting method, such as Figure 1 As shown; The system noise is , the crack length measurement error is , and All obey a normal distribution with a mean of 0, that is, , ,in, is the covariance matrix; The state transfer equation is: ; The observation equation is: ; S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed at 10-minute intervals; The specific process of determining the annual average wind speed under the standard wind turbine category is as follows: According to the latest version of the wind power design standard IEC61400-1Y2019 issued by the International Electrotechnical Commission (IEC), the annual average wind speed under the standard wind turbine category is divided into three levels: 7.5m / s, 8.5m / s and 10m / s; The process of determining the cut-in wind speed and cut-out wind speed of the wind turbine is as follows: According to the National Renewable Energy Laboratory (NREL) of the United States, the cut-in wind speed of a 5MW wind turbine is 3m / s and the cut-out wind speed is 25m / s; The probability density function of the Rayleigh distribution is: (2); In the formula, is the average wind speed, unit: m / s; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6 and continue iterative operation.

[0019] Example 3 The remaining life prediction method of a wind turbine blade based on unscented Kalman filtering proposed in this embodiment includes the following steps: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; The equation for blade crack growth is: (1); In the formula, a is the current blade crack length; R is the stress ratio, i.e., the ratio of the root mean square value of the minimum stress to the maximum stress in the stress cycle; A , m All are blade material related parameters; is the stress intensity factor, The calculation of requires blade cycle stress, which is obtained by rain flow counting method, such as Figure 1 As shown; The system noise is , the crack length measurement error is , and All obey a normal distribution with a mean of 0, that is, , ,in, is the covariance matrix; The state transfer equation is: ; The observation equation is: ; S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed at 10-minute intervals; The specific process of determining the annual average wind speed under the standard wind turbine category is as follows: According to the latest version of the wind power design standard IEC61400-1Y2019 issued by the International Electrotechnical Commission (IEC), the annual average wind speed under the standard wind turbine category is divided into three levels: 7.5m / s, 8.5m / s and 10m / s; The process of determining the cut-in wind speed and cut-out wind speed of the wind turbine is as follows: According to the National Renewable Energy Laboratory (NREL) of the United States, the cut-in wind speed of a 5MW wind turbine is 3m / s and the cut-out wind speed is 25m / s; The probability density function of the Rayleigh distribution is: (2); In the formula, is the average wind speed, unit: m / s; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; The specific process is: S4.1, if Figure 3 As shown in the figure, the turbulent wind speed data conforming to the Rayleigh distribution is updated every 10 minutes, and the time series load of the blade is calculated according to the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rain flow counting method is used to perform statistical analysis on these time series load data to obtain the number of cycles and the corresponding time series load amplitude range; S4.2, input the data in S4.1 into the blade crack growth model to calculate the blade crack growth rate; S4.3, multiply the rate in S4.2 by the time of the cycle to obtain the crack extension length of the blade in this cycle; S4.4, adding the extension length in this cycle to the blade crack length in the previous cycle, to obtain the crack length at the current moment; S4.5, repeat S4.1-S4.4, and finally obtain the observation data of the blade crack length in the whole time period; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6, continue iterative operation, such as Figure 4 shown.

[0020] Example 4 The remaining life prediction method of a wind turbine blade based on unscented Kalman filtering proposed in this embodiment includes the following steps: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; The equation for blade crack growth is: (1); In the formula, a is the current blade crack length; R is the stress ratio, i.e., the ratio of the root mean square value of the minimum stress to the maximum stress in the stress cycle; A , m All are blade material related parameters; is the stress intensity factor, The calculation of requires blade cycle stress, which is obtained by rain flow counting method, such as Figure 1 As shown; The system noise is , the crack length measurement error is , and All obey a normal distribution with a mean of 0, that is, , ,in, is the covariance matrix; The state transfer equation is: ; The observation equation is: ; S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed at 10-minute intervals; The specific process of determining the annual average wind speed under the standard wind turbine category is as follows: According to the latest version of the wind power design standard IEC61400-1Y2019 issued by the International Electrotechnical Commission (IEC), the annual average wind speed under the standard wind turbine category is divided into three levels: 7.5m / s, 8.5m / s and 10m / s; The process of determining the cut-in wind speed and cut-out wind speed of the wind turbine is as follows: According to the National Renewable Energy Laboratory (NREL) of the United States, the cut-in wind speed of a 5MW wind turbine is 3m / s and the cut-out wind speed is 25m / s; The probability density function of the Rayleigh distribution is: (2); In the formula, is the average wind speed, unit: m / s; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; The specific process is: S4.1, if Figure 3 As shown in the figure, the turbulent wind speed data conforming to the Rayleigh distribution is updated every 10 minutes, and the time series load of the blade is calculated according to the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rain flow counting method is used to perform statistical analysis on these time series load data to obtain the number of cycles and the corresponding time series load amplitude range; S4.2, input the data in S4.1 into the blade crack growth model to calculate the blade crack growth rate; S4.3, multiply the rate in S4.2 by the time of the cycle to obtain the crack extension length of the blade in this cycle; S4.4, adding the extension length in this cycle to the blade crack length in the previous cycle, to obtain the crack length at the current moment; S4.5, repeat S4.1-S4.4, and finally obtain the observation data of the blade crack length in the whole time period; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; Bayesian prediction assumes that the posterior probability density function of a given state is , then the state prior distribution at the current k moment is for: (3); In the formula, are independent observed variables; is an unknown variable of the system and obeys a first-order Markov process; is the state transition probability density function of the system; is the posterior probability density obtained at the previous moment; Bayesian updating is to obtain the observation value at the current k moment through the system observation model , update and adjust the system's prior probability density function , realize the transformation from the prior probability density function to the posterior probability density function at time k Derivation of: (4); Where: is the likelihood function, which describes the possibility of the unknown variable appearing under given observations; is the normalization constant; S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6 and continue iterative operation.

[0021] Example 5 The remaining life prediction method of a wind turbine blade based on unscented Kalman filtering proposed in this embodiment includes the following steps: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; The equation for blade crack growth is: (1); In the formula, a is the current blade crack length; R is the stress ratio, i.e., the ratio of the root mean square value of the minimum stress to the maximum stress in the stress cycle; A , m All are blade material related parameters; is the stress intensity factor, The calculation of requires blade cycle stress, which is obtained by rain flow counting method, such as Figure 1 As shown; The system noise is , the crack length measurement error is , and All obey a normal distribution with a mean of 0, that is, , ,in, is the covariance matrix; The state transfer equation is: ; The observation equation is: ; S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed at 10-minute intervals; The specific process of determining the annual average wind speed under the standard wind turbine category is as follows: According to the latest version of the wind power design standard IEC61400-1Y2019 issued by the International Electrotechnical Commission (IEC), the annual average wind speed under the standard wind turbine category is divided into three levels: 7.5m / s, 8.5m / s and 10m / s; The process of determining the cut-in wind speed and cut-out wind speed of the wind turbine is as follows: According to the National Renewable Energy Laboratory (NREL) of the United States, the cut-in wind speed of a 5MW wind turbine is 3m / s and the cut-out wind speed is 25m / s; The probability density function of the Rayleigh distribution is: (2); In the formula, is the average wind speed, unit: m / s; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; The specific process is: S4.1, if Figure 3 As shown in the figure, the turbulent wind speed data conforming to the Rayleigh distribution is updated every 10 minutes, and the time series load of the blade is calculated according to the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rain flow counting method is used to perform statistical analysis on these time series load data to obtain the number of cycles and the corresponding time series load amplitude range; S4.2, input the data in S4.1 into the blade crack growth model to calculate the blade crack growth rate; S4.3, multiply the rate in S4.2 by the time of the cycle to obtain the crack extension length of the blade in this cycle; S4.4, adding the extension length in this cycle to the blade crack length in the previous cycle, to obtain the crack length at the current moment; S4.5, repeat S4.1-S4.4, and finally obtain the observation data of the blade crack length in the whole time period; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; Bayesian prediction assumes that the posterior probability density function of a given state is , then the state prior distribution at the current k moment is for: (3); In the formula, are independent observed variables; is an unknown variable of the system and obeys a first-order Markov process; is the state transition probability density function of the system; is the posterior probability density obtained at the previous moment; Bayesian updating is to obtain the observation value at the current k moment through the system observation model , update and adjust the system's prior probability density function , realize the transformation from the prior probability density function to the posterior probability density function at time k Derivation of: (4); Where: is the likelihood function, which describes the possibility of the unknown variable appearing under given observations; is the normalization constant; The specific process of Bayesian prediction using the unscented Kalman filter is as follows: Step 1: Obtain the Sigma point set and its corresponding weights according to formula (5); (5); Where: is a scaling function used to reduce the total prediction error; is an estimated value; Step 2: Calculate the one-step prediction of 2n+1 Sigma point sets according to formula (6); (6); Step 3: Calculate the one-step prediction and covariance matrix of the system state quantity according to equations (7) and (8); (7); (8); Where: is the weight; Step 4: Based on the one-step prediction value, use the untraceable transformation again according to formula (9) to generate a new Sigma point set; (9); Step 5: Substitute the predicted Sigma point set in S4 into the observation equation according to formula (10) to obtain the predicted value of the Sigma point set; (10); Step 6: According to equations (11), (12) and (13), the predicted value of the Sigma point set is obtained from S5, and the mean and covariance of the system prediction are obtained by weighted summation; (11); (12); (13); Step 7: Calculate the Kalman gain matrix according to formula (14); (14); Step 8: Calculate the state update and covariance update of the system according to equations (15) and (16); (15); (16); S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6, continue iterative operation, such as Figure 4 shown.

[0022] Example 6 The remaining life prediction method of a wind turbine blade based on unscented Kalman filtering proposed in this embodiment includes the following steps: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; The equation for blade crack growth is: (1); In the formula, a is the current blade crack length; R is the stress ratio, i.e., the ratio of the root mean square value of the minimum stress to the maximum stress in the stress cycle; A , m All are blade material related parameters; is the stress intensity factor, The calculation of requires blade cycle stress, which is obtained by rain flow counting method, such as Figure 1 As shown; The system noise is , the crack length measurement error is , and All obey a normal distribution with a mean of 0, that is, , ,in, is the covariance matrix; The state transfer equation is: ; The observation equation is: ; S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed at 10-minute intervals; The specific process of determining the annual average wind speed under the standard wind turbine category is as follows: According to the latest version of the wind power design standard IEC61400-1Y2019 issued by the International Electrotechnical Commission (IEC), the annual average wind speed under the standard wind turbine category is divided into three levels: 7.5m / s, 8.5m / s and 10m / s; The process of determining the cut-in wind speed and cut-out wind speed of the wind turbine is as follows: According to the National Renewable Energy Laboratory (NREL) of the United States, the cut-in wind speed of a 5MW wind turbine is 3m / s and the cut-out wind speed is 25m / s; The probability density function of the Rayleigh distribution is: (2); In the formula, is the average wind speed, unit: m / s; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; The specific process is: S4.1, if Figure 3 As shown in the figure, the turbulent wind speed data conforming to the Rayleigh distribution is updated every 10 minutes, and the time series load of the blade is calculated according to the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rain flow counting method is used to perform statistical analysis on these time series load data to obtain the number of cycles and the corresponding time series load amplitude range; S4.2, input the data in S4.1 into the blade crack growth model to calculate the blade crack growth rate; S4.3, multiply the rate in S4.2 by the time of the cycle to obtain the crack extension length of the blade in this cycle; S4.4, adding the extension length in this cycle to the blade crack length in the previous cycle, to obtain the crack length at the current moment; S4.5, repeat S4.1-S4.4, and finally obtain the observation data of the blade crack length in the whole time period; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; Bayesian prediction assumes that the posterior probability density function of a given state is , then the state prior distribution at the current k moment is for: (3); In the formula, are independent observed variables; is an unknown variable of the system and obeys a first-order Markov process; is the state transition probability density function of the system; is the posterior probability density obtained at the previous moment; Bayesian updating is to obtain the observation value at the current k moment through the system observation model , update and adjust the system's prior probability density function , realize the transformation from the prior probability density function to the posterior probability density function at time k Derivation of: (4); Where: is the likelihood function, which describes the possibility of the unknown variable appearing under given observations; is the normalization constant; The specific process of Bayesian prediction using the unscented Kalman filter is as follows: Step 1: Obtain the Sigma point set and its corresponding weights according to formula (5); (5); Where: is a scaling function used to reduce the total prediction error; is an estimated value; Step 2: Calculate the one-step prediction of 2n+1 Sigma point sets according to formula (6); (6); Step 3: Calculate the one-step prediction and covariance matrix of the system state quantity according to equations (7) and (8); (7); (8); Where: is the weight; Step 4: Based on the one-step prediction value, use the untraceable transformation again according to formula (9) to generate a new Sigma point set; (9); Step 5: Substitute the predicted Sigma point set in S4 into the observation equation according to formula (10) to obtain the predicted value of the Sigma point set; (10); Step 6: According to equations (11), (12) and (13), the predicted value of the Sigma point set is obtained from S5, and the mean and covariance of the system prediction are obtained by weighted summation; (11); (12); (13); Step 7: Calculate the Kalman gain matrix according to formula (14); (14); Step 8: Calculate the state update and covariance update of the system according to equations (15) and (16); (15); (16); S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6, continue iterative operation, such as Figure 4 As shown; The wind speed was set according to the Rayleigh distribution of the IEC standard. In order to simulate the operation of the NREL 5MW wind turbine in the full wind speed range, the simulation was performed at intervals of 1m / s, and the turbulence intensity was set to Class C. Below and above the rated wind speed, the wind turbine operated according to the maximum wind energy tracking control and constant power control strategies respectively; Parameters related to blade materials , , , , wind speed v = 12m / s, turbulence intensity is C level, maximum wind energy tracking and constant power control strategy are implemented for NREL 5MW wind turbine to calculate the time series load of blade flapping moment, and the statistical results of rain flow counting method are as follows Figure 5 The blade crack length observations at annual average wind speeds of 10m / s, 8.5m / s, 7.5m / s and constant wind speed under the IEC standard are simulated, as shown in Figure 1. Figure 6 As shown in the figure, the damage threshold is set to 0.2m and the initial crack length is 0.03m. As can be seen from the figure, the higher the annual average wind speed, the faster the crack expansion speed of the wind turbine blade. When the annual average wind speed is 10m / s, the wind turbine reaches the damage threshold after 1.5 years of operation. The wind turbine model in this paper adopts a three-blade wind turbine, so the cracks at the three blade nodes 1 and 6 are observed, as shown in the figure. Figure 7 As shown, it can be seen that the crack propagation speed at the root of the blade is faster than that at 30% of the distance from the root of the blade. The cracks at the root of the blade are concentrated in 1.5 years and reach the damage threshold, while the cracks at 30% of the distance from the root of the blade are concentrated in 2.5 years and reach the damage threshold; Figure 8 As shown in Figure 1, the crack propagation of blades is predicted under an annual average wind speed of 10 m / s. Figure 8 It can be seen that when the measurement noise is large, the crack prediction results are still very close to the observed values. Fig. 9 The distribution of blade life. For different initial crack lengths a=0.01m, a=0.03m, a=0.05m, a=0.07m, a=0.09m, under the same load conditions and at the same location, the crack extension is also predicted under an annual average wind speed of 10m / s. Fig.10 The blade crack growth prediction curves for multiple initial crack lengths are shown. It can be seen that when the initial crack length is longer, the estimated remaining service life is lower.

Claims

1. A method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering, characterized in that: The following steps are involved: S1. Establish a blade damage model, run the wind turbine, and cracks begin to appear on the blades and gradually expand; S2, transform the blade crack extension into a discrete accumulation process, introduce system noise and crack length measurement error, and obtain the state transfer equation and observation equation of blade damage; S3, determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the wind speed at the hub of the wind turbine; S4. Obtaining observation data of blade crack length according to Rayleigh distribution; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transfer equation and the observation equation. Set the covariance matrix of the system noise and the crack length measurement error. and ; S6, use the unscented Kalman filter method to perform Bayesian prediction and combine it with the observation equation for Bayesian update; S7, set the blade crack extension damage threshold, determine whether the updated state exceeds the set damage threshold, if it exceeds, calculate the remaining service life, if not, return to S6 and continue iterative operation.

2. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The equation for blade crack growth described in S2 is: （1）； In the formula, a is the current blade crack length; R is the stress ratio, i.e., the ratio of the root mean square value of the minimum stress to the maximum stress in the stress cycle; A , m All are blade material related parameters; is the stress intensity factor, The calculation of requires blade cycle stress, which is obtained by the rainflow counting method.

3. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The system noise described in S2 is , the crack length measurement error is , and All obey a normal distribution with a mean of 0, that is, , ,in, is the covariance matrix.

4. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The state transfer equation described in S2 is: ; The observation equation is: .

5. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The specific process of determining the annual average wind speed under the standard wind turbine category described in S3 is as follows: According to the latest version of the wind power design standard IEC61400-1Y2019 issued by the International Electrotechnical Commission (IEC), the annual average wind speed under the standard wind turbine category is divided into three levels: 7.5m / s, 8.5m / s and 10m / s.

6. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The process of determining the cut-in wind speed and cut-out wind speed of the wind turbine generator set described in S3 is as follows: according to the regulations of the National Renewable Energy Laboratory (NREL) of the United States, the cut-in wind speed of a 5MW wind turbine generator set is 3m / s, and the cut-out wind speed is 25m / s.

7. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The probability density function of the Rayleigh distribution described in S3 is: （2）； In the formula, is the average wind speed, unit: m / s.

8. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The specific process of S4 is as follows: S4.

1. Update the turbulent wind speed data that conforms to the Rayleigh distribution every 10 minutes, and calculate the time series load of the blades according to the maximum wind energy tracking and constant power control strategy. After a complete cycle, use the rain flow counting method to perform statistical analysis on these time series load data to obtain the number of cycles and the corresponding time series load amplitude range; S4.2, input the data in S4.1 into the blade crack growth model to calculate the blade crack growth rate; S4.3, multiply the rate in S4.2 by the time of the cycle to obtain the crack extension length of the blade in this cycle; S4.4, adding the extension length in this cycle to the blade crack length in the previous cycle, to obtain the crack length at the current moment; S4.

5. Repeat S4.1-S4.4 to finally obtain the observation data of the blade crack length in the entire time period.

9. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The Bayesian prediction in S6 assumes that the posterior probability density function of a given state is , then the state prior distribution at the current k moment is for: （3）； In the formula, are independent observed variables; is an unknown variable of the system and obeys a first-order Markov process; is the state transition probability density function of the system; is the posterior probability density obtained at the previous moment; The Bayesian update is to obtain the observation value at the current k moment through the system observation model , update and adjust the system's prior probability density function , realize the transformation from the prior probability density function to the posterior probability density function at time k Derivation of: （4）； Where: is the likelihood function, which describes the possibility of the unknown variable appearing under given observations; is the normalization constant.

10. The method for predicting the remaining life of a wind turbine blade based on unscented Kalman filtering according to claim 1, characterized in that: The specific process of using the unscented Kalman filter method to perform Bayesian prediction is as follows: Step 1: Obtain the Sigma point set and its corresponding weights according to formula (5); （5）； Where: is a scaling function used to reduce the total prediction error; is an estimated value; Step 2: Calculate the one-step prediction of 2n+1 Sigma point sets according to formula (6); （6）； Step 3: Calculate the one-step prediction and covariance matrix of the system state quantity according to equations (7) and (8); （7）； （8）； Where: is the weight; Step 4: Based on the one-step prediction value, use the untraceable transformation again according to formula (9) to generate a new Sigma point set; （9）； Step 5: Substitute the predicted Sigma point set in S4 into the observation equation according to formula (10) to obtain the predicted value of the Sigma point set; （10）； Step 6: According to equations (11), (12) and (13), the predicted value of the Sigma point set is obtained from S5, and the mean and covariance of the system prediction are obtained by weighted summation; （11）； （12）； （13）； Step 7: Calculate the Kalman gain matrix according to formula (14); （14）； Step 8: Calculate the state update and covariance update of the system according to equations (15) and (16); （15）； （16）。

Images