Fatigue reliability self-adaptive evaluation method for wind power gear box under real-time wind condition

Abstract

The invention relates to the field of reliability design analysis of a transmission system, and discloses a fatigue reliability adaptive evaluation method for a wind power gear box under a real-time wind condition, and the method comprises the steps: obtaining a working condition curve of a wind turbine generator according to a wind energy utilization coefficient Cp; according to the real-time wind speed history, fitting to obtain average wind speed distribution and turbulent flow standard deviation characteristic values of the real-time wind speed history; according to the cut-in wind speed and the cut-out wind speed, generating an input rotating speed and an input torque history according to the working condition curve of the wind turbine generator, and obtaining a fatigue stress history of each failure mode according to external excitation; according to the fatigue stress distribution and the stress cycle index of the part under each working condition, obtaining the fatigue stress distribution of the part under the full wind condition; generating a corresponding equivalent fatigue load spectrum according to the fatigue stress distribution and the stress cycle index of the part under each working condition; obtaining the fatigue damage and residual strength history of the fatigue sample of each part; obtaining a function value of each failure mode of each fatigue sample; and obtaining the fatigue reliability of the wind power gear box according to the reliability of each failure mode.

Description

Technical Field

[0001] The present invention relates to the field of transmission system reliability design and analysis, and in particular to a fatigue reliability adaptive evaluation method for a wind power gearbox under real-time wind conditions. Background Art

[0002] Wind turbine gearbox is an important part of the wind turbine transmission chain, responsible for converting the low-speed rotation of the wind rotor into the high-speed rotation required by the generator. Gearboxes operate in high altitude and complex wind conditions, and are prone to fatigue failure due to variable loads. With the large-scale and scaled-up wind turbines, gearboxes need to have higher stability and reliability. The structural fatigue reliability assessment of wind turbine gearboxes is of great significance for improving structural design and maintaining the safe operation of the system.

[0003] The fatigue failure of wind turbine gearbox transmission system is mainly concentrated on gears and bearings. The failure mode of gears is pitting corrosion on tooth surface and fracture on tooth root, and the failure mode of bearings is subsurface fatigue pitting corrosion on inner ring, outer ring and rolling element. The fatigue performance of gears and bearings is usually affected by many factors such as material properties, operating load, geometric dimensions, etc. However, due to the instability of processing technology, random distribution of internal defects of materials and fluctuation of operating conditions, the above parameters show uncertainty, resulting in the dispersion of fatigue life of parts.

[0004] Existing structural fatigue reliability analysis methods usually only consider fatigue strength and fatigue stress as independent static variables, and solve static reliability through numerical analysis or simulation. However, under actual working conditions, various parts of wind turbine gearboxes bear the load source together, resulting in fatigue damage and common cause failure. Static reliability lacks the prediction of gearbox performance status and fatigue failure under variable working conditions. Summary of the invention

[0005] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a method for adaptively evaluating the fatigue reliability of a wind turbine gearbox under real-time wind conditions, comprising the following steps:

[0006] Step 1: Obtain the wind energy utilization coefficient C of the wind turbine blades p , and according to the wind energy utilization coefficient C p Obtain the operating curve of the wind turbine;

[0007] Step 2: According to the set time period, the average wind speed per unit time and the turbulence standard deviation of the real-time wind speed history are obtained, and the average wind speed distribution and turbulence standard deviation characteristic values ​​of the real-time wind speed history are obtained by fitting;

[0008] Step 3: extract the main working condition according to the cut-in wind speed and the cut-out wind speed with the set average wind speed increment value, the main working condition includes the extracted average wind speed, the corresponding turbulence standard deviation characteristic value and the instantaneous wind speed history, and generate the corresponding input speed and input torque history according to the working condition curve of the wind turbine generator set, the input speed and input torque history are the external excitation of the main wind condition;

[0009] Step 4: Input the external excitation of the main wind conditions into the gearbox dynamics model, use the Newmark-β time domain integration method to obtain the dynamic load history of the parts, and simultaneously calculate the gear tooth surface contact stress, tooth root bending stress and bearing subsurface shear stress to obtain the fatigue stress history of each failure mode;

[0010] Step 5, obtaining the fatigue stress and stress cycle number of the parts under the main wind conditions, and using the generalized extreme value distribution to fit the fatigue stress distribution of the parts under the main wind conditions, and using the piecewise function to fit the stress cycle number function of the parts to obtain the fatigue stress distribution of the parts under all wind conditions;

[0011] Step 6: Establish the PSN fatigue curve expression of the part, obtain the statistical data of failure probability and life at different stress levels, and fit the statistical data to obtain the PSN fatigue curves of different probabilities of the part; use Latin hypercube sampling to generate equal probability intervals from 0% to 100%. value, Substitute into the PSN fatigue curve formula to generate the PSN fatigue curve family;

[0012] Step 7, calculating the average wind speed as the current wind condition with unit time as the time period, interpolating the fatigue stress distribution and stress cycle number of the current wind condition, and generating the corresponding equivalent fatigue load spectrum;

[0013] Step 8: Monte Carlo simulation is performed on fatigue samples of each component, a load spectrum of the current wind condition is established per unit time, and the Pavlou nonlinear fatigue damage accumulation model is used to solve the fatigue damage and residual strength history of fatigue samples of each component;

[0014] Step nine, using adaptive bandwidth kernel density estimation method to fit the residual strength distribution of parts in real time; using the fitted residual strength distribution and fatigue stress distribution of parts to obtain the fatigue reliability of each failure mode of each part; according to the residual strength and fatigue stress of each fatigue sample of each part obtained by Latin hypercube sampling, calculate the functional function value of each failure mode of each fatigue sample;

[0015] Step 10: Based on the hierarchical correlation of system failures, a nested Copula structure is used to describe the failure correlation between various failure modes of the wind turbine gearbox, and a correlation network of each failure mode is established;

[0016] Step 11, determining the marginal distribution function value according to the obtained functional function of each component fatigue sample, and then using the Copula function to couple the marginal distribution to generate the joint distribution of the reliability of each failure mode, and obtain the reliability of the common cause failure;

[0017] Step 12: Use the AIC information criterion to measure the goodness of fit of the Copula function, select the Copula function with the highest goodness of fit, input the reliability of each failure mode into the correlation network, select the optimal Copula function layer by layer and connect them to obtain the fatigue reliability of the wind turbine gearbox at the top level.

[0018] Furthermore, the wind energy utilization coefficient C of the wind turbine blades is obtained. p , and according to the wind energy utilization coefficient C p The operating curve of the wind turbine is obtained, including:

[0019] Wind energy utilization coefficient C p for:

[0020]

[0021] The optimal tip speed ratio λ opt The input speed curve is obtained:

[0022]

[0023] The optimal tip speed ratio λ opt and the maximum wind energy utilization coefficient C pmax The input torque curve is obtained:

[0024]

[0025] Where β is the pitch angle, λ is the tip speed ratio, v is the incoming wind speed, ω is the rotor angular velocity, R is the blade radius; A, B, C, D, E, F, G, H are constants; n in is the input speed, n low is the minimum speed, n rate is the rated speed, v is the wind speed, v cut-in is the cut-in wind speed, v cut-off is the cut-out wind speed, v nrate T is the wind speed corresponding to the rated speed of the wind turbine; in is the input torque, T rate is the rated torque, v is the wind speed, v rate is the rated wind speed.

[0026] Furthermore, the method of obtaining the average wind speed per unit time and the turbulence standard deviation of the real-time wind speed history according to the set time period, and fitting the average wind speed distribution and the turbulence standard deviation characteristic value of the real-time wind speed history includes:

[0027] Probability density function of the Weibull distribution of mean wind speed:

[0028]

[0029] Turbulence standard deviation function of instantaneous wind speed:

[0030] σ1(V)=I 15 (0.75V+5.6)

[0031] Where V is the average wind speed, k is the shape parameter of the Weibull distribution, c is the scale parameter of the Weibull distribution, and I 15 It is the characteristic value of turbulence intensity at an average wind speed of 15m / s.

[0032] Furthermore, the PSN fatigue curve expression of the part is established to obtain the statistical data of failure probability and life at different stress levels, and the PSN fatigue curves of the part with different probabilities are obtained by fitting the statistical data; Latin hypercube sampling is used to generate equal probability intervals from 0% to 100%. value, Substitute the PSN fatigue curve formula to generate the PSN fatigue curve family, including:

[0033] The PSN fatigue curve of the gear tooth surface contact is:

[0034]

[0035] The PSN fatigue curve of the gear tooth root bending is:

[0036]

[0037] The PSN fatigue curve of the bearing subsurface shear stress is:

[0038]

[0039] Where N H is the number of stress cycles to contact failure, σ H is the contact stress, C H is the contact fatigue stress constant, A H is the contact stress life coefficient, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, Z R is the roughness coefficient, Z v is the speed coefficient, Z L is the lubricant coefficient, Z W is the working hardening coefficient of the tooth surface, Z X is the contact strength size factor, N F is the number of stress cycles to bending failure, σ Fis the bending stress, C F is the bending fatigue stress constant, A F is the bending stress life factor, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, Y ST is the stress correction factor, Y δrelT is the relative tooth root fillet sensitivity coefficient, Y RrelT is the relative tooth root surface condition coefficient, Y X is the bending strength size factor, N o is the number of cycles of bearing subsurface spalling, τ o is the subsurface shear stress, C o is the subsurface fatigue stress constant, A o is the subsurface shear stress life coefficient, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, K τ is the effective stress concentration factor, ε τ is the size coefficient, and β is the surface processing coefficient.

[0040] Furthermore, the Monte Carlo simulation is performed on the fatigue samples of each component, a load spectrum of the current wind condition is established per unit time, and the Pavlou nonlinear fatigue damage accumulation model is used to solve the fatigue damage and residual strength history of the fatigue samples of each component, including:

[0041] Fatigue damage accumulation D of parts k The function is as follows:

[0042]

[0043] The residual strength of the part material under multi-level loading is as follows:

[0044] r=r0-(r0-σ k )D k

[0045] Where D k is the k-level fatigue damage accumulation, σ k is the kth level fatigue stress amplitude, n k is the number of stress cycles at level k, N k is the number of failure cycles at level k, q(σ)=(σσ u ) -0.75 is the nonlinear damage index, σ u is the ultimate bearing stress of the part, r is the residual strength of the part, r0=σ u is the initial static strength of the part.

[0046] Furthermore, the adaptive bandwidth kernel density estimation method is used to fit the residual strength distribution of parts in real time; the fatigue reliability of each failure mode of each part is obtained by using the fitted residual strength distribution and fatigue stress distribution of the parts; and the functional function value of each failure mode of each fatigue sample is calculated according to the residual strength and fatigue stress of each fatigue sample of each part obtained by Latin hypercube sampling, including:

[0047] The mean of the equivalent normalized variable is:

[0048]

[0049] The standard deviation of the equivalent normalized variable is:

[0050]

[0051] Where f(x * ) is the probability density function of the design point, is the standard normal probability density function.

[0052] The reliability index and dynamic reliability are as follows:

[0053]

[0054] R(t)=Φ[β(t)]

[0055] where x * is the design point, F(x * ) is the cumulative distribution function of the design point, Φ- 1 is the standard normal inverse cumulative distribution function,

[0056] μ g is the mean value of the performance function, s g is the standard deviation of the performance function, μ rn is the normalized residual intensity mean, μ σn is the normalized fatigue stress mean, s rn is the normalized residual intensity standard deviation, s σn is the normalized fatigue stress standard deviation.

[0057] Furthermore, the AIC information criterion is used to measure the goodness of fit of the Copula function, the Copula function with the highest goodness of fit is selected, the reliability of each failure mode is input into the correlation network, the optimal Copula function is selected layer by layer and connected, and the fatigue reliability of the wind turbine gearbox at the top layer is obtained, including:

[0058] The expression of AIC information criterion is:

[0059] AIC=-2lnL(θ)+2K

[0060] Where L(θ) is the maximum likelihood function value of the Copula function, K is the number of parameters, the reliability of each failure mode is input into the correlation network, the optimal Copula function is selected layer by layer and connected, and finally the fatigue reliability of the wind turbine gearbox at the top layer is obtained.

[0061] The beneficial effect of the present invention is that through the technical solution provided by the present invention, fatigue damage and common cause failure caused by the common load source borne by various parts of the wind turbine gearbox can be achieved, and the performance status and fatigue failure of the gearbox under variable working conditions can be predicted. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a flow chart of a method for adaptively evaluating fatigue reliability of wind turbine gearboxes under real-time wind conditions;

[0063] Figure 2 is an overall flow chart of an embodiment of the present invention;

[0064] Figure 3 A surface diagram of wind energy utilization coefficient of the wind turbine generator system according to the present invention;

[0065] Figure 4 is an operating condition curve diagram of the wind turbine generator set according to the present invention;

[0066] Figure 5 is an operating condition curve diagram of the wind turbine generator set according to the present invention;

[0067] Figure 6 It is a statistical diagram of fatigue stress of gears and bearings according to the present invention;

[0068] Figure 7 A fatigue stress distribution estimation diagram of the parts described in the present invention;

[0069] Figure 8 A sample diagram of the PSN fatigue curve of the parts described in the present invention;

[0070] Fig. 9 A residual strength distribution estimation diagram of the part according to the present invention;

[0071] Fig.10 It is a fatigue reliability curve diagram of failure modes of various parts under real-time wind conditions described in the present invention;

[0072] Fig.11 A failure mode correlation network diagram of the wind power gearbox according to the present invention;

[0073] Fig.12 This is a fatigue reliability curve diagram of the wind turbine gearbox under real-time wind conditions described in the present invention. DETAILED DESCRIPTION

[0074] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings, but the protection scope of the present invention is not limited to the following.

[0075] The features and performance of the present invention are further described in detail below in conjunction with the embodiments.

[0076] like Figure 1 As shown, a fatigue reliability adaptive evaluation method for a wind turbine gearbox under real-time wind conditions comprises the following steps:

[0077] Step 1: Obtain the wind energy utilization coefficient C of the wind turbine blades p , and according to the wind energy utilization coefficient C p Obtain the operating curve of the wind turbine;

[0078] Step 2: According to the set time period, the average wind speed per unit time and the turbulence standard deviation of the real-time wind speed history are obtained, and the average wind speed distribution and turbulence standard deviation characteristic values ​​of the real-time wind speed history are obtained by fitting;

[0079] Step 3: extract the main working condition according to the cut-in wind speed and the cut-out wind speed with the set average wind speed increment value, the main working condition includes the extracted average wind speed, the corresponding turbulence standard deviation characteristic value and the instantaneous wind speed history, and generate the corresponding input speed and input torque history according to the working condition curve of the wind turbine generator set, the input speed and input torque history are the external excitation of the main wind condition;

[0080] Step 4: Input the external excitation of the main wind conditions into the gearbox dynamics model, use the Newmark-β time domain integration method to obtain the dynamic load history of the parts, and simultaneously calculate the gear tooth surface contact stress, tooth root bending stress and bearing subsurface shear stress to obtain the fatigue stress history of each failure mode;

[0081] Step 5, obtaining the fatigue stress and stress cycle number of the parts under the main wind conditions, and using the generalized extreme value distribution to fit the fatigue stress distribution of the parts under the main wind conditions, and using the piecewise function to fit the stress cycle number function of the parts to obtain the fatigue stress distribution of the parts under all wind conditions;

[0082] Step 6: Establish the PSN fatigue curve expression of the part, obtain the statistical data of failure probability and life at different stress levels, and fit the statistical data to obtain the PSN fatigue curves of different probabilities of the part; use Latin hypercube sampling to generate equal probability intervals from 0% to 100%. value, Substitute into the PSN fatigue curve formula to generate the PSN fatigue curve family;

[0083] Step 7, calculating the average wind speed as the current wind condition with unit time as the time period, interpolating the fatigue stress distribution and stress cycle number of the current wind condition, and generating the corresponding equivalent fatigue load spectrum;

[0084] Step 8: Monte Carlo simulation is performed on fatigue samples of each component, a load spectrum of the current wind condition is established per unit time, and the Pavlou nonlinear fatigue damage accumulation model is used to solve the fatigue damage and residual strength history of fatigue samples of each component;

[0085] Step nine, using adaptive bandwidth kernel density estimation method to fit the residual strength distribution of parts in real time; using the fitted residual strength distribution and fatigue stress distribution of parts to obtain the fatigue reliability of each failure mode of each part; according to the residual strength and fatigue stress of each fatigue sample of each part obtained by Latin hypercube sampling, calculate the functional function value of each failure mode of each fatigue sample;

[0086] Step 10: Based on the hierarchical correlation of system failures, a nested Copula structure is used to describe the failure correlation between various failure modes of the wind turbine gearbox, and a correlation network of each failure mode is established;

[0087] Step 11, determine the marginal distribution function value according to the obtained functional function of each part fatigue sample, and then use the Copula function to couple the marginal distribution to generate the joint distribution of the reliability of each failure mode, so as to obtain the reliability of the common cause failure;

[0088] Step 12: Use the AIC information criterion to measure the goodness of fit of the Copula function, select the Copula function with the highest goodness of fit, input the reliability of each failure mode into the correlation network, select the optimal Copula function layer by layer and connect them to obtain the fatigue reliability of the wind turbine gearbox at the top level.

[0089] The wind energy utilization coefficient C of the wind turbine blades is obtained. p , and according to the wind energy utilization coefficient C p The operating curve of the wind turbine is obtained, including:

[0090] Wind energy utilization coefficient C p for:

[0091]

[0092] The optimal tip speed ratio λ opt The input speed curve is obtained:

[0093]

[0094] The optimal tip speed ratio λ opt and the maximum wind energy utilization coefficient C pmax The input torque curve is obtained:

[0095]

[0096] Where β is the pitch angle, λ is the tip speed ratio, v is the incoming wind speed, ω is the rotor angular velocity, R is the blade radius; A, B, C, D, E, F, G, H are constants; n in is the input speed, n low is the minimum speed, n rate is the rated speed, v is the wind speed, v cut-in is the cut-in wind speed, v cut-off is the cut-out wind speed, v nrate T is the wind speed corresponding to the rated speed of the wind turbine; in is the input torque, T rate is the rated torque, v is the wind speed, v rate is the rated wind speed.

[0097] The method of obtaining the average wind speed per unit time and the turbulence standard deviation of the real-time wind speed history according to the set time period, and fitting the average wind speed distribution and the turbulence standard deviation characteristic value of the real-time wind speed history includes:

[0098] Probability density function of the Weibull distribution of mean wind speed:

[0099]

[0100] Turbulence standard deviation function of instantaneous wind speed:

[0101] σ1(V)=I 15 (0.75V+5.6)

[0102] Where V is the average wind speed, k is the shape parameter of the Weibull distribution, c is the scale parameter of the Weibull distribution, and I 15 It is the characteristic value of turbulence intensity at an average wind speed of 15m / s.

[0103] The PSN fatigue curve expression of the part is established, and the statistical data of failure probability and life at different stress levels are obtained. The PSN fatigue curves of different probabilities of the part are obtained by fitting the statistical data; Latin hypercube sampling is used to generate equal probability intervals from 0% to 100%. value, Substitute the PSN fatigue curve formula to generate the PSN fatigue curve family, including:

[0104] The PSN fatigue curve of the gear tooth surface contact is:

[0105]

[0106] The PSN fatigue curve of the gear tooth root bending is:

[0107]

[0108] The PSN fatigue curve of the bearing subsurface shear stress is:

[0109]

[0110] Where N H is the number of stress cycles to contact failure, σ H is the contact stress, C H is the contact fatigue stress constant, A H is the contact stress life coefficient, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, Z R is the roughness coefficient, Z v is the speed coefficient, Z L is the lubricant coefficient, Z W is the working hardening coefficient of the tooth surface, Z X is the contact strength size factor, N F is the number of stress cycles to bending failure, σ F is the bending stress, C F is the bending fatigue stress constant, A F is the bending stress life factor, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, Y ST is the stress correction factor, Y δrelT is the relative tooth root fillet sensitivity coefficient, Y RrelT is the relative tooth root surface condition coefficient, Y X is the bending strength size factor, N o is the number of cycles of bearing subsurface spalling, τ o is the subsurface shear stress, C o is the subsurface fatigue stress constant, A o is the subsurface shear stress life coefficient, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, K τ is the effective stress concentration factor, ε τ is the size coefficient and β is the surface processing coefficient.

[0111] The Monte Carlo simulation is performed on the fatigue samples of each component, a load spectrum of the current wind condition is established per unit time, and the Pavlou nonlinear fatigue damage accumulation model is used to solve the fatigue damage and residual strength history of the fatigue samples of each component, including:

[0112] Fatigue damage accumulation D of parts k The function is as follows:

[0113]

[0114] The residual strength of the part material under multi-level loading is as follows:

[0115] r=r0-(r0-σ k )D k

[0116] Where D k is the k-level fatigue damage accumulation, σ k is the kth level fatigue stress amplitude, n k is the number of stress cycles at level k, N k is the number of failure cycles at level k, q(σ)=(σσ u ) -0.75 is the nonlinear damage index, σ u is the ultimate bearing stress of the part, r is the residual strength of the part, r0=σ u is the initial static strength of the part.

[0117] The method adopts the adaptive bandwidth kernel density estimation method to fit the residual strength distribution of parts in real time; uses the fitted residual strength distribution and fatigue stress distribution of parts to obtain the fatigue reliability of each failure mode of each part; and calculates the functional function value of each failure mode of each fatigue sample based on the residual strength and fatigue stress of each fatigue sample obtained by Latin hypercube sampling, including:

[0118] The mean of the equivalent normalized variable is:

[0119] μ xn =x * -Φ -1 [F(x * )]s xn

[0120] The standard deviation of the equivalent normalized variable is:

[0121]

[0122] Where f(x * ) is the probability density function of the design point, is the standard normal probability density function.

[0123] The reliability index and dynamic reliability are as follows:

[0124]

[0125] R(t)=Φ[β(t)]

[0126] where x * is the design point, F(x * ) is the cumulative distribution function of the design point, Φ- 1 is the standard normal inverse cumulative distribution function,

[0127] μ g is the mean value of the performance function, s g is the standard deviation of the performance function, μ rn is the normalized residual intensity mean, μ σn is the normalized fatigue stress mean, s rn is the normalized residual intensity standard deviation, s σn is the normalized fatigue stress standard deviation.

[0128] The AIC information criterion is used to measure the goodness of fit of the Copula function, the Copula function with the highest goodness of fit is selected, the reliability of each failure mode is input into the correlation network, the optimal Copula function is selected layer by layer and connected, and the fatigue reliability of the wind turbine gearbox at the top layer is obtained, including:

[0129] The expression of AIC information criterion is:

[0130] AIC=-2lnL(θ)+2K

[0131] Where L(θ) is the maximum likelihood function value of the Copula function, K is the number of parameters, the reliability of each failure mode is input into the correlation network, the optimal Copula function is selected layer by layer and connected, and finally the fatigue reliability of the wind turbine gearbox at the top layer is obtained.

[0132] Specifically, the present invention aims to evaluate the fatigue reliability of wind turbine gearboxes in real time according to real-time wind conditions, and provide reference suggestions for health monitoring and maintenance of wind turbine gearboxes. Figure 2 Schematic diagram of the process, this embodiment is divided into three modules: fatigue stress solution, part fatigue reliability, and failure correlation to evaluate the fatigue reliability of the wind power gearbox in real time, specifically including the following steps:

[0133] Step 1: Obtain the wind energy utilization coefficient C of the wind turbine blades p for:

[0134]

[0135] Among them, β is the pitch angle, λ is the tip speed ratio, v is the incoming wind speed, ω is the angular velocity of the wind wheel, and R is the blade radius; A, B, C, D, E, F, G, and H are constants, and the parameters are provided by the blade manufacturer. Taking a 2MW wind turbine as an example, a large amount of data fitting results in A=0.5176, B=116, C=0.08, D=0.035, E=0.4, F=5, G=21, and H=0.0068. It can be seen from the wind energy utilization coefficient surface that when the wind speed is lower than the rated wind speed, the pitch angle β=0°, and the rotor speed ω is controlled to ensure that the tip speed ratio λ remains at the optimal tip speed ratio λ. opt , so that the wind turbine can achieve the maximum wind energy utilization coefficient Cpmax Continue to run, such as Figure 3 shown.

[0136] Step 2: Use wind energy utilization coefficient C p Calculate the operating curve of the wind turbine, including the input speed curve and input torque curve, such as Figure 4 As shown in (a) and (b).

[0137] The optimal tip speed ratio λ opt The input speed curve is obtained:

[0138]

[0139] where n in is the input speed, n low is the minimum speed, n rate is the rated speed, v is the wind speed, v cut-in is the cut-in wind speed, v cut-off To cut out the wind speed, It is the wind speed corresponding to when the wind turbine reaches the rated speed.

[0140] The optimal tip speed ratio λ opt and the maximum wind energy utilization coefficient C pmax The input torque curve is obtained:

[0141]

[0142] Where T in is the input torque, T rate is the rated torque, v is the wind speed, v cut-in is the cut-in wind speed, v rate is the rated wind speed, v cut-off To cut out wind speed.

[0143] Step 3: According to IEC 61400, the ten-minute average wind speed and turbulence standard deviation of the real-time wind speed history are counted in a ten-minute period, and the average wind speed distribution and turbulence standard deviation characteristic values ​​of the real-time wind speed history are fitted. The average wind speed distribution and turbulence standard deviation characteristic values ​​will change with the real-time wind speed history, thereby changing the fatigue stress distribution.

[0144] Probability density function of the Weibull distribution of mean wind speed:

[0145]

[0146] Where V is the average wind speed, k is the shape parameter of the Weibull distribution, and c is the scale parameter of the Weibull distribution.

[0147] Turbulence standard deviation function of instantaneous wind speed:

[0148] σ1(V)=I15 (0.75V+5.6) (6)

[0149] Where V is the average wind speed, I 15 It is the characteristic value of turbulence intensity at an average wind speed of 15m / s.

[0150] Step 4: Since the actual wind speed is complex and changeable, the main working conditions are extracted from the cut-in wind speed of 3m / s to the cut-out wind speed of 25m / s in increments of 1m / s for the simulation of wind speed history and external excitation. The main working conditions include the extracted average wind speed, the corresponding turbulence standard deviation characteristic value and the instantaneous wind speed history. The instantaneous wind speed needs to be generated using the Kaimal spectrum. Further, according to the working condition curves (3) and (4), the input speed and input torque history are generated as external excitations, such as Figure 5 As shown in (a) and (b).

[0151] The instantaneous wind speed generated by the Kaimal spectrum is:

[0152]

[0153] Where t is the time, V is the average wind speed, N is the number of samples, and T is the duration. is the resolution frequency, is a random phase.

[0154] The power spectral density S of the Kaimal spectrum is:

[0155]

[0156] Where f is the frequency, L is the turbulence scale, and V is the average wind speed.

[0157] Step 5: Wind turbine gearboxes mainly include five failure modes: pitting of tooth surface, tooth root fracture, and pitting of bearing inner ring, roller, and outer ring. The external excitation of the main wind conditions is input into the gearbox dynamics model, and the newmark-β time domain integration method is used to solve the dynamic load history of the parts. At the same time, the gear tooth surface contact stress, tooth root bending stress, and bearing subsurface shear stress are calculated to obtain the fatigue stress history of each failure mode.

[0158] The formula for tooth surface contact stress is as follows:

[0159]

[0160] where F t is the nominal tangential force, B is the tooth width, d is the pitch circle diameter of the pinion, u is the gear ratio, Z H is the node area coefficient, Z E is the elastic modulus, Z ε is the overlap coefficient, Z β is the helix angle coefficient, K Hβ and KHα are the tooth load distribution coefficient and inter-tooth load distribution coefficient of contact stress, K A is the service factor, K V is the dynamic load coefficient.

[0161] The formula for tooth root bending stress is as follows:

[0162]

[0163] Among them, F t is the nominal tangential force, B is the tooth width, m n is the normal modulus, Y Fα is the tooth form factor, Y Sα is the stress correction factor, Y ε is the overlap coefficient, Y β is the helix angle coefficient, K Fβ and K Fα are the tooth load distribution coefficient and inter-tooth load distribution coefficient of bending stress, K A is the service factor, K V is the dynamic load coefficient.

[0164] The formula for the bearing subsurface shear stress is as follows:

[0165]

[0166] Where y is the coordinate of the contact area raceway direction, l is the roller length, is the half width of the contact surface, is the curvature of the contact point and, is the raceway normal load, ψ is the azimuth angle, F r is the radial load, J r is the integral of the radial load.

[0167] Step 6: According to IEC 61400-4, the fatigue stress and stress cycle number of parts under the main wind conditions are counted. The fatigue stress of gears is calculated by using the load duration distribution method, which discretizes the continuous stress history into stress blocks of different ranges for statistics. The fatigue stress of bearings is obtained by statistically analyzing the peak value of stress pulsation, such as Figure 6 As shown in (a) and (b).

[0168] Step 7: Use generalized extreme value distribution to fit the fatigue stress distribution of parts under main wind conditions, use piecewise function to fit the stress cycle function of parts, and further solve the fatigue stress distribution of parts under all wind conditions, such as Figure 7 As shown in (a), (b), (c), (d), (e), and (f).

[0169] The probability density function of fatigue stress distribution of parts under main wind conditions is:

[0170]

[0171] where σ V is the fatigue stress of parts under various wind conditions, k is the shape parameter of the generalized extreme value distribution, c is the scale parameter of the generalized extreme value distribution, and μ is the location parameter of the generalized extreme value distribution.

[0172] The stress cycle function of the part is:

[0173]

[0174] Where V is the average wind speed, V cutin is the cut-in wind speed, V cutoff To cut out the wind speed, is the average wind speed corresponding to the maximum number of cycles, n m i n is the minimum number of cycles, n max is the maximum number of cycles.

[0175] The probability density function of fatigue stress distribution of parts under full wind conditions is:

[0176]

[0177] Where σ is the fatigue stress of the part, is the probability density function of fatigue stress generated by parts under various wind conditions, f(V) is the probability density function of average wind speed distribution, n(V) is the stress cycle number function of parts, f(σ V ) is the probability density function of fatigue stress distribution of parts under various wind conditions.

[0178] Step 8: Refer to the Basquin formula to establish the PSN fatigue curve expression of the parts. Taking the 2MW wind power gearbox as an example, the inner ring gear uses 42CrMoA tempered surface hardened soft tooth surface gear, the cylindrical outer gears are all 18CrNiMo7-6 carburized hardened hard tooth surface gears, and the bearings are made of GCr15 bearing steel. The group test method is used to perform fatigue tests on standard gear and bearing samples to obtain statistical data on the failure probability and life of several groups of different stress levels, and the data is fitted to obtain the PSN fatigue curves of different probabilities of parts.

[0179] The PSN fatigue curve of gear tooth contact is:

[0180]

[0181] Where N H is the number of stress cycles to contact failure, σ H is the contact stress, C H is the contact fatigue stress constant, A His the contact stress life coefficient, E δ is the coefficient of variation of logarithmic life span, is the standard normal deviation of the failure probability. The above parameters are related to the gear material and can be obtained based on the group test of standard gears. R is the roughness coefficient, Z v is the speed coefficient, Z L is the lubricant coefficient, Z W is the working hardening coefficient of the tooth surface, Z X is the contact strength size coefficient. The above fatigue strength influence coefficient can be calculated according to GB / T 3480.

[0182] The PSN fatigue curve of the gear tooth root bending is:

[0183]

[0184] Where N F is the number of stress cycles to bending failure, σ F is the bending stress, C F is the bending fatigue stress constant, A F is the bending stress life factor, E δ is the coefficient of variation of logarithmic life span, is the standard normal deviation of the failure probability. The above parameters are related to the gear material and can be obtained based on the group test of standard gears. ST is the stress correction factor, Y δrelT is the relative tooth root fillet sensitivity coefficient, Y RrelT is the relative tooth root surface condition coefficient, Y X is the bending strength size coefficient. The above fatigue strength influence coefficient can be calculated according to GB / T 3480.

[0185] The PSN fatigue curve of the bearing subsurface shear stress is:

[0186]

[0187] Where N o is the number of cycles of bearing subsurface spalling, τ o is the subsurface shear stress, C o is the subsurface fatigue stress constant, A o is the subsurface shear stress life coefficient, E δ is the coefficient of variation of logarithmic life span, is the standard normal deviation of the failure probability. The above parameters are related to the bearing material and can be obtained based on group tests of standard bearing steel. τ is the effective stress concentration factor, ε τ is the size coefficient, β is the surface processing coefficient, and the above fatigue strength influence coefficients can be obtained according to the mechanical design manual.

[0188] Step 9: Use Latin hypercube sampling to generate equally probable intervals from 0% to 100%. value, Substitute the PSN fatigue curve formula to generate the PSN fatigue curve family, such as Figure 8 (a), (b), and (c) are shown. Latin hypercube sampling is a uniform stratified sampling in multidimensional space, and minimizing the correlation between each part improves the sampling efficiency. The fatigue samples of each part generated can better represent the overall characteristics of the sample, and Monte Carlo simulation can be performed with a smaller amount of calculation.

[0189] Step 10: According to IEC 61400, the average wind speed is calculated as the current wind condition with a time period of ten minutes, and the fatigue stress distribution and stress cycle number of the current wind condition are interpolated to generate the corresponding equivalent fatigue load spectrum. The equivalent stress expression is:

[0190]

[0191] where σ V is the fatigue stress of the component under the current wind condition, A is the stress life coefficient, which is obtained according to the PSN fatigue curve family of the component, and f(σ V ) is the probability density function of fatigue stress distribution under current wind conditions.

[0192] Step 11: Perform Monte Carlo simulation on the fatigue samples of each part, establish the load spectrum (including equivalent stress and number of stress cycles) of the current wind conditions every ten minutes, and use the Pavlou nonlinear fatigue damage accumulation model to solve the fatigue damage and residual strength history of the fatigue samples of each part.

[0193] Fatigue damage accumulation D of parts k The function is as follows:

[0194]

[0195] Where D k is the k-level fatigue damage accumulation, σ k is the kth level fatigue stress amplitude, n k is the number of stress cycles at level k, N k is the number of failure cycles at level k, q(σ)=(σσ u ) -0.75 is the nonlinear damage index, σ u is the ultimate bearing stress of the part.

[0196] The residual strength of the part material under multi-level loading is as follows:

[0197] r=r0-(r0-σ k )D k (20)

[0198] Where r is the residual strength of the part, r0 = σ u is the initial static strength of the part.

[0199] Step 12: Since the residual strength distribution of the part is a typical long-tail distribution, the adaptive bandwidth kernel density estimation method is used to fit the residual strength distribution of the part in real time, such as Fig. 9 As shown. This method can adaptively correct the bandwidth parameters to fit the residual strength distribution in real time according to the residual strength statistics of the fatigue samples of the parts, thereby improving the fitting efficiency and accuracy. The probability density function of the residual strength distribution of the parts is:

[0200]

[0201] Where S is the number of fatigue samples of each part, r i is the residual strength of fatigue samples of each part, is the Gaussian kernel function for kernel density estimation, To modify the adaptive bandwidth obtained by fixing the fixed bandwidth h, f(r i ) is the probability density estimate of a fixed bandwidth.

[0202] Step 13: Use the fitted residual strength distribution and fatigue stress distribution of the parts to solve the fatigue reliability of each failure mode of the parts, including five failure modes: pitting of the tooth surface, tooth root fracture, and pitting of the inner ring, roller, and outer ring of the bearing. The fatigue stress distribution and residual strength distribution of the parts during the fatigue process are non-normal distributions, so the JC method is used to perform equivalent normalization on the non-normal variables, search for the point with the maximum failure probability on the failure boundary as the design point, and solve the reliability index and reliability through tail approximate normalization. The fatigue reliability of each failure mode of the parts under real-time wind conditions is obtained as follows: Fig.10 shown.

[0203] The mean of the equivalent normalized variable is:

[0204]

[0205] where x * is the design point, F(x * ) is the cumulative distribution function of the design point, Φ -1 is the standard inverse normal cumulative distribution function.

[0206] The standard deviation of the equivalent normalized variable is:

[0207]

[0208] Where f(x * ) is the probability density function of the design point, is the standard normal probability density function.

[0209] The reliability index and dynamic reliability are as follows:

[0210]

[0211] where μ g is the mean value of the performance function, s g is the standard deviation of the performance function, μ rn is the normalized residual intensity mean, μ σn is the normalized fatigue stress mean, s rn is the normalized residual intensity standard deviation, s σn is the normalized fatigue stress standard deviation.

[0212] Step 14: Based on the residual strength and fatigue stress of each fatigue sample of each part obtained by Latin hypercube sampling, calculate the functional function value of each failure mode of each fatigue sample. The functional function of each failure mode of the part can be expressed as

[0213] g(t)=r(t)-σ(t) (26)

[0214] Step 15: Based on the hierarchical association of system failures, a nested Copula structure is used to describe the failure correlation between various failure modes of the wind turbine gearbox, and a correlation network of each failure mode is established, such as Fig.11 As shown. The reliability of a single failure mode of a part in the figure is R n (n=H,F,o,m,i), respectively representing tooth surface contact, tooth root bending and pitting of the inner ring, roller and outer ring of the bearing. The reliability of common cause failure of gear is C g , the reliability of common cause failure of bearing is C b , the reliability of common cause failure of parts is C i (i=1,2···14), where the top layer C 14 is the system reliability of the gearbox.

[0215] Step 16: Determine the marginal distribution function value based on the obtained functional function of each part fatigue sample, and then use the Copula function to couple the marginal distribution to generate the joint distribution of the reliability of each failure mode. The reliability of common cause failure can be obtained as:

[0216] C 12 =P(g1>0I g2>0)=C(R1,R2) (27)

[0217] Wherein g1 and g2 are the functional functions of a single failure mode, obtained in step 14; R1 and R2 are the reliability of a single failure mode, obtained in step 13.

[0218] Commonly used Copula functions include Gumbel, Clayton and Frank Archimedean Copula functions. These three functions are used to describe the failure correlation of functional functions. The expressions are as follows:

[0219]

[0220] Among them, u1 and u2 are two related variables, and θ is the correlation parameter.

[0221] Step 17: Use the AIC information criterion to measure the goodness of fit of the Copula function. The smaller the AIC value, the higher the goodness of fit of the Copula function. Select the Copula function with the highest goodness of fit from the three Copula functions. The expression of the AIC information criterion is:

[0222] AIC=-2lnL(θ)+2K (31)

[0223] Where L(θ) is the maximum likelihood function value of the Copula function, and K is the number of parameters.

[0224] The reliability of each failure mode is input into the correlation network, and the optimal Copula function is selected layer by layer and connected, and finally the fatigue reliability of the wind turbine gearbox at the top layer is obtained, such as Fig.12 shown.

[0225] The above is only a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the form disclosed herein, and should not be regarded as excluding other embodiments, but can be used in various other combinations, modifications and environments, and can be modified within the scope of the concept described herein through the above teachings or the technology or knowledge of the relevant field. The changes and modifications made by those skilled in the art shall not deviate from the spirit and scope of the present invention, and shall be within the scope of protection of the claims attached to the present invention.

Claims

1. A method for adaptively evaluating fatigue reliability of wind turbine gearboxes under real-time wind conditions, characterized in that: The steps include: Step 1: Obtain the wind energy utilization coefficient C of the wind turbine blades p , and according to the wind energy utilization coefficient C p Obtain the operating curve of the wind turbine; Step 2: According to the set time period, the average wind speed per unit time and the turbulence standard deviation of the real-time wind speed history are obtained, and the average wind speed distribution and turbulence standard deviation characteristic values ​​of the real-time wind speed history are obtained by fitting; Step 3: extract the main working condition according to the cut-in wind speed and the cut-out wind speed with the set average wind speed increment value, the main working condition includes the extracted average wind speed, the corresponding turbulence standard deviation characteristic value and the instantaneous wind speed history, and generate the corresponding input speed and input torque history according to the working condition curve of the wind turbine generator set, the input speed and input torque history are the external excitation of the main wind condition; Step 4: Input the external excitation of the main wind conditions into the gearbox dynamics model, use the Newmark-β time domain integration method to obtain the dynamic load history of the parts, and simultaneously calculate the gear tooth surface contact stress, tooth root bending stress and bearing subsurface shear stress to obtain the fatigue stress history of each failure mode; Step 5, obtaining the fatigue stress and stress cycle number of the parts under the main wind conditions, and using the generalized extreme value distribution to fit the fatigue stress distribution of the parts under the main wind conditions, and using the piecewise function to fit the stress cycle number function of the parts to obtain the fatigue stress distribution of the parts under all wind conditions; Step 6: Establish the PSN fatigue curve expression of the part, obtain the statistical data of failure probability and life at different stress levels, and fit the statistical data to obtain the PSN fatigue curves of different probabilities of the part; use Latin hypercube sampling to generate equal probability intervals from 0% to 100%. value, Substitute into the PSN fatigue curve formula to generate the PSN fatigue curve family; Step 7, calculating the average wind speed as the current wind condition with unit time as the time period, interpolating the fatigue stress distribution and stress cycle number of the current wind condition, and generating the corresponding equivalent fatigue load spectrum; Step 8: Monte Carlo simulation is performed on fatigue samples of each component, a load spectrum of the current wind condition is established per unit time, and the Pavlou nonlinear fatigue damage accumulation model is used to solve the fatigue damage and residual strength history of fatigue samples of each component; Step nine, using adaptive bandwidth kernel density estimation method to fit the residual strength distribution of parts in real time; using the fitted residual strength distribution and fatigue stress distribution of parts to obtain the fatigue reliability of each failure mode of each part; according to the residual strength and fatigue stress of each fatigue sample of each part obtained by Latin hypercube sampling, calculate the functional function value of each failure mode of each fatigue sample; Step 10: Based on the hierarchical correlation of system failures, a nested Copula structure is used to describe the failure correlation between various failure modes of the wind turbine gearbox, and a correlation network of each failure mode is established; Step 11, determining the marginal distribution function value according to the obtained functional function of each component fatigue sample, and then using the Copula function to couple the marginal distribution to generate the joint distribution of the reliability of each failure mode, and obtain the reliability of the common cause failure; Step 12: Use the AIC information criterion to measure the goodness of fit of the Copula function, select the Copula function with the highest goodness of fit, input the reliability of each failure mode into the correlation network, select the optimal Copula function layer by layer and connect them to obtain the fatigue reliability of the wind turbine gearbox at the top level.

2. The method for adaptively evaluating fatigue reliability of a wind turbine gearbox under real-time wind conditions according to claim 1 is characterized in that: The wind energy utilization coefficient C of the wind turbine blades is obtained. p , and according to the wind energy utilization coefficient C p The operating curve of the wind turbine is obtained, including: Wind energy utilization coefficient C p for: The optimal tip speed ratio λ opt The input speed curve is obtained: The optimal tip speed ratio λ opt and the maximum wind energy utilization coefficient C pmax The input torque curve is obtained: Where β is the pitch angle, λ is the tip speed ratio, v is the incoming wind speed, ω is the rotor angular velocity, R is the blade radius; A, B, C, D, E, F, G, H are constants; n in is the input speed, n low is the minimum speed, n rate is the rated speed, v is the wind speed, v cut-in is the cut-in wind speed, v cut-off To cut out the wind speed, T is the wind speed corresponding to the rated speed of the wind turbine; in is the input torque, T rate is the rated torque, v is the wind speed, v rate is the rated wind speed.

3. The fatigue reliability adaptive evaluation method of a wind turbine gearbox under real-time wind conditions according to claim 1 is characterized in that: The method of obtaining the average wind speed per unit time and the turbulence standard deviation of the real-time wind speed history according to the set time period, and fitting the average wind speed distribution and the turbulence standard deviation characteristic value of the real-time wind speed history includes: Probability density function of the Weibull distribution of mean wind speed: Turbulence standard deviation function of instantaneous wind speed: σ1(V)=I 15 (0.75V+5.6) Where V is the average wind speed, k is the shape parameter of the Weibull distribution, c is the scale parameter of the Weibull distribution, and I 15 It is the characteristic value of turbulence intensity at an average wind speed of 15m / s.

4. The method for adaptively evaluating fatigue reliability of a wind turbine gearbox under real-time wind conditions according to claim 1 is characterized in that: The PSN fatigue curve expression of the part is established, and the statistical data of failure probability and life at different stress levels are obtained. The PSN fatigue curves of different probabilities of the part are obtained by fitting the statistical data; Latin hypercube sampling is used to generate equal probability intervals from 0% to 100%. value, Substitute the PSN fatigue curve formula to generate the PSN fatigue curve family, including: The PSN fatigue curve of the gear tooth surface contact is: The PSN fatigue curve of the gear tooth root bending is: The PSN fatigue curve of the bearing subsurface shear stress is: Where N H is the number of stress cycles to contact failure, σ H is the contact stress, C H is the contact fatigue stress constant, A H is the contact stress life coefficient, E δ is the coefficient of variation of logarithmic life span, Φ p -1 is the standard normal skewness of the failure probability, Z R is the roughness coefficient, Z v is the speed coefficient, Z L is the lubricant coefficient, Z W is the working hardening coefficient of the tooth surface, Z X is the contact strength size factor, N F is the number of stress cycles to bending failure, σ F is the bending stress, C F is the bending fatigue stress constant, A F is the bending stress life factor, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, Y ST is the stress correction factor, Y δrelT is the relative tooth root fillet sensitivity coefficient, Y RrelT is the relative tooth root surface condition coefficient, Y X is the bending strength size factor, N o is the number of cycles of bearing subsurface spalling, τ o is the subsurface shear stress, C o is the subsurface fatigue stress constant, A o is the subsurface shear stress life coefficient, E δ is the coefficient of variation of logarithmic life span, is the standard normal skewness of the failure probability, K τ is the effective stress concentration factor, ε τ is the size coefficient and β is the surface processing coefficient.

5. The method for adaptively evaluating fatigue reliability of a wind turbine gearbox under real-time wind conditions according to claim 1, characterized in that: The Monte Carlo simulation is performed on the fatigue samples of each component, a load spectrum of the current wind condition is established per unit time, and the Pavlou nonlinear fatigue damage accumulation model is used to solve the fatigue damage and residual strength history of the fatigue samples of each component, including: Fatigue damage accumulation D of parts k The function is as follows: The residual strength of the part material under multi-level loading is as follows: rr0-(r0-σ k )D k Where D k is the k-level fatigue damage accumulation, σ k is the kth level fatigue stress amplitude, n k is the number of stress cycles at level k, N k is the number of failure cycles at level k, q(σ)=(σσ u ) -0.75 is the nonlinear damage index, σ u is the ultimate bearing stress of the part, r is the residual strength of the part, r0=σ u is the initial static strength of the part.

6. The method for adaptively evaluating fatigue reliability of a wind turbine gearbox under real-time wind conditions according to claim 1, characterized in that: The adaptive bandwidth kernel density estimation method is used to fit the residual strength distribution of parts in real time; the fatigue reliability of each part failure mode is obtained by using the fitted residual strength distribution and fatigue stress distribution of the parts; According to the residual strength and fatigue stress of each fatigue sample of each part obtained by Latin hypercube sampling, the functional function value of each failure mode of each fatigue sample is calculated, including: The mean of the equivalent normalized variable is: m xn =x * -F -1 [F(x * )]s xn The standard deviation of the equivalent normalized variable is: Where f(x * ) is the probability density function of the design point, is the standard normal probability density function. The reliability index and dynamic reliability are as follows: R(t)=Φ[β(t)] where x * is the design point, F(x * ) is the cumulative distribution function of the design point, Φ -1 is the standard normal inverse cumulative distribution function, μ g is the mean value of the performance function, s g is the standard deviation of the performance function, μ rn is the normalized residual intensity mean, μ σn is the normalized fatigue stress mean, s rn is the normalized residual intensity standard deviation, s σn is the normalized fatigue stress standard deviation.

7. The method for adaptively evaluating fatigue reliability of a wind turbine gearbox under real-time wind conditions according to claim 1, characterized in that: The AIC information criterion is used to measure the goodness of fit of the Copula function, the Copula function with the highest goodness of fit is selected, the reliability of each failure mode is input into the correlation network, the optimal Copula function is selected layer by layer and connected, and the fatigue reliability of the wind turbine gearbox at the top layer is obtained, including: The expression of AIC information criterion is: AIC=-2lnL(θ)+2K Where L(θ) is the maximum likelihood function value of the Copula function, K is the number of parameters, the reliability of each failure mode is input into the correlation network, the optimal Copula function is selected layer by layer and connected, and finally the fatigue reliability of the wind turbine gearbox at the top layer is obtained.

Images