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

Abstract

The wind turbine blade residual life prediction method based on the unscented Kalman filtering disclosed by the application comprises the following steps: S1, a blade damage model is established, a wind turbine is operated, and cracks begin to appear on the blade and gradually expand; S2, the blade crack expansion is converted into a discrete cumulative process to obtain a state transition equation and an observation equation; S3, the annual average wind speed, the cut-in wind speed, the cut-out wind speed and the Rayleigh distribution are determined; S4, the observation data of the blade crack length are obtained; S5, the damage state process is defined, and the covariance matrix of the system noise and the crack length measurement error is set; S6, the Bayesian prediction is performed by using the unscented Kalman filtering method, and the Bayesian update is performed; S7, the blade crack expansion damage threshold is set, whether the updated state exceeds the set damage threshold is judged, and the residual service life is calculated. The prediction method solves the problem that it is difficult to coordinate the error reduction and the calculation complexity in the prior art.

Description

Technical Field

[0001] This invention belongs to the field of wind power generation technology, specifically relating to a method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering. Background Technology

[0002] Wind turbines consist of multiple components, and their blades are affected by various factors such as wind shear, tower shadow effect, turbulent wind, gravity load, centrifugal load, and surface erosion. Therefore, blade failure accounts for approximately 23% of all wind turbine failures, and replacing blades leads to prolonged downtime and increased maintenance costs. Current wind turbines are equipped with numerous sensors to collect operational data. This data is not only used for monitoring the turbine itself but also for assessing the performance and health of various components based on predictive and health management systems. Accurate prediction of the remaining service life of the turbine blades is particularly crucial, playing a vital role in ensuring normal turbine operation and preventing accidents.

[0003] The key to predicting the remaining service life of wind turbine blades lies in predicting the evolution of fatigue damage and addressing the nonlinearity and uncertainty of crack propagation. Currently, commonly used methods include the Extended Kalman Filter (EKF) and the Particle Filter (PFF). While the EKF has shown some effectiveness in estimating nonlinear systems, it also has limitations. Specifically, when the system's nonlinearity is low, the error introduced by the EKF is not significant; however, for strongly nonlinear systems, the EKF's linearization process, which directly ignores higher-order terms, leads to larger errors. In contrast, the PFF is also applicable to nonlinear systems and can effectively address the problem of large errors. However, the PFF has a relatively high computational cost and a complex implementation process. A feasible solution to reconcile the contradiction between reducing error and computational complexity remains elusive. Summary of the Invention

[0004] The purpose of this invention is to provide a method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering, which solves the problem of difficulty in coordinating the reduction of error and computational complexity in the prior art.

[0005] The technical solution adopted in this invention is a method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering, comprising the following steps:

[0006] S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand.

[0007] S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage.

[0008] S3. Determine the Rayleigh distribution of annual average wind speed, wind turbine cut-in wind speed, wind turbine cut-out wind speed, and wind speed at the hub.

[0009] S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length;

[0010] S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ;

[0011] S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update.

[0012] S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it exceeds the threshold, calculate the remaining service life. If it does not exceed the threshold, return to S6 and continue iterative operation.

[0013] The invention is further characterized by:

[0014] The equation for the crack propagation in blade S2 is:

[0015] (1);

[0016] In the formula, a This represents the current blade crack length; R The stress ratio is the ratio of the root mean square value of the minimum stress to the maximum stress in a stress cycle. A , m All parameters are related to blade materials; Stress intensity factor The calculation requires the cyclic stress of the blade, which is obtained through the rainflow counting method.

[0017] The system noise in S2 is The crack length measurement error is , and They all follow a normal distribution with a mean of 0, that is... , ,in, Let be the covariance matrix.

[0018] The state transition equation in S2 is: The observation equation is: .

[0019] The specific process for 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 published 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.

[0020] The process for determining the cut-in and cut-out wind speeds of the S3 wind turbine 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.

[0021] The probability density function of the Rayleigh distribution in S3 is:

[0022] (2);

[0023] In the formula, Average wind speed, unit: m / s.

[0024] The specific process of S4 is as follows:

[0025] S4.1 Update the turbulent wind speed data conforming to the Rayleigh distribution every 10 minutes, and calculate the time-series load of the blades based on the maximum wind energy tracking and constant power control strategy. After a complete cycle, use the rainflow 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.

[0026] S4.2 Input the data from S4.1 into the blade crack propagation model and calculate the blade crack propagation rate;

[0027] S4.3 Multiply the rate in S4.2 by the time of the cycle to obtain the blade crack propagation length in this cycle;

[0028] S4.4 Add the length of the crack extension in the current cycle to the length of the blade crack in the previous cycle to get the length of the crack at the current moment.

[0029] S4.5. Repeat S4.1-S4.4 to obtain the observation data of the blade crack length over the entire time period.

[0030] Bayesian prediction in S6 assumes that the posterior probability density function of a given state is... Then the prior distribution of the state at time k is... for:

[0031] (3);

[0032] In the formula, These are mutually independent observation variables; The system is an unknown variable and follows a first-order Markov process; Let be the state transition probability density function of the system; This is the posterior probability density obtained at the previous time step;

[0033] Bayesian update obtains the observation value at time k through a system observation model. Update and adjust the prior probability density function of the system. To realize the transformation from the prior probability density function at time k to the posterior probability density function. Derivation:

[0034] (4);

[0035] In the formula: The likelihood function describes the probability of an unknown variable occurring given observations. This is the normalization constant.

[0036] The specific process of Bayesian prediction using the unscented Kalman filter method is as follows:

[0037] Step 1: Obtain the Sigma point set and its corresponding weights according to equation (5);

[0038] (5);

[0039] In the formula: This is a scaling function used to reduce the total prediction error; This is an estimated value;

[0040] Step 2: Calculate the one-step prediction of the 2n+1 Sigma point set according to equation (6);

[0041] (6);

[0042] Step 3: Calculate the one-step prediction and covariance matrix of the system state variables according to equations (7) and (8);

[0043] (7);

[0044] (8);

[0045] In the formula: For weights;

[0046] Step 4: Based on the predicted value in the first step, use the unscented transformation again according to equation (9) to generate a new Sigma point set;

[0047] (9);

[0048] Step 5: Substitute the predicted Sigma point set in S4 into the observation equation according to equation (10) to obtain the predicted value of the Sigma point set;

[0049] (10);

[0050] Step 6: Obtain the predicted value of the Sigma point set from S5 according to equations (11), (12) and (13), and obtain the mean and covariance of the system prediction by weighted summation;

[0051] (11);

[0052] (12);

[0053] (13);

[0054] Step 7: Calculate the Kalman gain matrix according to equation (14);

[0055] (14);

[0056] Step 8: Calculate the system state update and covariance update according to equations (15) and (16);

[0057] (15);

[0058] (16).

[0059] The beneficial effects of this invention are:

[0060] The present invention provides a method for predicting the remaining service life of wind turbine blades based on unscented Kalman filtering. This method collects blade flapping moment data across the entire wind speed range and 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. It incorporates various actual operating conditions of wind turbine units, making the crack length observation more realistic. The method uses unscented Kalman filtering to predict the evolution of fatigue damage, addressing the nonlinearity and uncertainty issues of crack propagation. Compared with other methods, it has lower computational complexity, smaller errors, and more accurate predicted remaining service life. Attached Figure Description

[0061] Figure 1 This is a simplified diagram of the rainflow counting method in the wind turbine blade remaining life prediction method based on unscented Kalman filtering of this invention.

[0062] Figure 2 This is a Rayleigh distribution map of wind speed under different annual average wind speeds according to the present invention;

[0063] Figure 3This is a flowchart of the blade crack observation process of the present invention;

[0064] Figure 4 This is a flowchart of the blade crack prediction process of the present invention;

[0065] Figure 5 This is a stress range distribution diagram of the rainflow counting method in Embodiment 6 of the present invention;

[0066] Figure 6 This is a comparison diagram of crack observation in constant wind speed range and variable wind speed range in Embodiment 6 of the present invention;

[0067] Figure 7 This is a comparison image of the crack observation at the critical node of the three-blade blade in Embodiment 6 of the present invention;

[0068] Figure 8 This is the blade crack prediction diagram based on unscented Kalman filtering in Embodiment 6 of the present invention;

[0069] Figure 9 This is a blade lifetime distribution diagram based on unscented Kalman filtering in Embodiment 6 of the present invention;

[0070] Figure 10 This is a blade crack prediction diagram with different initial lengths in Embodiment 6 of the present invention. Detailed Implementation

[0071] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.

[0072] Example 1

[0073] The wind turbine blade remaining life prediction method based on unscented Kalman filtering proposed in this embodiment includes the following steps:

[0074] S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand.

[0075] S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage.

[0076] S3. Determine the Rayleigh distribution of annual average wind speed, wind turbine cut-in wind speed, wind turbine cut-out wind speed, and wind speed at the hub.

[0077] S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length;

[0078] S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ;

[0079] S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update.

[0080] S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it exceeds the threshold, calculate the remaining service life. If it does not exceed the threshold, return to S6 and continue iterative operation.

[0081] Example 2

[0082] The wind turbine blade remaining life prediction method based on unscented Kalman filtering proposed in this embodiment includes the following steps:

[0083] S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand.

[0084] S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage.

[0085] The equation for blade crack propagation is:

[0086] (1);

[0087] In the formula, a This represents the current blade crack length; R The stress ratio is the ratio of the root mean square value of the minimum stress to the maximum stress in a stress cycle. A , m All parameters are related to blade materials; Stress intensity factor The calculation requires the cyclic stress of the blade, which is obtained through the rainflow counting method, such as... Figure 1 As shown;

[0088] System noise is The crack length measurement error is , and They all follow a normal distribution with a mean of 0, that is... , ,in, It is the covariance matrix;

[0089] The state transition equation is: The observation equation is: ;

[0090] S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed of the wind turbine, the cut-out wind speed, and the wind speed at the hub, such as Figure 2As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed over a 10-minute interval;

[0091] The specific process for 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 published 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.

[0092] The process for determining the cut-in wind speed and cut-out wind speed of wind turbines is as follows: According to the regulations of the National Renewable Energy Laboratory (NREL) in 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.

[0093] The probability density function of the Rayleigh distribution is:

[0094] (2);

[0095] In the formula, Average wind speed, unit: m / s;

[0096] S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length;

[0097] S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ;

[0098] S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update.

[0099] S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it exceeds the threshold, calculate the remaining service life. If it does not exceed the threshold, return to S6 and continue iterative operation.

[0100] Example 3

[0101] The wind turbine blade remaining life prediction method based on unscented Kalman filtering proposed in this embodiment includes the following steps:

[0102] S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand.

[0103] S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage.

[0104] The equation for blade crack propagation is:

[0105] (1);

[0106] In the formula, a This represents the current blade crack length; R The stress ratio is the ratio of the root mean square value of the minimum stress to the maximum stress in a stress cycle. A , m All parameters are related to blade materials; Stress intensity factor The calculation requires the cyclic stress of the blade, which is obtained through the rainflow counting method, such as... Figure 1 As shown;

[0107] System noise is The crack length measurement error is , and They all follow a normal distribution with a mean of 0, that is... , ,in, It is the covariance matrix;

[0108] The state transition equation is: The observation equation is: ;

[0109] S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed of the wind turbine, the cut-out wind speed, and the wind speed at the hub, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed over a 10-minute interval;

[0110] The specific process for 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 published 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.

[0111] The process for determining the cut-in wind speed and cut-out wind speed of wind turbines is as follows: According to the regulations of the National Renewable Energy Laboratory (NREL) in 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.

[0112] The probability density function of the Rayleigh distribution is:

[0113] (2);

[0114] In the formula, Average wind speed, unit: m / s;

[0115] S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length;

[0116] The specific process is as follows:

[0117] S4.1, such as Figure 3 As shown, turbulent wind speed data conforming to Rayleigh distribution is updated every 10 minutes, and the time-series load of the blade is calculated based on the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rainflow counting method is used to statistically analyze these time-series load data to obtain the number of cycles and the corresponding time-series load amplitude range.

[0118] S4.2 Input the data from S4.1 into the blade crack propagation model and calculate the blade crack propagation rate;

[0119] S4.3 Multiply the rate in S4.2 by the time of the cycle to obtain the blade crack propagation length in this cycle;

[0120] S4.4 Add the length of the crack extension in the current cycle to the length of the blade crack in the previous cycle to get the length of the crack at the current moment.

[0121] S4.5 Repeat S4.1-S4.4 to finally obtain the observation data of the blade crack length over the entire time period;

[0122] S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ;

[0123] S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update.

[0124] S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it does, calculate the remaining service life; otherwise, return to S6 and continue iteratively. Figure 4 As shown.

[0125] Example 4

[0126] The wind turbine blade remaining life prediction method based on unscented Kalman filtering proposed in this embodiment includes the following steps:

[0127] S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand.

[0128] S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage.

[0129] The equation for blade crack propagation is:

[0130] (1);

[0131] In the formula, a This represents the current blade crack length; R The stress ratio is the ratio of the root mean square value of the minimum stress to the maximum stress in a stress cycle. A , m All parameters are related to blade materials; Stress intensity factor The calculation requires the cyclic stress of the blade, which is obtained through the rainflow counting method, such as... Figure 1 As shown;

[0132] System noise is The crack length measurement error is , and They all follow a normal distribution with a mean of 0, that is... , ,in, It is the covariance matrix;

[0133] The state transition equation is: The observation equation is: ;

[0134] S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed of the wind turbine, the cut-out wind speed, and the wind speed at the hub, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed over a 10-minute interval;

[0135] The specific process for 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 published 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.

[0136] The process for determining the cut-in wind speed and cut-out wind speed of wind turbines is as follows: According to the regulations of the National Renewable Energy Laboratory (NREL) in 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.

[0137] The probability density function of the Rayleigh distribution is:

[0138] (2);

[0139] In the formula, Average wind speed, unit: m / s;

[0140] S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length;

[0141] The specific process is as follows:

[0142] S4.1, such as Figure 3 As shown, turbulent wind speed data conforming to Rayleigh distribution is updated every 10 minutes, and the time-series load of the blade is calculated based on the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rainflow counting method is used to statistically analyze these time-series load data to obtain the number of cycles and the corresponding time-series load amplitude range.

[0143] S4.2 Input the data from S4.1 into the blade crack propagation model and calculate the blade crack propagation rate;

[0144] S4.3 Multiply the rate in S4.2 by the time of the cycle to obtain the blade crack propagation length in this cycle;

[0145] S4.4 Add the length of the crack extension in the current cycle to the length of the blade crack in the previous cycle to get the length of the crack at the current moment.

[0146] S4.5 Repeat S4.1-S4.4 to finally obtain the observation data of the blade crack length over the entire time period;

[0147] S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ;

[0148] S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update.

[0149] Bayesian prediction assumes that the posterior probability density function of a given state is... Then the prior distribution of the state at time k is... for:

[0150] (3);

[0151] In the formula, These are mutually independent observation variables; The system is an unknown variable and follows a first-order Markov process; Let be the state transition probability density function of the system; This is the posterior probability density obtained at the previous time step;

[0152] Bayesian update obtains the observation value at time k through a system observation model. Update and adjust the prior probability density function of the system. To realize the transformation from the prior probability density function at time k to the posterior probability density function. Derivation:

[0153] (4);

[0154] In the formula: The likelihood function describes the probability of an unknown variable occurring given observations. This is the normalization constant;

[0155] S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it exceeds the threshold, calculate the remaining service life. If it does not exceed the threshold, return to S6 and continue iterative operation.

[0156] Example 5

[0157] The wind turbine blade remaining life prediction method based on unscented Kalman filtering proposed in this embodiment includes the following steps:

[0158] S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand.

[0159] S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage.

[0160] The equation for blade crack propagation is:

[0161] (1);

[0162] In the formula, a This represents the current blade crack length; R The stress ratio is the ratio of the root mean square value of the minimum stress to the maximum stress in a stress cycle. A , m All parameters are related to blade materials; Stress intensity factor The calculation requires the cyclic stress of the blade, which is obtained through the rainflow counting method, such as... Figure 1 As shown;

[0163] System noise is The crack length measurement error is , and They all follow a normal distribution with a mean of 0, that is... , ,in, It is the covariance matrix;

[0164] The state transition equation is: The observation equation is: ;

[0165] S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed of the wind turbine, the cut-out wind speed, and the wind speed at the hub, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed over a 10-minute interval;

[0166] The specific process for 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 published 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.

[0167] The process for determining the cut-in wind speed and cut-out wind speed of wind turbines is as follows: According to the regulations of the National Renewable Energy Laboratory (NREL) in 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.

[0168] The probability density function of the Rayleigh distribution is:

[0169] (2);

[0170] In the formula, Average wind speed, unit: m / s;

[0171] S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length;

[0172] The specific process is as follows:

[0173] S4.1, such as Figure 3 As shown, turbulent wind speed data conforming to Rayleigh distribution is updated every 10 minutes, and the time-series load of the blade is calculated based on the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rainflow counting method is used to statistically analyze these time-series load data to obtain the number of cycles and the corresponding time-series load amplitude range.

[0174] S4.2 Input the data from S4.1 into the blade crack propagation model and calculate the blade crack propagation rate;

[0175] S4.3 Multiply the rate in S4.2 by the time of the cycle to obtain the blade crack propagation length in this cycle;

[0176] S4.4 Add the length of the crack extension in the current cycle to the length of the blade crack in the previous cycle to get the length of the crack at the current moment.

[0177] S4.5 Repeat S4.1-S4.4 to finally obtain the observation data of the blade crack length over the entire time period;

[0178] S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ;

[0179] S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update.

[0180] Bayesian prediction assumes that the posterior probability density function of a given state is... Then the prior distribution of the state at time k is... for:

[0181] (3);

[0182] In the formula, These are mutually independent observation variables; The system is an unknown variable and follows a first-order Markov process; Let be the state transition probability density function of the system; This is the posterior probability density obtained at the previous time step;

[0183] Bayesian update obtains the observation value at time k through a system observation model. Update and adjust the prior probability density function of the system. To realize the transformation from the prior probability density function at time k to the posterior probability density function. Derivation:

[0184] (4);

[0185] In the formula: The likelihood function describes the probability of an unknown variable occurring given observations. This is the normalization constant;

[0186] The specific process of Bayesian prediction using the unscented Kalman filter method is as follows:

[0187] Step 1: Obtain the Sigma point set and its corresponding weights according to equation (5);

[0188] (5);

[0189] In the formula: This is a scaling function used to reduce the total prediction error; This is an estimated value;

[0190] Step 2: Calculate the one-step prediction of the 2n+1 Sigma point set according to equation (6);

[0191] (6);

[0192] Step 3: Calculate the one-step prediction and covariance matrix of the system state variables according to equations (7) and (8);

[0193] (7);

[0194] (8);

[0195] In the formula: For weights;

[0196] Step 4: Based on the predicted value in the first step, use the unscented transformation again according to equation (9) to generate a new Sigma point set;

[0197] (9);

[0198] Step 5: Substitute the predicted Sigma point set in S4 into the observation equation according to equation (10) to obtain the predicted value of the Sigma point set;

[0199] (10);

[0200] Step 6: Obtain the predicted value of the Sigma point set from S5 according to equations (11), (12) and (13), and obtain the mean and covariance of the system prediction by weighted summation;

[0201] (11);

[0202] (12);

[0203] (13);

[0204] Step 7: Calculate the Kalman gain matrix according to equation (14);

[0205] (14);

[0206] Step 8: Calculate the system state update and covariance update according to equations (15) and (16);

[0207] (15);

[0208] (16);

[0209] S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it does, calculate the remaining service life; otherwise, return to S6 and continue iteratively. Figure 4As shown.

[0210] Example 6

[0211] The wind turbine blade remaining life prediction method based on unscented Kalman filtering proposed in this embodiment includes the following steps:

[0212] S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand.

[0213] S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage.

[0214] The equation for blade crack propagation is:

[0215] (1);

[0216] In the formula, a This represents the current blade crack length; R The stress ratio is the ratio of the root mean square value of the minimum stress to the maximum stress in a stress cycle. A , m All parameters are related to blade materials; Stress intensity factor The calculation requires the cyclic stress of the blade, which is obtained through the rainflow counting method, such as... Figure 1 As shown;

[0217] System noise is The crack length measurement error is , and They all follow a normal distribution with a mean of 0, that is... , ,in, It is the covariance matrix;

[0218] The state transition equation is: The observation equation is: ;

[0219] S3. Determine the Rayleigh distribution of the annual average wind speed, the cut-in wind speed of the wind turbine, the cut-out wind speed, and the wind speed at the hub, such as Figure 2 As shown, the Rayleigh distribution determines the possible values ​​of the average wind speed over a 10-minute interval;

[0220] The specific process for 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 published 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.

[0221] The process for determining the cut-in wind speed and cut-out wind speed of wind turbines is as follows: According to the regulations of the National Renewable Energy Laboratory (NREL) in 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.

[0222] The probability density function of the Rayleigh distribution is:

[0223] (2);

[0224] In the formula, Average wind speed, unit: m / s;

[0225] S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length;

[0226] The specific process is as follows:

[0227] S4.1, such as Figure 3 As shown, turbulent wind speed data conforming to Rayleigh distribution is updated every 10 minutes, and the time-series load of the blade is calculated based on the maximum wind energy tracking and constant power control strategy. After a complete cycle, the rainflow counting method is used to statistically analyze these time-series load data to obtain the number of cycles and the corresponding time-series load amplitude range.

[0228] S4.2 Input the data from S4.1 into the blade crack propagation model and calculate the blade crack propagation rate;

[0229] S4.3 Multiply the rate in S4.2 by the time of the cycle to obtain the blade crack propagation length in this cycle;

[0230] S4.4 Add the length of the crack extension in the current cycle to the length of the blade crack in the previous cycle to get the length of the crack at the current moment.

[0231] S4.5 Repeat S4.1-S4.4 to finally obtain the observation data of the blade crack length over the entire time period;

[0232] S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ;

[0233] S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update.

[0234] Bayesian prediction assumes that the posterior probability density function of a given state is... Then the prior distribution of the state at time k is... for:

[0235] (3);

[0236] In the formula, These are mutually independent observation variables; The system is an unknown variable and follows a first-order Markov process; Let be the state transition probability density function of the system; This is the posterior probability density obtained at the previous time step;

[0237] Bayesian update obtains the observation value at time k through a system observation model. Update and adjust the prior probability density function of the system. To realize the transformation from the prior probability density function at time k to the posterior probability density function. Derivation:

[0238] (4);

[0239] In the formula: The likelihood function describes the probability of an unknown variable occurring given observations. This is the normalization constant;

[0240] The specific process of Bayesian prediction using the unscented Kalman filter method is as follows:

[0241] Step 1: Obtain the Sigma point set and its corresponding weights according to equation (5);

[0242] (5);

[0243] In the formula: This is a scaling function used to reduce the total prediction error; This is an estimated value;

[0244] Step 2: Calculate the one-step prediction of the 2n+1 Sigma point set according to equation (6);

[0245] (6);

[0246] Step 3: Calculate the one-step prediction and covariance matrix of the system state variables according to equations (7) and (8);

[0247] (7);

[0248] (8);

[0249] In the formula: For weights;

[0250] Step 4: Based on the predicted value in the first step, use the unscented transformation again according to equation (9) to generate a new Sigma point set;

[0251] (9);

[0252] Step 5: Substitute the predicted Sigma point set in S4 into the observation equation according to equation (10) to obtain the predicted value of the Sigma point set;

[0253] (10);

[0254] Step 6: Obtain the predicted value of the Sigma point set from S5 according to equations (11), (12) and (13), and obtain the mean and covariance of the system prediction by weighted summation;

[0255] (11);

[0256] (12);

[0257] (13);

[0258] Step 7: Calculate the Kalman gain matrix according to equation (14);

[0259] (14);

[0260] Step 8: Calculate the system state update and covariance update according to equations (15) and (16);

[0261] (15);

[0262] (16);

[0263] S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it does, calculate the remaining service life; otherwise, return to S6 and continue iteratively. Figure 4 As shown;

[0264] Based on the Rayleigh distribution of the IEC standard, wind speeds were set. To simulate the operation of the NREL 5MW wind turbine across the entire wind speed range, simulations were conducted at 1 m / s intervals, with turbulence intensity set to level C. Below and above the rated wind speed, the wind turbine operated under maximum wind energy tracking control and constant power control strategies, respectively.

[0265] Blade material related parameters , , , With a wind speed of v=12m / s and turbulence intensity of C, maximum wind energy tracking and constant power control strategies were implemented on the NREL 5MW wind turbine to calculate the time-series load of blade flapping moment. Statistical results using the rainflow counting method are as follows: Figure 5 As shown, simulations were performed to observe the blade crack length at annual average wind speeds of 10 m / s, 8.5 m / s, 7.5 m / s, and a constant wind speed under the IEC standard. Figure 6 As shown, the damage threshold was set to 0.2m, and the initial crack length was 0.03m. The figure shows that the higher the average annual wind speed, the faster the crack propagation rate of the wind turbine blades. Specifically, at an average annual wind speed of 10m / s, the wind turbine reaches the damage threshold after 1.5 years of operation. This paper uses a three-bladed wind turbine model, so cracks at nodes 1 and 6 of the three blades are observed. Figure 7 As shown, the crack propagation rate at the leaf root is faster than that at 30% of the leaf root. Cracks at the leaf root reach the damage threshold in 1.5 years, while cracks at 30% of the leaf root reach the damage threshold in 2.5 years. Figure 8 As shown, the propagation of blade cracks under an annual average wind speed of 10 m / s is predicted. From... Figure 8 As can be seen, even when the measurement noise is large, the crack prediction results are still very close to the observed values. Figure 9 The distribution of blade lifetime is shown. For different initial crack lengths a=0.01m, a=0.03m, a=0.05m, a=0.07m, and a=0.09m, under the same load conditions and at the same location, the crack propagation under the same annual average wind speed of 10m / s is predicted. Figure 10 The blade crack propagation prediction curves for multiple initial crack lengths are shown. It can be seen that when the initial crack length is long, the estimated remaining service life is low.

Claims

1. A method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering, characterized in that, Includes the following steps: S1. Establish a blade damage model, run the wind turbine, and the blades begin to show cracks that gradually expand. S2. The blade crack propagation is transformed into a discrete cumulative process, and system noise and crack length measurement error are introduced to obtain the state transition equation and observation equation for blade damage. S3. Determine the Rayleigh distribution of annual average wind speed, wind turbine cut-in wind speed, wind turbine cut-out wind speed, and wind speed at the hub. S4. Based on the Rayleigh distribution, obtain the observation data of the blade crack length; S5. Define the damage state process as a function of the previous damage state. The damage state process includes the state transition equation and the observation equation. Set the covariance matrix of system noise and crack length measurement error. and ; S6. Use the unscented Kalman filter to perform Bayesian prediction, and combine it with the observation equation to perform Bayesian update. S7. Set the blade crack propagation damage threshold, and determine whether the updated state exceeds the set damage threshold. If it exceeds the threshold, calculate the remaining service life. If it does not exceed the threshold, return to S6 and continue iterative operation. The equation for blade crack propagation described in S2 is: （1）； In the formula, a This represents the current blade crack length; R The stress ratio is the ratio of the root mean square value of the minimum stress to the maximum stress in a stress cycle. A , m All parameters are related to blade materials; Stress intensity factor The calculation requires the cyclic stress of the blade, which is obtained through the rainflow counting method; The specific process of S4 is as follows: S4.1 Update the turbulent wind speed data conforming to the Rayleigh distribution every 10 minutes, and calculate the time-series load of the blades based on the maximum wind energy tracking and constant power control strategy. After a complete cycle, use the rainflow 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 from S4.1 into the blade crack propagation model and calculate the blade crack propagation rate; S4.3 Multiply the rate in S4.2 by the time of the cycle to obtain the blade crack propagation length in this cycle; S4.4 Add the length of the crack extension in the current cycle to the length of the blade crack in the previous cycle to get the length of the crack at the current moment. S4.

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

2. The method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering according to claim 1, characterized in that, The system noise mentioned in S2 is The crack length measurement error is The and They all follow a normal distribution with a mean of 0, that is... , ,in, Let be the covariance matrix.

3. The method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering according to claim 1, characterized in that, The state transition equation described in S2 is: The observation equation is: .

4. The method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering according to claim 1, characterized in that, The specific process for 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 published by the International Electrotechnical Commission, 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.

5. The method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering according to claim 1, characterized in that, The process for determining the cut-in wind speed and cut-out wind speed of the wind turbine 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 is 3m / s and the cut-out wind speed is 25m / s.

6. The method for predicting the remaining life of wind turbine blades 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, Average wind speed, unit: m / s.

7. The method for predicting the remaining life of wind turbine blades 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 prior distribution of the state at time k is... for: （3）； In the formula, These are mutually independent observation variables; The system is an unknown variable and follows a first-order Markov process; Let be the state transition probability density function of the system; This is the posterior probability density obtained at the previous time step; The Bayesian update obtains the observation value at time k through a system observation model. Update and adjust the prior probability density function of the system. To realize the transformation from the prior probability density function at time k to the posterior probability density function. Derivation: （4）； In the formula: The likelihood function describes the probability of an unknown variable occurring given observations. This is the normalization constant.

8. The method for predicting the remaining life of wind turbine blades based on unscented Kalman filtering according to claim 1, characterized in that, The specific process of Bayesian prediction using the unscented Kalman filter method is as follows: Step 1: Obtain the Sigma point set and its corresponding weights according to equation (5); （5）； In the formula: This is a scaling function used to reduce the total prediction error; This is an estimated value; Step 2: Calculate the one-step prediction of the 2n+1 Sigma point set according to equation (6); （6）； Step 3: Calculate the one-step prediction and covariance matrix of the system state variables according to equations (7) and (8); （7）； （8）； In the formula: For weights; Step 4: Based on the predicted value in the first step, use the unscented transformation again according to equation (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 equation (10) to obtain the predicted value of the Sigma point set; （10）； Step 6: Obtain the predicted value of the Sigma point set from S5 according to equations (11), (12) and (13), and obtain the mean and covariance of the system prediction by weighted summation; （11）； （12）； （13）； Step 7: Calculate the Kalman gain matrix according to equation (14); （14）； Step 8: Calculate the system state update and covariance update according to equations (15) and (16); （15）； （16）。

Images