Fan blade stress evaluation system based on digital twinning

Abstract

The invention relates to the technical field of digital twinning, in particular to a fan blade stress evaluation system based on digital twinning, which comprises a twinning model sensing mapping module, a global stress field resolving module, a stress gradient anomaly diagnosis module and a dangerous section comprehensive evaluation module. According to the method, strain vibration time sequence data associated with wind field load-load intensity are screened, a measured value is taken as a reference, and blade global stress field distribution is calculated by using a spatial interpolation model, so that accurate mapping from a limited sensing point to an overall stress state is realized, and then potential stress concentration points are rapidly positioned according to stress space gradient changes; and then automatically identifying and delineating a dangerous section range by using a density clustering algorithm, and performing weighted summation on gradient values to form a set of comprehensive evaluation indexes which are more objective to the dangerous section.

Description

TECHNICAL FIELD

[0001] The application relates to the technical field of digital twinning, in particular to a fan blade stress evaluation system based on digital twinning. BACKGROUND

[0002] The technical field of digital twinning relates to real-time mapping and interaction between a physical entity and its digital model, and core matters include three-dimensional modeling of a physical object, real-time data acquisition, virtual simulation, dynamic behavior analysis and prediction optimization, etc. Through the mapping system of the physical entity in the digital space, combined with multi-source data and physical rule models, the monitoring, analysis and prediction of the entity state are realized, and the technology is widely applied to intelligent manufacturing, equipment health management, industrial process optimization and other scenes. This technical field not only depends on accurate modeling of the physical equipment state, but also requires the establishment of a data and behavior mapping relationship to support the continuous reflection and response of virtual objects to the running state and behavior of the entity. Among them, the traditional fan blade stress evaluation system refers to a technical scheme for calculating and analyzing the stress state of the fan blade under complex wind load and operating conditions during the operation of the wind energy equipment. The main problem it aims to solve is the fatigue damage risk evaluation of the blade due to long-term stress. The traditional method usually collects stress on site by laying strain gauges, and then establishes a blade model by combining a blade structure simulation analysis system based on the finite element method. The blade surface load distribution is derived through the wind speed distribution parameters, and then the stress calculation of the static load and the periodic load is carried out according to the material mechanics parameters, and the stress response characteristics of the blade service state are derived combined with the fatigue life prediction criterion.

[0003] The traditional fan blade stress evaluation system excessively relies on limited physical measuring points and idealized simulation models. Due to the limitations of the number and position of strain gauges, it is difficult to capture the real stress distribution of the blade in the whole domain, especially the stress concentration phenomenon in the non-measuring point area is easily ignored. At the same time, the load distribution derived based on the theoretical wind field model deviates from the actual complex and changeable operating environment, which limits the accuracy of the stress calculation results. For example, sudden turbulence or extreme weather may cause local stress anomalies that the model cannot predict, thereby causing misjudgment of the fatigue damage risk of the blade. SUMMARY

[0004] The purpose of the present application is to solve the shortcomings in the prior art, and a fan blade stress evaluation system based on digital twinning is proposed.

[0005] In order to achieve the above-mentioned purpose, the technical scheme adopted by the present application is as follows: a fan blade stress evaluation system based on digital twinning comprises: The twinning model sensing mapping module screens stress evaluation sensitive parameters, collects strain and vibration time series data, calculates the cross-correlation coefficient with the wind field load parameters, screens the associated parameter set, and transmits it to the global stress field solving module. A global stress field solving module is configured to solve the global stress of the blade, call the real-time numerical values in the associated parameter set as known point stress inputs, perform stress spatial interpolation calculation on unknown nodes in the finite element grid by using a Kriging interpolation model, generate global stress distribution data, and deliver the global stress distribution data to a stress gradient anomaly diagnosis module; The stress gradient anomaly diagnosis module is configured to diagnose stress concentration areas, calculate stress spatial gradient values of each node according to the global stress distribution data, identify abnormal gradient nodes by comparing the stress spatial gradient values with local average gradient difference values, generate an abnormal gradient node set, and deliver the abnormal gradient node set to a dangerous section comprehensive evaluation module. The dangerous section comprehensive evaluation module is configured to evaluate dangerous sections, call the abnormal gradient node set to define dangerous sections by density clustering, extract stress values of nodes in the sections from the global stress distribution data, and output a section stress evaluation index according to weighted summation of the stress spatial gradient values.

[0006] As a further scheme of the present application, the associated parameter set includes a parameter identifier, an associated coefficient value, and a timestamp sequence, the global stress distribution data specifically refers to a grid node index, a three-dimensional stress component, and an interpolation confidence degree, the abnormal gradient node set includes a node spatial coordinate, a gradient deviation amplitude, and an abnormal state identifier, and the section stress evaluation index specifically refers to a dangerous section number, a comprehensive risk weight value, and a stress concentration degree.

[0007] As a further scheme of the present application, the twin model sensing mapping module includes: A time series data acquisition sub-module is configured to acquire strain and vibration time series data, perform timestamp alignment operation under a unified time reference on the strain and vibration time series data, format and integrate multi-channel data streams, and generate a multi-source synchronous time series data stream. A cross-correlation coefficient solving sub-module is configured to acquire a wind field load parameter, perform point-by-point time domain cross-correlation operation on multiple data sequences in the multi-source synchronous time series data stream and the wind field load parameter, and generate a load correlation quantitative coefficient value representing coupling strength between the multiple data sequences in the multi-source synchronous time series data stream and the wind field load parameter. An associated parameter screening sub-module is configured to set a sensitivity discrimination threshold, perform numerical comparison between the load correlation quantitative coefficient value and the sensitivity discrimination threshold, screen parameters with numerical values higher than the sensitivity discrimination threshold, and construct an associated parameter set.

[0008] As a further scheme of the present application, the global stress field solving module includes: The known point stress input submodule calls real-time values in the correlation parameter set, extracts known point node coordinates of the finite element grid, binds and matches the real-time values and the known point node coordinates, establishes a mapping relationship between known point stress input and node position, and generates a stress measurement point space coordinate set; The unknown node interpolation solver submodule calculates the spatial distance between known nodes and the stress difference value fitting spatial variation function based on the stress measurement point space coordinate set, and performs covariance matrix operation on unknown nodes in the finite element grid using the spatial variation function to obtain an interpolated node stress estimation value; The global data fusion generation submodule integrates the stress measurement point space coordinate set and the interpolated node stress estimation value to the node database of the finite element grid, sorts and formats the data structure according to the node number index, and establishes global stress distribution data.

[0009] As a further scheme of the present application, the stress gradient anomaly diagnosis module comprises: The stress field data analysis submodule indexes all discrete grid nodes in the global stress distribution data, extracts the three-dimensional space coordinates and stress tensor components of each node, and performs structured packaging to establish a node stress feature vector; The gradient anomaly recognition submodule calculates the stress spatial gradient of each node based on the node stress feature vector, defines the arithmetic mean of the stress spatial gradients of all nodes within the statistical range of the local spatial range, obtains the gradient deviation amount by subtracting the stress spatial gradient of a single node from the arithmetic mean, and simultaneously determines the gradient deviation amount according to a preset gradient deviation threshold to generate a node gradient anomaly identifier; The diagnosis set construction submodule calls the node gradient anomaly identifier, traverses all the nodes, collects the node identifiers identified as abnormal states, and generates an abnormal gradient node set.

[0010] As a further scheme of the present application, the dangerous cross section comprehensive evaluation module comprises: The dangerous cross section definition submodule calls the Euclidean distance between multiple nodes in the abnormal gradient node set, performs density clustering based on a neighborhood radius threshold and a core point minimum number benchmark, filters node clusters whose density reaches the standard, and generates a dangerous cross section core node set; The node stress extraction submodule retrieves global stress distribution data according to the node index in the dangerous cross section core node set, extracts stress state components corresponding to multiple node indexes, constructs an ordered numerical sequence, and generates a cross section node stress sequence; The cross-section index weighted summation submodule uses the stress space gradient value as a weight factor for each stress value in the node stress sequence within the cross-section, performs weighted product calculation, and then accumulates and sums all the weighted products in the sequence to generate a cross-section stress assessment index.

[0011] As a further solution of the present invention, the unknown node interpolation operator module specifically adopts a Gaussian model to fit the spatial variation function; The unknown node interpolation operator module calculates the model parameters of the Gaussian model using the least squares method based on the Euclidean distance and stress value difference between multiple known nodes in the stress measurement point spatial coordinate set. The calculation formula of the Gaussian model is: ; in, The spatial separation distance is The stress semivariance value of two known nodes is, represents the three-dimensional Euclidean distance between any two known nodes, represents the nugget effect value reflecting measurement error and micro-scale spatial variability, represents the sill value corresponding to the total spatial variation variance of the stress regionalization variable, It indicates the maximum spatial autocorrelation distance between nodes when the spatial variation function reaches the base value, that is, the range; The unknown node interpolation operator module uses the Gaussian model whose model parameters have been obtained to construct a covariance matrix between known nodes and a covariance vector between an unknown node and all known nodes, obtains an unbiased stress estimate of multiple unknown nodes by solving the Kriging equations, and establishes the interpolation node stress estimate.

[0012] As a further solution of the present invention, the gradient anomaly identification submodule calls the node stress characteristic vector to calculate the equivalent stress value of each node; The gradient anomaly identification submodule uses the central difference method to calculate the stress change rate of each node in the three coordinate axis directions and calculates the modulus of the stress space gradient. The calculation formula of the gradient deviation is: ; in, Indicates that the index is The gradient deviation of the node, Indicates that the index is The modulus of the spatial gradient of the equivalent stress at the node, Indicates that the index is The total number of nodes in the neighborhood determined by the local spatial range with the node as the center, representing the sum of the modulus values of the equivalent stress space gradient of all nodes in the neighborhood; The gradient anomaly identification submodule compares the calculated gradient deviation amount with the preset gradient deviation threshold value point by point, and when the gradient deviation amount is greater than the preset gradient deviation threshold value, sets the node gradient anomaly identification of the corresponding node to be abnormal.

[0013] As a further scheme of the present application, the cross-section exponential weighted summation submodule calls the dangerous cross-section core node set, and for each core node in the set, extracts the equivalent stress value from the global stress distribution data, and obtains the stress space gradient value from the calculation result of the gradient anomaly identification submodule; The cross-section exponential weighted summation submodule calculates the gradient weight coefficient of each core node, which is the ratio of the stress space gradient value of a single core node to the sum of the stress space gradient values of all core nodes in the dangerous cross-section core node set; The cross-section exponential weighted summation submodule multiplies the equivalent stress value of each core node by the corresponding gradient weight coefficient to obtain the weighted stress contribution value of the core node; The cross-section exponential weighted summation submodule accumulates and sums the weighted stress contribution values of all core nodes in the dangerous cross-section core node set by the following formula: ; generates the cross-section stress assessment index; wherein, represents the cross-section stress assessment index, represents the dangerous cross-section core node set, and are indexes traversing the dangerous cross-section core node set, represents the summation operation on all core nodes in the dangerous cross-section core node set, represents the equivalent stress value of the core node with index , represents the stress space gradient value of the core node with index , represents the cumulative summation of the stress space gradient values of all core nodes in the dangerous cross-section core node set.

[0014] Compared with the prior art, the present application has the following advantages and positive effects: In the application, by screening strain vibration time series data with strong correlation with wind field load, using spatial interpolation model to calculate global stress field distribution based on measured values, accurate mapping from limited sensing points to overall stress state is realized, and then potential stress concentration points are quickly located according to stress space gradient change, and dangerous section range is automatically identified and circled by using density clustering algorithm, and through gradient value weighted summation, a more objective comprehensive evaluation index of dangerous section is formed. BRIEF DESCRIPTION OF DRAWINGS

[0015] Figure 1 The stress evaluation system overall technical scheme flowchart of the application is shown in the figure. Figure 2 The sensitive parameter screening and correlation parameter set construction flowchart of the application is shown in the figure. Figure 3 The global stress field solving flowchart based on Kriging interpolation of the application is shown in the figure. Figure 4 The stress gradient anomaly diagnosis flowchart of the application is shown in the figure. Figure 5 The dangerous section clustering and comprehensive evaluation flowchart of the application is shown in the figure. DETAILED DESCRIPTION

[0016] In order to make the purpose, technical scheme and advantages of the application clearer, the technical scheme realized based on software is described in detail below in combination with system architecture diagram and embodiments. It should be understood that the specific embodiments described herein are only used to explain the technical scheme of the application, and do not constitute a limitation on the protection scope.

[0017] In the description of the application, the system architecture relationship or data processing flow indicated by the terms "level", "module", "interface", "data flow", "client", "server" and the like are defined based on the corresponding architecture diagram or flowchart of the embodiments, and this description method is only used to clearly explain the logical relationship of each element in the technical scheme, and is not limited to the physical deployment form. The "multiple" contains two or more technical units, including but not limited to multiple data nodes, processing threads, service instances or functional components, and the specific number is determined according to the actual business scenario.

[0018] Please refer to Figure 1 and Figure 2 , the application provides a technical scheme: a fan blade stress evaluation system based on digital twinning includes: The twinning model sensing mapping module screens the stress evaluation sensitive parameters, collects strain and vibration time series data, calculates the cross-correlation coefficient with the wind field load parameter, screens the correlation parameter set, and transmits it to the global stress field solving module. The correlation parameter set includes parameter identification, correlation coefficient value and time stamp sequence. The twin model sensor mapping module comprises: The time series data acquisition submodule acquires strain and vibration time series data, performs timestamp alignment operation under a unified time reference on the strain and vibration time series data, and formats and integrates the multi-channel data stream to generate a multi-source synchronous time series data stream; The cross-correlation coefficient solving submodule obtains the wind field load parameter, performs point-by-point time domain cross-correlation operation on the multiple data sequences in the multi-source synchronous time series data stream and the wind field load parameter, and generates a load correlation quantization coefficient value representing the coupling strength between the multiple data sequences in the multi-source synchronous time series data stream and the wind field load parameter; The correlation parameter screening submodule sets a sensitivity discrimination threshold, compares the load correlation quantization coefficient value with the sensitivity discrimination threshold, screens parameters with a value higher than the sensitivity discrimination threshold, and constructs a correlation parameter set.

[0019] The specific implementation process of the twin model sensor mapping module is as follows: first, the function of the time series data acquisition submodule is executed, and the submodule acquires data for No. 2 blade of a horizontal axis wind turbine with a model number of FD160, a single machine capacity of 5.0 MW, and a blade length of 78 meters. In the key area of 50% to 80% of the blade span, 12 strain monitoring points (labeled as SG01 to SG12) and 4 three-axis vibration monitoring points (labeled as ACC01 to ACC04) are arranged along the chord. The strain monitoring points use distributed fiber Bragg grating sensors, and the vibration monitoring points use MEMS acceleration sensors. The data acquisition system synchronously acquires data of all sensor channels at a sampling frequency of 100 Hz. In a certain evaluation period, for example, from 10:00:00 to 10:00:10 on September 5, 2025, 1000 time points of data are collected. The acquired raw data stream is first subjected to timestamp alignment operation, and the system uses the GPS clock signal as the unified time reference to correct the microsecond-level time deviation that may be caused by data transmission delay, so as to ensure that the data points of all channels strictly correspond in time sequence. Subsequently, formatting and integration are performed, and the data (wavelength drift of fiber Bragg grating sensor and voltage signal of MEMS acceleration sensor) from different sensor sources are uniformly converted into standard physical quantity units, strain data are converted into micro-strain (με), vibration data are converted into acceleration (m / s²), and a data matrix containing time stamp, sensor identifier and physical quantity value is integrated to generate a multi-source synchronous time series data stream.

[0020] Next, the function of the cross-correlation coefficient solver module is executed. This module first acquires the wind farm load parameters in the same time period (10:00:00-10:00:10 on September 5, 2025) from the supervisory control and data acquisition (SCADA) system. In this embodiment, the average wind speed at the hub center is selected as the key wind farm load parameter. The acquired wind speed time series data also has 1000 data points at a sampling frequency of 100 Hz. Subsequently, the 16 data sequences (12 strain sequences and 4 vibration sequences) in the multi-source synchronous time series data stream are respectively subjected to point-by-point time domain cross-correlation operation with the wind speed time series data. Taking the time series data of strain monitoring point SG01 and the wind speed time series data as an example, the cross-correlation coefficient is calculated. First, the average values of the two sequences are calculated, assuming that the average value of SG01 strain is με, and the average value of the wind speed is m / s. Then, the sum of the deviation products at each time point , and the respective deviation squares and are calculated. By dividing the sum of the deviation products by the square root of the product of the two deviation squares, the cross-correlation coefficient of the SG01 sequence and the wind speed sequence, i.e., the load correlation quantification coefficient value, is obtained. For example, the coefficient value of SG01 is calculated to be 0.91, the coefficient value of SG02 is 0.88, the coefficient value of ACC01 is 0.75, and the coefficient value of ACC04 is 0.62. Repeat this process for all 16 channels to generate a set of 16 load correlation quantification coefficient values.

[0021] The specific calculation formula of the above cross-correlation coefficient is: wherein represents the cross-correlation coefficient between two time series x and y, represents the total number of data points in the time series, represents the index of the data point in the time series, represents the value of the i-th data point in the time series x, represents the value of the i-th data point in the time series y, represents the average value of the time series x, represents the average value of the time series y, represents the summation operation on all data points with index from 1 to n.

[0022] ​Finally, the function of the correlation parameter screening submodule is executed. The submodule first sets a sensitivity discrimination threshold. The setting of the threshold is based on the following: in a series of controllable wind tunnel experiments, the same type of blade is subjected to different levels of stable wind load from 5 m / s to 25 m / s, and the response data of all sensors are recorded. By analyzing the experimental data, it is found that when the cross-correlation coefficient of the sensor response and the wind load is less than 0.80, the error introduced by the data in the subsequent stress field calculation exceeds 5%, and the stress distribution of the key area of the blade cannot be accurately reflected. To ensure the accuracy of subsequent calculation, 100 sets of repetitive experiments are performed; the statistical results show that when the cross-correlation coefficient is higher than 0.80, the parameters and external loads show strong correlation. Therefore, the sensitivity discrimination threshold in this embodiment is set to 0.80. The 16 load correlation quantification coefficient values calculated in the previous step are compared with the threshold value 0.80 point by point. For example, the coefficient value of SG01 is 0.91, which is higher than 0.80 and is screened and retained; The coefficient value of SG02 is 0.88, which is higher than 0.80 and is screened and retained; The coefficient value of ACC01 is 0.75, which is lower than 0.80 and is excluded; The coefficient value of ACC04 is 0.62, which is lower than 0.80 and is excluded. After traversing all 16 coefficient values, all values higher than 0.80 are screened out. In this example, it is assumed that SG01, SG02, SG03, SG05, SG06, SG08, SG09, SG11, a total of 8 strain parameters and ACC02, a total of 1 vibration parameter are finally screened out. Based on these screened parameters, a correlation parameter set is constructed, and the data structure of the set is: [{‘parameter identification’:‘SG01’,‘correlation coefficient value’:0.91,‘time stamp sequence’:[…]},{‘parameter identification’:‘SG02’,‘correlation coefficient value’:0.88,‘time stamp sequence’:[…]},…,{‘parameter identification’:‘ACC02’,‘correlation coefficient value’:0.82,‘time stamp sequence’:[…]}].

[0023] Please refer to Figure 1 and Figure 3 , the global stress field calculation module calculates the global stress of the blade, calls the real-time values in the correlation parameter set as known point stress input, uses the Kriging interpolation model to perform stress space interpolation calculation on the unknown nodes in the finite element grid, generates global stress distribution data, and transmits the data to the stress gradient anomaly diagnosis module; The global stress distribution data is specifically the grid node index, the three-dimensional stress component, and the interpolation confidence; The global stress field calculation module includes: The known point stress input submodule calls the real-time numerical value in the associated parameter set, extracts the known point node coordinates of the finite element mesh, binds and matches the real-time numerical value with the known point node coordinates, establishes a mapping relationship between the known point stress input and the node position, and generates a stress measurement point spatial coordinate set; The unknown node interpolation operator module calculates the spatial distance between known nodes and the stress difference based on the spatial coordinate set of stress measurement points, fits the spatial variation function, and performs covariance matrix operations on unknown nodes in the finite element mesh using the spatial variation function to obtain the estimated stress value of the interpolated node; The unknown node interpolation operator module specifically uses the Gaussian model to fit the spatial variation function; The unknown node interpolation operator module uses the least squares method to calculate the model parameters of the Gaussian model based on the Euclidean distance and stress value difference between multiple known nodes in the stress measurement point spatial coordinate set. The calculation formula of the Gaussian model is: ; in, The spatial separation distance is The stress semivariance value of two known nodes is, represents the three-dimensional Euclidean distance between any two known nodes, represents the nugget effect value reflecting measurement error and micro-scale spatial variability, represents the sill value corresponding to the total spatial variation variance of the stress regionalization variable, It indicates the maximum spatial autocorrelation distance between nodes when the spatial variation function reaches the base value, that is, the range; The unknown node interpolation operator module uses the Gaussian model with obtained model parameters to construct the covariance matrix between known nodes and the covariance vector between the unknown node and all known nodes. By solving the Kriging equations, it obtains unbiased stress estimates for multiple unknown nodes and establishes interpolated node stress estimates. The global data fusion generation submodule integrates the two sets of data into the node database of the finite element mesh based on the spatial coordinate set of the stress measurement points and the interpolated node stress estimates, sorts them by node number index and formats the data structure to establish the global stress distribution data.

[0024] The specific implementation process of the global stress field solving module is as follows: first, the function of the known point stress input submodule is executed, and the associated parameter set generated by the previous module is called. Taking the instantaneous data with a time stamp of "10:00:10.000" as an example, the real-time values in the associated parameter set at this time are extracted. For example, the real-time strain values of SG01, SG02 and SG03 are 350με, 380με and 365με respectively, and the real-time data of the other 6 sensors are obtained. These strain data are converted into stress values through the elastic modulus of the material (for glass fiber reinforced composite material, the value is 45GPa), for example, the stress value corresponding to SG01 is MPa. At the same time, the vibration data are converted into equivalent stress values through the transfer function relationship model with the strain data. The submodule then extracts the precise three-dimensional spatial coordinates of the 9 sensors from the pre-established blade finite element model database. For example, the node coordinates of SG01 are (40.5, 1.2, 0.3) meters, and the node coordinates of SG02 are (40.8, 1.5, 0.3) meters. Subsequently, the calculated real-time stress values are matched with the corresponding node coordinates one by one, for example, 15.75MPa is bound with (40.5, 1.2, 0.3) meters. After the binding of all 9 known points is completed, a mapping relationship between the known point stress input and the node position is established, and a stress measurement point spatial coordinate set containing 9 data points is generated, and the data structure is as follows: [{coordinate:(40.5, 1.2, 0.3), stress value: 15.75MPa}, {coordinate:(40.8, 1.5, 0.3), stress value: 17.10MPa}, …].

[0025] Next, the function of the unknown node interpolation solving submodule is executed. Based on the stress measurement point spatial coordinate set, the spatial distance and stress difference between the known nodes are calculated to fit the spatial variation function. Taking the two known nodes SG01 and SG02 as an example, first, the three-dimensional Euclidean distance between them is calculated , that is meters. Then the semi-variance of the stress values of the two nodes is calculated , that is MPa². All 9 known nodes are combined in pairs (a total of pairs), and 36 groups of data points are calculated. Subsequently, the least squares method is used to fit the 36 data points to determine the model parameters , and of the Gaussian model.

[0026] In formula , the explanations of the parameters and the operation logic are as follows: represents the spatial semi-variance, which measures the spatial distance of a measure of the degree of difference in stress values between two points; represents the three-dimensional Euclidean distance between two known nodes; is the nugget value, representing the semivariance when the distance between two points tends to zero, reflecting the randomness caused by measurement error and spatial variability smaller than the sampling scale; is the sill value, representing the total spatial variability variance of the regionalized variable (stress), which is the stable value that the semivariance function tends to as the distance increases; is the range, referring to the distance between nodes when the semivariance function reaches the sill value, representing the spatial autocorrelation range of stress values. The entire formula describes the process of stress difference gradually increasing with distance and eventually reaching a stable level through an exponential function, which conforms to the general law in the physical world that similar points have similar properties and distant points have no correlation. Among them, is the core of the Gaussian model, determining the shape of the function curve, and the coefficient 3 ensures that the range a is obtained when the function value reaches about 95% of the sill value. By introducing the parameter , , that reflects the spatial structure characteristics, the model can fit the spatial variability of the stress field, thereby improving the accuracy of interpolation calculation.

[0027] Parameter acquisition and calculation process: , the nugget value, is determined by repeatedly measuring the same monitoring point under stable load. In the blade static test, a constant load is applied to the SG01 position, and 100 stress data are continuously collected, and the variance of these data is calculated to be 0.04 MPa². This variance is the embodiment of measurement error, so it is set to MPa². , the sill value, reflects the overall variance of stress in the observation area. According to the 9 stress values (e.g. 15.75, 17.10, …, 16.50 MPa) in the stress point spatial coordinate set, the sample variance is calculated to be 2.85 MPa². Considering that the sample variance may be slightly lower than the overall variance, and combining historical data, the sill value is set to MPa². , the range, is determined by observing and fitting the experimental variogram graph. The 36 data points are plotted, and it is observed that the semivariance value tends to be stable after a distance of about 5 meters. Through least squares fitting, the optimal range value is meters. Substituting these parameters into the Gaussian model formula, the specific spatial variability function is obtained: The advantage of this formula is that it not only considers the distance between data points, but also 、 、 The three parameters precisely describe the spatial structural characteristics of the stress field under a specific physical environment, making the interpolation result closer to the actual physical distribution. Compared with simple interpolation methods such as inverse distance, it has higher physical fidelity and prediction accuracy. Then, for an unknown node in the finite element mesh, such as the U01 node with coordinates of (40.6, 1.3, 0.3) meters, the Gaussian model with the obtained model parameters is used to construct the known nodes. Covariance matrix and the unknown node U01 with all 9 known nodes Covariance vector. By solving the Kriging equations, an unbiased stress estimate of 16.28 MPa is obtained for node U01. This process is repeated for all unknown nodes to establish interpolated node stress estimates. The result of 16.28 MPa is a linear unbiased estimate calculated based on the stress values ​​of the surrounding known points and spatial autocorrelation. It shows that under the current load, the stress level at node U01 is between its neighbors SG01 (15.75 MPa) and SG02 (17.10 MPa), and is closer to SG01, which is consistent with its spatial location.

[0028] Finally, the global data fusion generation submodule is executed. This submodule integrates the spatial coordinates of the stress measurement points (nine known points) with the interpolated nodal stress estimates (for example, 150,000 unknown nodes). These two sets of data are imported into the node database of the finite element mesh and indexed and sorted by node number. The data structure is formatted to include fields such as node index, 3D coordinates, three principal stress components (obtained through supplementary calculations), von Mises equivalent stress, and interpolation confidence. Ultimately, high-resolution global stress distribution data covering the entire blade surface is established.

[0029] See also Figure 1 and Figure 4 ,Stress gradient anomaly diagnosis module, diagnoses stress concentration areas, calculates the stress spatial gradient value of each node based on the global stress distribution data, identifies abnormal gradient nodes by comparing with the local average gradient difference, generates an abnormal gradient node set, and transmits it to the dangerous section comprehensive assessment module; The abnormal gradient node set includes node space coordinates, gradient deviation amplitude, and abnormal state identification; The stress gradient anomaly diagnosis module includes: The stress field data analysis submodule indexes all discrete grid nodes in the global stress distribution data, extracts the three-dimensional spatial coordinates and stress tensor components of each node, performs structured packaging, and establishes the node stress characteristic vector; The gradient anomaly identification submodule calculates the stress space gradient of each node based on the node stress feature vector, and defines the arithmetic mean of the stress space gradient of all nodes within the local space range. The stress space gradient of a single node is subtracted from the arithmetic mean to obtain the gradient deviation amount, and the gradient deviation amount is determined according to the preset gradient deviation threshold to generate a node gradient anomaly identifier; The gradient anomaly identification submodule calls the node stress feature vector to calculate the equivalent stress value of each node; The gradient anomaly identification submodule calculates the stress change rate of each node in three coordinate axis directions using the central difference method, and calculates the modulus value of the stress space gradient. The calculation formula of the gradient deviation amount is: ; Wherein, represents the gradient deviation amount of the node with index , represents the modulus value of the equivalent stress space gradient of the node with index , represents the total number of nodes in the neighborhood determined according to the local space range with the node with index as the center, represents the sum of the modulus values of the equivalent stress space gradients of all nodes in the neighborhood; The gradient anomaly identification submodule compares the calculated gradient deviation amount with the preset gradient deviation threshold point by point, and when the gradient deviation amount is greater than the preset gradient deviation threshold, the node gradient anomaly identifier of the corresponding node is set to abnormal; The diagnostic set construction submodule calls the node gradient anomaly identifier, traverses all nodes, and collects the node identifiers identified as abnormal state to generate an abnormal gradient node set.

[0030] The specific implementation process of the stress gradient anomaly diagnosis module is as follows: first, execute the function of the stress field data analysis submodule. This submodule processes the global stress distribution data generated by the previous module. The global stress distribution data contains the stress information of all (for example, 150009) discrete grid nodes in the blade finite element model. The submodule indexes these nodes one by one, for example, starting from node 1, extracts its three-dimensional space coordinates, such as (40.5, 1.2, 0.3) meters, and extracts the stress tensor component of the node. These data are structured and packaged to form the stress feature vector of the node, in the format [node index, x, y, z, , , …]. After all nodes are processed, a complete node stress feature vector database is established.

[0031] Next, the function of the gradient anomaly identification submodule is executed. Based on the node stress feature vector, the von Mises equivalent stress value of each node is first calculated . Taking node i with index i as an example, its stress tensor components are , , , , , (unit: MPa), and its value is MPa. The von Mises equivalent stress value of all nodes is calculated. Subsequently, the central difference method is used to calculate the stress change rate of each node in three coordinate axis directions, and then the modulus value of the stress space gradient is obtained. Taking node i as an example, the stress change rate of node i in the x direction is calculated by the stress values and coordinate differences of its adjacent nodes (i+1, i-1), i.e. . Similarly, the change rates in the y and z directions are calculated. Assuming that the three components of node i are calculated as 50.0 MPa / m, 60.0 MPa / m, and 25.0 MPa / m, respectively, then the modulus value of the stress space gradient of node i is MPa / m.

[0032] In formula , the parameters and operation logic are as follows: is the gradient deviation of node i, representing the prominence of the stress gradient of the node relative to its local environment; is the modulus value of the von Mises equivalent stress space gradient of node i, reflecting the severity of the stress change of the node; is the total number of nodes in the local spatial neighborhood centered on node i; is the sum of the stress space gradient modulus values of all nodes in the neighborhood. The calculation logic of the entire formula is as follows: first, the stress gradient size of node i itself is calculated, then the average stress gradient size of its local neighborhood is calculated, and finally the two are subtracted. This difference reveals whether node i is a stress gradient “hot spot”. This formula uses a relative comparison method to identify anomalies; in different regions of the blade, the normal stress gradient level differs, and using the local average value as a reference can adapt to such differences to locate abnormal stress concentration points caused by material defects or minor damage.

[0033] Parameter acquisition and calculation process: define the local spatial range, in this embodiment, set the neighborhood radius to 0.1 meters, that is, all nodes with a Euclidean distance less than or equal to 0.1 meters from node i belong to its neighborhood. After searching, there are 10 nodes (including itself) in the neighborhood of node i, that is . , the gradient modulus of node i has been calculated as 82.0 MPa / m. , the gradient modulus of all 10 nodes in the neighborhood is summed up. Assuming that the gradient modulus of these 10 nodes is [82.0, 35.5, 40.1, 38.2, 41.5, 36.8, 39.0, 42.3, 37.6, 40.0] (MPa / m), then the sum is 433.0 MPa / m. MPa / m. Substituting the formula to calculate the gradient deviation: MPa / m. Next, the calculated gradient deviation is compared with the preset gradient deviation threshold. The threshold is set based on a large amount of simulation analysis and fatigue test data. By performing finite element analysis on blade models containing known micro defects (such as pores, inclusions), the gradient deviation distribution of nodes around the defects is statistically analyzed. The results show that 98% of the defect areas have a gradient deviation exceeding 30.0 MPa / m. Combined with the data verification of the early crack initiation position in the accelerated fatigue experiment, the gradient deviation threshold is finally determined to be 30.0 MPa / m. The gradient deviation of node i is 38.7 MPa / m, which is greater than the threshold of 30.0 MPa / m, so the node gradient anomaly identification of node i is set to "abnormal" (or marked as 1). Repeat this determination process for all nodes.

[0034] Finally, the function of the diagnostic set construction submodule is executed. This submodule traverses the gradient anomaly identification of all 150009 nodes. The identifiers of all nodes identified as "abnormal" (including node numbers, spatial coordinates), gradient deviation amplitudes, and abnormal state identifiers (for example, 1 represents abnormal, 0 represents normal) are collected. Finally, the abnormal gradient node set is generated, for example: [{node coordinates: (x_i, y_i, z_i), gradient deviation amplitude: 38.7 MPa / m, abnormal state identifier: 1}, {node coordinates: (x_k, y_k, z_k), gradient deviation amplitude: 35.2 MPa / m, abnormal state identifier: 1}, …].

[0035] Please refer to Figure 1 and Figure 5, a dangerous section comprehensive evaluation module, evaluating a dangerous section, calling an abnormal gradient node set to perform density clustering to define the dangerous section, extracting node stress values in the section from global stress distribution data, and performing weighted summation according to stress space gradient values to output a section stress evaluation index; The section stress evaluation index specifically refers to a dangerous section number, a comprehensive risk weight, and a stress concentration degree. The dangerous section comprehensive evaluation module includes: The dangerous section definition submodule calls the Euclidean distance between multiple nodes in the abnormal gradient node set, performs density clustering based on a neighborhood radius threshold and a minimum number of core points reference, filters node clusters that meet the density standard, and generates a dangerous section core node set. The node stress extraction submodule retrieves global stress distribution data according to node indexes in the dangerous section core node set, extracts stress state components corresponding to multiple node indexes, constructs an ordered numerical sequence, and generates a node stress sequence in the section. The section index weighted summation submodule calls the stress space gradient value as a weight factor for each stress value in the node stress sequence in the section, performs weighted product calculation, and then accumulates and sums all weighted products in the sequence to generate a section stress evaluation index. The section index weighted summation submodule calls the dangerous section core node set and extracts equivalent stress values from global stress distribution data for each core node in the set, and obtains stress space gradient values from the calculation results of the gradient anomaly identification submodule. The section index weighted summation submodule calculates the gradient weight coefficient of each core node, which is the ratio of the stress space gradient value of a single core node to the sum of the stress space gradient values of all core nodes in the dangerous section core node set. The section index weighted summation submodule multiplies the equivalent stress value of each core node by the corresponding gradient weight coefficient to obtain the weighted stress contribution value of the core node. The section index weighted summation submodule accumulates and sums the weighted stress contribution values of all core nodes in the dangerous section core node set according to the following formula: ; to generate a section stress evaluation index. wherein, represents the section stress evaluation index, represents the dangerous section core node set, and are indexes traversing the dangerous section core node set, represents a summation operation on all core nodes in the dangerous section core node set, represents the equivalent stress value of the core node with index , the stress space gradient value of the core node with index the stress space gradient value of the core node with index

[0036] The specific implementation process of the dangerous section comprehensive evaluation module is as follows: first, the function of the dangerous section definition submodule is executed, and the abnormal gradient node set generated by the previous module is called. It is assumed that the set contains 50 abnormal gradient nodes. In order to identify the dangerous area with physical meaning from these discrete abnormal points, a density-based spatial clustering algorithm (DBSCAN) is used. The algorithm needs to set two core parameters: neighborhood radius threshold (Eps) and core point minimum number benchmark (MinPts). The setting of the two parameters is as follows: the setting of the neighborhood radius threshold (Eps) refers to the average mesh size of the blade finite element model and the typical micro-crack propagation scale. In this embodiment, the average mesh size is 0.05 meters, and considering that a dangerous area should cover at least several grid units, Eps is set to 0.15 meters. The setting of the core point minimum number benchmark (MinPts) is based on the definition of the critical crack size in material fracture mechanics, which considers that a damage area that can continue to expand should at least contain a certain number of continuous high stress gradient points. Through statistical analysis of historical blade damage cases, it is found that most structural damage starts from an area where at least 5 or more abnormal gradient nodes are gathered. Therefore, MinPts is set to 5. The clustering process is as follows: from the abnormal gradient node set, an unvisited node P is selected, and the number of nodes in the Eps neighborhood (0.15 meter radius) of P is calculated. If the number is greater than or equal to MinPts (5), P is a core point, and a new cluster is created. Then all nodes in the neighborhood of P are added to the cluster, and the neighborhood check process is repeated for the newly added nodes to expand the cluster range, until all density-reachable nodes are included. If the number of nodes in the neighborhood of P is less than MinPts, P is temporarily marked as a noise point. After traversing all nodes, the clustering is completed. It is assumed that 3 node clusters are identified from the 50 abnormal nodes, containing 12, 8 and 6 nodes respectively, and 24 nodes are identified as noise points. The 3 node clusters are defined as dangerous sections, and a dangerous section core node set is generated.

[0037] Next, the function of the node stress extraction submodule is executed. According to the node index in the dangerous section core node set, for example, for the first dangerous section (marked as C1) containing 12 nodes, the global stress distribution data is retrieved. The von Mises equivalent stress values corresponding to the 12 node indices are extracted, and an ordered numerical sequence is constructed.

[0038] Table 1 Dangerous section C1 core node data table​

[0039] The von Mises equivalent stress values and stress space gradient values of the 12 core nodes in the dangerous section C1 are listed in Table 1. These data will be used for subsequent section index weighted summation calculation.

[0040] Finally, the function of the section index weighted summation submodule is executed. This submodule calls the corresponding stress space gradient value as the weight factor for each stress value in the stress sequence of the nodes in the section, and performs weighted summation.

[0041] In the formula , the parameters and operation logic are as follows: is the final calculated section stress assessment index, which integrates the stress level and stress concentration degree in the section; represents the core node set of the dangerous section, which is the C1 set in this example; Index and are iteration variables that traverse all nodes in the set C; indicates the summation operation on all core nodes in the set C; is the von Mises equivalent stress value of the core node with index k; is the stress space gradient value of the core node with index k. The fractional term calculates the proportion of the stress gradient of a single node k in the total gradient of the entire section, and this ratio is the gradient weight coefficient. The calculation logic of the entire formula is that for each core node in the section, its own stress value is multiplied by its gradient weight coefficient to obtain the weighted stress contribution value of the node; Finally, the weighted stress contribution values of all core nodes in the section are added up to obtain the final assessment index. This formula introduces stress gradient as the weight, so that the assessment index can comprehensively reflect the average stress size and the severity of stress concentration in the section; the greater the stress gradient of a node, the higher the weight it is given in the comprehensive evaluation, making the dangerous indication of high stress gradient area more sensitive.

[0042] The parameter acquisition and calculation process are based on the data in Table 1: ={1001,1002,…,1012} and are directly obtained from Table 1. First, calculate the denominator i.e. the stress space gradient values of all core nodes are accumulated and summed up: MPa / m. Then, the weighted stress contribution value of each core node is calculated one by one and accumulated: . Finally, the cross-section stress evaluation index is 257.96.

[0043] The result is used for risk level determination; the risk level standard is formulated based on statistical analysis of the full life cycle data of the blade, and is specifically divided into: <150 is "low risk"; 150≤ <250 is "medium risk"; ≥250 is "high risk".

[0044] The calculation result in this embodiment is 257.96, which indicates that the risk level of the dangerous cross-section C1 is "high risk". The relationship between this numerical result and the step result is direct: the calculated value is directly mapped to the risk state identification of the cross-section by comparing it with the preset risk interval. For operation and maintenance, this means that high-precision non-destructive detection needs to be immediately performed on the physical location corresponding to the dangerous cross-section C1, and limiting operation or preventive maintenance measures may need to be taken. Finally, the generated evaluation result includes the dangerous cross-section number "C1", the comprehensive risk weight value (i.e. the evaluation index) "257.96", and the stress concentration degree "high".

[0045] The above embodiments demonstrate the preferred embodiments of the present application, and any equivalent adjustment of the technical solutions based on software engineering methods is within the protection scope, including but not limited to: implementing algorithm logic using different programming languages, service reconstruction of functional modules, adjustment of data interaction protocols, optimization of resource scheduling strategies, etc. Any implementation scheme derived through reasonable modification of the data processing flow, service calling link or system architecture level without deviating from the technical core of the present application should be considered within the protection scope defined by the claims of the present application.

Claims

1. A wind turbine blade stress assessment system based on digital twinning, characterized by, The system comprises: A twin model sensor mapping module filters stress evaluation sensitive parameters, collects strain and vibration time series data, calculates cross-correlation coefficients with wind field load parameters, filters a set of associated parameters, and delivers them to a global stress field solving module; The global stress field solving module solves blade global stress, calls real-time values in the set of associated parameters as known point stress input, performs stress spatial interpolation calculation on unknown nodes in the finite element grid using the Kriging interpolation model, generates global stress distribution data, and delivers them to a stress gradient anomaly diagnosis module; The stress gradient anomaly diagnosis module diagnoses stress concentration areas, calculates stress spatial gradient values for each node according to the global stress distribution data, identifies abnormal gradient nodes by comparing with local average gradient difference, generates an abnormal gradient node set, and delivers it to a dangerous section comprehensive evaluation module; The dangerous section comprehensive evaluation module evaluates dangerous sections, calls the abnormal gradient node set to define dangerous sections by density clustering, extracts node stress values within the section from the global stress distribution data, and outputs a section stress evaluation index according to stress spatial gradient value weighted summation.

2. The digital-twin-based fan blade stress assessment system of claim 1, wherein, The set of associated parameters includes parameter identification, correlation coefficient value, and time stamp sequence, the global stress distribution data specifically refers to grid node index, three-dimensional stress component, and interpolation confidence, the abnormal gradient node set includes node spatial coordinates, gradient deviation amplitude, and abnormal state identification, and the section stress evaluation index specifically refers to dangerous section number, comprehensive risk weight, and stress concentration degree.

3. The digital-twin-based fan blade stress assessment system of claim 2, wherein, The twin model sensor mapping module comprises: A time series data collection sub-module collects strain and vibration time series data, performs timestamp alignment operation under a unified time reference on the strain and vibration time series data, formats and integrates multi-channel data streams, and generates multi-source synchronous time series data streams; A cross-correlation coefficient solving sub-module obtains wind field load parameters, performs point-by-point time domain cross-correlation operation on multiple data sequences in the multi-source synchronous time series data streams and the wind field load parameters, and generates load correlation quantitative coefficient values representing coupling strength between the multiple data sequences in the multi-source synchronous time series data streams and the wind field load parameters; An associated parameter filtering sub-module sets a sensitivity discrimination threshold, compares the load correlation quantitative coefficient values with the sensitivity discrimination threshold, filters parameters with values higher than the sensitivity discrimination threshold, and constructs a set of associated parameters.

4. The digital-twin-based fan blade stress assessment system of claim 3, wherein, The global stress field solving module comprises: A known point stress input sub-module calls real-time values in the set of associated parameters, extracts known point node coordinates of the finite element grid, binds and matches the real-time values with the known point node coordinates, establishes a mapping relationship between known point stress input and node position, and generates a stress measurement point spatial coordinate set; An unknown node interpolation solving sub-module calculates spatial distance and stress difference fitting spatial variation functions between known nodes based on the stress measurement point spatial coordinate set, and performs covariance matrix operation on unknown nodes in the finite element grid using the spatial variation functions to obtain interpolation node stress estimation values; The global data fusion generation submodule integrates two groups of data to the node database of the finite element grid according to the stress measurement point spatial coordinate set and the interpolation node stress estimation value, sorts and formats the data structure according to the node number index, and establishes the global stress distribution data.

5. The digital-twin-based fan blade stress assessment system of claim 4, wherein, The stress gradient anomaly diagnosis module comprises: The stress field data analysis submodule indexes all the discrete grid nodes in the global stress distribution data, extracts the three-dimensional spatial coordinates and stress tensor components of each node, and performs structured packaging to establish a node stress feature vector; The gradient anomaly identification submodule calculates the stress spatial gradient of each node based on the node stress feature vector, defines the arithmetic mean of the stress spatial gradient of all nodes within the statistical range of the local spatial range, and obtains the gradient deviation amount by subtracting the arithmetic mean from the stress spatial gradient of a single node. Meanwhile, the gradient deviation amount is determined according to a preset gradient deviation threshold to generate a node gradient anomaly identifier; The diagnosis set construction submodule calls the node gradient anomaly identifier, traverses all the nodes, and collects the node identifiers identified as abnormal states to generate an abnormal gradient node set.

6. The digital-twin-based fan blade stress assessment system of claim 5, wherein, The dangerous section comprehensive evaluation module comprises: The dangerous section definition submodule calls the Euclidean distance between multiple nodes in the abnormal gradient node set, performs density clustering based on a neighborhood radius threshold and a core point minimum number benchmark, filters node clusters with a density reaching a standard, and generates a dangerous section core node set; The node stress extraction submodule retrieves the global stress distribution data according to the node index in the dangerous section core node set, extracts stress state components corresponding to multiple node indexes, constructs an ordered numerical sequence, and generates a section node stress sequence; The section index weighted summation submodule calls the stress spatial gradient value as a weight factor for each stress value in the section node stress sequence, performs a weighted product calculation, and then accumulates and sums all the weighted products in the sequence to generate a section stress evaluation index.

7. The digital-twin-based fan blade stress assessment system of claim 6, wherein, The unknown node interpolation solver module specifically adopts a Gaussian model to fit the spatial variation function; The unknown node interpolation solver module calculates the model parameters of the Gaussian model using the least squares method based on the Euclidean distance between multiple known nodes in the stress measurement point spatial coordinate set and the stress value difference, and the calculation formula of the Gaussian model is: ； where, denotes the stress semi-variance value of two known nodes with spatial separation distance , denotes the three-dimensional Euclidean distance between any two known nodes, denotes the nugget value reflecting the measurement error and micro-scale spatial variability, denotes the sill value corresponding to the total spatial variability variance of the stress regionalization variable, denotes the maximum spatial autocorrelation distance between nodes when the spatial variability function reaches the sill value, i.e., the range. The unknown node interpolation solver module uses the Gaussian model with the obtained model parameters to construct a covariance matrix between known nodes and a covariance vector between unknown nodes and all known nodes, and obtains stress unbiased estimates of multiple unknown nodes by solving the Kriging equation set to establish the interpolation node stress estimation value.

8. The digital-twin-based fan blade stress assessment system of claim 7, wherein, The gradient anomaly identification submodule calls the node stress feature vector to calculate the equivalent stress value of each node; The gradient anomaly identification submodule calculates the stress change rate of each node in three coordinate axis directions using the central difference method, and calculates the modulus value of the stress spatial gradient, and the calculation formula of the gradient deviation amount is: ； wherein, denotes the gradient deviation of the node with index , denotes the modulus of the stress space gradient of the node with index , denotes the total number of nodes within the neighborhood determined according to the local spatial range centered at the node with index , denotes the sum of the modulus of the stress space gradient of all nodes within the neighborhood. The gradient anomaly identification submodule calculates the gradient deviation A point-by-point comparison is performed with the preset gradient deviation threshold, and when the gradient deviation is greater than the preset gradient deviation threshold, the node gradient anomaly flag of the corresponding node is set to abnormal.

9. The digital-twin-based fan blade stress assessment system of claim 8, wherein, The cross-section exponential weighted summation submodule calls the dangerous cross-section core node set, and for each core node in the set, extracts the equivalent stress value from the global stress distribution data, and obtains the stress space gradient value from the calculation result of the gradient anomaly identification submodule; The cross-section exponential weighted summation submodule calculates the gradient weight coefficient of each core node, which is the ratio of the stress space gradient value of a single core node to the sum of the stress space gradient values of all core nodes in the dangerous cross-section core node set; The cross-section exponential weighted summation submodule multiplies the equivalent stress value of each core node by the corresponding gradient weight coefficient to obtain the weighted stress contribution value of the core node; The cross-section exponential weighted summation submodule accumulates and sums the weighted stress contribution values of all core nodes in the dangerous cross-section core node set, through the following formula: ； Generate the cross-section stress assessment index; wherein, represents a cross-sectional stress evaluation index, represents a dangerous cross-sectional core node set, wherein, are both indexes of the dangerous cross-sectional core node set, represents a summation operation on all core nodes in the dangerous cross-sectional core node set, represents an equivalent stress value of a core node with index , represents a stress space gradient value of a core node with index , represents an accumulated summation of stress space gradient values of all core nodes in the dangerous cross-sectional core node set.

Images