Initial damage prediction method for interlaminar bonding defects of wind turbine blades under extreme wind conditions

Abstract

The present application proposes a method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions, which belongs to the technical field of wind turbine status monitoring, including: inputting sample points of wind characteristic parameters into a sensitivity analysis model to obtain the sensitivity of wind characteristic parameters to interlayer bonding defects at the main beam of the blade, the sensitivity analysis model is established using a neural network and the Sobol method, and the KL divergence of the strain residual energy density at the maximum damage point and the non-maximum damage point is used as training data; the first N wind characteristic parameters with the highest first-order Sobol effect index and total effect Sobol index are used as highly sensitive factors; the actual measurement range of the highly sensitive factors is input into the cohesive force finite element model of the interlayer bonding defect of the blade to obtain the strain residual energy density, and the strain residual energy density is used as the damage degree of the blade. The present application can identify high-risk areas for early damage caused by interlayer bonding defects of wind turbine blades under extreme wind conditions.

Description

Technical Field

[0001] The present invention belongs to the technical field of wind turbine generator set status monitoring, and in particular relates to a method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions. Background Art

[0002] Blades are key components of wind turbines. With the increasing trend of large-scale and offshore wind turbines, the maintenance cost and failure rate of blades are also increasing. After the performance of blades deteriorates, they need to be repaired, hoisted, transported, and replaced, which brings huge downtime losses to wind farms and greatly reduces the power generation income. The main bearing structure of wind turbine blades is mostly glass fiber reinforced plastic (GFRP). The most commonly used molding process is vacuum infusion. It is very easy to produce manufacturing defects due to molding process factors. Among them, due to process factors such as contamination of the adhesive interface between glass fiber cloth layers, uneven coating of adhesives, and aging of adhesives, it is easy to produce bonding defects between layers of blade composite materials. This is a very common manufacturing defect of wind turbine blades, resulting in the blades after molding and curing are not intact. During the service process, the stress concentration phenomenon will appear at the interlayer bonding defects of the blades and gradually evolve into fatigue damage. Especially when the wind turbine stops running under extreme wind conditions in the field environment, the mechanical properties of the blades will also deteriorate and degenerate rapidly, making it difficult to exert the proper mechanical properties in a long-term and stable manner.

[0003] According to the IEC61400-13 standard, the wind measurement parameters involved in the description of wind conditions include five wind measurement parameters: wind speed, wind direction, wind shear, turbulence intensity and wind frequency. Among them, wind speed plays a significant role in the changes in the aerodynamic loads of various components of the blades. At the same time, the wind direction, wind shear, turbulence intensity and wind frequency cause the wind turbine tower to cause the average wind speed of the wind turbine blades to be non-constant, which also causes the aerodynamic loads to fluctuate frequently. In recent years, most scholars have focused on studying the aerodynamic response of blades in a rotating state, while there is little research on the wind-induced aerodynamic response of blades in a stopped state. When a wind turbine is in a stopped state, the centrifugal stiffening effect of its blades disappears. Therefore, the bending stiffness of the wind turbine blades is low in the stopped state, and it is more susceptible to the drastic changes in strong wind speed and wind direction, resulting in a significant change in the wind attack angle of the blade airfoil, thereby causing aerodynamic instability, and then causing damage or even destruction to the blade structure. The aerodynamic load of the blades seriously exceeds the load-bearing load under normal wind conditions, and the local structure is subjected to large deflection and bending, which is prone to sudden fracture failure, resulting in the collapse of the entire wind turbine, causing the collapse or damage of the wind turbine unit. If the highly sensitive wind characteristic parameters are clear under extreme wind conditions, timely regulation can be made in the shutdown state to adapt to extreme wind conditions to avoid blade damage, or the most dangerous parts of wind turbine blades with bonding defects can be predicted in advance to take safety measures, that is, preventive measures can be taken in the initial period of risk to make initial damage predictions.

[0004] However, there are many factors that affect the evolution of interlayer bonding defects into fatigue damage. The service conditions of wind turbines are also very complex. The blade airfoil structures and ply designs of different manufacturers, in previous blade fracture accidents caused by extreme weather such as typhoons, still lack a quantitative analysis method for the degree of blade damage based on wind parameters. If the sensitive factors of blade damage containing interlayer bonding defects under extreme service wind conditions can be evaluated, it will help provide a new method for establishing intelligent and healthy operation and maintenance control strategies. At present, the shutdown load control methods of wind turbines, brake systems, yaw systems, pitch control systems, etc. are mostly based on experience to set threshold parameters such as extreme wind speed and wind direction, which often causes the internal bonding defects of the blades to evolve into initial damage, greatly reducing the service life of the blades. In order to prevent the interlayer bonding defects that are common inside the blades from causing high stress concentration and resulting in damage evolving into local cracking or buckling, on the one hand, the wind turbine shutdown and load reduction control strategy needs to clarify the dominant factor among the five wind characteristic parameters that is most sensitive to damage evolution; on the other hand, wind turbine blades with quality defects need to quantitatively evaluate the degree of damage under extreme wind conditions and predict the risk of failure, so as to coordinate with other components to formulate a reasonable operation and inspection plan.

[0005] Therefore, the damage evolution law of wind turbine blades under extreme wind conditions was studied, and the highly sensitive wind characteristic parameters for the damage evolution of interlayer bonding defects in blades under extreme wind conditions were screened out. The influence of wind speed, wind direction, wind shear and turbulence intensity on blade damage was quantitatively evaluated. Early health monitoring of initial delamination damage of blades can be achieved by monitoring wind speed and wind direction data. Summary of the invention

[0006] In view of the shortcomings of the prior art, this application proposes a method for predicting damage of interlayer bonding defects of wind turbine blades under extreme wind conditions, including:

[0007] Step S1: obtaining sample points of wind characteristic parameters according to the extreme wind conditions in the wind farm;

[0008] Step S2: inputting the sample points of the wind characteristic parameters into a pre-established sensitivity analysis model to obtain the sensitivity of the wind characteristic parameters to the interlayer bonding defects at the main beam of the blade, wherein the pre-established sensitivity analysis model is established by using a BP neural network and a Sobol method, and the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is used as training data;

[0009] Step S3: Calculate the first-order Sobol effect index and the total effect Sobol index according to the wind characteristic parameters and the sensitivity, so that the first N wind characteristic parameters with the highest first-order Sobol effect index and / or total effect Sobol index are used as high-sensitivity factors;

[0010] Step S4: Input the actual measurement range of the highly sensitive factor into the pre-established cohesive force finite element model of the interlayer bonding defect of the blade to obtain the strain complementary energy density of the damaged unit in the delaminated damage area. The strain complementary energy density is used as the damage degree of the blade, and the ratio between the damage degree of the blade and the design load of the defect-free blade is calculated. If the ratio is greater than a preset threshold, the prediction result is that initial damage occurs at the corresponding position of the blade.

[0011] The method of obtaining the sample points of the wind characteristic parameters according to the extreme wind conditions in the wind farm includes:

[0012] Determine wind characteristic parameters of wind turbine blades and parameter ranges of wind characteristic parameters according to the extreme wind conditions in the wind farm, wherein the wind characteristic parameters include wind speed, wind direction, wind shear, turbulence intensity and wind frequency;

[0013] Performing orthogonal design on the wind characteristic parameters and parameter intervals of the wind characteristic parameters, and selecting representative parameter combinations and corresponding parameter intervals;

[0014] According to the representative parameter combination and the corresponding parameter interval, the Sobol sequence sampling method is used to obtain the sample points of the wind characteristic parameters.

[0015] The training data acquisition process includes:

[0016] According to the sample points of the wind characteristic parameters, under the action of different wind characteristic parameter combinations, the time domain load data of each airfoil center of the blade is calculated;

[0017] Converting the time domain load data into the distribution of equivalent stress at different cross-sectional positions of a full-size blade;

[0018] The equivalent stress is input as an external load into the pre-established cohesive finite element model of the interlaminar bonding defect of the blade to perform delamination damage evolution and obtain the stress and strain response data of delamination damage.

[0019] According to the stress and strain response data of the delamination damage, the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is calculated, and the KL divergence is used as training data for the pre-established sensitivity analysis model.

[0020] The method of calculating the time domain load data of each airfoil center of the blade based on the sample points of the wind characteristic parameters under the action of different wind characteristic parameter combinations includes: sequentially importing the sample points of the wind characteristic parameters into the load calculation software to establish a wind model, inputting the design parameters of the wind turbine and the actual extreme service wind conditions of the blades, and calculating the time domain load data of the blade flapping and swing directions at the center of each airfoil.

[0021] The step of converting the time domain load data into the distribution of equivalent stress at different cross-sectional positions of a full-size blade includes:

[0022] According to the time domain load data of each airfoil center of the blade, the normal stress and shear stress of the blade section are calculated; the normal stress σ0 of the blade section is calculated as follows:

[0023]

[0024] Among them, M x is the aerodynamic bending moment in the blade flapping direction, M y is the aerodynamic bending moment in the swing direction, W n is the bending section coefficient;

[0025] The calculation formula of blade section shear stress τ0 is as follows:

[0026]

[0027] Among them, F x is the shear force in the blade swinging direction, F y is the shear force in the swing direction, A is the cross-sectional area of ​​the blade,

[0028] According to the normal stress of the blade section and the shear stress of the blade section, the equivalent stress at different cross-sectional positions of the full-size blade is calculated, and the calculation formula is as follows:

[0029]

[0030] Where σ is the equivalent stress at different cross-sectional positions of the full-size blade;

[0031] According to the equivalent stress at different cross-sectional positions of the full-size blade, an equivalent stress load curve of the full-size blade is drawn, wherein the position of the highest stress area in the full-size blade equivalent stress load curve is the maximum equivalent stress of the blade;

[0032] The ultimate strength of the blade is determined according to the maximum equivalent stress of the blade.

[0033] The equivalent stress is input as an external load into the pre-established cohesive force finite element model of the interlaminar bonding defect position of the blade to perform delamination damage evolution and obtain stress and strain response data of delamination damage, including:

[0034] The equivalent stress is input as an external load into the pre-established cohesive force finite element model of the interlayer bonding defect of the blade to perform delamination damage evolution, wherein the interface unit introduces an energy-based bilinear cohesive force model to describe the mechanical behavior of the evolution of the interlayer bonding defect, and the association between the stress and relative displacement of the damage unit is established to describe the evolution process of the delamination damage;

[0035] In the evolution of delamination damage, the stress and strain response data of delamination damage are collected, and the failure proportion of damaged units at different times is calculated. The time of initial delamination damage, damage area and stress distribution are determined according to the proportion of failed units, and the stress and displacement difference curve is plotted. The stress and displacement difference curve is integrated to obtain the strain energy density of the damaged failure unit.

[0036] The step of calculating the KL divergence of the strain complementary energy density at the maximum damage location and the non-maximum damage location according to the stress and strain response data of the layer damage includes:

[0037] The strain complementary energy density at the maximum damage is calculated based on the energy release rate at the maximum damage and the critical fracture energy at the maximum damage. The strain complementary energy density at the non-maximum damage is calculated based on the energy release rate at the non-maximum damage and the critical fracture energy at the non-maximum damage.

[0038] The probability distribution function of strain complementary energy at the maximum damage is obtained according to the strain complementary energy density at the maximum damage, and the probability distribution function of strain complementary energy at the non-maximum damage is obtained according to the strain complementary energy density at the non-maximum damage;

[0039] According to the probability distribution function of strain complementary energy at the maximum damage point and the probability distribution function of strain complementary energy at the non-maximum damage point, the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is calculated.

[0040] The pre-established sensitivity analysis model is established using BP neural network and Sobol method, including:

[0041] Determine the BP neural network topology structure, including: determine the number of neurons in the input layer, the number of neurons in the output layer, the number of hidden layers and the number of neurons in each hidden layer, wherein the number of hidden layers is not equal to the number of neurons in each hidden layer;

[0042] The Sobol sequence method is used to randomly extract sample points of wind characteristic parameters as the input of the BP neural network. The KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is used as the output of the BP neural network for fitting, and multiple sensitivity analysis models are obtained.

[0043] By comparing the output error statistical parameters of multiple sensitivity analysis models, the one with the smallest error is selected as the pre-established sensitivity analysis model.

[0044] The first-order Sobol effect index is calculated as follows:

[0045]

[0046] Among them, S iis the first-order Sobol effect index, which indicates the wind characteristic parameter X i The influence of a single action on the initial damage evolution of bonding defects, D Xi (E X～i (Y|X i )) is the wind characteristic parameter X i The variance caused by, D(Y) is the variance of sensitivity Y, E{·} is the mathematical expectation, X ～i to consider all random values ​​of variable i while keeping variable i fixed.

[0047] The total effect Sobol index is calculated as follows:

[0048]

[0049] Among them, S Ti is the total effect Sobol index, E X～i (D Xi (Y|X ～i )) is to determine the exclusion of X ～i The mean variance of sensitivity Y assuming the true values ​​of all variables except .

[0050] Beneficial effects:

[0051] The present application proposes a method for predicting damage caused by interlayer bonding defects in wind turbine blades under extreme wind conditions. The method of the present application is based on blade design parameters and actual wind conditions in local wind farms. It has the ability to dynamically detect blade quality defects based on design and service conditions. It can identify high-risk areas prone to early damage based on blade design parameters, defect types, and geometric morphology, and guide detection equipment to perform precise measurements, fault tracing, and optimized design. It also helps to provide control thresholds for wind characteristic parameters when formulating shutdown and load reduction control strategies under extreme wind conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 Flow chart of a method for predicting damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to an embodiment of the present application;

[0053] Figure 2 Schematic diagram of a method for predicting damage of interlayer bonding defects of a wind turbine blade under extreme wind conditions according to an embodiment of the present application;

[0054] Figure 3 Coordinate system and external force diagram of a full-size blade in an embodiment of the present application;

[0055] Figure 4 The flapping bending moment of the sixth airfoil section of the blade in the embodiment of the present application;

[0056] Figure 5Equivalent stress distribution of each airfoil section of a full-size blade in an embodiment of the present application;

[0057] Figure 6 Finite element model of a blade laminate containing interlayer bonding defects according to an embodiment of the present application;

[0058] Figure 7 The process in which the interlayer bonding defect of the embodiment of the present application evolves into the initial delamination damage;

[0059] Figure 8 Schematic diagram of the difference in initial damage distribution of the sixth airfoil section of the blade main beam in the embodiment of the present application;

[0060] Fig. 9 Schematic diagram of the BP neural network topology structure of the embodiment of the present application;

[0061] Fig.10 The first-order effect index and total effect index of the wind characteristic factor of the embodiment of the present application;

[0062] Fig.11 The first-order effect index and the total effect index mean of the wind characteristic factors of different cross sections in the embodiment of the present application;

[0063] Fig.12 Initial delamination damage prediction of the main airfoil section of the blade main beam part in the embodiment of the present application. DETAILED DESCRIPTION

[0064] The specific implementation of the present application is further described in detail below in conjunction with the drawings and examples.

[0065] In response to the uncertainty of under what extreme wind conditions wind turbine blades, which generally contain interlayer bonding defects, will evolve into initial damage, this application proposes a damage prediction method for interlayer bonding defects of wind turbine blades under extreme wind conditions. A sensitivity analysis method is used to screen out highly sensitive factors in the input wind characteristic parameters from the output deterministic indicators, and then the key inducements for the evolution of bonding defects into delamination damage under extreme wind conditions are obtained, thereby reducing the uncertainty and complexity of the quantitative assessment of initial damage, and then the damage degree of the blade is quantitatively assessed within the wind measurement parameter range of the actual wind field and the maximum risk position is predicted.

[0066] Firstly, the distribution characteristics of wind characteristic parameters of extreme wind conditions in the local wind field were collected to obtain sample points. The time domain load data of each airfoil center of the full-size blade under the combined action of different sampling points were calculated using wind load calculation software and converted into cross-sectional equivalent stress distribution. Secondly, a cohesive zone model (CZM, Cohesive Zone Model) of the blade bonding defect area was established, and the maximum cross-sectional equivalent stress was loaded on the finite element model to simulate the initial damage evolution process and obtain the stress and strain response data of the damage evolution process. Then, a sensitivity analysis model combining the BP neural network (Back Propagation Neural Network) and the Sobol method was established to calculate the KL divergence (KLD, Kullback–Leibler divergence) is used as the sensitivity analysis response index of initial damage to screen the highly sensitive factors that affect the evolution of wind turbine blade bonding defects into initial damage; finally, based on the actual extreme wind conditions of the local wind farm, the actual measurement range of the highly sensitive factors is input on the basis of the cohesive zone (CZM) finite element model of the blade bonding defect, and the strain residual energy density of the damaged unit in the damaged area is quantitatively calculated to evaluate the damage degree and predict the location of the maximum risk.

[0067] Example:

[0068] This application proposes a method for predicting the damage of interlayer bonding defects of wind turbine blades under extreme wind conditions, such as Figure 1 , Figure 2 As shown, including:

[0069] Step S1: obtaining sample points of wind characteristic parameters according to the extreme wind conditions in the wind farm;

[0070] In this embodiment, the local extreme wind conditions of the wind farm, the wind characteristic influencing factors and distribution information of wind turbine blade damage are first determined, and sample points are obtained accordingly. According to the sample points of wind characteristic parameters, the load calculation software is used to calculate the time domain load data of each airfoil center of the blade under the combined action of extreme wind characteristic factors, and converted into equivalent stress distribution at different cross-sectional positions of the full-size blade;

[0071] The method of obtaining the sample points of the wind characteristic parameters according to the extreme wind conditions in the wind farm includes:

[0072] Step S1.1: determining wind characteristic parameters of wind turbine blades and parameter ranges of wind characteristic parameters according to the extreme wind conditions in the wind farm, wherein the wind characteristic parameters include wind speed, wind direction, wind shear, turbulence intensity and wind frequency;

[0073] This application takes a 2MW wind turbine blade in service and in extreme wind conditions as an example. The air density ρ of the local wind farm is 1.225kg / m3 The main design parameters of the wind turbine are shown in Table 1:

[0074] Table 1 Main design parameters of wind turbines

[0075]

[0076]

[0077] The distribution characteristics of the local wind characteristic parameters (wind speed, wind direction, wind shear, turbulence intensity, and wind frequency) are determined. The extreme changes in these wind characteristic parameters cause drastic changes in the aerodynamic load on the blade surface. In order to prevent the common interlayer bonding defects inside the blade from causing high stress concentration and damage evolution into local cracking or buckling, on the one hand, the wind turbine shutdown and load reduction control strategy needs to clarify the dominant factor that is most sensitive to damage evolution among the five wind characteristic parameters. On the other hand, wind turbine blades with quality defects need to quantitatively evaluate the degree of damage under extreme wind conditions and predict the risk of failure, so as to coordinate with other components to formulate a reasonable operation and inspection plan.

[0078] In this embodiment, the wind speed is the 10-minute average maximum wind speed at the hub height that occurs once every 50 years. Since the wind speed of the wind turbine varies with the position, considering the extreme wind shear effect and the tower shadow effect, the wind speed referred to here is the wind speed at the center of the rotor. The extreme wind speed is a Gumbel distribution with an average value of 38m / s and a standard deviation of 7.65m / s. The wind speed model is as follows:

[0079]

[0080] Among them, v0 is the initial wind speed, v is the wind speed, the parameters required by the wind speed model are R, the rotor radius, h, the tower height, and v h is the speed at the tower height, C0 is the tower radius, κ is the roughness coefficient of the local wind farm ground, q is the turbulent kinetic energy dissipation rate, and u0 is the distance between the blade rotation plane and the tower centerline. The wind direction azimuth θ is set to 0 degrees in the due north direction.

[0081] Wind frequency reflects the frequency of change of uncertain wind speed, which is the cut-in wind speed V in and cut-out wind speed V out The wind speed probability density function P(V hub ), V hub is the average wind speed at the hub center. In order to apply fatigue cyclic loads to the defective parts of the blade, the Fast Fourier Transform (FFT) is used to integrate the time domain into the frequency domain over the entire time axis for the uncertain wind speed frequency, and the time domain components of the wind stress are converted into frequencies in the spectrum through the transfer equation:

[0082]

[0083] Among them, P(u) is the frequency domain function, u is the spatial frequency, and f(x) is the spatial domain function. Usually, within the measurement time, the frequency of uncertain wind speed follows a normal distribution.

[0084] Turbulence intensity I reflects the intensity of wind speed changes and is generally considered to be log-normally distributed. v is the standard deviation of horizontal wind speed relative to the average wind speed v. The turbulence intensity I is expressed as:

[0085]

[0086] The wind shear α is generally expressed using an empirical wind shear index, which is adjusted by increasing or decreasing the simple power law exponent U(H), where U(H) is the average wind speed at the reference height. Here, z is the height of the calculated section, and the reference height H is the hub center height, which is 100 meters in this embodiment, which is basically the same as the tower height. The wind shear has the following functional relationship:

[0087]

[0088] Step S1.2: performing orthogonal design on the wind characteristic parameters and parameter intervals of the wind characteristic parameters, and selecting representative parameter combinations and corresponding parameter intervals;

[0089] In this embodiment, when the wind turbine blades are shut down in extreme wind conditions, they are assumed to be rigid structures, the blade pitch is kept unchanged at 0 degrees, and the aeroelastic effect is not considered. At this time, the aerodynamic load on the wind turbine blades is determined by the wind field in which they are located. The parameter ranges of the five wind characteristic parameters are further determined in combination with the actual wind conditions of the local wind field, and representative parameter combinations are selected based on this orthogonal design, as shown in Table 2.

[0090] Considering the randomness and dispersion of the sample points of each wind characteristic parameter, the Gaussian probability distribution function is used to represent the distribution information h(x i ):

[0091]

[0092] Among them, x i represents the wind measurement parameters of each actual wind field (i=1,2,…5), μ i is the mean of the distribution interval of each wind measurement parameter, ω i is the standard deviation of the distribution interval of each wind measurement parameter. According to the mean μ of the distribution function of each wind measurement parameter i and standard deviation ω i , determine the sampling interval of each wind characteristic parameter.

[0093] Table 2 Sampling interval of wind characteristic parameters

[0094]

[0095]

[0096] Step S1.3: According to the representative parameter combination and the corresponding parameter interval, the sample points of the wind characteristic parameters are obtained by using the Sobol sequence sampling method.

[0097] Step S2: inputting the sample points of the wind characteristic parameters into a pre-established sensitivity analysis model to obtain the sensitivity of the wind characteristic parameters to the interlayer bonding defects at the main beam of the blade, wherein the pre-established sensitivity analysis model is established by using a BP neural network and a Sobol method, and the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is used as training data;

[0098] The training data acquisition process includes:

[0099] Step S2.11: Calculating the time domain load data of each airfoil center of the blade under the action of different wind characteristic parameter combinations according to the sample points of the wind characteristic parameters;

[0100] The method of calculating the time domain load data of each airfoil center of the blade according to the sample points of the wind characteristic parameters under the action of different wind characteristic parameter combinations includes: importing the sample points of the wind characteristic parameters into the GH-Bladed software in sequence to establish a wind model, inputting the design parameters of the wind turbine and the actual extreme service wind conditions of the blades, and calculating the time domain load data of the blade flapping and swing directions at the center of each airfoil.

[0101] In this embodiment, according to Table 2, the orthogonal design test method was used to generate 25 groups of optimal parameter combinations with 5 parameter combinations, which were successively imported into the calculation tool GH-Bladed based on the blade element-momentum theory to establish 25 wind models, and the dynamic aerodynamic loads at the center of each airfoil of the blade were calculated respectively. A time-varying wind field across the entire rotor swept area was established using the calculation tool for random time and space simulation, and the wind vibration coefficient β z Take 1, the body shape coefficient μ s Take it as 0.6, input the parameters of the wind turbine set to establish the horizontal axis wind turbine blade model and the load analysis model of the set, which mainly includes: (1) blades, hub, tower, transmission system, control system, etc.; (2) external environmental parameters such as land and sea environment and wind conditions; (3) parameters such as blade length, cross-sectional airfoil, cross-sectional distance, chord length, torsion angle and torsion center position.

[0102] Under the action of airflow, the blade will undergo alternating motion in two directions, namely the flapping direction and the swinging direction. The aerodynamic load is mainly composed of loads in these two directions. Figure 3 As shown, the loads in the blade flapping and swinging directions are calculated as:

[0103] (1) Shear force F in the blade flapping and swinging direction x and F y for:

[0104]

[0105] (2) Aerodynamic shear force:

[0106]

[0107] (3) Aerodynamic bending moment:

[0108]

[0109] Where x, y are the distances from blade r to the hub center; C L is the lift characteristic coefficient; C is the drag characteristic coefficient, which is usually determined by a semi-empirical curve; is the relative flow velocity, m / s; ρ is the air density, kg / m 3 , which is generally assumed to be a constant value; v is the incoming wind speed, m / s; ω is the speed of the wind turbine, rad / s; is the inflow angle, rad; r is the distance from the blade root to the cross section, m; C(r) is the blade chord length at position r, m, R is the rotor radius, C D is the characteristic coefficient of resistance.

[0110] The fatigue load of wind turbine blades is usually simulated by combining aerodynamics and structural dynamics with a 10-minute variable wind load to simulate the fatigue cycle load value within the cycle. The simulation time of each extreme wind condition is 600s, and the time step is 0.4×10 -3 In order to avoid the instability of blade response dynamic load caused by changes in wind speed and to prevent the existence of adopted data, the first 60s is the warm-up time.

[0111] Step S2.12: converting the time domain load data into the distribution of equivalent stress at different cross-sectional positions of the full-size blade;

[0112] In this embodiment, the GH-Blade tool is used to obtain the central aerodynamic loads of the airfoil section for the 25 wind characteristic combinations in Table 1, calculate the equivalent stress of each airfoil section of the full-size blade, and determine the maximum equivalent stress position of the blade and the ultimate strength of the blade.

[0113] The step of converting the time domain load data into the distribution of equivalent stress at different cross-sectional positions of a full-size blade includes:

[0114] Step S2.12.1: Calculate the normal stress and shear stress of the blade section according to the time domain load data of each airfoil center of the blade;

[0115] Generally, since the input external force of the aerodynamic model can only be based on experience, the discrete resultant force and moment at the aerodynamic center of each airfoil of the blade are provided for application in finite element analysis (FEA). This may cause two problems: first, the concentrated force ignores the real pressure distribution on the blade surface; second, local stress concentration may occur when using concentrated force. For this reason, this embodiment converts the resultant force and moment into an equivalent stress synthesized by loads in the bending, torsion and shear directions for fatigue damage analysis.

[0116] Combined with the blade coordinate system, the force in the blade swinging direction can be decomposed into the shear force F x and bending moment M y , the force in the swing direction is decomposed into M y and M x , where M y and M x Determines the normal stress in the blade section, F x and F y Determine the shear stress of the section. According to the load at the center of each airfoil section of the blade, the normal stress σ0 of the blade section is calculated as:

[0117]

[0118] The shear stress τ0 of the blade section is:

[0119]

[0120] Where Wn represents the bending section coefficient; A represents the blade cross-sectional area,

[0121] Step S2.12.2: Calculate the equivalent stress at different cross-sectional positions of the full-size blade according to the normal stress of the blade cross section and the shear stress of the blade cross section;

[0122] In this embodiment, after the stress condition of the blade cross section is obtained, the stress synthesis of the blade root cross section is performed according to the fourth strength theory of material mechanics, and the equivalent stress expression of each cross section is:

[0123]

[0124] Where σ is the equivalent stress of the cross section.

[0125] Step S2.12.1.3: Draw an equivalent stress load curve of the full-size blade according to the equivalent stress at different cross-sectional positions of the full-size blade. In the equivalent stress load curve of the full-size blade, the position where the highest stress area is located is the maximum equivalent stress of the blade;

[0126] Step S2.12.1.4: Determine the ultimate strength of the blade based on the maximum equivalent stress of the blade.

[0127] The wind turbine blade of this embodiment is designed with 17 airfoil sections, the span from the blade root to the airfoil is 5.1m, and the 6th section is the maximum chord length. Among them, the flapping direction bending moment M at the 6th airfoil section within 1s is x (calculated according to a combination of one set of wind characteristic parameters), such as Figure 4 The loads of 17 airfoil sections are calculated and substituted into formula (14) to obtain the equivalent stress mean σ of each section, and the full-scale equivalent stress load curve is plotted as shown in the following figure: Figure 5 shown.

[0128] from Figure 5 It can be seen that the 4th to 9th sections of the blade are local high stress areas, reaching 125Mpa at the 6th section, which is the maximum equivalent stress position. The 6th airfoil section of the blade is the maximum chord length of the blade. The geometry and material changes at this position are the greatest. The suction side will undergo large nonlinear bending, causing high stress concentration. With the dynamic changes of aerodynamic wind loads, it will lead to instability phenomena such as local buckling, delamination and debonding, which are the initial damage evolution. If this position already contains native interlayer bonding defects, it is most likely to exceed the ultimate strength, so Figure 5 The maximum equivalent stress position of the full-size blade and the ultimate strength of the blade can be determined.

[0129] Step S2.13: Input the equivalent stress as an external load into the pre-established cohesive finite element model of the interlaminar bonding defect of the blade to perform delamination damage evolution and obtain stress and strain response data of delamination damage, including:

[0130] Step S2.13.1: the equivalent stress is input as an external load into the pre-established cohesive force finite element model of the interlayer bonding defect of the blade to perform delamination damage evolution, wherein the interface unit introduces an energy-based bilinear cohesive force model to describe the mechanical behavior of the evolution of the interlayer bonding defect, and the association between the stress and relative displacement of the damage unit is established to describe the evolution process of the delamination damage;

[0131] In this embodiment, the geometric dimensions of the blade interlayer bonding defect cohesion (CZM) finite element model are: length 100mm, thickness 20mm, width 50mm, the glass fiber ply angle is 0 degrees, the model adopts uniaxial compression loading, the left end of the blade main beam laminate is fixed, and the displacement load is applied to the right end. The loading amount along the X direction is 0.01mm / min. The material properties of each layer are defined according to Table 3, Table 4, and Table 5. The eight-node hexahedral unit (C3D8) is used to mesh the fibers, matrix and homogenized area in the laminate, and a layer of eight-node cohesion unit (COH3D8) with a thickness of 10μm is arranged on the preset layering path in the middle plane of the laminate to simulate the interface debonding. The prefabricated size of the bonding defect is 5*5mm 2 , and define the automatic contact model at the same time. Figure 6 shown.

[0132] Table 3 Blade material parameters

[0133]

[0134] Table 4 Strength parameters of blade main beam laminate

[0135]

[0136] Table 5 Cohesive interface parameters of laminates

[0137]

[0138] The interface unit introduces an energy-based bilinear cohesion model to describe the mechanical behavior of the evolution of interlaminar bonding defects, and establishes the relationship between the stress and relative displacement of the damage unit to describe the evolution of delamination damage. The model is calculated in the explicit solver ABAQUS / Explicit, and a zero-thickness quasi-static cohesion unit is constructed based on the user-defined unit subroutine (UEL). The equivalent load stress of the sixth airfoil section in step 4 is selected as the model external load to simulate the initial delamination damage progressive evolution process of the interlaminar bonding defect of the blade.

[0139] When the interlayer bonding defect evolves into initial damage, the damage variable matrix is ​​expressed as follows: The initial stage of the interface is the linear elastic stage, and the constitutive relationship is:

[0140] σ=Dδ (15)

[0141] Among them, σ = {σ1, σ2, σ3} is the normal and tangential stress; δ = {δ1, δ2, δ3} is the normal and tangential relative displacement; D is the unit stiffness matrix.

[0142] Among them, the relative displacement δ m The expression is:

[0143]

[0144] Among them, δ shear is the interface tangential relative displacement; 〈〉 is the Macaulay operator.

[0145] Interface damage initial displacement for:

[0146]

[0147] Where N is the interface normal strength, S and T are the interface shear strength.

[0148] Interface complete failure displacement for:

[0149]

[0150] Among them,:G Ⅰ , G Ⅱ , G Ⅲ - are the energy release rates of mode I, mode II and mode III cracks respectively; G ⅠC , G ⅡC , G ⅢC The critical strain energy of mode I, mode II and mode III cracks, respectively, and the damage variable after the initiation of interlaminar interface damage are expressed as:

[0151]

[0152] Among them, the interface damage is considered to be irreversible, and the maximum relative displacement of the interface is defined as The maximum relative displacement that occurs during loading is defined as

[0153] Step S2.13.2: During the evolution of delamination damage, the stress and strain response data of delamination damage are collected, and the failure proportion of damaged units at different times is calculated. The time of initial delamination damage, damage area and stress distribution are determined according to the proportion of failed units, and the stress and displacement difference curve is plotted. The stress and displacement difference curve is integrated to obtain the strain energy density of the damaged failure unit.

[0154] In this embodiment, during the evolution of the delamination damage of the sixth airfoil section, the interface unit fails continuously and the strength of the blade defect position decreases continuously. The energy release rate G of the interface damage unit reaches the critical fracture energy G c , the cohesive unit fails completely and the damage spreads. Figure 7 The figure shows the process of interlayer bonding defects evolving into initial delamination damage.

[0155] In order to screen out the main contributing factors of extreme wind conditions to damage evolution, it is necessary to quantify the sensitivity response index to accurately quantify the influence of each wind characteristic parameter on the initial stratified damage, so as to solve the uncertainty and complexity of the initial damage evolution. According to the strain energy density theory, for the cohesion model, the damage degree of the damaged unit is represented by the strain complementary energy density (energy release rate G), which can be obtained by integrating the cohesion-displacement difference curve. Figure 6 The established finite element model collects stress σ and strain δ response data of the initial damage process, conducts extreme wind characteristic sensitivity analysis, calculates the failure proportion of damaged units at different times, determines the time of initial delamination damage, damage area and stress distribution according to the proportion of failed units, and plots the integral of the cohesion-displacement difference curve to obtain the strain energy density of the damaged failure unit, thereby realizing the prediction of the maximum damage failure risk of blades under extreme wind conditions.

[0156] Step S2.14: Calculate the KL divergence of the strain complementary energy density at the maximum damage location and the non-maximum damage location based on the stress and strain response data of the layer damage, and use the KL divergence as training data for the pre-established sensitivity analysis model, including:

[0157] Step S2.14.1: Calculate the strain complementary energy density at the maximum damage location based on the energy release rate at the maximum damage location and the critical fracture energy at the maximum damage location; calculate the strain complementary energy density at the non-maximum damage location based on the energy release rate at the non-maximum damage location and the critical fracture energy at the non-maximum damage location;

[0158] Step S2.14.2: Obtain the strain complementary energy probability distribution function at the maximum damage location according to the strain complementary energy density at the maximum damage location, and obtain the strain complementary energy probability distribution function at the non-maximum damage location according to the strain complementary energy density at the non-maximum damage location;

[0159] Step S2.14.3: Calculate the KL divergence of the strain complementary energy density at the maximum damage location and the strain complementary energy density at the non-maximum damage location based on the strain complementary energy probability distribution function at the maximum damage location and the strain complementary energy probability distribution function at the non-maximum damage location.

[0160] The form and order of initial damage are guided by the weakest point, that is, the maximum damage point. The degree of initial damage of the bonding defect can be expressed by the distribution difference of the strain complementary energy density between the maximum damage point and the non-maximum damage point. Therefore, quantifying the distribution difference of the strain complementary energy density between the maximum damage point and the non-maximum damage point is an effective way to obtain the sensitivity response index.

[0161] (1) The KL divergence is used to calculate the relative entropy of the strain complementary energy density of the maximum and non-maximum damaged units as a sensitivity index to measure the stress distribution difference rate in the key damage area, providing a new judgment index for quantitative analysis of the uncertainty of the initial damage mode. C , the cohesive unit fails completely, satisfying formula (20)

[0162] G≥G C (20)

[0163] Where G = G Ⅰ +G Ⅱ +G Ⅲ is the energy release rate of the hybrid model, G C is the critical fracture energy.

[0164] Usually G C The BK criterion proposed by Benzeggaph and Kenane is adopted. This criterion takes into account the coupling relationship of strain energy release rate under different loading modes, and its expression is as follows:

[0165]

[0166] Among them, G shear =G Ⅱ +G Ⅲ is the shear strain energy, G T =G Ⅰ +G shear， The power exponent of the η-BK criterion is 2.

[0167] (2) When the energy release rate G of the interface damage unit reaches the critical fracture energy Gc, the energy release rate of the damage unit exceeds the cohesive strength, the cohesive unit fails completely, and the damage expands.

[0168] The degree of damage of interlayer bonding defect is represented by the distribution difference of strain complementary energy density between the maximum damage site and the non-maximum damage site, which is an effective way to measure the sensitivity response index. The strain complementary energy density at the damage site is expressed as,

[0169] ΔG=G C -G (22)

[0170] Assume that the strain complementary energy density at the maximum damage point of delamination damage is ΔG max , the strain complementary energy density at the non-maximum damage point is ΔG r , ΔG max and ΔG r The distribution difference of can be obtained using the difference measurement method KL divergence.

[0171] ΔG max =max(ΔG) (23)

[0172]

[0173] (3) In the damage simulation process, the symmetric KL divergence in discrete form is used:

[0174]

[0175] Where p and q are ΔG max and ΔG r The probability distribution function of . KL(p∥q) represents ΔG max and ΔG r The distribution difference of is used as the sensitivity response index W of layered damage.

[0176] W=KL(p||q) (26)

[0177] Calculated from step 4 Figure 5 In the figure, high stress concentration occurs at the 4th to 9th airfoil sections of the blade, which are the locations where initial damage of bonding defects is prone to occur. The stress and strain response data of the initial delamination damage stage of the bonding defect under 25 groups of wind characteristic factor combinations of these 6 sections are obtained respectively, and the data are substituted into formulas (22-26) to calculate the distribution difference of the initial damage KL divergence of the 4th to 9th airfoil sections of the blade main beam, as shown in Figure 8 shown.

[0178] Figure 8 It shows that the distribution of the initial damage degree of different sections of the blade main beam is different. The median and mean of the KL divergence at the 5th and 6th sections are the highest, indicating that the damage stratification at these two locations is the most serious. The distance difference between the upper quartile and the lower quartile of the box plot at the 6th section is the largest, and the range of the 1.5IQR quartile difference is also the largest, indicating that the degree of stratification damage here is the largest. Therefore, it can be determined that the initial position of the stratification damage has the highest probability of occurring at the 6th section of the blade main beam. Therefore, this method can predict the risk location based on the design parameters of this type of blade and the extreme conditions of the wind field.

[0179] The pre-established sensitivity analysis model is established using a BP neural network and a variance-based GSA-Sobol method, including:

[0180] Step S2.21: determining the BP neural network topology structure, including: determining the number of neurons in the input layer, the number of neurons in the output layer, the number of hidden layers, and the number of neurons in each hidden layer, wherein the number of hidden layers is not equal to the number of neurons in each hidden layer;

[0181] Step S2.22: Using the Sobol sequence method to randomly extract sample points of wind characteristic parameters as the input of the BP neural network, and using the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point as the output of the BP neural network, fitting is performed to obtain multiple sensitivity analysis models;

[0182] Step S2.23: By comparing the output error statistical parameters of multiple sensitivity analysis models, the one with the smallest error is selected as the pre-established sensitivity analysis model.

[0183] Since this embodiment takes into account the design parameters of the blades and the service wind conditions, it is difficult to take into account the coupling effects between complex factors during finite element calculations. Based on the error gradient descent method, the BP neural network can establish a mathematical model of the mapping relationship between the input and output of a multivariate nonlinear function, and estimate the discrete or continuous unknown nonlinear response of any input, which is widely used in multi-parameter, nonlinear model prediction. Based on the training data set, the input signal is transferred to the output layer layer by layer through the forward propagation algorithm, and the error between the output value and the target value is calculated. Then, the back propagation algorithm is used to reversely transfer the error back to the input layer layer by layer, and the error is used to continuously adjust the weights and biases until an acceptable model prediction accuracy is reached, so as to achieve the purpose of continuously optimizing the network fitting effect, and is used to study the highly nonlinear relationship between design, wind parameters and the evolution of blade structure fatigue damage.

[0184] This embodiment uses a BP neural network to estimate the response surface model of the extreme wind characteristic parameter X and the KL divergence Y of the initial delamination damage of the interlayer bonding defect of the blade main beam. The wind characteristic parameter x obtained in step S1 i The sampling points and the sensitivity response index W corresponding to each sampling point obtained in step 5 are used as a training data set. The specific steps of establishing a sensitivity analysis model of wind turbine blade damage based on the BP neural network are as follows:

[0185] (1) Determine the topological structure of the BP neural network as follows Fig. 9 As shown, the number of neurons in the input layer and the number of neurons in the output layer, the number of hidden layers and the number of neurons in each hidden layer are determined. The 25 orthogonal design combination factors of extreme wind characteristic parameters (Table 2) are used as input, the sensitivity response index KL divergence is used as output, the BP neural network substitution model is trained, and the sensitivity analysis model is established. According to the input and output distribution information, the original sample space is divided and sampled.

[0186] (2) To ensure the randomness and independence of the sample points, the sample points in each interval are randomly sorted to ensure that the sampling points are evenly covered in each sample space dimension to avoid sampling bias. This embodiment adopts the quasi-random sampling technology of the Sobol sequence method (Sobol 1967) so that the data points can be evenly distributed in the sampling interval of the wind characteristic parameters in Table 2.

[0187] (3) Since increasing the number of hidden layers and neurons may improve the performance of the training data set and enhance the prediction accuracy of the neural network, more parameters are added during model training to better fit the response of the input data set. Several neural network models with different numbers of hidden layers and neurons are established. Fig. 9 The network topology structure is fitted and predicted, and the best sensitivity analysis model is screened by comparing the output statistical parameters of the network model. The three statistical parameters include mean absolute error MAE, root mean square error RMSE, and mean absolute percentage error MAPE. The smaller the error, the better the consistency with the data and the stronger the predictability. The error calculation formula is as follows:

[0188]

[0189] Where N is the total number of data points in the training data set or the test data set; P i and O i Represent the predicted value and observed value of a data point respectively.

[0190] This case uses a data set with a total of 5000 samples. To train the network and evaluate the prediction effect, the sample data points are randomly divided into a training set (90%, 4500 data points) and a test set (10%, 500 data points) in proportion. According to Tables 2, 3 and 4, the Sobol sequence method is used to sample 4500 of each of the five wind characteristic parameters and sensitivity response indicators. After preprocessing, the sample data points are trained using Matlab software to implement the BP neural network model. It has been verified that the BP neural model with 2 hidden layers and 24 neurons in each hidden layer has the best prediction performance for the test data, as shown in Table 6.

[0191] Table 6 Performance comparison of BP neural network on training and test data sets

[0192]

[0193] Step S3: Calculate the first-order Sobol effect index and the total effect Sobol index according to the wind characteristic parameters and the sensitivity, so that the first N wind characteristic parameters with the highest first-order Sobol effect index and / or total effect Sobol index are used as high-sensitivity factors;

[0194] In this embodiment, according to the sensitivity analysis model in step S2, highly sensitive factors of damage evolution of interlayer bonding defects of wind turbine blades under extreme wind conditions are screened out, and the influence of extreme wind characteristic parameters on the initial damage evolution of bonding defects is quantitatively evaluated. The variance analysis method can decompose the influence of each input and the interaction between inputs on the output. This method uses the total variance of the output to normalize the decomposition terms and defines it as the Sobol index.

[0195] Among them, the first-order Sobol effect index represents the factor x i The degree of influence of a single action on the response index. The first-order Sobol index (S i ) is expressed as:

[0196]

[0197] Among them, D Xi (E X～i (Y|X i )), is due to the single parameter (X i ) is the variance caused by the main effect, and D(Y) is the variance of the output Y. E{·} represents the mathematical expectation, X ～i Consider all random values ​​while keeping the variable i fixed.

[0198] The total effect index takes into account the total contribution of a specific input to the output change, which is the sum of the first-order effect index and the parameter interaction effect index. The higher the index, the greater the contribution of the parameter to the output uncertainty. Ti ) is expressed as:

[0199]

[0200] E X～i (D Xi (Y|X ～i )) is to determine the exclusion of X ～i The variance mean of Y under the premise of the true value of all variables except Ti The highest factor will be evaluated as the most influential factor, which is a highly sensitive factor. On the contrary, it is an insensitive factor.

[0201] Fig.10 The sensitivity effect index of the initial damage stage of the interlayer bonding defect at the 6th section is calculated. The results show that wind speed is the most sensitive factor, followed by wind direction, which is a relatively sensitive factor. Wind frequency is in the middle, and turbulence intensity and wind shear are at the end, which are relatively insensitive factors. It can be screened that wind speed and wind direction are the most critical factors for the initial damage evolution of the blade main beam bonding defect under extreme wind conditions, and are also important monitoring parameters when serving in extreme weather. Once the design value is exceeded, an alarm will be issued.

[0202] The blade main beam laminate undergoes nonlinear bending due to extreme swinging loads, causing different degrees of initial delamination damage to the bonding defects within the layer. The highly sensitive factors of damage evolution in different sections may not be exactly the same. In order to verify whether the highly sensitive factors of the maximum risk section (the 6th airfoil section) are correct, this embodiment calculates the first-order effect index and the average value of the total effect index of the wind characteristic factors of the 4th to 9th airfoil sections respectively, and compares them with the calculation results of the maximum risk position (the 6th airfoil section). Fig.11 shown.

[0203] By comparison, it was found that the most sensitive factors for other interfaces are wind speed and wind direction. However, the difference is that the wind frequency is more sensitive to the 9th section than other sections. The main reason is that under extreme wind speed and direction conditions, as the section is farther away from the blade root, the degree of blade bending is higher and the amplitude is more obvious, causing the section on the suction side of the main beam to be subjected to higher pressure stress, aggravating the interface stratification of bonding defects. Therefore, this section is more sensitive to changes in wind frequency.

[0204] According to the level range of wind characteristic factors in Table 2, the low-sensitivity parameters of wind frequency, wind shear and turbulence intensity are fixed and taken as the mean value of Table 2. At the same time, the value range of sensitive factors wind speed and wind direction is adjusted as input, and substituted into the BP neural network model trained in step S2 to predict the initial stratification damage degree, and obtain the data point Y of KL divergence. P , and compared with the original KL divergence data point Y under the full factor control condition. To facilitate comparative observation, the original KL divergence data point Y of the full factor is used as the horizontal coordinate, and an orthogonal solid line is connected to divide the coordinate area into two, highlighting the predicted value Y P The regression degree with the original data point Y, the comparison results are as follows Fig.12 .

[0205] pass Fig.12 The results show that although the damage prediction result data points do not always cross the original data solid line, they are relatively close. The three statistical parameter indicators are shown in Table 7. It can be seen that the data consistency is high, indicating the accuracy of the initial delamination damage prediction results. This shows that the health monitoring of the initial delamination damage of the main beam can be achieved when only monitoring the wind speed and wind direction data, and it also proves the correctness of the sensitivity analysis to screen the damage causes.

[0206] Table 7 Prediction accuracy of initial damage of six different sections of the blade

[0207]

[0208] Step S4: Input the actual measurement range of the highly sensitive factor into the pre-established cohesive force finite element model of the interlayer bonding defect of the blade to obtain the strain complementary energy density of the damaged unit in the delaminated damage area. The strain complementary energy density is used as the damage degree of the blade, and the ratio between the damage degree of the blade and the design load of the defect-free blade is calculated. If the ratio is greater than a preset threshold, the prediction result is that initial damage occurs at the corresponding position of the blade.

[0209] In this embodiment, targeting the actual extreme wind conditions in the local wind farm, the actual measurement range of the highly sensitive factors is input based on the cohesive force (CZM) finite element model of the blade bonding defect, and the stress and strain data of the delaminated damage area are obtained. The data are compared with the design load of the defect-free blade, thereby quantitatively evaluating the initial damage evolution degree of the interlayer bonding defect in the wind turbine blade.

[0210] By evaluating the sensitive factors of defects turning into damage, real-time monitoring of wind measurement parameter sensitive factors, and formulating control strategies to prevent wind turbine damage, we can avoid the evolution of destructive damage to the blades in extreme weather, and reduce the frequency and cost of blade maintenance. Repeat steps S2.12 to S2.13. The damage degree of the damaged unit can be expressed by the strain residual energy density, that is, Figure 4 The area enclosed by the middle curve and the coordinates can be used to quantitatively evaluate the damage degree of the interlayer bonding defects of the blade under the shutdown state under extreme wind conditions by integrating the cohesion-displacement difference curve.

[0211] This method is based on blade design parameters and actual wind conditions in local wind farms. It has the ability to dynamically detect blade quality defects based on design and service conditions. It can identify high-risk areas that are prone to early damage based on blade design parameters, defect types, and geometric morphology, and guide detection equipment to perform precise measurements, fault tracing, and optimized design. It also helps provide control thresholds for wind characteristic parameters when formulating shutdown and load reduction control strategies under extreme wind conditions.

[0212] The various embodiments in the present application are described in a progressive manner, and the same or similar parts between the various embodiments can be referenced to each other, and each embodiment focuses on the differences from other embodiments.

[0213] The protection scope of the present application is not limited to the above-mentioned embodiments. Obviously, those skilled in the art can make various changes and modifications to the present disclosure without departing from the scope and spirit of the present disclosure. If these changes and modifications fall within the scope of the claims of the present disclosure and their equivalents, the intention of the present disclosure also includes these changes and modifications.

Claims

1. A method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions, characterized in that: include: According to the extreme wind conditions in the wind farm, sample points of wind characteristic parameters are obtained; Inputting the sample points of the wind characteristic parameters into a pre-established sensitivity analysis model to obtain the sensitivity of the wind characteristic parameters to the interlayer bonding defects at the main beam of the blade, wherein the pre-established sensitivity analysis model is established by using a BP neural network and a Sobol method, and the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is used as training data; According to the wind characteristic parameters and the sensitivity, a first-order Sobol effect index and a total effect Sobol index are calculated, so that the first N wind characteristic parameters with the highest first-order Sobol effect index and / or total effect Sobol index are used as highly sensitive factors; The actual measurement range of the highly sensitive factor is input into the pre-established cohesive force finite element model of the interlayer bonding defect of the blade to obtain the strain complementary energy density of the damaged unit in the delaminated damage area. The strain complementary energy density is used as the damage degree of the blade, and the ratio between the damage degree of the blade and the design load of the defect-free blade is calculated. If the ratio is greater than the preset threshold, the prediction result is that initial damage occurs at the corresponding position of the blade.

2. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 1 is characterized in that: The method of obtaining the sample points of the wind characteristic parameters according to the extreme wind conditions in the wind farm includes: Determine wind characteristic parameters of wind turbine blades and parameter ranges of wind characteristic parameters according to the extreme wind conditions in the wind farm, wherein the wind characteristic parameters include wind speed, wind direction, wind shear, turbulence intensity and wind frequency; Performing orthogonal design on the wind characteristic parameters and parameter intervals of the wind characteristic parameters, and selecting representative parameter combinations and corresponding parameter intervals; According to the representative parameter combination and the corresponding parameter interval, the Sobol sequence sampling method is used to obtain the sample points of the wind characteristic parameters.

3. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 1 is characterized in that: The training data acquisition process includes: According to the sample points of the wind characteristic parameters, under the action of different wind characteristic parameter combinations, the time domain load data of each airfoil center of the blade is calculated; Converting the time domain load data into the distribution of equivalent stress at different cross-sectional positions of a full-size blade; The equivalent stress is input as an external load into the pre-established cohesive finite element model of the interlaminar bonding defect of the blade to perform delamination damage evolution and obtain the stress and strain response data of delamination damage. According to the stress and strain response data of the delamination damage, the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is calculated, and the KL divergence is used as training data for the pre-established sensitivity analysis model.

4. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 3 is characterized in that: The method of calculating the time domain load data of each airfoil center of the blade based on the sample points of the wind characteristic parameters under the action of different wind characteristic parameter combinations includes: sequentially importing the sample points of the wind characteristic parameters into the load calculation software to establish a wind model, inputting the design parameters of the wind turbine and the actual extreme service wind conditions of the blades, and calculating the time domain load data of the blade flapping and swing directions at the center of each airfoil.

5. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 3 is characterized in that: The step of converting the time domain load data into the distribution of equivalent stress at different cross-sectional positions of a full-size blade includes: According to the time domain load data of each airfoil center of the blade, the normal stress and shear stress of the blade section are calculated; Among them, the calculation formula of the normal stress σ0 of the blade section is as follows: Among them, M x is the aerodynamic bending moment in the blade flapping direction, M y is the aerodynamic bending moment in the swing direction, W n is the bending section coefficient; The calculation formula of blade section shear stress τ0 is as follows: Among them, F x is the shear force in the blade swinging direction, F y is the shear force in the swing direction, A is the cross-sectional area of ​​the blade, According to the normal stress of the blade section and the shear stress of the blade section, the equivalent stress at different cross-sectional positions of the full-size blade is calculated, and the calculation formula is as follows: Where σ is the equivalent stress at different cross-sectional positions of the full-size blade; According to the equivalent stress at different cross-sectional positions of the full-size blade, an equivalent stress load curve of the full-size blade is drawn, wherein the position of the highest stress area in the full-size blade equivalent stress load curve is the maximum equivalent stress of the blade; The ultimate strength of the blade is determined according to the maximum equivalent stress of the blade.

6. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 3 is characterized in that: The equivalent stress is input as an external load into the pre-established cohesive force finite element model of the interlaminar bonding defect position of the blade to perform delamination damage evolution and obtain stress and strain response data of delamination damage, including: The equivalent stress is input as an external load into the pre-established cohesive force finite element model of the interlayer bonding defect of the blade to perform delamination damage evolution, wherein the interface unit introduces an energy-based bilinear cohesive force model to describe the mechanical behavior of the evolution of the interlayer bonding defect, and the association between the stress and relative displacement of the damage unit is established to describe the evolution process of the delamination damage; In the evolution of delamination damage, the stress and strain response data of delamination damage are collected, and the failure proportion of damaged units at different times is calculated. The time of initial delamination damage, damage area and stress distribution are determined according to the proportion of failed units, and the stress and displacement difference curve is plotted. The stress and displacement difference curve is integrated to obtain the strain energy density of the damaged failure unit.

7. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 3 is characterized in that: The step of calculating the KL divergence of the strain complementary energy density at the maximum damage location and the non-maximum damage location according to the stress and strain response data of the layer damage includes: The strain complementary energy density at the maximum damage is calculated based on the energy release rate at the maximum damage and the critical fracture energy at the maximum damage. The strain complementary energy density at the non-maximum damage is calculated based on the energy release rate at the non-maximum damage and the critical fracture energy at the non-maximum damage. The probability distribution function of strain complementary energy at the maximum damage is obtained according to the strain complementary energy density at the maximum damage, and the probability distribution function of strain complementary energy at the non-maximum damage is obtained according to the strain complementary energy density at the non-maximum damage; According to the probability distribution function of strain complementary energy at the maximum damage point and the probability distribution function of strain complementary energy at the non-maximum damage point, the KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is calculated.

8. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 1 is characterized in that: The pre-established sensitivity analysis model is established using BP neural network and Sobol method, including: Determine the BP neural network topology structure, including: determine the number of neurons in the input layer, the number of neurons in the output layer, the number of hidden layers and the number of neurons in each hidden layer, wherein the number of hidden layers is not equal to the number of neurons in each hidden layer; The Sobol sequence method is used to randomly extract sample points of wind characteristic parameters as the input of the BP neural network. The KL divergence of the strain complementary energy density at the maximum damage point and the non-maximum damage point is used as the output of the BP neural network for fitting, and multiple sensitivity analysis models are obtained. By comparing the output error statistical parameters of multiple sensitivity analysis models, the one with the smallest error is selected as the pre-established sensitivity analysis model.

9. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 8, characterized in that: The first-order Sobol effect index is calculated as follows: Among them, S i is the first-order Sobol effect index, which indicates the wind characteristic parameter X i The influence of a single action on the initial damage evolution of bonding defects, D Xi (E X～i (Y|X i )) is the wind characteristic parameter X i The variance caused by, D(Y) is the variance of sensitivity Y, E{·} is the mathematical expectation, X ～i to consider all random values ​​of variable i while keeping variable i fixed.

10. The method for predicting initial damage of interlayer bonding defects of wind turbine blades under extreme wind conditions according to claim 8, characterized in that: The total effect Sobol index is calculated as follows: Among them, S Ti is the total effect Sobol index, E X～i (D Xi (Y|X ～i )) is to determine the exclusion of X ～i The mean variance of sensitivity Y assuming the true values ​​of all variables except .

Images