Fan fault early warning method based on number and object fusion

Abstract

The invention relates to the technical field of fan fault early warning, and provides a fan fault early warning method based on number-object fusion. The fan fault early warning method comprises the steps of collecting multi-source data; calculating the maximum contact stress of the bearing according to the high-precision mechanism model of the double-row self-aligning roller bearing; carrying out pre-processing on the multi-source data; establishing a PEAKR model, and inputting the processed multi-source data and the maximum contact stress data into the PEAKR model for training; using the PEAKR model to predict the operation data to be monitored, generating a prediction value and calculating a residual sequence, and determining an MSE index according to the residual sequence; adopting an MSE index based on a moving window strategy to analyze the output residual error of the PEAKR model; and when the MSE index exceeds an MSE alarm threshold, indicating that the health state of the equipment is abnormal, and triggering an alarm. According to the invention, accurate identification of potential faults in equipment state monitoring can be realized, and safe, stable and efficient operation of machinery equipment is ensured to the greatest extent.

Description

Technical Field

[0001] The present invention relates to the technical field of wind turbine fault warning, and in particular to a wind turbine fault warning method based on digital-physical fusion. Background Art

[0002] Wind turbines, as essential renewable energy devices, are widely used worldwide. With the growing demand for green energy, wind power generation has gradually become a key alternative to traditional energy sources. However, the complex structure and long-term continuous operation of wind turbines make them prone to failure during actual operation, which in turn causes unplanned downtime. Unplanned downtime of wind turbines not only results in huge economic losses, but can also cause serious safety accidents, endangering the safety of personnel and equipment. Studies have shown that losses caused by equipment downtime in the wind power industry have accounted for a considerable portion of the total cost of wind power generation. Therefore, how to reduce these risks through effective condition monitoring and intelligent maintenance has become a key issue that needs to be addressed in the wind power field.

[0003] In recent years, data-driven approaches have become a mainstream focus in fault early warning research due to their high flexibility and predictive accuracy. By deeply mining wind turbine operating data, data-driven approaches can identify abnormal conditions in key components and predict failure trends. However, these approaches rely heavily on the availability and completeness of high-quality data. In the absence of sufficient historical fault samples or when faced with unknown failure modes, the generalization and robustness of the models remain significant challenges, potentially leading to the risk of false positives or missed detections.

[0004] Although data-driven approaches have shown great potential in wind turbine fault early warning, their application is still subject to multiple limitations. First, the operating environment of wind turbines is complex and changeable, significantly affected by environmental factors such as wind speed, temperature, and humidity, which often leads to operating data with strong non-stationary and nonlinear characteristics. Second, because data collection is constrained by conditions such as sensor layout and network transmission, there are problems such as incomplete data collection, severe noise pollution, and uneven sample distribution. More importantly, while relying solely on data-driven models can extract fault characteristics, they lack the ability to physically explain the causes of faults, making it difficult to effectively guide engineering practice. Summary of the Invention

[0005] The present invention mainly solves the technical problems existing in the above-mentioned existing technology of using data-driven methods for wind turbine fault early warning, and proposes a wind turbine fault early warning method based on data-physics fusion. It adopts an enhanced autoassociative kernel regression method combined with physical mechanism and data-driven (PEAKR), introduces the physical mechanism model of the wind turbine main shaft bearing, and combines it with multi-source data for deep fusion. It adopts the mean square error (MSE) evaluation method based on the moving window strategy to achieve accurate identification of potential faults in equipment status monitoring, and ensure the safe, stable and efficient operation of machinery and equipment to the greatest extent.

[0006] The present invention provides a fan fault early warning method based on digital-physical fusion, comprising:

[0007] Step 1: Collect multi-source data;

[0008] Step 2: Calculate the maximum contact stress of the bearing based on the high-precision mechanism model of the double-row spherical roller bearing;

[0009] Step 3: Pre-process the multi-source data;

[0010] Step 4: Establish a PEAKR model, and input the processed multi-source data and maximum contact stress data into the PEAKR model for training to obtain a trained EAKR model;

[0011] Step 5: Use the PEAKR model to predict the operating data to be monitored, generate predicted values ​​and calculate the residual sequence, and determine the MSE index based on the residual sequence for subsequent equipment status evaluation;

[0012] In step 6, the MSE indicator based on the moving window strategy is used to analyze the residual output of the PEAKR model. When the MSE indicator exceeds the MSE alarm threshold, it indicates that the health status of the device is abnormal and an alarm is triggered.

[0013] Furthermore, in step 2, the high-precision mechanism model of the double-row spherical roller bearing includes a load calculation module, an SRBs dynamic mechanism model, and a Hertz contact theory module; first, the load calculation module calculates the bearing load, and the bearing load is input into the SRBs dynamic mechanism model to obtain the rolling element support load, and then the maximum contact stress is calculated according to the Hertz contact theory.

[0014] Furthermore, the load calculation module includes hub load and bearing load;

[0015] The hub load includes M x 、M y 、M z 、F x 、F y and F z six independent variables;

[0016] The thrust and drag elements of the wind rotor can be expressed as:

[0017]

[0018] where L and D are the lift, W is the resultant velocity, c and C L are the airfoil chord length and airfoil lift coefficient, C D is the drag coefficient;

[0019] According to the momentum theorem, the thrust and torque acting on the element of the rotor (r ，r +dr) are:

[0020]

[0021] The energy equation is derived by combining the blade element theory and the momentum theory:

[0022]

[0023] The induced factors a and a' are obtained by solving the equations (1)-(3) simultaneously, and the hub load can be calculated.

[0024] The bearing load is equivalent to three mutually perpendicular concentrated loads of F B_x , F B_y , and F B_z ;

[0025] The axial load is entirely borne by the main bearing, i.e., F x =F B_z ;

[0026]

[0027] Further, the SRBs dynamic mechanism model includes:

[0028] The geometric parameters of the double-row self-aligning roller bearing, R o represents the profile radius of the outer ring of the bearing, R i represents the profile radius of the inner ring of the bearing, D p is the displacement of the center of the rolling element to the center of the bearing, and a is the nominal contact angle of the rolling element, D o represents the diameter of the rolling element; the inner ring of the bearing has 3 degrees of freedom, and the outer ring is fixed; the total kinetic energy can be calculated as:

[0029]

[0030] where x in , y in , and z in are the displacements of the inner ring / axle in the x, y, and z directions, respectively, z j is the radial displacement of the jth rolling element, and yj is the angular position of the jth roller, m j is the mass of the jth roller, I j is the moment of inertia of the jth roller, y j is the angular displacement of the jth roller, I s and y s are the moment of inertia and angular displacement of the shaft, N r is the number of rollers, α n is the nominal contact angle;

[0031] Introduced operator (-1) j , the total potential energy is calculated as:

[0032]

[0033] In formulas (8)-(10), δ in+ and δ out+ They represent the contact deformation between the jth roller and the inner and outer rings respectively; the contact type of the spherical roller bearing is point contact, and the contact deformation is calculated as follows:

[0034]

[0035] In the above formula, θ j represents the angular position of the jth roller, R ct_inner and R ct_outer Respectively represent the contour radius of the inner and outer rings, α n is the nominal contact angle, j is the roller number; θ j Calculated as:

[0036]

[0037] The inner radius can be expressed as:

[0038] R i =k i0 +κ ip sin(N i θ ij ) (14)

[0039] The inner ring contact angle θ ij It can be expressed as:

[0040]

[0041] The bearing outer ring raceway radius can be expressed as:

[0042] R o =κ o0 +κ op sin(N o θ oj) (16)

[0043] θ oj It can be written as:

[0044]

[0045] According to the Langrange equations, the motion equations of the rolling element in the contact normal direction and the motion equations of the bearing inner ring in the x, y and z directions are:

[0046]

[0047] Where, is the total damping force borne by the jth roller, which is the sum of the damping forces of the inner and outer rings;

[0048] F ax , F ay , F az It is the external force applied to the inner ring of the bearing in three directions, which can be static or dynamic; is the damping force generated in the x, y and z directions when the bearing inner ring contacts the rolling element; R ct_inner ,R ct_outer and R Roller are the contour radii of the inner and outer rings and rolling elements of the bearing respectively; Formulas (18) to (21) are second-order differential equations with a total degree of freedom of N r +3, where N r is the number of rolling elements.

[0049] Furthermore, the Hertz contact theory module includes:

[0050] The distribution of contact pressure on the contact ellipse is as follows:

[0051]

[0052] In Hertz theory, the maximum contact stress occurs at the geometric center of the contact surface, and its formula is:

[0053]

[0054] Where E1, E2, v1, and v2 represent the elastic modulus and Poisson's ratio of the two elastic bodies, respectively. The ∑ρ formula is:

[0055]

[0056] Parameter a * and b * The calculation formulas for the variables are:

[0057]

[0058] where, κ = a / b, e = 1-1 / κ 2 , K(e) and E(e) are the first and second kind complete elliptic integrals, respectively. In the equation, κ is determined by the following equation:

[0059]

[0060] where F(ρ) is the curvature difference of the two contacting elasticities:

[0061]

[0062] According to the Hertz theory, the relative displacement of the two elasticities is:

[0063]

[0064] For the two elasticities in point contact, the relationship between the load and the deformation is as follows:

[0065] Q = K p δ 3 / 2 (32)

[0066] Finally, the contact stiffness of the two elasticities can be derived from the formula (30) and (32) as follows:

[0067]

[0068] Further, the step 2 includes steps 201 to 207:

[0069] Step 201, hub load data calculation: the hub load data is calculated from the wind speed monitored on site of the wind farm wind turbine, to calculate the bearing load required subsequently;

[0070] Step 202, bearing load calculation: the load applied on the bearing under the action of the hub load is calculated through the static equilibrium condition of the simply supported beam, as the input of the bearing dynamics model;

[0071] Step 203, contact stiffness calculation: the contact stiffness between the rolling body and the inner and outer rings inside the double-row self-aligning roller bearing is calculated according to the Hertz contact theory, and the load-deformation relationship is introduced into the SRBs dynamics mechanism model;

[0072] Step 204, parameter analysis: the bearing load data and the contact stiffness obtained from the steps 201 to 203 are substituted into the SRBs dynamics mechanism model;

[0073] Step 205, support load calculation: the rolling body support load under the given bearing load is calculated based on the bearing dynamics model, for the subsequent mechanical analysis;

[0074] Step 206, contact stress calculation: Based on the support load calculation results, the contact stress is calculated using the Hertz contact formula as the final output of the model;

[0075] Step 207, model verification: Calculate the maximum contact stress under 16 extreme working conditions and compare it with the SKF calculation report results to verify the accuracy of the model;

[0076] The maximum contact stress of the bearing is calculated according to the method from step 201 to step 206.

[0077] Furthermore, step 4 includes the following steps 401 to 405:

[0078] Step 401: Divide the multi-source data and the maximum contact stress data calculated in step 2 into a training set, an optimization set, and a test set;

[0079] Step 402, establishing a PEAKR model;

[0080] Extract the historical health data of the equipment operation and create the memory matrix X of the PEAKR model, where X i，j represents the i-th vector value of the j-th key variable; for n memory vectors, the memory matrix X of p variables can be expressed as:

[0081]

[0082] The monitoring vector is represented by a 1×p matrix v:

[0083] v=[v1 v2 … v p ] (35)

[0084] The predicted value of the PEAKR model can be obtained by taking a weighted average of each memory vector of the memory matrix X, where the weighted average parameter is estimated using the unit health data;

[0085] Step 403: training the established EAKR model.

[0086] Furthermore, step 4 includes the following steps 4031 to 4034:

[0087] Step 4031, using the wavelet packet Bayesian noise reduction method to process the training set, optimization set and test set;

[0088] Step 4032, use the Manhattan method to calculate the monitoring vector v and each memory vector X i The distance between them, we get an n×1 distance vector d i ;

[0089] Step 4033: Calculate the weight w using the obtained distance matrix d and the Gaussian kernel function. Each element is calculated using the following formula:

[0090]

[0091] Among them, the weight w is an n×1 vector matrix; the parameter h is the bandwidth of the kernel function;

[0092] Step 4034 , automatically optimizing the kernel bandwidth in the model using a new set of health data and the Nelder-Mead optimization algorithm, and obtaining the optimal bandwidth by minimizing the mean square error during model training;

[0093] The mean square error (MSE) is used to detect the residual between the model prediction value and the true value. The formula for calculating the average mean square error of N test samples is as follows:

[0094]

[0095] Where N is the total number of test samples, p is the total number of model variables, is the true value of the jth variable of the i-th sample, is the predicted value of the jth variable of the model corresponding to the i-th sample.

[0096] Furthermore, after step 403, the method further includes:

[0097] Step 404, evaluating the trained EAKR model;

[0098] Calculate the performance indicators of the PEAKR model using the test set, including R 2 , MSE and MAE, the calculation formula is as follows

[0099]

[0100] Furthermore, in step 5, the predicted value of the monitoring vector v is predicted by the calculated weight w By each memory vector X i The weighted average of is calculated as follows:

[0101]

[0102] The present invention provides a wind turbine fault warning method based on digital-physical fusion. Digital-physical fusion refers to combining data-driven methods with physical mechanism models, making full use of the characterization capabilities of data and the constraints of physical laws to improve the accuracy, generalization and interpretability of modeling. The enhanced auto-associative kernel regression (EAKR) method is combined with physical mechanism and data-driven (PEAKR) to solve the key problems in wind turbine fault prediction. Based on the physical mechanism model of the main shaft bearing of the wind turbine generator, the present invention integrates multi-source data (such as vibration signals and temperature data, etc.), and realizes comprehensive characterization and monitoring of the equipment operating status by constructing a mean square error (MSE) evaluation method based on a moving window strategy. The PEAKR model not only inherits the ability of the data-driven method to process multidimensional nonlinear data, but also improves the robustness of the model and the accuracy of fault warning by embedding physical mechanisms. Through experimental analysis, it is found that compared with the traditional data-driven method, the PEAKR model has significant advantages under complex working conditions such as incomplete data and system nonlinearity, and can effectively identify potential faults of wind turbine generator sets. BRIEF DESCRIPTION OF THE DRAWINGS

[0103] Figure 1 This is a flow chart of the implementation of the fan fault early warning method based on digital-physical fusion provided by the present invention;

[0104] Figure 2a This is a schematic diagram of a fan drive shaft with a high-precision mechanism model of a double-row spherical roller bearing;

[0105] Figure 2b It is a schematic diagram of the high-precision mechanism model of a double-row spherical roller bearing;

[0106] Figure 3a It is a simplified diagram of the fan drive chain load in the xy plane;

[0107] Figure 3b It is a simplified diagram of the fan drive chain load in the xz plane;

[0108] Figure 4 It is a schematic diagram of the geometric shape of the spherical roller bearing;

[0109] Figure 5a It is a schematic diagram of elastic contact in Hertz contact;

[0110] Figure 5b is a schematic diagram of the contact stress ellipse in Hertzian contact. DETAILED DESCRIPTION

[0111] To make the technical problems solved, the technical solutions adopted, and the technical effects achieved by the present invention more clearly apparent, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the present invention and are not intended to limit the present invention. It should also be noted that, for ease of description, the accompanying drawings only illustrate portions relevant to the present invention, rather than all of the contents.

[0112] like Figure 1 As shown, an embodiment of the present invention provides a wind turbine fault early warning method based on digital-physical fusion, including:

[0113] Step 1: Collect multi-source data.

[0114] The multi-source data in the present invention includes two different types of data.

[0115] The multi-source data includes SCADA data and CMS data; the SCADA data is real-time monitoring data collected by the SCADA system, and the SCADA data includes wind speed; the CMS data is vibration data provided by the CMS system.

[0116] Step 2: Calculate the maximum contact stress of the bearing based on the high-precision mechanism model of the double-row spherical roller bearing.

[0117] The high-precision mechanism model for double-row spherical roller bearings includes a load calculation module, an SRBs dynamic mechanism model, and a Hertz contact theory module. The load calculation module first calculates the bearing load, which is then input into the SRBs dynamic mechanism model to determine the rolling element support load. The maximum contact stress is then calculated based on Hertz contact theory.

[0118] 1. Load calculation module

[0119] In the dynamic mechanism model of the double-row spherical roller bearing, the input load cannot usually be obtained directly and needs to be indirectly derived through information such as the hub load and the equipment's own gravity. Based on wind speed and other data monitored by the wind field and combined with blade element momentum theory, the corresponding hub load can be calculated. The hub load contains six independent variables, namely, M x 、M y 、M z 、F x 、F y and F z (like Figure 2a These hub loads and the weight of the equipment are transmitted to the main shaft bearings via the wind turbine main shaft.

[0120] 1. Hub load

[0121] The hub load is calculated using the AeroDyn module of FAST, which primarily uses the blade element momentum theory. The blade element theory was first proposed by Richard Fronde in 1889 and subsequently improved by S. Drzewiecki. This theory uses the blade element as the research object, analyzes the aerodynamic lift and torque on the blade element airfoil, and integrates them along the span to ultimately calculate the rotor torque or wind turbine power. The blade is divided into multiple microelements along the span direction, and the theory analyzes the blade forces through the flow around the blade element. The thrust element and drag element of the wind rotor can be expressed as

[0122]

[0123] Where L and D are lift (N), W is the resultant velocity (m / s), c and C L is the airfoil chord length (m) and the airfoil lift coefficient, C D is the drag coefficient.

[0124] The momentum theory was proposed by William Rankime in 1865. This theory differs from the blade element theory in that it proposes that the force acting on the blade element is related to the momentum change of the airflow through the blade element ring.

[0125] According to the momentum theorem, the thrust and torque acting on the blade ring (r, r+dr) are

[0126]

[0127] The energy equation is derived by combining blade element theory and momentum theory.

[0128]

[0129] It illustrates the relationship between the incoming wind speed and the thrust and torque of the wind turbine. The induction factors a and a′ are obtained by simultaneous solution, and the hub load can be calculated by applying them back to formulas (1)-(3).

[0130] 2. Bearing load

[0131] Since double row spherical roller bearings (SRBs) cannot withstand bending moments, the bearing load can be equivalent to three mutually perpendicular concentrated loads, namely F B_x 、F B_y , and F B_z (like Figure 2b The fan main shaft transmits the hub load and equipment weight to the main shaft bearings.

[0132] The fan drive chain adopts a single bearing and three-point support structure. The left end of the main shaft is constrained by the main shaft, the right end is elastically supported, and the gearbox torque arm is rigidly connected to the right end of the main shaft. The overall structure is statically determinate. In the horizontal and vertical planes, this structure can be regarded as a simply supported beam bearing loads, such as Figure 3a -b as shown.

[0133] Among them, L1 and L5 represent the length between the parts of the main axis, G and G 齿轮箱 Represent the gravity of the main shaft and gearbox respectively. The fan hub load and the main bearing load meet the static equilibrium condition, see formula (5). By solving this set of equations, the bearing load can be obtained. The axial load is all borne by the main shaft, that is, F x =F B_z .

[0134]

[0135] 2. Dynamic mechanism model of SRBs (Spherical Roller Bearings, double-row spherical roller bearings).

[0136] This model is derived using the Lagrange equations. It assumes the outer ring is fixed and the rollers do not slip on the contact surfaces with the inner and outer rings. This paper defines the kinetic and potential energies of a conservative system. The kinetic energy expression takes into account the velocity of the inner ring and all rollers. The potential energy of this system consists of two main components: the first is the gravitational potential energy of the inner ring and rollers; the second is the contact potential energy caused by the contact deformation between the rollers and the inner and outer rings.

[0137] Double row spherical roller bearing geometric parameters see Figure 4 , where R o Indicates the outer ring radius of the bearing, R i Indicates the contour radius of the bearing inner ring, D p is the displacement from the center of the rolling element to the center of the bearing, α is the nominal contact angle of the rolling element, D o represents the rolling element diameter. The inner ring of the bearing has 3 degrees of freedom (assuming that rotation about the horizontal and vertical axes is negligible), while the outer ring is fixed. Therefore, the total kinetic energy can be calculated as:

[0138]

[0139] Among them, x in 、y in 、z in are the displacements of the inner ring / shaft in the x, y and z directions respectively, j is the radial displacement of the jth roller, y j is the angular position of the jth roller, m j is the mass of the jth roller, I j is the moment of inertia of the jth roller, y j is the angular displacement of the jth roller, I s and y s are the moment of inertia and angular displacement of the shaft, N r is the number of rollers, αn is the nominal contact angle.

[0140] As shown in formula (7), the contact angle is now a function of the axial displacement of the inner ring and is no longer considered a constant. Unlike single-row deep groove bearings, the operator (-1) is introduced. j , to reflect that there is a pair of rollers at each evenly distributed cage position, and their relative displacements relative to the inner ring are always opposite. The total potential energy is then calculated as:

[0141]

[0142] In formulas (8)-(10), δ in+ and δ out+ Represents the contact deformation between the jth roller and the inner ring and outer ring respectively. N c Different values ​​are determined according to the contact type (point contact, line contact). The contact type of the spherical roller bearing studied in this paper is point contact, and the contact deformation calculation method is:

[0143]

[0144] In the above formula, θ j represents the angular position of the jth roller, R ct_inner and R ct_outer Represents the contour radius of the inner and outer rings, α n is the nominal contact angle, j is the roller number. j This can be calculated as (note that there are two rollers at each position):

[0145]

[0146] It is worth noting that R i , γ r and R o The value is determined by adding the nominal value to the surface waviness (from the inner ring, rollers and outer ring respectively).

[0147] Assume that the surface waviness of the inner / outer rings and rollers can be expressed as a sinusoidal function. The amplitude and number of waves will vary depending on the manufacturing process. Therefore, the inner radius can be expressed as:

[0148] R i =κ i0 +κ ip sin(N i θ ij ) (14)

[0149] The inner ring contact angle θ ij It can be expressed as:

[0150]

[0151] Assuming the surface profile of the two rows of raceways of the inner ring is the same, where N r represents the number of undulations of the raceway of the inner ring. Similar to the inner ring, the radius of the raceway of the outer ring can be represented as:

[0152] R o = κ o0 + κ op sin(N o θ oj ) (16)

[0153] θ oj can be written as:

[0154]

[0155] According to the Lagrange equation, the motion equation of the rolling element in the direction of the contact normal and the motion equation of the inner ring in the x, y, and z directions are respectively:

[0156]

[0157]

[0158] where, is the total damping force borne by the jth rolling element, and is the sum of the damping forces of the inner and outer rings.

[0159] The motion of the rolling element is represented by the elastic force and the damping force generated when it contacts the inner and outer rings of the bearing. F ax , F ay , F az are external forces applied to the inner ring in three directions, which can be static forces or dynamic forces; are the damping forces generated in the x, y, and z directions when the inner ring of the bearing contacts the rolling element; R ct_inner , R ct_outer , and R Roller are the profile radii of the inner and outer rings and the rolling element. Equations (18) to (21) are second-order differential equations with a total degree of freedom of N r + 3, where N r is the number of rolling elements.

[0160] III. Hertz Contact Theory Module

[0161] In 1882, Hertz proposed the Hertz contact theory for studying the interaction between two elastic bodies. As shown in Figure 5a , there are two elastic bodies I and II, and their radii in two orthogonal planes 1 and 2 are ρ I1 , ρ I2 , ρ II1 , and ρII2 Hertz believed that whether it is point contact or line contact, the applied load will be distributed over a certain area to prevent infinite stress. Figure 5b As shown, σ is the contact stress, the contact area is elliptical, and the maximum contact pressure occurs at the center of the ellipse. Hertz contact theory is based on the following assumptions:

[0162] 1. The deformation is within the elastic range and the material will not exceed the proportional limit;

[0163] 2. The load is perpendicular to the contact surface, and the influence of surface shear stress is not considered;

[0164] 3. The contact area is very small relative to the curvature radius of the elastic body.

[0165] The bearing contact problem conforms to the Hertz contact hypothesis mentioned above. The contact between the bearing rolling element and the raceway is a typical Hertz contact problem. The Hertz formula is a classic solution derived from the analysis of elastic mechanics theory and is often used in finite element analysis of elastic bodies. The simulation results of the fatigue life prediction method of the wind turbine main bearing based on the surrogate model are compared. This section calculates the contact stiffness and stress distribution between the bearing rolling element and the inner and outer ring raceways of the bearing based on the Hertz theoretical formula.

[0166] like Figure 5b As shown, the distribution of contact pressure on the contact elliptical surface is as follows:

[0167]

[0168] In Hertz theory, the maximum contact stress occurs at the geometric center of the contact surface, and its formula is:

[0169]

[0170] Where E1, E2, v1, and v2 represent the elastic modulus and Poisson's ratio of the two elastic bodies, respectively. The ∑ρ formula is:

[0171]

[0172] Parameter a * and b * The calculation formulas for the variables are:

[0173]

[0174] Among them, κ=a / b, e=1-1 / κ 2 , K(e) and E(e) are the complete elliptic integrals of the first and second kinds, respectively. Where κ is determined by the following equation:

[0175]

[0176] Where F(ρ) is the difference in curvature between the two contact elastic bodies:

[0177]

[0178] According to Hertz's theory, the relative displacement of the two elastic bodies is:

[0179]

[0180] For two elastic bodies in point contact, the relationship between load and deformation is as follows:

[0181] Q=K p δ 3 / 2(32)(40)

[0182] Finally, the contact stiffness of the two elastic bodies can be derived from formulas (30) and (32):

[0183]

[0184] Step 2 includes steps 201 to 207:

[0185] Step 201 , hub load data calculation: hub load data is calculated from the wind speed monitored on-site at the wind turbine generator system in the wind farm, so as to calculate the bearing load required subsequently.

[0186] Among them, the wind speed is given by SCADA data.

[0187] Step 202, bearing load calculation: The load applied to the bearing under the action of the hub load is calculated by using the static equilibrium condition of the simply supported beam, and is used as the input of the bearing dynamics model.

[0188] Step 203, contact stiffness calculation: Based on Hertz contact theory, the contact stiffness between the inner rolling elements and the inner and outer rings of the double-row spherical roller bearing is calculated, and the load-deformation relationship is introduced into the SRBs dynamic mechanism model.

[0189] Step 204, parameter analysis: Substitute the bearing load data and contact stiffness obtained from steps 201 to 203 into the SRBs dynamic mechanism model.

[0190] Step 205, support load calculation: Based on the bearing dynamics model, the rolling element support load under a given bearing load is calculated for subsequent mechanical analysis.

[0191] Step 206, contact stress calculation: Based on the support load calculation results, the contact stress is calculated using the Hertz contact formula as the final output of the model.

[0192] Step 207, model verification: Calculate the maximum contact stress under 16 extreme working conditions and compare it with the SKF calculation report results to verify the accuracy of the model.

[0193] The maximum contact stress of the bearing is calculated according to the method from step 201 to step 206.

[0194] The 16 extreme working conditions are the load components in the hub coordinate system (such as M x 、M y 、M z 、M yz 、F x 、F y 、F z and F yz ) maximum and minimum values. Among them, M yz and F yz Represents M y 、M z With F y 、F z Combined bending moment and combined load.

[0195] The SKF calculation report is a bearing life calculation report.

[0196] Step 3: Pre-process the multi-source data.

[0197] Pre-processing of SCADA data:

[0198] First, the missing values ​​of the data were processed and the median of the reliable data before and after the outlier was used for interpolation.

[0199] Subsequently, the Bayesian wavelet packet threshold denoising method is applied to remove the noise.

[0200] Next, the downsampled data were normalized to zero mean to eliminate the dimension difference.

[0201] Afterwards, the processed data is downsampled to once an hour to reduce the temporal resolution of the data and reduce the computational burden.

[0202] Finally, principal component analysis (PCA) is used to reduce the dimensionality of the denoised data to extract the main features and reduce data redundancy.

[0203] Pre-processing of CMS data:

[0204] Firstly, the Bayesian wavelet packet threshold denoising method is used to denoise the original vibration data.

[0205] Next, the denoised data is downsampled to once every two hours to reduce the temporal resolution of the data and reduce computational complexity.

[0206] Subsequently, key features are extracted based on the time domain feature extraction method, and variance analysis is further performed to select features that contribute more to the system.

[0207] Finally, missing value processing and standardization were performed to ensure the consistency and comparability of the data.

[0208] Step 4: Establish a PEAKR model and input the processed multi-source data and maximum contact stress data into the PEAKR model for training to obtain a trained EAKR model. Step 4 includes the following steps 401 to 405:

[0209] Step 401 : Divide the multi-source data and the maximum contact stress data calculated in step 2 into a training set, an optimization set, and a test set.

[0210] Divide the health data into three parts: training set, optimization set, and test set. The monitoring data is used as the prediction set. The health data should cover the normal operating status of the unit as much as possible.

[0211] The training set is used to generate the PEAKR model memory matrix and should encompass as many normal operating conditions as possible. The optimization set is used to optimize the model kernel bandwidth. After generating the memory matrix, the kernel bandwidth h is optimized using the optimization set and the simplex method to obtain the optimized model. The test set is used to test the model's predictive performance.

[0212] Step 402: Establish a PEAKR model.

[0213] The PEAKR (Enhanced Auto-Associative Kernel Regression) model is a similarity-based, nonparametric empirical modeling technique that uses historical, fault-free data to predict unit performance. This method is independent of the equipment and fault type and is applicable to multivariate operational monitoring and fault early warning for a wide range of equipment.

[0214] Extract the historical health data of the equipment operation and create the memory matrix X of the PEAKR model, where X i，j Represents the i-th vector value of the j-th key variable. For n memory vectors, the memory matrix X of p variables can be expressed as:

[0215]

[0216] The monitoring vector is represented by a 1×p matrix v:

[0217] v=[v1 v2 … v p ] (35)

[0218] The predicted value of the PEAKR model can be obtained by taking a weighted average of each memory vector of the memory matrix X, where the weighted average parameter is estimated using the unit health data.

[0219] Step 403: train the established EAKR model. Step 4 includes the following steps 4031 to 4034:

[0220] Step 4031: Process the training set, optimization set, and test set using the wavelet packet Bayesian denoising method.

[0221] The data collected by the sensor will be interfered by complex nonlinear noise. This study uses the wavelet packet Bayesian denoising method to process the signal. This method can avoid the loss of real information caused by excessive noise reduction and the hiding of fault information caused by insufficient noise reduction, and can improve the accuracy of the original AAKR fault warning method.

[0222] Step 4032, use the Manhattan method to calculate the monitoring vector v and each memory vector X i The distance between them, we get an n×1 distance vector d i .

[0223] Step 4033: Calculate the weight w using the obtained distance matrix d and the Gaussian kernel function. Each element is calculated using the following formula:

[0224]

[0225] The weight w is an n×1 vector matrix; the parameter h is the bandwidth of the kernel function, which determines the smoothness of the function. A smaller h reveals more detail but does not produce a smoother tail, while a larger h loses detail in system identification, leading to more false positives. Therefore, in industrial equipment condition monitoring, the selection of bandwidth h plays a crucial role in the accuracy of fault warnings. To improve the accuracy of model predictions, the simplex method is used in the PEAKR method to optimize h.

[0226] In step 4034, the kernel function bandwidth in the model is automatically optimized using a new set of health data and the Nelder-Mead optimization algorithm. The optimal bandwidth is obtained by minimizing the mean square error (MSE) during model training, thereby reducing the possibility of false alarms in equipment status monitoring.

[0227] The mean square error (MSE) is used to detect the residual between the model prediction value and the true value. The formula for calculating the average mean square error of N test samples is as follows:

[0228]

[0229] where N is the total number of test samples, p is the total number of model variables, is the true value of the jth variable of the ith sample, is the predicted value of the jth variable of the ith sample by the model.

[0230] Step 404, evaluate the trained EAKR model.

[0231] Calculate the performance indicators of the PEAKR model, including R 2 , MSE and MAE, through the test set to evaluate the prediction ability of the optimized model. The calculation formula is as follows

[0232]

[0233] Step 5, use the PEAKR model to predict the monitored operation data, generate predicted values and calculate residual sequences, and determine the MSE index according to the residual sequences for subsequent equipment state evaluation.

[0234] The predicted value of the monitoring vector v is calculated by the calculated weight w is calculated by the weighted average of each memory vector X, and the calculation formula is as follows:

[0235]

[0236] In the process of establishing the PEAKR model, the memory matrix should contain as much unit health data of various scenes as possible.

[0237] For the monitored operation data, use the pre-processing procedure mentioned in the foregoing to operate and denoise. The processed data is input into the established PEAKR model to generate predicted values and calculate residual sequences, and the MSE index is determined according to the residual sequences for subsequent equipment state evaluation. The trained PEAKR model is used for equipment operation monitoring and fault alarm.

[0238] Step 6, use the MSE index based on the moving window strategy to analyze the residual of the PEAKR model output; when the MSE index exceeds the MSE alarm threshold, it indicates that the equipment health state is abnormal, triggering an alarm.

[0239] The MSE (mean square error) indicator is a more convenient way to measure the average error. The MSE can measure the degree of change in the data. The smaller the MSE value, the healthier the unit operation status. The present invention uses the MSE indicator based on the moving window strategy to analyze the residuals of each model output. It should be noted that N is the length of the moving window at this time. When the trained model is used for equipment operation monitoring and fault alarm, if the MSE exceeds a given MSE alarm threshold, it indicates that the health status of the equipment is abnormal and triggers an alarm. The MSE alarm threshold at the time of alarm is usually determined by the model testing phase. In the model testing phase, the test data is the unit health operation data. The calculated MSE value represents the unit health operation status. The alarm threshold setting should not be too small to avoid false alarms. The value calculated in the model testing phase should be Set as fault alarm threshold.

[0240] Taking actual wind turbine operation data as an example, the method proposed in the present invention is verified:

[0241] 1. The data is derived from actual monitoring data from a wind turbine that experienced a main bearing failure on March 27, 2024. The CMS data collection point is located at the vertical vibration position of the front bearing of the main bearing. The sampling frequency is 2560 Hz, and the collection period is from December 7, 2023, to April 11, 2024, with a sampling interval of 2 hours. SCADA data is collected from November 1, 2023, to April 30, 2024, with a sampling interval of 30 seconds. Data for a total of 21 variables are collected.

[0242] 2. The bearing is a double-row spherical roller bearing for the main shaft of the wind turbine. Its model is SKF 240 / 750ECA / W33. The bearing material is 42CrMoA with a density of 7.85×10 3 kg / m 3 , the elastic modulus is 2.1×10 5 Mpa, Poisson's ratio is 0.3. The main geometric parameters of the bearing are shown in Table 1.

[0243] Table 1 Main parameters of SKF 240 / 750ECA / W33 bearings

[0244]

[0245] 3. Using a high-precision mechanism model of a double-row spherical roller bearing, the maximum contact stress corresponding to the wheel hub's extreme operating load was calculated and compared with the results reported by SKF. The results are shown in Table 2. As can be seen, the error between the calculated results of the mechanism model established in this paper and the calculated values ​​reported by SKF is small, with the maximum error being less than 3%. This demonstrates the high reliability and accuracy of the mechanism model developed in this paper.

[0246] Table 2 Comparison of maximum contact pressure results of bearings under extreme working conditions

[0247]

[0248]

[0249] 4. Based on the above verified mechanism model, the SCADA data is calculated to obtain the corresponding maximum contact stress value.

[0250] 5. Pre-process SCADA data and CMS data separately;

[0251] 6. Multi-source data is divided into health data and monitoring data in a ratio of 7:3. Among them, the first 70% of the health data is used to build the memory matrix to train the model, and the remaining health data is used for bandwidth optimization and model testing. The model fitting accuracy indicators corresponding to each data set include mean square error (MSE), mean absolute error (MAE) and coefficient of determination (R 2 ), summarized in Table 3. The results show that the PEAKR model has higher prediction accuracy than the EAKR model based solely on data-driven.

[0252] Table 3 Model fitting accuracy

[0253]

[0254]

[0255] 7. Use the trained model to make predictions and obtain the residual sequence between the predicted value and the actual value.

[0256] 8. For the generated residual sequence, use the MSE alarm indicator for fault warning.

[0257] Early Warning Results: The prediction accuracy and alarm results of the EAKR model using multi-source data are summarized in Table 4. This table shows that the EAKR model, which integrates multi-source data, can identify faults early during wind turbine operation. The model combines historical wind turbine operating data with the main bearing dynamic model response to construct a memory matrix for the wind turbine's healthy operating state. The multi-source data enhances the credibility of the memory matrix, enabling earlier identification of wind turbine faults during the prediction phase. The EAKR model, which integrates multi-source data, outperforms models based on single-source data across various early warning metrics. The EAKR model based on CMS+Physical Mechanism achieves a relative error reduction of over 90% in prediction accuracy and an 80% improvement in early warning. Furthermore, a comparison with an LSTM model based on CMS+Physical Mechanism was also performed. The results demonstrate that the EAKR model outperforms the LSTM model in both fault prediction accuracy and early warning capability, further validating the effectiveness and advantages of the multi-source data fusion approach.

[0258] Table 4 Multi-source data prediction accuracy and alarm results

[0259]

[0260] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications to the technical solutions described in the above embodiments, or equivalent replacement of some or all of the technical features therein, do not deviate the essence of the corresponding technical solutions from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A wind turbine fault early warning method based on digital-physical fusion, characterized in that: include: Step 1: Collect multi-source data; Step 2: Calculate the maximum contact stress of the bearing based on the high-precision mechanism model of the double-row spherical roller bearing; Step 3: Pre-process the multi-source data; Step 4: Establish a PEAKR model, and input the processed multi-source data and maximum contact stress data into the PEAKR model for training to obtain a trained EAKR model; Step 5: Use the PEAKR model to predict the operating data to be monitored, generate predicted values ​​and calculate the residual sequence, and determine the MSE index based on the residual sequence for subsequent equipment status evaluation; In step 6, the MSE indicator based on the moving window strategy is used to analyze the residual output of the PEAKR model. When the MSE indicator exceeds the MSE alarm threshold, it indicates that the health status of the device is abnormal and an alarm is triggered.

2. The wind turbine fault early warning method based on digital-physical fusion according to claim 1 is characterized in that: In step 2, the high-precision mechanism model of the double-row spherical roller bearing includes a load calculation module, an SRBs dynamic mechanism model, and a Hertz contact theory module; first, the load calculation module calculates the bearing load, and the bearing load is input into the SRBs dynamic mechanism model to obtain the rolling element support load, and then the maximum contact stress is calculated according to the Hertz contact theory.

3. The wind turbine fault early warning method based on digital-physical fusion according to claim 2 is characterized in that: The load calculation module includes hub load and bearing load; The hub load includes M x 、M y 、M z 、F x 、F y and F z six independent variables; The thrust element and resistance element of the wind wheel can be expressed as: Where L and D are lift forces, W is the resultant velocity, c and C L is the airfoil chord length and airfoil lift coefficient, C D is the drag coefficient; According to the momentum theorem, acting on the blade ring (r ，r +dr) is: The energy equation is derived by combining blade element theory and momentum theory: The induction factors a and a′ are obtained by simultaneous solution, and the hub load can be obtained by using formulas (1)-(3); The bearing load is equivalent to F B_x 、F B_y , and F B_z Three mutually perpendicular concentrated loads; The axial load is entirely borne by the main shaft, that is, F x =F B_z ; 4. The wind turbine fault early warning method based on digital-physical fusion according to claim 2 is characterized in that: The SRBs kinetic mechanism model includes: Double row spherical roller bearing geometric parameters, R o Indicates the outer ring radius of the bearing, R i Indicates the contour radius of the bearing inner ring, D p is the displacement from the center of the rolling element to the center of the bearing, α is the nominal contact angle of the rolling element, D o represents the rolling element diameter; the inner ring of the bearing has 3 degrees of freedom, while the outer ring is fixed; the total kinetic energy can be calculated as: Among them, x in 、y in 、z in are the displacements of the inner ring / shaft in the x, y and z directions respectively, j is the radial displacement of the jth roller, y j is the angular position of the jth roller, m j is the mass of the jth roller, I j is the moment of inertia of the jth roller, y j is the angular displacement of the jth roller, I s and y s are the moment of inertia and angular displacement of the shaft, N r is the number of rollers, α n is the nominal contact angle; Introduced operator (-1) j , the total potential energy is calculated as: In formulas (8)-(10), δ in+ and δ out+ They represent the contact deformation between the jth roller and the inner and outer rings respectively; the contact type of the spherical roller bearing is point contact, and the contact deformation is calculated as follows: In the above formula, θ j represents the angular position of the jth roller, R ct_inner and R ct_outer Respectively represent the contour radius of the inner and outer rings, α n is the nominal contact angle, j is the roller number; θ j Calculated as: The inner radius can be expressed as: R i =k i0 +k ip sin(N i i ij ) (14) The inner ring contact angle θ ij It can be expressed as: The bearing outer ring raceway radius can be expressed as: R o =k o0 +k op sin(N o i oj ) (16) θ oj It can be written as: According to the Langrange equations, the motion equations of the rolling element in the contact normal direction and the motion equations of the bearing inner ring in the x, y and z directions are: Where, is the total damping force borne by the jth roller, which is the sum of the damping forces of the inner and outer rings; F ax , F ay , F az It is the external force applied to the inner ring of the bearing in three directions, which can be static or dynamic; is the damping force generated in the x, y and z directions when the bearing inner ring contacts the rolling element; R ct_inner ,R ct_outer and R Roller are the contour radii of the inner and outer rings and rolling elements of the bearing respectively; Formulas (18) to (21) are second-order differential equations with a total degree of freedom of N r +3, where N r is the number of rolling elements.

5. The wind turbine fault early warning method based on digital-physical fusion according to claim 2 is characterized in that: The Hertz contact theory module includes: The distribution of contact pressure on the contact ellipse is as follows: In Hertz theory, the maximum contact stress occurs at the geometric center of the contact surface, and its formula is: Where E1, E2, v1, and v2 represent the elastic modulus and Poisson's ratio of the two elastic bodies, respectively. The ∑ρ formula is: Parameter a * and b * The calculation formulas for the variables are: Among them, κ=a / b, e=1-1 / κ 2 , K(e) and E(e) are the complete elliptic integrals of the first and second kinds, respectively; where κ is determined by the following equation: Where F(ρ) is the difference in curvature between the two contact elastic bodies: According to Hertz's theory, the relative displacement of the two elastic bodies is: For two elastic bodies in point contact, the relationship between load and deformation is as follows: Q=K p δ 3 / 2 (32) Finally, the contact stiffness of the two elastic bodies can be derived from formulas (30) and (32):

6. The wind turbine fault early warning method based on digital-physical fusion according to claim 2 is characterized in that: Step 2 includes steps 201 to 207: Step 201, hub load data calculation: hub load data is calculated from the wind speed monitored on-site at the wind farm wind turbines to calculate the subsequent required bearing load; Step 202, bearing load calculation: using the static equilibrium condition of a simply supported beam, calculate the load applied to the bearing under the hub load as input to the bearing dynamics model; Step 203, contact stiffness calculation: Based on Hertz contact theory, the contact stiffness between the rolling elements and the inner and outer rings of the double-row spherical roller bearing is calculated, and the load-deformation relationship is introduced into the SRBs dynamic mechanism model; Step 204, parameter analysis: Substitute the bearing load data and contact stiffness obtained from steps 201 to 203 into the SRBs dynamic mechanism model; Step 205, support load calculation: Based on the bearing dynamics model, the rolling element support load under a given bearing load is calculated for subsequent mechanical analysis; Step 206, contact stress calculation: Based on the support load calculation results, the contact stress is calculated using the Hertz contact formula as the final output of the model; Step 207, model verification: Calculate the maximum contact stress under 16 extreme working conditions and compare it with the SKF calculation report results to verify the accuracy of the model; The maximum contact stress of the bearing is calculated according to the method from step 201 to step 206.

7. The wind turbine fault early warning method based on digital-physical fusion according to claim 1 is characterized in that: Step 4 includes the following steps 401 to 405: Step 401: Divide the multi-source data and the maximum contact stress data calculated in step 2 into a training set, an optimization set, and a test set; Step 402, establishing a PEAKR model; Extract the historical health data of the equipment operation and create the memory matrix X of the PEAKR model, where X i，j represents the i-th vector value of the j-th key variable; for n memory vectors, the memory matrix X of p variables can be expressed as: The monitoring vector is represented by a 1×p matrix v: v=[v1 v2…v p ] (35) The predicted value of the PEAKR model can be obtained by taking a weighted average of each memory vector of the memory matrix X, where the weighted average parameter is estimated using the unit health data; Step 403: training the established EAKR model.

8. The wind turbine fault early warning method based on digital-physical fusion according to claim 7 is characterized in that: Step 4 includes the following steps 4031 to 4034: Step 4031, using the wavelet packet Bayesian noise reduction method to process the training set, optimization set and test set; Step 4032, use the Manhattan method to calculate the monitoring vector v and each memory vector X i The distance between them, we get an n×1 distance vector d i ; Step 4033: Calculate the weight w using the obtained distance matrix d and the Gaussian kernel function. Each element is calculated using the following formula: Among them, the weight w is an n×1 vector matrix; the parameter h is the bandwidth of the kernel function; Step 4034 , automatically optimizing the kernel bandwidth in the model using a new set of health data and the Nelder-Mead optimization algorithm, and obtaining the optimal bandwidth by minimizing the mean square error during model training; The mean square error (MSE) is used to detect the residual between the model prediction value and the true value. The formula for calculating the average mean square error of N test samples is as follows: Where N is the total number of test samples, p is the total number of model variables, is the true value of the jth variable of the i-th sample, is the predicted value of the jth variable of the model corresponding to the i-th sample.

9. The wind turbine fault early warning method based on digital-physical fusion according to claim 8 is characterized in that: After step 403, the method further includes: Step 404, evaluating the trained EAKR model; Calculate the performance indicators of the PEAKR model using the test set, including R 2 , MSE and MAE, the calculation formula is as follows 10. The wind turbine fault early warning method based on digital-physical fusion according to claim 1 is characterized in that: In step 5, the predicted value of the monitoring vector v is predicted by the calculated weight w By each memory vector X i The weighted average of is calculated as follows:

Images