Gearbox fatigue analysis method and system
By performing abnormal frequency domain analysis and fatigue load quantification on the gearbox vibration data and combining it with a decision tree algorithm to build a model, the problem of inaccurate analysis in traditional methods is solved, and more accurate fatigue prediction and equipment health monitoring are achieved.
Patent Information
- Application Number
- CN202510890066.0
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-06-30
- Publication Date
- 2025-09-16
- Estimated Expiration
- 2045-06-30
AI Technical Summary
Traditional gearbox fatigue analysis methods are inaccurate in analyzing abnormal vibration of gearbox gears, resulting in large fatigue analysis errors.
Vibration sensors are used to monitor the vibration operating status data of the gearbox, extract abnormal vibration states and perform frequency domain structure analysis. Combined with fractional frequency mutation intensity analysis and decision tree algorithm, a gearbox fatigue analysis model is constructed to achieve accurate fatigue load quantification and prediction.
It improves the accuracy of gearbox fatigue analysis, reduces errors, provides a scientific basis and predictive capability for equipment health monitoring, and supports equipment maintenance and preventive management.
Smart Images

Figure CN120426189B_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of gearbox fatigue analysis, and in particular to a gearbox fatigue analysis method and system. Background Art
[0002] Intelligent fatigue analysis methods based on vibration monitoring, sensor technology, and data analysis have gradually gained widespread application. By using vibration sensors to monitor the operating status of gearboxes in real time, it is possible to obtain vibration data during operation. This data can then be deeply analyzed to extract potential abnormal signals and identify the fatigue state of the gearbox. Abnormal vibration signals are often a precursor to gearbox failure. Using frequency domain and time domain analysis techniques, the internal structural characteristics of the gearbox can be extracted from these signals, revealing its health status. However, a traditional gearbox fatigue analysis method suffers from inaccurate analysis of abnormal vibrations in gearbox gears, resulting in large errors in gearbox gear fatigue analysis. Summary of the Invention
[0003] Based on this, it is necessary to provide a gearbox fatigue analysis method and system to solve at least one of the above technical problems.
[0004] To achieve the above object, a gearbox fatigue analysis method is provided, the method comprising the following steps:
[0005] Step S1: monitoring the vibration operating state data of the wind turbine gearbox through a vibration sensor to obtain gear vibration operating state data; extracting abnormal vibration states from the gear vibration operating state data to obtain abnormal vibration state data, and then performing abnormal vibration frequency domain structure analysis to obtain abnormal vibration frequency domain structure data;
[0006] Step S2: performing fractional frequency multiplication mutation intensity analysis on the abnormal vibration frequency domain structure data to obtain sideband mutation intensity data; performing gearbox fatigue load simulation and quantification based on the sideband mutation intensity data to obtain gearbox fatigue load quantification data;
[0007] Step S3: performing fatigue load feature learning on the gearbox fatigue load quantification data to obtain fatigue load feature learning data; constructing a gearbox fatigue analysis model on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; and sending the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0008] Preferably, step S1 includes the following steps:
[0009] Step S11: monitoring the vibration operating state data of the wind turbine gearbox through a vibration sensor to obtain gear vibration operating state data;
[0010] Step S12: Filling missing values in the gear vibration operation state data to obtain vibration operation state filled data;
[0011] Step S13: extracting abnormal vibration state from the vibration operation state filling data to obtain abnormal vibration state data;
[0012] Step S14: performing abnormal vibration frequency domain structure analysis on the abnormal vibration state data to obtain abnormal vibration frequency domain structure data.
[0013] Preferably, step S2 includes the following steps:
[0014] Step S21: performing sideband phase intensity change analysis on the abnormal vibration frequency domain structure data to obtain sideband phase intensity change data;
[0015] Step S22: performing fractional frequency mutation intensity analysis on the sideband phase intensity change data to obtain sideband mutation intensity data;
[0016] Step S23: Deducing the gear meshing offset impact load increment index based on the sideband mutation intensity data to obtain the meshing offset impact load increment index;
[0017] Step S24: performing gearbox fatigue load simulation and quantification based on the meshing offset impact load increment index to obtain gearbox fatigue load quantification data.
[0018] Preferably, step S22 includes the following steps:
[0019] Step S221: calculating the adjacent mean difference of sideband phase frequency mutations on the sideband phase intensity change data to obtain the adjacent mean difference of phase frequency mutations;
[0020] Step S222: performing abnormal narrowband peak nonlinear regression analysis on the sideband phase intensity change data according to the adjacent mean difference of phase frequency mutations to obtain abnormal narrowband peak nonlinear regression data;
[0021] Step S223: performing kurtosis increment distribution analysis on the abnormal narrow-band peak nonlinear regression data to obtain abnormal kurtosis increment distribution data;
[0022] Step S224: performing sideband abnormal frequency kurtosis fractional factorial difference calculation on the abnormal kurtosis increment distribution data to obtain the abnormal frequency kurtosis fractional factorial difference;
[0023] Step S225: performing fractional octave mutation intensity analysis based on the abnormal frequency kurtosis fractional factorial difference and the abnormal narrowband peak nonlinear regression data to obtain sideband mutation intensity data.
[0024] Preferably, step S23 includes the following steps:
[0025] Step S231: performing gear meshing torsional vibration frequency instability mapping based on the sideband mutation intensity data to obtain gear meshing torsional vibration frequency instability mapping data;
[0026] Step S232: Deducing the progressive offset of the torsional vibration axial angular displacement according to the gear meshing torsional vibration frequency instability mapping data to obtain progressive offset data of the torsional vibration axial angular displacement;
[0027] Step S233: performing offset impact energy distribution calculation based on the progressive offset data of the torsional vibration axial angular displacement and the gear meshing torsional vibration frequency instability mapping data to obtain offset impact energy distribution calculation data;
[0028] Step S234: performing gear meshing offset impact load increment index deduction on the offset impact energy distribution calculation data to obtain the meshing offset impact load increment index.
[0029] Preferably, step S24 includes the following steps:
[0030] Step S241: acquiring basic attribute data of the gears in the gearbox, wherein the basic attribute data of the gears include gear pressure angle, gear center distance and gear material hardness data;
[0031] Step S242: performing pressure angle and center distance offset numerical mapping measurement on the gear pressure angle and gear center distance in the basic attribute data of the gear according to the meshing offset impact load increment index, and performing simulation and calculation of the variance of the local contact stress concentration distribution of the impact load to obtain the variance of the local contact stress concentration distribution of the gear;
[0032] Step S243: performing gear temperature rise equivalent mapping matching based on the meshing offset impact load increment index to obtain gear temperature rise equivalent mapping data;
[0033] Step S244: Deducing the gear deformation stiffness loss ratio from the gear material hardness data based on the gear local contact stress concentration distribution variance and gear temperature rise mapping data to obtain the gear deformation stiffness loss ratio;
[0034] Step S245: performing gearbox fatigue load simulation and quantification based on the gear deformation stiffness loss ratio to obtain gearbox fatigue load quantification data.
[0035] Preferably, step S244 includes the following steps:
[0036] The gear local contact stress concentration distribution variance is identified by stress multi-peak skewness to obtain local contact stress multi-peak skewness data;
[0037] Based on the multi-peak deflection data of local contact stress and the gear temperature rise equivalent mapping data, the thermal stress cyclic strength geometric increment analysis is performed to obtain the thermal stress cyclic strength geometric increment data;
[0038] The thermal stress intensity increment time series variance regression analysis is performed on the thermal stress cyclic intensity geometric increment data to obtain the intensity increment time series variance regression data;
[0039] Based on the strength increment time series variance regression data, the gear material hardness data is simulated for material toughness equal loss to obtain material toughness equal loss data;
[0040] The gear deformation stiffness loss ratio is deduced based on the strength increment time series variance regression data and the material toughness equivalent loss data to obtain the gear deformation stiffness loss ratio.
[0041] Preferably, step S3 includes the following steps:
[0042] Step S31: performing convolution processing on the gearbox fatigue load quantization data to obtain gearbox fatigue load convolution data;
[0043] Step S32: performing fatigue load feature learning on the gearbox fatigue load convolution data to obtain fatigue load feature learning data;
[0044] Step S33: constructing a gearbox fatigue analysis model based on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model;
[0045] Step S34: Send the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0046] Preferably, the present invention further provides a gearbox fatigue analysis system for executing the gearbox fatigue analysis method described above, the gearbox fatigue analysis system comprising:
[0047] The abnormal vibration frequency domain analysis module is used to monitor the vibration operation status data of the wind turbine gearbox through a vibration sensor to obtain the gear vibration operation status data; perform abnormal vibration state extraction on the gear vibration operation status data to obtain abnormal vibration state data, and then perform abnormal vibration frequency domain structure analysis to obtain abnormal vibration frequency domain structure data;
[0048] The fatigue load simulation and quantification module is used to perform fractional frequency mutation intensity analysis on abnormal vibration frequency domain structure data to obtain sideband mutation intensity data; based on the sideband mutation intensity data, the gearbox fatigue load is simulated and quantified to obtain gearbox fatigue load quantification data;
[0049] The fatigue analysis model construction module is used to perform fatigue load feature learning on the gearbox fatigue load quantification data to obtain fatigue load feature learning data; construct a gearbox fatigue analysis model on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; and send the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0050] The beneficial effect of the present invention is that by monitoring the vibration of the wind turbine gearbox through a vibration sensor, the gear vibration operating status data can be obtained in real time, providing important basic information for subsequent fault diagnosis. By extracting abnormal vibration status data, the abnormal operating conditions of the equipment can be effectively identified. Further abnormal vibration frequency domain structure analysis not only helps to gain a deep understanding of the operating status of the gearbox, but also reveals potential fault signs, such as gear wear, bearing failure, etc., providing a more accurate basis for equipment health monitoring. The fractional frequency multiplication mutation intensity analysis of abnormal vibration frequency domain structure data can extract more detailed frequency domain features from the vibration data, especially the mutation intensity analysis of the sideband, which can effectively discover potential fatigue damage in the gearbox. This frequency feature analysis helps to more accurately quantify the fatigue load of the gearbox and provide scientific input data for subsequent fatigue analysis. By combining frequency domain analysis with fatigue load simulation and quantification, the fatigue condition of the gearbox can be better predicted, thereby providing data support for equipment maintenance and preventive management. By using fatigue load quantification data to learn fatigue load characteristics, we can further explore the fatigue characteristics and service life of the gearbox, reveal the specific impact of the load on the health of the gearbox, and build a fatigue analysis model based on the decision tree algorithm. It can accurately predict the fatigue damage condition of the gearbox based on historical data and learned characteristics. This model not only has high prediction accuracy, but also can adapt to changes in different working conditions and different equipment. After sending the fatigue analysis model to the terminal, it can realize real-time monitoring and automated analysis, making the fatigue analysis of the gearbox more intelligent and accurate, and significantly improving the efficiency and safety of equipment management. Therefore, the present invention is an optimization of a traditional gearbox fatigue analysis method, which solves the problem of inaccurate abnormal vibration analysis of gearbox gears in a traditional gearbox fatigue analysis method, thereby causing large errors in gearbox gear fatigue analysis, improves the accuracy of abnormal vibration analysis of gearbox gears, and reduces the error of gearbox gear fatigue analysis. BRIEF DESCRIPTION OF THE DRAWINGS
[0051] Figure 1 A schematic flow chart of the steps of a gearbox fatigue analysis method is shown;
[0052] Figure 2 for Figure 1 Detailed implementation steps of step S2 in FIG.
[0053] Figure 3 for Figure 1 Detailed implementation steps of step S3 in FIG. DETAILED DESCRIPTION
[0054] See also Figures 1 to 3 , a gearbox fatigue analysis method, the method comprising the following steps:
[0055] Step S1: monitoring the vibration operating state data of the wind turbine gearbox through a vibration sensor to obtain gear vibration operating state data; extracting abnormal vibration states from the gear vibration operating state data to obtain abnormal vibration state data, and then performing abnormal vibration frequency domain structure analysis to obtain abnormal vibration frequency domain structure data;
[0056] Step S2: performing fractional frequency multiplication mutation intensity analysis on the abnormal vibration frequency domain structure data to obtain sideband mutation intensity data; performing gearbox fatigue load simulation and quantification based on the sideband mutation intensity data to obtain gearbox fatigue load quantification data;
[0057] Step S3: performing fatigue load feature learning on the gearbox fatigue load quantification data to obtain fatigue load feature learning data; constructing a gearbox fatigue analysis model on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; and sending the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0058] In the embodiment of the present invention, reference Figure 1 The above is a schematic flow chart of the steps of a gearbox fatigue analysis method of the present invention. In this example, the gearbox fatigue analysis method includes the following steps:
[0059] Step S1: monitoring the vibration operating state data of the wind turbine gearbox through a vibration sensor to obtain gear vibration operating state data; extracting abnormal vibration states from the gear vibration operating state data to obtain abnormal vibration state data, and then performing abnormal vibration frequency domain structure analysis to obtain abnormal vibration frequency domain structure data;
[0060] In this embodiment of the present invention, triaxial piezoelectric acceleration vibration sensors are installed above the wind turbine gearbox housing, near the input and output shafts, to collect gearbox vibration acceleration signals during operation. The sampling frequency is set to 25.6 kHz, and the sampling time window is set to 5 seconds to ensure that the high-frequency dynamic changes during gear meshing are fully captured. The collected raw vibration operating status data is bandpass filtered with a filter cutoff frequency set between 1 kHz and 10 kHz to remove low-frequency environmental interference and high-frequency background noise. The filtered signal is then Hilbert transformed and its envelope signal extracted. Time-frequency analysis is performed using a short-time Fourier transform (STFT) to determine whether the vibration signal contains impulsive pulse characteristics. A multidimensional threshold judgment condition is then set using the root mean square value, kurtosis, and crest factor to extract abnormal pulse signal segments, generating abnormal vibration status data. Fast Fourier transform (FFT) is performed on the abnormal signal segment to extract frequency domain features, identify the main frequency, harmonics, sidebands and their amplitude distribution to constitute the abnormal vibration frequency domain structure data, which is used to reflect the abnormal change characteristics of the gear meshing state under potential fatigue failure.
[0061] Step S2: performing fractional frequency multiplication mutation intensity analysis on the abnormal vibration frequency domain structure data to obtain sideband mutation intensity data; performing gearbox fatigue load simulation and quantification based on the sideband mutation intensity data to obtain gearbox fatigue load quantification data;
[0062] In this embodiment of the present invention, frequency domain data is sliced using a sliding window in steps of 0.1 kHz. The energy ratio between the main frequency signal and its sidebands in each frequency slice is calculated. Subsequently, the signal is reconstructed in the fractional domain using the fractional Fourier transform method to extract local energy discontinuity points within each order of the fractional frequency components. Sideband discontinuity points are detected on the reconstructed spectrum using a continuous wavelet transform. Using the Morlet wavelet as a basis function, local extreme points of frequency amplitude variations in the high-order sideband regions are extracted. The density variation of discontinuity points at the sideband locations is statistically analyzed to generate sideband discontinuity intensity data. Next, based on the discontinuity intensity data, the deconvolution method is used to infer the excitation function characteristics caused by the gear mesh impact. Combined with the energy density distribution of the frequency components corresponding to each discontinuity point, the impact load energy integral per unit time is deduced. Combined with the speed data, the periodic impact response load is then estimated. Using the impeller moment of inertia and gear mesh stiffness from structural dynamics as input parameters, the fatigue load amplitude under the corresponding load is calculated using an equivalent mechanical transformation model, generating quantitative data on gearbox fatigue loads.
[0063] Step S3: performing fatigue load feature learning on the gearbox fatigue load quantification data to obtain fatigue load feature learning data; constructing a gearbox fatigue analysis model on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; and sending the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0064] In an embodiment of the present invention, after obtaining the above-mentioned fatigue load quantification data, a time series sliding window process is adopted, and the sliding window length is set to 1024 data points and the step size is 128 points. The characteristic parameters such as energy density, extreme amplitude, load average amplitude change rate, zero crossing rate, maximum gradient change, etc. in each window are extracted, and these parameter sequences are normalized. The normalized sequence data is subjected to dimensionality reduction processing using principal component analysis (PCA), and the first three principal components are selected as fatigue load feature learning data, and a sample set is constructed in the form of a feature vector. In the modeling process, the ID3 algorithm is used to construct a decision tree model, and the average amplitude, maximum fluctuation gradient and energy density in the fatigue load quantification characteristics are used as partitioning attributes. The best partitioning node is selected by calculating the information gain value, and the tree structure is repeatedly constructed in the training set. The corresponding characteristic conditions and fatigue level classification labels are recorded at each node, and finally the gearbox fatigue analysis model is constructed. After the model is constructed, the model data is uploaded to the embedded industrial control terminal through the TCP / IP communication protocol. The terminal receives the model parameters and combines them with the real-time vibration data to perform the intelligent diagnosis task of the gearbox fatigue status, realizing the online fatigue status assessment and early warning function.
[0065] Step S1 includes the following steps:
[0066] Step S11: monitoring the vibration operating state data of the wind turbine gearbox through a vibration sensor to obtain gear vibration operating state data;
[0067] Step S12: Filling missing values in the gear vibration operation state data to obtain vibration operation state filled data;
[0068] Step S13: extracting abnormal vibration state from the vibration operation state filling data to obtain abnormal vibration state data;
[0069] Step S14: performing abnormal vibration frequency domain structure analysis on the abnormal vibration state data to obtain abnormal vibration frequency domain structure data.
[0070] In this embodiment of the present invention, a triaxial piezoelectric accelerometer is installed symmetrically on the wind turbine gearbox housing, near the input shaft bearing cap and the high-speed output shaft, to ensure a complete meshing state response. The sensor is secured using bolt pre-tightening to ensure close contact and improve signal response accuracy. The data acquisition system utilizes an NI cDAQ-9178 data acquisition chassis and an NI 9234 vibration acquisition module. The sampling frequency is set to 25.6 kHz, the recording time is 5 seconds, and each acquisition generates 128,000 sample points of raw vibration operating state data. Wind speed and rotational speed data are recorded simultaneously during the acquisition process. The rotational speed is obtained by a Hall-effect speed sensor mounted on the main shaft at a sampling frequency of 1 kHz and is used for synchronous compensation processing in subsequent analysis. The vibration signal is filtered through a low-pass filter to remove high-frequency noise above 10 kHz to prevent aliasing. A Butterworth filter with a cutoff frequency of 10 kHz and a filter order of 4 is used for filtering. The filtered signal is stored in CSV format. The gear vibration operating status data was tested for integrity and missing value analysis. Sampling segments were reconstructed using a sliding time window method with a window width of 1000 sampling points and a step size of 200 points. Window-by-window analysis was performed to determine if there were any abnormal time periods or sampling signal interruptions. For detected sampling discontinuities, missing values were filled using a spline interpolation method. This method employed cubic spline interpolation, constructing a cubic polynomial interpolation function for the 100 valid sample points before and after. The interpolation function was fitted using the least squares method, forcing continuous derivatives and function values to ensure smooth data flow. The filled signal was stored in MAT format for subsequent analysis. The data filling process was accompanied by monitoring the maximum, minimum, and mean squared value of the original signal to determine if the filled segment exhibited abrupt or discontinuous deviations. Segments with deviations exceeding ±10% of the mean squared value were re-interpolated to ensure data structure stability. Ultimately, complete, missing-free filled vibration operating status data was obtained. Abnormal vibration states were extracted from the filled vibration operating status data, and impact vibration signal segments were identified using continuous wavelet transform combined with envelope demodulation. First, the envelope of the vibration signal is calculated using the Hilbert transform, and an energy envelope curve is constructed for the envelope signal. Statistical indicators such as kurtosis, crest factor, and root mean square amplitude within each 1000-point sliding window are used to determine whether the signal exhibits shocks or periodic mutations. A kurtosis greater than 4.5 and a crest factor greater than 3.0 are used as anomaly thresholds, while a root mean square amplitude greater than 1.5 times the mean of the entire sample sequence is used as a linkage criterion. Anomalous data is extracted within each time window that meets these three criteria. Abnormal segments are labeled and extracted, accompanied by their time index and original numerical sequence. The extracted data ultimately constitutes the abnormal vibration state data. This processing method utilizes wavelet packet functions and time-domain statistical functions implemented in MATLAB to perform computations, combining a looping and iterative approach to extract complete abnormal state segments from continuous time periods.The extracted abnormal vibration state data was analyzed in the frequency domain. The fast Fourier transform (FFT) method was used to convert the spectrum of each abnormal data segment. A 2048-point Hamming window function was used to avoid spectral leakage. After conversion, the frequency domain range was 0 to 12.8 kHz, with a frequency resolution of 12.5 Hz. The main peak frequency, secondary peak frequency, and their corresponding amplitudes were extracted from the spectrum results. The symmetry of the sideband distribution was calculated, and the ratio of the main frequency peak to the upper and lower sideband amplitudes was further calculated to analyze the presence of meshing frequency or gear fault characteristic frequency components. Kurtosis spectrum analysis was used to extract the kurtosis spectrum of sharp narrowband frequency segments in the frequency domain to identify abnormal narrowband frequency components. The extracted spectral features, sideband symmetry index, and narrowband abnormal frequency segments were combined to construct a frequency domain structure description matrix. This matrix was saved in a four-dimensional matrix format of frequency component, amplitude, frequency band index, and local peak density. This ultimately formed the abnormal vibration frequency domain structure data for subsequent fatigue load analysis modeling. The processing flow is completely based on the implementation functions of FFT, wavelet transform and spectrum feature extraction algorithm in MATLAB signal processing toolbox.
[0071] Step S2 includes the following steps:
[0072] Step S21: performing sideband phase intensity change analysis on the abnormal vibration frequency domain structure data to obtain sideband phase intensity change data;
[0073] Step S22: performing fractional frequency mutation intensity analysis on the sideband phase intensity change data to obtain sideband mutation intensity data;
[0074] Step S23: Deducing the gear meshing offset impact load increment index based on the sideband mutation intensity data to obtain the meshing offset impact load increment index;
[0075] Step S24: performing gearbox fatigue load simulation and quantification based on the meshing offset impact load increment index to obtain gearbox fatigue load quantification data.
[0076] As an example of the present invention, refer to Figure 2 As shown, in this example, step S2 includes:
[0077] Step S21: performing sideband phase intensity change analysis on the abnormal vibration frequency domain structure data to obtain sideband phase intensity change data;
[0078] In this embodiment of the present invention, sideband phase intensity variation analysis is performed on the acquired abnormal vibration frequency domain structure data. First, each abnormal spectral segment in the frequency domain structure data is subjected to a short-time Fourier transform (SFT). A 512-point Hanning window with a sliding step size of 128 points is used to ensure a balance between sufficient frequency domain resolution and time resolution. For the upper and lower sideband regions adjacent to the main meshing frequency and its multiples in each short-time spectrum matrix, amplitude and phase information within the sideband frequency range is extracted, with a frequency window width set to ±120 Hz. Phase difference mean and standard deviation analysis is used to calculate the phase change rate of each sideband position within adjacent time windows, and a phase intensity variation curve is constructed. The angular standard deviation method is used to convert the phase intensity variation into a stability index. Phase mutation segments in each spectrum segment are extracted, and frequency points with a phase change slope greater than 20 degrees / second and a standard deviation greater than 15 degrees are selected as valid phase mutation frequencies. All valid mutation frequencies are calibrated in the spectrum, and a two-dimensional frequency-phase variation intensity matrix is constructed. The data is saved in JSON format and named "sideband phase intensity variation data."
[0079] Step S22: performing fractional frequency mutation intensity analysis on the sideband phase intensity change data to obtain sideband mutation intensity data;
[0080] In this embodiment of the present invention, the sideband phase intensity variation data is analyzed for fractional octave mutation intensity. Based on fractional octave theory, the ratio of the sideband frequency to the main meshing frequency is set between 0.25 and 4.0, with a step size of 0.05, and the frequency ratio distribution is discretely sampled. For each ratio node, a determination is made as to whether there is a concentrated distribution of effective phase octave frequency points. If there are three or more phase intensity octave mutation frequency points within the ratio range of ±0.03, and the mutation intensity is greater than two standard deviations of the mean, then the octave point is considered a fractional octave outlier. The average, maximum, and distribution density of the mutation intensity of each fractional octave point are calculated to construct a three-dimensional distribution map of fractional octave-intensity. The local maximum region is extracted to form an intensity mutation envelope. The envelope boundary is limited to mutation intensity greater than the mean + 3σ. The data within the envelope is the sideband mutation intensity data, which is saved as a triple array, representing the octave coefficient, mutation frequency, and mutation intensity value, respectively.
[0081] Step S23: Deducing the gear meshing offset impact load increment index based on the sideband mutation intensity data to obtain the meshing offset impact load increment index;
[0082] In an embodiment of the present invention, the gear meshing offset impact load incremental index is deduced based on the sideband mutation intensity data. First, the frequency multiplication coefficient and frequency position corresponding to the peak frequency of the mutation intensity are extracted, and compared with the typical meshing frequency multiplication frequency in the gear transmission theory to identify the position where there is a significant offset. The frequency offset between the offset frequency and the theoretical frequency multiplication frequency is calculated, and the offset impact factor function is constructed in combination with the mutation intensity value, which is defined as the normalized product of the frequency offset and the mutation intensity. The corresponding offset impact factor is calculated for all mutation frequency points to form an offset impact factor vector sequence. Subsequently, the mean, variance, and maximum value of the factor sequence are statistically calculated, and the mutation factors that are more than twice the mean are integrated, and the integral result is extracted as the meshing offset impact load incremental index. The index is expressed as a single-valued number, reflecting the upward trend of the additional impact load intensity caused by meshing mismatch under the current abnormal vibration state. The deduction process is implemented using the signal processing and numerical integration tool functions in MATLAB.
[0083] Step S24: performing gearbox fatigue load simulation and quantification based on the meshing offset impact load increment index to obtain gearbox fatigue load quantification data.
[0084] In this embodiment of the present invention, gearbox fatigue load simulation and quantification are performed based on the mesh offset impact load increment index. First, a basic model of the gear mesh force is established. The initial load reference value is determined based on the rotor input torque and the main meshing frequency, with a typical torque value of 6500 Nm. The main meshing frequency is calculated from the gear ratio and the main shaft speed. Based on the deduced load increment index, the basic meshing force is amplitude modulated. The increment index is multiplied by a modulation coefficient of 0.3 and used as an additional load factor. A nonlinear perturbation function is superimposed within the corresponding meshing cycle to construct a simulated impact load curve. The load curve is constructed based on the periodic meshing frequency, with the cycle length consistent with the main frequency. The simulation duration is set to 10 seconds, and the time resolution is 1 ms. The modulated impact load is combined with the basic static load to generate fatigue simulation load data containing irregular amplitude variations. The final constructed fatigue load quantification data is stored in an array structure, which records the time series, instantaneous load value, impact modulation term value, and meshing cycle label.
[0085] Step S22 includes the following steps:
[0086] Step S221: calculating the adjacent mean difference of sideband phase frequency mutations on the sideband phase intensity change data to obtain the adjacent mean difference of phase frequency mutations;
[0087] Step S222: performing abnormal narrowband peak nonlinear regression analysis on the sideband phase intensity change data according to the adjacent mean difference of phase frequency mutations to obtain abnormal narrowband peak nonlinear regression data;
[0088] Step S223: performing kurtosis increment distribution analysis on the abnormal narrow-band peak nonlinear regression data to obtain abnormal kurtosis increment distribution data;
[0089] Step S224: performing sideband abnormal frequency kurtosis fractional factorial difference calculation on the abnormal kurtosis increment distribution data to obtain the abnormal frequency kurtosis fractional factorial difference;
[0090] Step S225: performing fractional octave mutation intensity analysis based on the abnormal frequency kurtosis fractional factorial difference and the abnormal narrowband peak nonlinear regression data to obtain sideband mutation intensity data.
[0091] In this embodiment of the present invention, the sideband phase intensity change data is subjected to calculation of the adjacent mean difference of sideband phase frequency mutations. First, all phase intensity change curves within the sideband frequency range are sorted in ascending frequency order to extract phase change mutation points, which are defined as frequency points with a phase difference greater than 20 degrees. With the mutation point as the center, the phase change values of the three frequency points immediately preceding and following it are taken. The average phase change value of these six points is calculated and subtracted from the phase change value of the mutation point itself to obtain the adjacent mean difference of the mutation point. To ensure processing consistency, all sideband mutation points are batch processed using this method, and the results are stored as a triplet of frequency point, mutation value, and adjacent mean difference. To eliminate the influence of background noise, all adjacent mean difference values are z-score normalized, and the abnormal threshold is set to points greater than the mean plus two standard deviations as significant mutation points. This highlights non-stationary mutation frequency segments and provides input indicators for subsequent nonlinear regression. Based on the phase frequency mutation adjacent mean difference data obtained in the previous step, the original sideband phase intensity change data is subjected to abnormal narrowband peak nonlinear regression analysis. First, the sideband frequency domain was divided into non-overlapping frequency bands of 50 Hz width. The number of breakpoints within each band and the average of their mean deviation amplitudes were counted. Nonlinear regression was performed on the areas where breakpoints were concentrated using a quintic polynomial fitting method, with frequency as the independent variable and phase intensity as the dependent variable. The fitting residuals were obtained through least squares fitting, and the standard deviation of the fitting residuals was calculated. After the fitting was completed, the data at the maximum peak of the fitting curve and within the 5 Hz frequency range to the left and right were extracted. The deviation between the local maximum and the fitting trend line was calculated and defined as an anomalous narrowband peak. If the deviation of the peak exceeded three times the standard deviation of the residual across the entire frequency band, it was marked as a nonlinear anomalous peak. The combination of the anomalous point frequency, fitting residual value, and deviation value was recorded as the anomalous narrowband peak nonlinear regression data. Kurtosis incremental distribution analysis was performed on the obtained anomalous narrowband peak nonlinear regression data. Specifically, the bandpass filtering results corresponding to the anomalous peak frequency were extracted from the original time-domain vibration signal. The filter used an 8th-order IIR bandpass filter with a center frequency set to the anomalous peak frequency and a bandwidth of ±10 Hz. Kurtosis was calculated for the filtered signal segments (each segment was 1024 points long and had a 50% overlap) using the fourth-order center distance divided by the square of the second-order center distance. The difference between the kurtosis value of each segment and the previous segment was extracted to form a kurtosis increment sequence. The kurtosis increments over the entire analysis period were statistically analyzed using a distribution histogram, divided into 30 intervals. Kernel density estimation was performed on the frequency density, and the increment values and probability densities at the distribution peaks were extracted to generate a kurtosis increment distribution map. In this distribution map, each abnormal peak frequency point corresponds to a corresponding kurtosis increment density distribution curve. The results are stored as a two-dimensional matrix, named abnormal kurtosis increment distribution data.The fractional factorial differences of the sideband abnormal frequency kurtosis are calculated for the abnormal kurtosis increment distribution data. The calculation steps are as follows: first, the kurtosis increment distribution is reconstructed into a sequence of kurtosis increments with equal frequency intervals according to the frequency distribution. The kurtosis increment values at each frequency point are then subjected to first-, second-, and third-order differences in sequence order. Fractional-order differences are then processed, with the order of the fractional differences set to 1.5. Discrete differences are calculated using the Grünwald–Letnikov approximation. The maximum and standard deviation of the resulting fractional-order differences are calculated. If the fractional-order difference of a frequency point is greater than twice the standard deviation of the sequence, it is marked as an abnormal frequency. Finally, all frequency points that meet the above conditions are the abnormal frequency kurtosis fractional-factorial difference results. The kurtosis increments and fractional-order differences of the corresponding frequency points are combined to form a triple record, forming the abnormal frequency kurtosis fractional-factorial difference data matrix. The abnormal frequency kurtosis fractional-factorial difference data are combined with the abnormal narrow-band peak nonlinear regression data to perform fractional-octave mutation intensity analysis. The frequency points in the two data sets were matched, with a matching tolerance of ±5 Hz. For each successfully matched frequency point, its corresponding harmonic coefficient (based on the main meshing frequency) was extracted, and its mutation intensity was calculated, defined as the product of the fractional factorial difference in kurtosis and the nonlinear regression deviation. All matching frequency points were sorted in ascending order by harmonic coefficient to construct a fractional harmonic mutation intensity curve. A local extreme value detection algorithm (with a local window length of 5 points) was applied to this curve to extract the peak point, which was the strong mutation frequency location. The corresponding mutation intensity value was then used as the sideband mutation intensity data. The final output format was a two-dimensional array, with each row representing a harmonic mutation point. The data included harmonic coefficient, mutation frequency, and mutation intensity value, and was used in the subsequent fatigue load increment calculation.
[0092] Step S23 includes the following steps:
[0093] Step S231: performing gear meshing torsional vibration frequency instability mapping based on the sideband mutation intensity data to obtain gear meshing torsional vibration frequency instability mapping data;
[0094] Step S232: Deducing the progressive offset of the torsional vibration axial angular displacement according to the gear meshing torsional vibration frequency instability mapping data to obtain progressive offset data of the torsional vibration axial angular displacement;
[0095] Step S233: performing offset impact energy distribution calculation based on the progressive offset data of the torsional vibration axial angular displacement and the gear meshing torsional vibration frequency instability mapping data to obtain offset impact energy distribution calculation data;
[0096] Step S234: performing gear meshing offset impact load increment index deduction on the offset impact energy distribution calculation data to obtain the meshing offset impact load increment index.
[0097] In an embodiment of the present invention, the gear meshing torsional vibration frequency instability mapping is performed based on the sideband mutation intensity data. First, the main meshing frequency and its low-order harmonic characteristics must be confirmed. During the processing, the rotation frequency of the wind turbine gearbox output shaft is used as the reference fundamental frequency, the meshing frequency and the sideband mutation frequency are calibrated, and the mutation point threshold of the sideband mutation intensity data is set to the mean plus two standard deviations. All mutation point frequencies are normalized to the main meshing frequency and mapped to the equivalent harmonic axis domain. Subsequently, in the normalized frequency domain, the harmonic sequence distribution density of the frequency points whose mutation intensity is higher than the threshold is counted, and the continuity of their occurrence frequency, jump amplitude and frequency domain offset are calculated, which is defined as the frequency instability index. On this basis, a two-dimensional mapping matrix is constructed, in which the horizontal axis is the time index or operating condition number, the vertical axis is the normalized frequency multiple, and the matrix unit value is the frequency instability index value of the corresponding point. To enhance recognition accuracy, the local sliding standard deviation of frequency variation (with a window length of five frequency points) is introduced as a frequency jitter factor and superimposed on the normalized value of the mutation intensity to form enhanced instability mapping data. The progressive offset of the torsional axial angular displacement of the gear mesh torsional vibration frequency instability mapping data is deduced. First, the frequency instability index of each time slice in the mapping matrix is analyzed. The time slice corresponding to the column with the largest frequency variation amplitude is selected as the key period. The instantaneous phase of the torsional vibration signal in this time slice is extracted using the Hilbert transform, and the instantaneous frequency is then determined by phase differentiation. The angular velocity variation is inferred from the instantaneous frequency change rate, and the axial angular displacement time series is obtained by integrating it. Due to the periodic nature of the meshing system, the main meshing frequency is used for period synchronization during the processing. The angular displacement variation within each main cycle is statistically analyzed, and the average angular displacement variation difference between adjacent cycles is analyzed and defined as the progressive offset value of the angular displacement. The progressive angular displacement values over multiple consecutive cycles are constructed as time series data, recording the displacement direction (positive or negative), displacement rate (degrees / ms), and displacement stability (standard deviation). This data is then output as torsional vibration axial angular displacement progressive offset data. Based on this torsional vibration axial angular displacement progressive offset data and the gear mesh torsional vibration frequency instability mapping data, the offset impact energy distribution calculation is performed. First, within the time window corresponding to each progressive angular displacement point, the torsional vibration acceleration signal is extracted and second-order integrated to obtain the displacement. The velocity of this segment of the signal is then calculated, and the impact kinetic energy is calculated using the known moment of inertia parameters (using the NREL 750kW gearbox sample as an example, the output shaft moment of inertia is 37.6 kg·m²). The angular displacement offset within each cycle is multiplied by the change in angular velocity and then by the moment of inertia to obtain the impact energy per unit time. This impact energy data is sorted by time to generate an energy sequence, and metrics such as the total impact energy value, maximum value, and peak frequency within each time window are statistically analyzed. This results in a three-dimensional energy distribution array with dimensions of time, frequency, and impact energy intensity.If the impact energy within multiple time windows shows a concentrated peak (defined as exceeding twice the overall mean value for three consecutive cycles), it is marked as a high-intensity impact segment. Its frequency position, duration, and cumulative energy intensity are also recorded and output as offset impact energy distribution calculation data. To deduce the gear mesh offset impact load increment index from the offset impact energy distribution calculation data, a time-accumulated integral is performed based on the obtained energy distribution sequence. For each high-intensity impact segment, a time window length is set (for example, 50 ms), and the ratio of the total impact energy to the angular displacement increment within the window is calculated as the impact unit angular load. The impact unit angular load values for all high-intensity impact segments are serialized, and the time interval distribution density is calculated to obtain the time-domain clustering characteristics of the impact load. An impact intensity factor, defined as the logarithm of the product of the square of the impact energy and the impact frequency, is introduced and combined with the above impact unit angular load values to construct an impact load index factor sequence. After noise reduction of the impact load index factor using an exponential smoothing filter (smoothing coefficient α is 0.35), local extreme values are extracted, and the difference between the maximum peak and the baseline is calculated as the load increment index. Finally, the timestamp, load increment index value, and corresponding frequency range triplet corresponding to each strong impact segment are summarized as the result, and the output is the meshing offset impact load increment index.
[0098] Step S24 includes the following steps:
[0099] Step S241: acquiring basic attribute data of the gears in the gearbox, wherein the basic attribute data of the gears include gear pressure angle, gear center distance and gear material hardness data;
[0100] Step S242: performing pressure angle and center distance offset numerical mapping measurement on the gear pressure angle and gear center distance in the basic attribute data of the gear according to the meshing offset impact load increment index, and performing simulation and calculation of the variance of the local contact stress concentration distribution of the impact load to obtain the variance of the local contact stress concentration distribution of the gear;
[0101] Step S243: performing gear temperature rise equivalent mapping matching based on the meshing offset impact load increment index to obtain gear temperature rise equivalent mapping data;
[0102] Step S244: Deducing the gear deformation stiffness loss ratio from the gear material hardness data based on the gear local contact stress concentration distribution variance and gear temperature rise mapping data to obtain the gear deformation stiffness loss ratio;
[0103] Step S245: performing gearbox fatigue load simulation and quantification based on the gear deformation stiffness loss ratio to obtain gearbox fatigue load quantification data.
[0104] In an embodiment of the present invention, when obtaining the basic attribute data of the gears in the gearbox, it is necessary to obtain the basic geometric and material parameters of each gear pair through design drawings, technical documents or measurements. Specifically, it includes the pressure angle, center distance and material hardness value of the gear. Taking a certain type of 3MW wind power gearbox as an example, the pressure angle of the first-stage transmission gear is 20°, the center distance is 215 mm, and 18CrNiMo7-6 carburized and quenched steel is used, with a Rockwell hardness of 61 HRC. During the data collection process, it is necessary to check the deviation between the theoretical design value in the drawing and the actual measured value on site, and at the same time ensure that the material hardness data source is the hardness test result of the standard metallographic specimen, and the test accuracy does not exceed ±1 HRC. All data should be recorded in a structured manner in a table, and the fields include gear number, module, number of teeth, pressure angle, center distance, tooth width, hardness value and heat treatment process number to ensure that they can be accurately indexed and matched in subsequent calculations. The gear pressure angle and center distance are numerically mapped according to the meshing offset impact load increment index. With the increment index as the horizontal axis, construct its relationship with the pressure angle change and center distance changes ΔaThe mapping relationship between the gear teeth and the contact stress was calculated using a numerical mapping formula, with the index value interval set to a step size of 0.05. The mapping formula was obtained by calibration with historical fatigue load test data. After numerical mapping, the gear contact stress was simulated and estimated using the contact force calculation formula. Combined with the gear tooth contact area (calculated based on the tooth width and instantaneous meshing line length), the stress distribution per unit contact area was derived using the Hertz contact theory formula. A two-dimensional mesh was established on the meshing tooth surface, with each mesh node corresponding to a local stress value. The variance of these node values was then calculated as an indicator of the stress concentration distribution variance. In the experiment, a gear sample with a tooth width of 90 mm was selected. The mesh was divided into 2 mm × 2 mm cells, forming a 45 × 45 node matrix. The standard deviation of the stress distribution was 172 MPa, corresponding to a local impact load mapping value of 32.7 kN. The simulation was repeated at each incremental index level to obtain stress concentration trend curves under different loads. Gear temperature rise equivalent mapping was performed. Using the incremental index of mesh offset impact load as an input variable, a mapping relationship between it and tooth surface temperature rise was established. Based on thermal-mechanical coupled fatigue test data, the gear surface temperature rise was calculated as the ratio of the friction work per unit area of the tooth surface to the heat conduction heat dissipation rate multiplied by the contact time. Taking a typical incremental index of 0.67 as an example, the corresponding contact impact energy is 68.3 J. For a gear tooth surface area of 85 × 90 mm², a heat conduction coefficient of 23 W / (m·K) under air cooling conditions, and a contact time of 0.15 s, the calculated surface temperature rise is approximately 21.8°C. The temperature rise values at each index level were tabulated to construct an equivalent mapping dataset. The fields include incremental index, unit impact energy, tooth surface temperature rise, heat conduction efficiency, and cooling conditions. Based on the variance of the local contact stress concentration distribution and the equivalent mapping data for temperature rise, the gear deformation stiffness loss ratio was deduced from the gear material hardness data. The process first evaluates the microscopic deformation caused by local contact stress variance, assuming the material's yield strength is a known value of 1100 MPa. When the local stress variance exceeds 20% of the material's yield strength, it is considered a potential source of microplastic deformation. The hardness value is then dynamically corrected based on the known softening rate (approximately 0.5 HRC per 10°C in the 20-100°C range) and the effect of temperature rise on material thermal softening. Finally, the gear's overall deformation stiffness loss ratio is calculated by comparing the corrected hardness change to the initial hardness, combined with gear geometry analysis (using the deformation influence factor in the mesh stiffness formula). For example, when the temperature rise is 25°C and the stress concentration distribution variance is 185 MPa, the corresponding corrected hardness decreases by 2.1 HRC, resulting in a calculated stiffness loss ratio of 0.17, representing a 17% decrease compared to the initial state. This gear deformation stiffness loss ratio is then used to quantify gearbox fatigue load simulation. First, the stiffness loss ratio of each gear group is input into the load distribution matrix of the whole machine as the local weakening parameter.Based on the multi-degree-of-freedom rigid body dynamics model of the transmission system, the torque transmission relationship of each level of transmission elements is constructed, and the stiffness coefficient of each node is corrected. The instantaneous torque, reaction force and bearing support reaction force changes are calculated in each operating cycle using a dynamic iteration method. Taking the main stage gear stiffness loss ratio of 0.12 as an example, the simulation shows that the decrease in the output shaft torque transmission capacity is about 6.8%, and the corresponding fatigue load coefficient increases by 1.21 times. The fatigue load simulation results of each level of gear are superimposed, and the output is a three-dimensional matrix data containing time series, load amplitude, and number of cycles. Finally, a gearbox fatigue load quantitative data set is formed, with fields including parameters such as gear stage, loss ratio, load peak, number of cycles and load correction coefficient.
[0105] Step S244 includes the following steps:
[0106] The gear local contact stress concentration distribution variance is identified by stress multi-peak skewness to obtain local contact stress multi-peak skewness data;
[0107] Based on the multi-peak deflection data of local contact stress and the gear temperature rise equivalent mapping data, the thermal stress cyclic strength geometric increment analysis is performed to obtain the thermal stress cyclic strength geometric increment data;
[0108] Perform thermal stress intensity increment time series variance regression analysis on the thermal stress cycle intensity geometric increment data to obtain the strength increment time series variance regression data;
[0109] Based on the strength increment time series variance regression data, the gear material hardness data is simulated for material toughness equal loss to obtain material toughness equal loss data;
[0110] The gear deformation stiffness loss ratio is deduced based on the strength increment time series variance regression data and the material toughness equivalent loss data to obtain the gear deformation stiffness loss ratio.
[0111] In an embodiment of the present invention, when identifying stress multimodal skewness based on the variance of the local contact stress concentration distribution of gears, local stress distribution data must first be extracted from the contact area of the meshing tooth surfaces. This data is derived from the two-dimensional nodal stress matrix established in step S242. In a specific implementation, taking a cylindrical gear with a tooth width of 90 mm and a meshing surface width of 85 mm as an example, the meshing surface is divided into 3825 rectangular cells of 2 mm × 1 mm, each corresponding to an average contact stress value, forming a one-dimensional stress vector. Using kurtosis, skewness, and a probability density function fitting method based on kernel density estimation, the presence of two or more main peaks in the vector is identified, and the main peak positions and the degree of distribution skewness are further evaluated. If the local minimum distance between the detected main peaks is less than 1.5 times the local stress standard deviation and the skewness is greater than 0.3, the distribution is determined to be multimodal. In the experimental samples, the main peaks were observed at 1120 MPa and 1370 MPa, respectively, with a skewness of 0.41, confirming a bimodal distribution. The corresponding positions and intensities were recorded as local contact stress multimodal skewness data. In the geometric incremental analysis of thermal stress cyclic strength based on the local contact stress multimodal skewness data and gear temperature rise equivalent mapping data, the multimodal principal values and their intermittent frequencies were used as inputs to establish a joint input space with the aforementioned meshing impact temperature rise data to analyze the fatigue strength variation under combined thermal-mechanical effects. The specific method was to set the contact cycle frequency corresponding to the principal stress peak (e.g., 1370 MPa) to 38 Hz, corresponding to a temperature rise of 27.3°C. A geometric function was introduced to construct a sequential increment factor, and the corresponding decrease in fatigue strength threshold was calculated at each incremental step. Using metal fatigue fracture test data, a geometric relationship curve was derived based on the cycle life decrease interval corresponding to each 10 MPa load increase. It was found that under the current load and temperature rise conditions, the fatigue limit geometrically decreases by approximately 0.87 for each doubling of the stress in the tooth surface thermal stress cyclic strength. Multiple sets of geometric increment data tables were constructed for different stress principal value and temperature rise matches. During the time-series variance regression analysis of the thermal stress cyclic strength geometric increment data, the geometric increment data was converted into a time series format, and the magnitude of thermal stress changes within different time periods was statistically analyzed. Specifically, the simulation run time was set to 30 minutes, divided into 1800 1-second cycles, and the actual number and magnitude of thermal stress cycles were interpolated within each cycle. The mean and variance of the strength increments were then statistically analyzed using a sliding window (60-second window width, 10-second step size), and a third-order polynomial regression was used to analyze the time-series variance trend. The sum of squared residuals in the regression function was used as the regression goodness-of-fit criterion, with an R² value exceeding 0.92 considered valid.In the experimental data, stress variance increased rapidly during the initial loading phase and reached saturation after 8 minutes. Fitting results showed a third-order coefficient of 1.2e-6, indicating cumulative and nonlinear stress variations. The final output, the strength increment time-series variance regression data, includes the fitting coefficient, residual variance, and mean deviation for each time period. When simulating the equivalent toughness loss of gear material hardness data based on the strength increment time-series variance regression data, the stress time-series variance regression results must first be mapped to a microstructural change model of the metal material to derive the damage effect of stress perturbations on material toughness. For example, 18CrNiMo7-6 steel, with a room-temperature impact energy absorption of 42 J, was observed after high-temperature, high-frequency loading. Based on the microstructural evolution after thermo-mechanical fatigue in the material database, a toughness reduction rate of 0.3% was set for each 10 MPa thermal shock load applied. Using the regression variance coefficient and the cyclic stress fluctuation intensity as mapping inputs, the toughness change for each stress increment interval was integrated and accumulated to form a complete equivalent toughness loss curve. Under conditions of a peak stress variance of 185 MPa and an average stress rise rate of 1.7 MPa / s, the measured toughness decreased by approximately 13% after 30 minutes, resulting in an equivalent toughness loss of 5.46 J. The corresponding loss ratio was recorded in the data set. In the process of deducing the gear deformation stiffness loss ratio based on the strength increment time-series variance regression data and the material toughness equivalent loss data, it is necessary to convert the toughness change and stress perturbation trend into a dynamic degradation factor of the material stiffness. Assuming a positive correlation between material stiffness and impact toughness, and referring to the fatigue microcrack initiation threshold and the local plastic zone extension criterion, a stiffness reduction coefficient function was constructed. The input variables were the material impact absorbed energy, the slope of the regression curve, and the main peak stress difference. By comparing the elastic recovery during single-tooth meshing with the unit load energy consumption, the reduction in the load required for unit deformation of the gear after continuous loading was evaluated. In the specific deduction, using an initial toughness of 42 J and degradation to 36.5 J as input, combined with a variance regression slope of 0.021 MPa / s², the overall stiffness reduction rate was calculated to be 0.163, indicating a gear deformation stiffness loss ratio of 0.163. This ratio will be used as a key adjustment factor in the transmission system's mechanical equilibrium equations in the subsequent fatigue load simulation.
[0112] Step S3 includes the following steps:
[0113] Step S31: performing convolution processing on the gearbox fatigue load quantization data to obtain gearbox fatigue load convolution data;
[0114] Step S32: performing fatigue load feature learning on the gearbox fatigue load convolution data to obtain fatigue load feature learning data;
[0115] Step S33: constructing a gearbox fatigue analysis model based on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model;
[0116] Step S34: Send the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0117] As an example of the present invention, refer to Figure 3 As shown, in this example, step S3 includes:
[0118] Step S31: performing convolution processing on the gearbox fatigue load quantization data to obtain gearbox fatigue load convolution data;
[0119] In this embodiment of the present invention, during the convolution process of the quantized gearbox fatigue load data, the quantized gearbox fatigue load data obtained in step S24 is first converted into a continuous time series format, with time as the horizontal axis and the unit fatigue load value as the vertical axis, forming a sequence of equally spaced samples. A sliding window with a fixed length of 128 is used to perform one-dimensional convolution on the load sequence. Each window corresponds to a sampling segment. The convolution kernel is a real vector with a length of 5 and a step size of 1. The convolution operation uses an unpadded boundary processing method, and the local load response intensity corresponding to each window is output. After the convolution kernel weights are initialized, the maximum response amplitude is used as the normalization factor, so that the convolution result reflects the strength comparison of the fatigue response in each interval. In the experiment, the total length of the sample sequence is 36,000 time points, each corresponding to 1 Hz sampling. Six different convolution kernels are used for parallel convolution to generate six feature channels, each with a length of 35,996, ultimately forming a fatigue load convolution data tensor with six-dimensional features.
[0120] Step S32: performing fatigue load feature learning on the gearbox fatigue load convolution data to obtain fatigue load feature learning data;
[0121] In an embodiment of the present invention, when fatigue load feature learning is performed on the gearbox fatigue load convolution data, a method based on Fourier feature extraction combined with multi-scale gradient analysis is adopted. First, the convolution response signal of each channel is subjected to a fast Fourier transform, and the amplitudes of the first 10 main frequency components are extracted as frequency domain fatigue features. Subsequently, the local fluctuation rate is calculated using a high-order difference method, and the third-order central difference is obtained within every 100 time points and the average value is taken to obtain the local gradient change trend on each channel. The cumulative frequency is then statistically calculated based on the maximum response point position information after convolution to generate a 22-dimensional feature vector set including frequency domain features (10 dimensions), time domain fluctuation features (6 dimensions) and response peak frequency features (6 dimensions). After all feature data are Z-score standardized, a fatigue load feature learning dataset is constructed. The total number of data samples is consistent with the number of convolution windows, and each sample contains 22-dimensional features within a corresponding time period.
[0122] Step S33: constructing a gearbox fatigue analysis model based on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model;
[0123] In an embodiment of the present invention, when constructing a gearbox fatigue analysis model based on fatigue load feature learning data using a decision tree algorithm, fatigue data samples are first labeled with fatigue levels. Based on fatigue test results, the fatigue response in each time window is divided into five levels (L1 to L5, from low to high). This labeling criteria is based on the tooth surface spalling depth and the load-bearing duration under unit load in actual gear fatigue life tests. In the experiments, samples with a depth exceeding 0.25 mm and a duration of less than 800 seconds were classified as L5. An ID3 decision tree construction method was employed, using entropy increase as the partitioning criterion. The features with the highest information gain were progressively selected from the 22-dimensional features for node partitioning, constructing a multi-branch tree structure with a depth of no more than 6. The training set consisted of 28,000 samples and the test set consisted of 8,000. After training, a decision tree model consisting of 43 decision nodes and 12 leaf nodes was output. The model structure included feature indexes for each node layer, partitioning thresholds, and sample distribution statistics. The constructed gearbox fatigue analysis model achieved a classification accuracy of 87.3% in the test set, with the primary classification error concentrated at the boundary between the L3 and L4 levels.
[0124] Step S34: Send the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0125] In an embodiment of the present invention, in the process of sending the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method, the model structure file and the standardized parameter file are packaged into a unified configuration file using a binary compression method. The model structure file includes a decision tree node structure table and a feature index mapping table, and the standardized parameter file contains the mean and standard deviation data of each feature. The packaged model data is sent to the on-site monitoring terminal device via the TCP communication protocol. The terminal configuration has the function of real-time collection of vibration signals and torque data, and can automatically execute the constructed decision tree model to judge the fatigue level of the current sampled data. In the actual industrial deployment scenario, the model file size is 1.12MB, the deployment cycle is less than 30 seconds, and in the on-site terminal verification, it has the processing capacity of processing no less than 10 sets of time series data per second, ensuring the stable operation of the gearbox fatigue analysis method in the online monitoring system.
[0126] The present invention also provides a gearbox fatigue analysis system for executing the gearbox fatigue analysis method described above. The gearbox fatigue analysis system comprises:
[0127] The abnormal vibration frequency domain analysis module is used to monitor the vibration operation status data of the wind turbine gearbox through a vibration sensor to obtain the gear vibration operation status data; perform abnormal vibration state extraction on the gear vibration operation status data to obtain abnormal vibration state data, and then perform abnormal vibration frequency domain structure analysis to obtain abnormal vibration frequency domain structure data;
[0128] The fatigue load simulation and quantification module is used to perform fractional frequency mutation intensity analysis on abnormal vibration frequency domain structure data to obtain sideband mutation intensity data; based on the sideband mutation intensity data, the gearbox fatigue load is simulated and quantified to obtain gearbox fatigue load quantification data;
[0129] The fatigue analysis model construction module is used to perform fatigue load feature learning on the gearbox fatigue load quantification data to obtain fatigue load feature learning data; construct a gearbox fatigue analysis model on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; and send the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
[0130] The foregoing description is intended only to provide specific embodiments of the present invention, which will enable those skilled in the art to understand and implement the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not intended to be limited to the embodiments shown herein, but is to be construed in the widest possible manner consistent with the principles and novel features disclosed herein.
Claims
1. A gearbox fatigue analysis method, characterized in that: The following steps are involved: Step S1: monitoring the vibration operating state data of the wind turbine gearbox through a vibration sensor to obtain gear vibration operating state data; extracting abnormal vibration states from the gear vibration operating state data to obtain abnormal vibration state data, and then performing abnormal vibration frequency domain structure analysis to obtain abnormal vibration frequency domain structure data; Step S2: performing fractional frequency multiplication mutation intensity analysis on the abnormal vibration frequency domain structure data to obtain sideband mutation intensity data; performing gearbox fatigue load simulation and quantification based on the sideband mutation intensity data to obtain gearbox fatigue load quantification data; Step S3: performing fatigue load feature learning on the gearbox fatigue load quantification data to obtain fatigue load feature learning data; A gearbox fatigue analysis model is constructed based on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; the gearbox fatigue analysis model is sent to a terminal to execute a gearbox fatigue analysis method.
2. The gearbox fatigue analysis method according to claim 1, characterized in that: Step S1 includes the following steps: Step S11: monitoring the vibration operating state data of the wind turbine gearbox through a vibration sensor to obtain gear vibration operating state data; Step S12: Filling missing values in the gear vibration operation state data to obtain vibration operation state filled data; Step S13: extracting abnormal vibration state from the vibration operation state filling data to obtain abnormal vibration state data; Step S14: performing abnormal vibration frequency domain structure analysis on the abnormal vibration state data to obtain abnormal vibration frequency domain structure data.
3. The gearbox fatigue analysis method according to claim 1, characterized in that: Step S2 includes the following steps: Step S21: performing sideband phase intensity change analysis on the abnormal vibration frequency domain structure data to obtain sideband phase intensity change data; Step S22: performing fractional frequency mutation intensity analysis on the sideband phase intensity change data to obtain sideband mutation intensity data; Step S23: Deducing the gear meshing offset impact load increment index based on the sideband mutation intensity data to obtain the meshing offset impact load increment index; Step S24: performing gearbox fatigue load simulation and quantification based on the meshing offset impact load increment index to obtain gearbox fatigue load quantification data.
4. The gearbox fatigue analysis method according to claim 3, characterized in that: Step S22 includes the following steps: Step S221: calculating the adjacent mean difference of sideband phase frequency mutations on the sideband phase intensity change data to obtain the adjacent mean difference of phase frequency mutations; Step S222: performing abnormal narrowband peak nonlinear regression analysis on the sideband phase intensity change data according to the adjacent mean difference of phase frequency mutations to obtain abnormal narrowband peak nonlinear regression data; Step S223: performing kurtosis increment distribution analysis on the abnormal narrow-band peak nonlinear regression data to obtain abnormal kurtosis increment distribution data; Step S224: performing sideband abnormal frequency kurtosis fractional factorial difference calculation on the abnormal kurtosis increment distribution data to obtain the abnormal frequency kurtosis fractional factorial difference; Step S225: performing fractional octave mutation intensity analysis based on the abnormal frequency kurtosis fractional factorial difference and the abnormal narrowband peak nonlinear regression data to obtain sideband mutation intensity data.
5. The gearbox fatigue analysis method according to claim 3, characterized in that: Step S23 includes the following steps: Step S231: performing gear meshing torsional vibration frequency instability mapping based on the sideband mutation intensity data to obtain gear meshing torsional vibration frequency instability mapping data; Step S232: Deducing the progressive offset of the torsional vibration axial angular displacement according to the gear meshing torsional vibration frequency instability mapping data to obtain progressive offset data of the torsional vibration axial angular displacement; Step S233: performing offset impact energy distribution calculation based on the progressive offset data of the torsional vibration axial angular displacement and the gear meshing torsional vibration frequency instability mapping data to obtain offset impact energy distribution calculation data; Step S234: performing gear meshing offset impact load increment index deduction on the offset impact energy distribution calculation data to obtain the meshing offset impact load increment index.
6. The gearbox fatigue analysis method according to claim 3, characterized in that: Step S24 includes the following steps: Step S241: acquiring basic attribute data of the gears in the gearbox, wherein the basic attribute data of the gears include gear pressure angle, gear center distance and gear material hardness data; Step S242: performing pressure angle and center distance offset numerical mapping measurement on the gear pressure angle and gear center distance in the basic attribute data of the gear according to the meshing offset impact load increment index, and performing simulation and calculation of the variance of the local contact stress concentration distribution of the impact load to obtain the variance of the local contact stress concentration distribution of the gear; Step S243: performing gear temperature rise equivalent mapping matching based on the meshing offset impact load increment index to obtain gear temperature rise equivalent mapping data; Step S244: Deducing the gear deformation stiffness loss ratio from the gear material hardness data based on the gear local contact stress concentration distribution variance and gear temperature rise mapping data to obtain the gear deformation stiffness loss ratio; Step S245: performing gearbox fatigue load simulation and quantification based on the gear deformation stiffness loss ratio to obtain gearbox fatigue load quantification data.
7. The gearbox fatigue analysis method according to claim 6, characterized in that: Step S244 includes the following steps: The gear local contact stress concentration distribution variance is identified by stress multi-peak skewness to obtain local contact stress multi-peak skewness data; Based on the multi-peak deflection data of local contact stress and the gear temperature rise equivalent mapping data, the thermal stress cyclic strength geometric increment analysis is performed to obtain the thermal stress cyclic strength geometric increment data; The thermal stress intensity increment time series variance regression analysis is performed on the thermal stress cyclic intensity geometric increment data to obtain the intensity increment time series variance regression data; Based on the strength increment time series variance regression data, the gear material hardness data is simulated for material toughness equal loss to obtain material toughness equal loss data; The gear deformation stiffness loss ratio is deduced based on the strength increment time series variance regression data and the material toughness equivalent loss data to obtain the gear deformation stiffness loss ratio.
8. The gearbox fatigue analysis method according to claim 1, characterized in that: Step S3 includes the following steps: Step S31: performing convolution processing on the gearbox fatigue load quantization data to obtain gearbox fatigue load convolution data; Step S32: performing fatigue load feature learning on the gearbox fatigue load convolution data to obtain fatigue load feature learning data; Step S33: constructing a gearbox fatigue analysis model based on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; Step S34: Send the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
9. A gearbox fatigue analysis system, characterized in that: Used to perform the gearbox fatigue analysis method according to any one of claims 1 to 8, the gearbox fatigue analysis system comprises: The abnormal vibration frequency domain analysis module is used to monitor the vibration operation status data of the wind turbine gearbox through a vibration sensor to obtain the gear vibration operation status data; perform abnormal vibration state extraction on the gear vibration operation status data to obtain abnormal vibration state data, and then perform abnormal vibration frequency domain structure analysis to obtain abnormal vibration frequency domain structure data; The fatigue load simulation and quantification module is used to perform fractional frequency mutation intensity analysis on abnormal vibration frequency domain structure data to obtain sideband mutation intensity data; based on the sideband mutation intensity data, the gearbox fatigue load is simulated and quantified to obtain gearbox fatigue load quantification data; The fatigue analysis model construction module is used to perform fatigue load feature learning on the gearbox fatigue load quantification data to obtain fatigue load feature learning data; construct a gearbox fatigue analysis model on the fatigue load feature learning data based on a decision tree algorithm to obtain a gearbox fatigue analysis model; and send the gearbox fatigue analysis model to the terminal to execute the gearbox fatigue analysis method.
Citation Information
Patent Citations
Vibration monitoring-based wind generator set automatic fault diagnosis method
CN101858778A
Composite fault diagnosis method and system of gear case
CN102937522A