A fitting blade rigidity and linear density construction method based on finite element analysis

Abstract

This invention relates to a method for constructing blade stiffness and linear density based on finite element analysis, belonging to the field of wind turbine blade technology. It solves one of the problems of lacking a universal model, which is detrimental to the early prediction and rapid evaluation of test bench schemes. The invention includes the following steps: Step S1: Calculate the linear density and stiffness of the wind turbine blade; Step S2: Normalize and preprocess the stiffness and linear density to obtain a discrete point data set; Step S3: Fit a double exponential function with 15 different initial values, use the Levenberg-Marquardt algorithm to solve for the optimal parameters, and select R... 2 The highest fitting result is taken as the initial optimal solution; Step S4: Adaptive outlier detection and removal, comparing the R-values ​​of several methods. 2 For the degree of improvement, select R. 2 The method that maximizes the improvement is taken as the final result; Step S5: Use the double exponential function model for subsequent prediction of the linear density and stiffness of ultra-long blades. This invention is used to improve the accuracy of blade physical performance curve fitting and provides a model basis for predicting the physical performance of ultra-long blades.

Description

Technical Field

[0001] This invention relates to the field of wind turbine blade technology, and in particular to a method for constructing the fitting blade stiffness and linear density based on finite element analysis. Background Technology

[0002] As the global energy structure shifts towards cleaner and lower-carbon energy, wind power, as a technologically mature and economically viable renewable energy source, has become an important component of my country's energy structure. The main development direction of the industry is towards larger and higher-power wind turbines, with single-unit capacity increasing from megawatts to tens of megawatts, and the length of wind turbine blades continuously increasing, with mainstream models reaching 120-150 meters in length. According to international and domestic practices, to ensure quality, wind turbine blades require full-scale structural testing.

[0003] However, the number and capacity of existing test rigs both domestically and internationally are significantly insufficient. The largest testing capability abroad is for blades up to 115 meters long, while my country's test rigs, starting in 2014 / 2015, could only test blades no longer than 100 meters, only gradually improving to the ability to test 150-160 meter long blades in recent years. The current goal is to break through the key technologies of test rigs for 160-180 meter blades. Therefore, it is necessary to estimate the parameters of 180-meter blades and design test rigs based on this, thereby enabling the testing of large-sized blades. However, the required blade data can only be obtained from the main engine manufacturer, making data acquisition difficult. The obtained data is also limited to specific models and sizes, lacking a universal model when extrapolating to larger blades (such as 180-meter blades). Furthermore, the existing complete blade modeling process is complex and involves a large workload, hindering the early prediction and rapid evaluation of test rig designs. Summary of the Invention

[0004] In view of the above analysis, the present invention aims to provide a method for constructing fitting blade stiffness and linear density based on finite element analysis, in order to solve one of the problems of lack of a general model when extrapolating to larger blades, the complexity and large amount of engineering work in the existing complete blade modeling process, and the unfavorability to the early prediction and rapid evaluation of test bench schemes.

[0005] This invention provides a method for constructing the fitting blade stiffness and linear density based on finite element analysis, comprising the following steps: Step S1: Based on finite element analysis, perform dimensionality reduction modeling, calculate the linear density of the wind turbine blades, and set a rigid plane to calculate the blade stiffness; Step S2: Normalize and preprocess the calculated blade stiffness and linear density to obtain a discrete point dataset; Step S3: Fit a double exponential function with 15 different initial values, use the Levenberg-Marquardt algorithm to solve for the optimal parameters, and select R... 2 The highest fitting result is taken as the preliminary optimal solution; Step S4: Adaptive outlier detection and removal, comparing the R-squared values ​​of the standard deviation method, interquartile range method, and percentile method. 2 For the degree of improvement, select R. 2 The method that maximizes improvement is taken as the final result; Step S5: Use the double exponential function model for subsequent prediction of the linear density and stiffness of ultra-long blades.

[0006] Further, step S1 includes the following steps: Step S1.1: Based on the multi-source wind turbine blade model data, establish a three-dimensional physical model of the wind turbine blade to be analyzed, and extract the wind turbine blade axis and the geometric dimensions and material parameters of each cross section distributed along the axis; Step S1.2: In the finite element analysis software, the wind turbine blade is equivalent to a variable cross-section beam element model arranged along the axis; the wind turbine blade axis is divided into several short beam elements according to a preset unit length. Step S1.3: Calculate the mass, linear density value and gravitational moment of each short beam unit based on the three-dimensional physical model of the wind turbine blade.

[0007] Furthermore, step S1 also includes the following steps: Step S1.4: Set the short beam element near the leaf root as the fixed end and the end away from the leaf root as the cantilever end. Apply rigid planar constraints to the cantilever end to form simplified cantilever beam boundary conditions. Step S1.5: Calculate the swing stiffness, sway stiffness and torsional stiffness of each short beam element.

[0008] Furthermore, in step S1.5, unit rotation angles around the swing direction and the oscillation direction are applied to the cantilever end respectively, and the corresponding fixed end bending moments Mx and My are solved by finite element software. Based on the pre-set short beam length, the section moment of inertia Ix and Iy, and the elastic modulus E, the swing stiffness and oscillation stiffness of each short beam unit are calculated.

[0009] Furthermore, in step S1.5, a unit rotation angle about the torsional direction is applied to the cantilever end, and the fixed end torque Mz is solved by finite element software. Based on the short beam length, the polar moment of inertia Ip of the cross section and the shear modulus G, the torsional stiffness of each short beam element is calculated.

[0010] Further, step S2 includes the following steps: Step S2.1: Normalize the axial position of each short beam unit obtained in step S1 according to the full length of the wind turbine blade, and map the axial length of the wind turbine blade to the interval [0,1] to obtain the corresponding dimensionless position coordinates; Step S2.2: Normalize the amplitude of the linear density and linear stiffness in each direction of each short beam element to obtain dimensionless data that is easy to fit numerically; export the normalized axial position, linear density and linear stiffness as a CSV file in tabular form. Step S2.3: Read the CSV file and related original parameters, and automatically remove invalid values ​​to ensure that the number of valid data points is not less than the preset threshold, so as to increase the degree of freedom for subsequent fitting. Step S2.4: Establish a correspondence between the normalized axial position and the normalized linear density and linear stiffness to form a discrete point dataset for fitting.

[0011] Further, step S3 includes the following steps: Step S3.1: Select the double exponential function as the fitting model for the distribution of wind turbine blade stiffness or linear density along the axial length; Step S3.2: Construct multiple sets of initial parameter sets; Step S3.3: For each set of initial parameters, the Levenberg-Marquardt nonlinear least squares algorithm is used to fit and solve the double exponential function model. Under the premise of satisfying the preset convergence conditions and the maximum number of iterations constraint, the corresponding fitting parameters are obtained.

[0012] Furthermore, step S3 also includes step S3.4: calculating the coefficient of determination R² and root mean square error RMSE for each group of fitting results, and using R² as the main evaluation index, selecting the group of parameters with the highest R² value as the preliminary optimal solution of the double exponential function model.

[0013] Further, step S5 includes the following steps: Step S5.1: Substitute the final fitting parameters obtained in step S4 into the double exponential function model to obtain the continuous distribution equation of the wind turbine blade linear density or stiffness as a function of normalized axial position. Step S5.2: Based on the actual length of the target wind turbine blade and the reference linear density and stiffness level, scale the distribution equation to restore or scale it, and map the normalized position to the true axis coordinates to obtain the predicted values ​​of linear density and stiffness of the target wind turbine blade at any axial position. Step S5.3: Use the prediction results for preliminary estimation of the structural strength and driving capability of the test bench, selection of design parameters for ultra-large blade structures, and related performance evaluation.

[0014] Furthermore, the parameters of the wind turbine blades are simplified to linear density, flapping stiffness, oscillation stiffness, and torsional stiffness.

[0015] Compared with the prior art, the present invention can achieve at least one of the following beneficial effects: (1) The present invention adopts a multi-initial value strategy with 15 sets of systematic design, covering various double exponential decay types, with strong adaptability, able to identify and remove outliers, and improve fitting accuracy; it can fit the physical performance curves of blades of different lengths in the field of blade design, identify abnormal designs, provide a reliable basis for production quality control, and provide a model basis for predicting the physical performance of ultra-long blades, solving the problem of lacking a general model when extrapolating to larger blades.

[0016] (2) The model of the present invention has strong versatility: the double exponential function model combined with normalization processing is used to unify the stiffness and linear density distribution of blades of different lengths into a dimensionless coordinate system, and the soft boundary constraint ensures that the fitting curve conforms to the physical characteristics of the blade, so that the model does not depend on specific models and sizes, and can be applied to the parameter prediction of blades of different lengths, solving the problem of lack of a universal model when extrapolating to larger blades.

[0017] (3) The present invention has high modeling efficiency: by simplifying the complex blade structure into the distribution fitting of four key parameters (linear density, flapping stiffness, oscillation stiffness and torsional stiffness), the blade parameter distribution model can be established without obtaining complete detailed design data or measured data, which significantly simplifies the modeling process, facilitates the early prediction and rapid evaluation of test bench schemes, and solves the problems that the existing complete blade modeling process is relatively complex, has a large workload, and is not conducive to the early prediction and rapid evaluation of test bench schemes.

[0018] (4) The present invention has high fitting accuracy and strong robustness: it adopts a multi-initial value strategy of 15 sets combined with the Levenberg-Marquardt algorithm for global optimization to avoid getting trapped in local optima; at the same time, it adopts multiple threshold judgment methods to automatically identify and remove abnormal data points, and improves fitting accuracy through secondary fitting optimization, providing a reliable basis for predicting the structural parameters of ultra-long blades.

[0019] In this invention, the above-described technical solutions can be combined with each other to achieve more preferred combinations. Other features and advantages of this invention will be set forth in the following description, and some advantages may become apparent from the description or be learned by practicing the invention. The objects and other advantages of this invention can be realized and obtained from what is particularly pointed out in the description and drawings. Attached Figure Description

[0020] The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. Throughout the drawings, the same reference numerals denote the same parts.

[0021] Figure 1 This is a schematic diagram of the process for fitting the blade stiffness and linear density based on finite element analysis in a specific embodiment. Figure 2This is a schematic diagram of the force analysis of the wind turbine blades in a specific embodiment; Figure 3 This is a schematic diagram of the residual distribution during the fitting process of the normalized curve of the linear density distribution with respect to the blade length in a specific embodiment. Figure 4 This is a schematic diagram of the residual histogram during the fitting process of the normalized curve of the linear density distribution with respect to the blade length in a specific embodiment. Figure 5 This is a schematic diagram of the double exponential fitting process in the normalized curve fitting process of the linear density distribution with blade length in a specific embodiment. Figure 6 This is a schematic diagram of the residual distribution during the fitting process of the normalized curve of the flapping stiffness as a function of blade length in a specific embodiment. Figure 7 This is a schematic diagram showing the comparison of residuals during the fitting process of the normalized curve of the flapping stiffness distribution with blade length in a specific embodiment. Figure 8 This is a schematic diagram showing the fitting comparison during the fitting process of the normalized curve of the flapping stiffness distribution with blade length in a specific embodiment. Figure 9 This is a schematic diagram of the final fitting result in the process of fitting the normalized curve of the flapping stiffness distribution with blade length in a specific embodiment. Figure 10 This is a schematic diagram of the residual histogram during the fitting process of the normalized curve of the distribution of shimmy stiffness with blade length in a specific embodiment. Figure 11 This is a schematic diagram of the residual distribution during the fitting process of the normalized curve of the oscillation stiffness as a function of blade length in a specific embodiment. Figure 12 This is a schematic diagram of the double exponential fitting process in the normalized curve fitting process of the oscillation stiffness distribution with blade length in a specific embodiment. Figure 13 This is a schematic diagram of the residual distribution during the fitting process of the normalized curve of the torsional stiffness distribution with respect to the blade length in a specific embodiment. Figure 14 This is a schematic diagram comparing the residuals during the fitting process of the normalized curve of the torsional stiffness distribution with blade length in a specific embodiment. Figure 15 This is a schematic diagram showing the fitting comparison during the normalized curve fitting process of the torsional stiffness distribution with blade length in a specific embodiment. Figure 16 This is a schematic diagram of the final fitting result in the process of fitting the normalized curve of the torsional stiffness distribution with respect to the blade length in a specific embodiment. Detailed Implementation

[0022] The preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings, which constitute a part of the present invention and are used together with the embodiments of the present invention to illustrate the principles of the present invention, but are not intended to limit the scope of the present invention.

[0023] A specific embodiment of the present invention, such as Figure 1 As shown, a method for constructing blade stiffness and linear density based on finite element analysis is disclosed. Specifically, it is a method for constructing the blade stiffness and linear density along the length distribution based on finite element analysis, including the following steps: Step S1: Based on finite element analysis, perform dimensionality reduction modeling, calculate the linear density of the wind turbine blades, and set a rigid plane to calculate the stiffness of the wind turbine blades.

[0024] Step S2: Normalize the calculated stiffness and linear density of the wind turbine blades to obtain a discrete point data set.

[0025] Step S3: Fit the double exponential function with 15 different initial values, and use the Levenberg-Marquardt (LM) algorithm to solve for the optimal parameters. Select R0... 2 The highest fitting result is taken as the preliminary optimal solution.

[0026] In step S3, under small sample conditions, the risk of getting trapped in local optima is reduced by using "multiple initial values ​​+ LM" to obtain a benchmark solution, providing a unified reference for subsequent residual and outlier judgment.

[0027] Step S4: Adaptive outlier detection and removal, comparing the R-squared values ​​of the standard deviation method, interquartile range method, and percentile method. 2 For the degree of improvement, select R. 2 The method that maximizes improvement is taken as the final result.

[0028] Step S5: Use the double exponential function model for subsequent prediction of the linear density and stiffness of ultra-long blades.

[0029] Compared with existing technologies, the method for constructing blade stiffness and linear density based on finite element analysis provided in this embodiment adopts a multi-initial value strategy with 15 sets of systematically designed systems, covering various double exponential decay types. It has strong adaptability, can identify and remove outliers, and improves fitting accuracy. It can fit the physical performance curves of blades of different lengths in the field of blade design, identify abnormal designs, provide a reliable basis for production quality control, and provide a model basis for predicting the physical performance of ultra-long blades, solving the problem of lacking a general model when extrapolating to larger blades.

[0030] For step S1: Based on finite element analysis, dimensionality reduction modeling is performed to calculate the linear density of the wind turbine blades and to calculate the stiffness of the wind turbine blades by setting a rigid plane; in other words, dimensionality reduction modeling and cross-sectional parameter calculation are performed. Step S1 specifically includes the following steps: Step S1.1: Based on the multi-source wind turbine blade model data, establish a three-dimensional physical model of the wind turbine blade to be analyzed, and extract the wind turbine blade axis and the geometric dimensions and material parameters of each section distributed along the axis. The geometric parameters include the overall cross-sectional dimensions, chord length, main beam and web dimensions, skin thickness, fiber cloth layup sequence and quantity, and the material parameters include the material elastic modulus and material density.

[0031] Step S1.2: In the finite element analysis software, the wind turbine blade is equivalent to a variable cross-section beam element model arranged along the axis. The wind turbine blade axis is divided into several short beam elements with a preset unit length and a length not exceeding 1m. It is assumed that the cross-section of each short beam element remains basically unchanged in the length direction.

[0032] Step S1.3: Based on the three-dimensional physical model of the wind turbine blade, calculate the mass of each short beam unit, and obtain the linear density value of each short beam unit by dividing the mass by the unit length. Calculate the corresponding gravity and center of gravity position through integration, and then calculate the gravitational torque.

[0033] Step S1.4: Set the short beam element near the leaf root as the fixed end and the end away from the leaf root as the cantilever end. Apply rigid planar constraints to the cantilever end to form simplified cantilever beam boundary conditions.

[0034] Step S1.5: Apply unit rotation angles around the swinging and oscillation directions to the cantilever end respectively, and solve the corresponding fixed end bending moments using finite element software. M x , M y Based on the pre-set short beam length and cross-sectional moment of inertia I x , I y and elastic modulus E The swing stiffness and oscillation stiffness of each short beam element were calculated, as shown in equations (I) and (II). Figure 2 This is a schematic diagram of the force analysis of the wind turbine blades.

[0035] (I) (II) in, M x The bending moment at the fixed end in the direction of the oscillation; M y Fixed end bending moment in the direction of waving; I x Wind turbine blade cross section x Moment of inertia of the axis; I y For the blade section yMoment of inertia of the axis; θ The unit rotation angle in the swinging or waving direction is always 1; E It is the elastic modulus; l This is the length of the short beam.

[0036] A unit rotation angle about the torsional direction is applied to the cantilever end, and the torque at the fixed end is solved using finite element software. M z Based on the short beam length and the polar moment of inertia of the cross section I p and shear modulus G The torsional stiffness of each short beam element is calculated, as shown in Equation (Ⅲ). This yields the linear stiffness data distributed along the blade axis.

[0037] (III) in, M z This refers to the torque of the wind turbine blades; I p Let be the polar moment of inertia of the wind turbine blade section; θ The unit angle of rotation in the direction of twisting is always 1; G It is the shear modulus; l This is the length of the short beam.

[0038] For step S2: The calculated stiffness and linear density of the wind turbine blades are normalized and preprocessed to obtain a discrete point data set; in other words, data normalization and discrete point construction. This solves the problem of lacking a general model when extrapolating to larger blade sizes. Step S2 specifically includes the following steps: Step S2.1: Normalize the axial position of each short beam unit obtained in step S1 according to the full length of the wind turbine blade, and map the axial length of the wind turbine blade to the interval [0,1] to obtain the corresponding dimensionless position coordinates.

[0039] Step S2.2: Normalize the amplitude of the linear density and linear stiffness in each direction of each short beam element to obtain dimensionless data that is easy to fit numerically; export the normalized axial position, linear density and linear stiffness as a CSV file in tabular form.

[0040] Step S2.3: Read the CSV file and automatically remove invalid values, including "null", "non-numeric", "infinity", etc., to ensure that the number of valid data points is not less than the preset threshold (e.g., not less than 6) to ensure the degree of freedom of subsequent fitting.

[0041] Step S2.4: Establish a correspondence between the normalized axial position and the normalized linear density and linear stiffness to form a discrete point dataset for fitting.

[0042] To address the lack of a universal model when extrapolating to larger blade sizes, step S3 of this embodiment involves fitting a double exponential function with 15 different initial values, using the Levenberg-Marquardt algorithm to solve for the optimal parameters, and selecting R... 2 The highest fitting result is taken as the preliminary optimal solution; in other words, the preliminary solution is obtained by fitting the double exponential function. Step S3 specifically includes the following steps: Step S3.1: Select a double exponential function as the fitting model for the distribution of wind turbine blade stiffness or linear density along the axial length. The fitting model is as follows: (IV) in, x To normalize the axial position, y This corresponds to the normalized linear density or stiffness value. a , b , c , d These are the model parameters to be determined.

[0043] Step S3.2: Construct multiple sets of initial parameters.

[0044] The initial parameter set covers various typical forms of the double exponential function, including the standard double-decaying type, the fast-slow combination type, the asymmetric type, the extreme decay rate type, and the single-exponential approximation type, to improve the convergence stability of the nonlinear fitting. Specifically, the standard double-decaying type has parameters such as a≈0.5, b≈-2, c≈0.5, and d≈-2; the fast-slow combination type has a large difference between b and d; the asymmetric type has a large difference between a and c; the extreme type has an extremely fast or extremely slow decay rate; and the single-exponential approximation type has a value close to 0 for either a or c.

[0045] Step S3.3: For each set of initial parameters, the Levenberg-Marquardt nonlinear least squares algorithm is used to fit and solve the double exponential function model. Under the premise of satisfying the preset convergence conditions and the maximum number of iterations constraint, the corresponding fitting parameters a, b, c, and d are obtained. During the fitting process, soft boundary constraints can be optionally introduced to make the function value at the normalized position x=0 close to 1 and the function value at x=1 close to 0, so as to reflect the relative changing trend of linear density or stiffness at the blade root and blade tip.

[0046] Step S3.4: Calculate the coefficient of determination R² and RMSE (Root Mean Squared Error) for each group of fitting results, and use R² as the main evaluation index to select the group of parameters with the highest R² value as the preliminary optimal solution of the double exponential function model.

[0047] For step S4: Adaptive outlier detection and removal, comparing the R-values ​​of the standard deviation method, interquartile range method, and percentile method.2 For the degree of improvement, select R. 2 The method that maximizes improvement is taken as the final result; in other words: adaptive outlier detection and result correction.

[0048] To improve the fitting accuracy, step S4 specifically includes the following steps: Step S4.1: Based on the preliminary optimal solution obtained in step S3, calculate the fitted value of each data point using the double exponential function model, and use the difference between the fitted value and the actual data value as the residual to form a residual sequence.

[0049] Step S4.2: For the residual sequence, various outlier detection methods are used to determine candidate outliers. These outlier detection methods include thresholding based on the mean and standard deviation, thresholding based on the interquartile range, and thresholding based on residual percentiles. Different threshold combinations are set for each method.

[0050] Step S4.3: Under each detection method and its threshold combination, data points with residuals greater than the corresponding threshold are marked as outliers. The proportion of the number of data points removed to the total number of data points is constrained to be within a preset range (e.g., 1% to 20%), and the number of data points retained is not less than the minimum amount of data required for fitting.

[0051] Step S4.4: Update the dataset with marked outliers, remove outliers, and re-execute the double exponential function fitting to obtain the corrected fitting parameters and new R² and RMSE indices for each detection method.

[0052] Step S4.5: Calculate the R² improvement of each detection method relative to the preliminary optimal solution. When the R² improvement of a certain detection method is not lower than the preset threshold, the fitting result corresponding to the method is recorded as an acceptable result. Among all acceptable results, the fitting result with the largest R² improvement is selected as the final fitting result. If the R² improvement of each method is not significant, the preliminary optimal solution obtained in step S3 is retained as the final fitting result.

[0053] In step S4, to address the differences in patterns caused by different wind turbine blades, a "multi-strategy outlier handling + adaptive selection based on ΔR²" approach is used to achieve robustness and adaptability to heterogeneous data, and the final fitted curve is output.

[0054] For step S5: the double exponential function model is used for subsequent prediction of the linear density and stiffness of the ultra-long blades; in other words, the distribution equation is applied. Step S5 specifically includes the following steps: Step S5.1: Substitute the final fitting parameters obtained in step S4 into the double exponential function model to obtain the continuous distribution equation of the wind turbine blade linear density or stiffness as a function of normalized axial position.

[0055] Step S5.2: Based on the actual length of the target wind turbine blade and the reference linear density and stiffness level, scale the distribution equation to restore or scale it, and map the normalized position to the true axis coordinates to obtain the predicted values ​​of linear density and stiffness of the target wind turbine blade at any axial position.

[0056] Step S5.3: Use the prediction results for preliminary estimation of the structural strength and driving capability of the test bench, selection of design parameters for ultra-large blade structures, and related performance evaluation, so as to provide a reference for engineering design.

[0057] Compared with existing technologies, this embodiment uses 15 systematically designed multi-initial-value strategies to cover various double-exponential decay types. It exhibits strong adaptability, can identify and remove outliers, and improves fitting accuracy. It can fit the physical performance curves of blades of different lengths in the field of blade design, identify abnormal designs, provide a reliable basis for production quality control, and provide a model foundation for predicting the physical performance of ultra-long blades.

[0058] The model in this embodiment has strong versatility: it adopts a double exponential function model combined with normalization processing to unify the stiffness and linear density distribution of blades of different lengths into a dimensionless coordinate system, and ensures that the fitted curve conforms to the physical characteristics of the blade through soft boundary constraints. This makes the model independent of specific models and sizes, and applicable to parameter prediction of blades of different lengths, solving the problem of lacking a universal model when extrapolating to larger blades.

[0059] This embodiment boasts high modeling efficiency: the complex blade structure is simplified to four key parameters—linear density, flapping stiffness, teeter stiffness, and torsional stiffness—with a certain inherent relationship between the three stiffnesses. By establishing the relationship between the distribution of these four key parameters and the blade length, a simplified design method based on beam elements for blade finite element models can be achieved, improving the efficiency of blade design and analysis and simplifying the modeling process. By simplifying the complex blade structure to the distribution fitting of four key parameters (linear density, flapping stiffness, teeter stiffness, and torsional stiffness), a blade parameter distribution model can be established without obtaining complete and detailed design data or measured data, significantly simplifying the modeling process and facilitating early prediction and rapid evaluation of test bench schemes.

[0060] This embodiment features high fitting accuracy and robustness: it employs a multi-initial-value strategy with 15 sets of initial values ​​combined with the Levenberg-Marquardt algorithm for global optimization, avoiding getting trapped in local optima; at the same time, it uses multiple threshold determination methods to automatically identify and remove abnormal data points, and improves fitting accuracy through secondary fitting optimization, providing a reliable basis for predicting the structural parameters of ultra-long blades.

[0061] This embodiment provides a solution based on existing data. The method can still predict the distribution of key structural parameters of blades relatively accurately even in the absence of complete and detailed design data or measured data, providing a reference for test bench design, tower design, blade structure design and performance evaluation.

[0062] Example 1 A specific embodiment of the present invention, such as Figure 1 , Figure 3 , Figure 4 and Figure 5 As shown, the normalized curve fitting of linear density distribution with wind turbine blade length is presented, including... Figure 1 The process specifically includes the following steps: Step S1: Dimensionality reduction modeling and line density calculation.

[0063] Based on the three-dimensional model of the multi-source wind turbine blades, the geometric dimensions and material parameters of the blade axis and the cross-section distributed along the axis were extracted. In the finite element software ANSYS, the six wind turbine blades were equivalent to variable cross-section beam elements, and a total of 384 short beam elements were divided along the axis by unit length, so that the cross-section of each short beam element remained basically unchanged within its length. The mass was calculated based on the material density and cross-sectional area of ​​each short beam element, and the corresponding linear density data was obtained by dividing by the element length, forming a linear density sequence distributed along the axis.

[0064] Step S2: Normalization and construction of discrete point dataset.

[0065] The blade axis length was normalized to the [0,1] interval, resulting in 384 dimensionless position coordinates x. The linear density values ​​were then normalized to distribute them within the [0,1] interval and exported as a CSV file in the form of "position x - normalized linear density y". The CSV file was then read, and non-existent values ​​(NaN, Inf, etc.) were removed. The number of valid data points was confirmed to be 384, which is greater than the preset minimum data point threshold, thus forming a discrete point dataset for fitting.

[0066] Step S3: Fitting a double exponential function with soft boundary constraints.

[0067] In this embodiment, a double exponential function is selected as the fitting model for the linear density normalized distribution: Where x is the normalized axial position, y is the corresponding normalized linear density value, and a, b, c, and d are the model parameters to be determined.

[0068] To reflect the relative variation trend of leaf linear density at the leaf root and leaf tip, a soft boundary constraint is introduced during the fitting process: at the normalized position x=0, the function value is expected to be close to 1, and at x=1, the function value is expected to be close to 0. This constraint is then added to the residual as a penalty term by setting the boundary constraint weight (e.g., 40).

[0069] An initial parameter set covering various typical morphologies was constructed, and the Levenberg-Marquardt nonlinear least squares algorithm was used to fit all 384 data points. After multiple solutions under different initial parameters, a set of parameters that satisfied the boundary constraints and had good fitting quality was obtained: a = 0.761853574127964; b = 160.442611549335737;c=0.238152643517669;d= 3.443626710358401.

[0070] At this point, the coefficient of determination R² for the normalized linear density data is 0.94337959, and the root mean square error RMSE is 0.0305025659. At the boundary, the function values ​​satisfy f(0) = 1.000006 and f(1) = 0.007609, with deviations from the target values ​​of 1.0 and 0.0 of 0.000006 and 0.007609 respectively, meeting the boundary accuracy requirements set in this embodiment. This result serves as the preliminary fitting result for the linear density double exponential model.

[0071] Step S4: Adaptive outlier detection and fitting correction.

[0072] Based on the preliminary results above, the residuals of each data point were calculated, and various methods such as standard deviation (3σ, 2.5σ), residual percentiles (95th percentile, 90th percentile), and interquartile range (IQR) were used to threshold the residuals, attempt to mark outliers and remove them before refitting.

[0073] Under the data conditions of this embodiment, if various threshold combinations are directly applied, all 384 data points will be judged as outliers, that is, the proportion of outliers reaches 100%, which does not meet the preset constraint that "the proportion of outliers is between 1% and 20%, and the number of retained data points is not less than the minimum amount of data required for fitting". Therefore, the above outlier detection results are all judged as invalid.

[0074] Therefore, in this embodiment, outlier removal is not performed on the original data, and the fitting result obtained in step S3 is directly retained as the final result.

[0075] Step S5: Construct the linear density distribution equation and its application.

[0076] Substituting the final fitting parameters from step S3 into the double exponential function model, we obtain the distribution equation of the linear density as a function of the normalized axial position: y ρ (x)=0.7618535741exp(-160.4426115493x)+0.2381526435exp(-3.4436267104x) Where x∈[0,1] is the normalized axial position of the blade, and y ρ (x) represents the corresponding normalized linear density.

[0077] In practical applications, the above distribution equation can be scaled back by combining the actual length of the target wind turbine blade and the reference linear density level, and the normalized position can be mapped to the true axis coordinates, thereby obtaining the predicted value of linear density at any axial position. Figure 3 , Figure 4 and Figure 5 The normalized fitting curve of linear density distribution with blade length in this embodiment is shown, which can provide a reference for mass distribution estimation in the design of ultra-large wind turbine blade test benches.

[0078] Example 2 A specific embodiment of the present invention, such as Figure 1 , Figure 6 , Figure 7 , Figure 8 and Figure 9 As shown, in this embodiment... Figure 1 Based on the method flow shown, the target physical quantity is set as flapping stiffness, which is used for fitting the normalized curve of the flapping stiffness distribution with blade length. The method still includes steps S1-S5, but only the specific implementation and fitting results related to flapping stiffness are explained. Specifically, it includes the following steps: Step S1: Dimensional reduction modeling and line swing stiffness calculation.

[0079] Based on the three-dimensional model of the multi-source wind turbine blades, the geometric dimensions and material parameters of the blade axis and the cross-sections distributed along the axis are extracted. In the finite element analysis software ANSYS, the nine wind turbine blades are equivalent to variable cross-section beam elements. The blades are divided into several short beam elements along the blade axis at a preset unit length. In this embodiment, 345 short beam elements are obtained.

[0080] The end closer to the blade root is designated as the fixed end, and the end farther from the blade root is designated as the cantilever end. Rigid planar constraints are applied to the cantilever end, and a unit rotation angle is applied around the flapping direction. The fixed end bending moment in the flapping direction of each short beam element is obtained through finite element analysis. The corresponding linear flapping stiffness is calculated by combining the short beam length, section moment of inertia, and elastic modulus, forming a flapping stiffness sequence distributed along the axis.

[0081] Step S2: Normalization and construction of discrete point dataset.

[0082] The blade axis length was normalized to the [0,1] interval, resulting in 345 dimensionless position coordinates x. The amplitude of the linear flapping stiffness was normalized so that the normalized stiffness y was distributed in the [0,1] interval. The "position x - normalized flapping stiffness y" data was exported as a CSV file, and the file was read to remove any invalid values ​​(NaN, Inf, etc.). The number of valid data points was confirmed to be 345, which is greater than the minimum number of data points required for fitting, thus forming a discrete point dataset for fitting.

[0083] Step S3: Fitting a double exponential function with soft boundary constraints (using all data).

[0084] In this embodiment, a double exponential function is selected as the fitting model for the normalized distribution of swing stiffness: Where x is the normalized axial position, y is the corresponding normalized swing stiffness value, and a, b, c, and d are the model parameters to be determined.

[0085] To reflect the relative change trend of flapping stiffness at the leaf root and leaf tip, a soft boundary constraint is introduced during the fitting process: at the normalized position x=0, the function value is expected to be close to 1, and at x=1, the function value is expected to be close to 0. The constraint is then superimposed on the residual function as a penalty term by setting the boundary constraint weight (40 in this embodiment).

[0086] Based on this, an initial parameter set covering various typical forms was constructed, and the Levenberg-Marquardt nonlinear least squares algorithm was used to fit all 345 data points. After solving multiple sets of initial values, a set of parameters that satisfied the boundary constraints and had good fitting quality was obtained: a = 0.34792480; b = 7.74732955;c=0.65209330;d= 87.76200466.

[0087] At this point, the coefficient of determination R² = 0.868113, and the root mean square error RMSE = 0.06915986. At the boundary, the function values ​​satisfy f(0) = 1.000018 and f(1) = 0.000150, with deviations from the target values ​​of 1.0 and 0.0 of 0.000018 and 0.000150, respectively, which meets the boundary accuracy requirements set in this embodiment. This result serves as the preliminary fitting result of the double exponential model of swing stiffness.

[0088] Step S4: Adaptive outlier detection and fitting correction.

[0089] Based on the preliminary results above, the residuals of each data point were calculated, and candidate outliers were identified using a variety of outlier detection methods, including: 3σ and 2.5σ thresholding methods based on standard deviation; 95th and 90th percentile thresholding methods based on residual percentiles; and thresholding methods based on interquartile range (IQR).

[0090] For each thresholding method, provided that the proportion of outliers is limited (e.g., between 1% and 20%) and the number of retained data points is not less than the minimum number required for fitting, data points with residuals exceeding the threshold are marked as outliers and temporarily deleted. The remaining data are then used to re-perform the double exponential function fitting with boundary soft constraints to obtain new parameters and their R² and RMSE indices.

[0091] In this embodiment, the representative results of each method are as follows: 3σ method: Remove 10 data points, accounting for 2.90% of the total, and after refitting, R²≈0.981104; 2.5σ method: Remove 11 data points, accounting for 3.19% of the total, and after refitting, R²≈0.981221; 95th percentile method: Remove 18 data points, accounting for 5.22% of the total, and after refitting, R²≈0.988723; 90th percentile method: Remove 35 data points, accounting for 10.14% of the total, and after refitting, R²≈0.995855; IQR method: Remove 60 data points, accounting for 17.39% of the total, and after refitting, R²≈0.998926.

[0092] Under the constraint of outlier ratio, the IQR method yielded the largest improvement in R², increasing from 0.868113 to 0.998926 compared to the result using all data, an improvement of approximately +0.1308 (+15.07%). The RMSE decreased from 0.06915986 to 0.00545575, a reduction of approximately 0.0637 (an improvement of 92.11%). Therefore, the fitting result after removing 60 outliers using the IQR method is taken as the final fitting result in this embodiment.

[0093] In this case, 285 data points were retained, and the final fitting parameters were: a = 0.298987194119906; b = 7.228965832556773;c=0.701012162295850;d= 63.103354418145003. Corresponding R²=0.99892555, RMSE=0.0054557495. At the boundary, the function values ​​are f(0)=0.999999 and f(1)=0.000217, with deviations from the target values ​​of 1.0 and 0.0 of 0.000001 and 0.000217 respectively, satisfying the boundary soft constraint requirements.

[0094] Step S5: Construction and application of the swing stiffness distribution equation.

[0095] Substituting the final fitting parameters determined in step S4 into the double exponential function model, we obtain the distribution equation of the swing stiffness as a function of the normalized axial position: y kf (x)=0.2989871941exp(-7.2289658326x)+0.7010121623exp(-63.1033544181x) Where x∈[0,1] is the normalized axial position of the blade, and y kf (x) represents the normalized line swing stiffness at the corresponding position.

[0096] In engineering applications, the above distribution equation can be scaled back by combining the actual length of the target blade and the reference flapping stiffness level, and the normalized position can be mapped to the real axis coordinates to obtain the predicted flapping stiffness value at any axial position. Figure 6 , Figure 7 , Figure 8 and Figure 9 The normalized fitting curve of the flapping stiffness distribution with blade length in this embodiment is shown, which can provide a reference for the load and deformation analysis of ultra-large wind turbine blades in the flapping direction.

[0097] Example 3 A specific embodiment of the present invention, such as Figure 1 , Figure 10 , Figure 11 and Figure 12 As shown, in this embodiment... Figure 1 Based on the method flow shown, the target physical quantity is set as the teeter stiffness, which is used for fitting the normalized curve of the teeter stiffness distribution with respect to blade length. The method also includes steps S1-S5. Only the specific implementation and fitting results related to the teeter stiffness are explained, which specifically include the following steps: Step S1: Dimensional reduction modeling and calculation of linear oscillation stiffness.

[0098] Based on the three-dimensional model of the multi-source wind turbine blades, the geometric dimensions and material parameters of the blade axis and the cross-sections distributed along the axis are extracted. In the finite element analysis software ANSYS, the seven wind turbine blades are equivalent to variable cross-section beam elements. The blades are divided into several short beam elements along the blade axis at a preset unit length. In this embodiment, 301 short beam elements are obtained.

[0099] The end closer to the blade root is designated as the fixed end, and the end farther from the blade root is designated as the cantilever end. A rigid planar constraint is applied to the cantilever end, and a unit rotation angle is applied around the oscillation direction. The fixed end bending moment in the oscillation direction of each short beam element is obtained through finite element analysis. The corresponding linear oscillation stiffness is calculated by combining the short beam length, section moment of inertia, and elastic modulus, forming an oscillation stiffness sequence distributed along the axis.

[0100] Step S2: Normalization and construction of discrete point dataset.

[0101] The blade axis length was normalized to the [0,1] interval, resulting in 301 dimensionless position coordinates x. The linear oscillation stiffness was normalized to ensure that the normalized stiffness y was distributed within the [0,1] interval. The "position x - normalized oscillation stiffness y" data was exported as a CSV file and read. Invalid values ​​(NaN, Inf, etc.) were removed, confirming that there were 301 valid data points, which is greater than the minimum number of data points required for fitting. This constitutes the discrete point dataset for fitting.

[0102] Step S3: Fitting a double exponential function with soft boundary constraints (using all data).

[0103] In this embodiment, a double exponential function is selected as the fitting model for the normalized distribution of oscillation stiffness: Where x is the normalized axial position, y is the corresponding normalized oscillation stiffness value, and a, b, c, and d are the model parameters to be determined.

[0104] To reflect the relative change trend of the oscillation stiffness at the blade root and blade tip, a soft boundary constraint is introduced during the fitting process: at the normalized position x=0, the function value is expected to be close to 1, and at x=1, the function value is expected to be close to 0. The constraint is then superimposed on the residual function as a penalty term by setting the boundary constraint weight (40 in this embodiment).

[0105] Based on this, an initial parameter set covering various typical forms was constructed, and the Levenberg-Marquardt nonlinear least squares algorithm was used to fit all 301 data points. After solving multiple sets of initial values, a set of parameters that satisfied the boundary constraints and had good fitting quality was obtained: a = 0.698393148763837; b = 261.688231917459063;c=0.301607346132748;d= 4.438066711268386.

[0106] At this point, the coefficient of determination R² = 0.97166546, and the root mean square error RMSE = 0.0262654362. At the boundary, the function values ​​satisfy f(0) = 1.000000 and f(1) = 0.003565, with deviations from the target values ​​of 1.0 and 0.0 of 0.000000 and 0.003565, respectively, meeting the boundary accuracy requirements set in this embodiment. This result serves as the preliminary fitting result of the bi-exponential model of oscillation stiffness.

[0107] Step S4: Adaptive outlier detection and result determination.

[0108] Based on the preliminary results above, the residuals of each data point were calculated, and various outlier detection methods, such as standard deviation (3σ, 2.5σ), residual percentiles (95th percentile, 90th percentile), and interquartile range (IQR), were used to threshold the residuals, attempt to mark outliers, remove them, and then refit.

[0109] Under the data conditions of this embodiment, if various threshold methods are applied directly, all 301 data points will be identified as outliers, that is, the proportion of outliers will reach 100%, which does not meet the preset constraint that "the proportion of outliers is between 1% and 20%, and the number of data points retained is not less than the minimum amount of data required for fitting". Therefore, the above outlier detection results are all determined to be invalid.

[0110] Therefore, in this embodiment, outlier removal is not performed on the original data, and the fitting result obtained in step S3 is directly retained as the final result.

[0111] Step S5: Construction and application of the oscillation stiffness distribution equation.

[0112] Substituting the final fitting parameters from step S3 into the double exponential function model, we obtain the distribution equation of the oscillation stiffness as a function of the normalized axial position: y kb (x)=0.6983931488exp(-261.6882319175x)+0.3016073461exp(-4.4380667113x) Where x∈[0,1] is the normalized axial position of the blade, and y kb (x) represents the normalized linear oscillation stiffness at the corresponding position.

[0113] In engineering applications, the above distribution equation can be scaled back by combining the actual length of the target blade and the reference shimmy stiffness level, and the normalized position can be mapped to the real axis coordinates to obtain the predicted value of shimmy stiffness at any axial position. Figure 10 , Figure 11 and Figure 12The normalized fitting curve of the shimmy stiffness distribution with blade length in this embodiment is shown, which can provide a reference for the load and deformation analysis of ultra-large wind turbine blades in the shimmy direction.

[0114] Example 4 A specific embodiment of the present invention, such as Figure 1 , Figure 13 , Figure 14 , Figure 15 and Figure 16 As shown, in this embodiment... Figure 1 Based on the method flow shown, the target physical quantity is set as torsional stiffness, which is used for fitting the normalized curve of the torsional stiffness distribution with blade length. The method also includes steps S1-S5. Only the specific implementation and fitting results related to the shimmy stiffness are explained, which specifically include the following steps: Step S1: Dimensional reduction modeling and linear torsional stiffness calculation.

[0115] Based on the three-dimensional model of the multi-source wind turbine blades, the geometric dimensions and material parameters of the blade axis and the cross-sections distributed along the axis are extracted. In the finite element analysis software ANSYS, the seven wind turbine blades are equivalent to variable cross-section beam elements. The blades are divided into several short beam elements along the blade axis at a preset unit length. In this embodiment, 300 short beam elements are obtained.

[0116] The end closer to the blade root is designated as the fixed end, and the end farther from the blade root is designated as the cantilever end. A rigid planar constraint is applied to the cantilever end, and a unit rotation angle is applied around the torsional direction. The fixed-end torque of each short beam element is obtained through finite element analysis. Combined with the short beam length, polar moment of inertia of the section, and shear modulus, the corresponding linear torsional stiffness is calculated, forming a torsional stiffness sequence distributed along the axis.

[0117] Step S2: Normalization and construction of discrete point dataset.

[0118] The blade axis length was normalized to the [0,1] interval, resulting in 300 dimensionless position coordinates x. The linear torsional stiffness was normalized to ensure that the normalized stiffness y was distributed within the [0,1] interval. The "position x - normalized torsional stiffness y" data was exported as a CSV file, and the file was read to remove any invalid values ​​(NaN, Inf, etc.). The number of valid data points was confirmed to be 300, which is greater than the minimum number of data points required for fitting, thus forming a discrete point dataset for fitting.

[0119] Step S3: Fitting a double exponential function with soft boundary constraints (using all data).

[0120] In this embodiment, a double exponential function is selected as the fitting model for the normalized distribution of torsional stiffness: Where x is the normalized axial position, y is the corresponding normalized torsional stiffness value, and a, b, c, and d are the model parameters to be determined.

[0121] To reflect the relative change trend of torsional stiffness at the blade root and blade tip, a soft boundary constraint is introduced during the fitting process: at the normalized position x=0, the function value is expected to be close to 1, and at x=1, the function value is expected to be close to 0. The constraint is then superimposed on the residual function as a penalty term by setting the boundary constraint weight (40 in this embodiment).

[0122] Based on this, an initial parameter set covering various typical forms was constructed, and the Levenberg-Marquardt nonlinear least squares algorithm was used to fit all 300 data points. After solving multiple sets of initial values, a set of parameters that satisfied the boundary constraints and had good fitting quality was obtained: a = 0.31518245; b = 9.74535862;c=0.68481755;d= 3210.88718535.

[0123] At this point, the coefficient of determination R² = 0.931386, and the root mean square error RMSE = 0.04075179. At the boundary, the function values ​​satisfy f(0) = 1.000000 and f(1) = 0.000018, with deviations from the target values ​​of 1.0 and 0.0 of 0.000000 and 0.000018, respectively, which meets the boundary accuracy requirements set in this embodiment. This result serves as the preliminary fitting result of the torsional stiffness double exponential model.

[0124] Step S4: Adaptive outlier detection and fitting correction.

[0125] Based on the above preliminary results, the residuals of each data point were calculated, and candidate outliers were identified using various outlier detection methods, including: 3σ and 2.5σ thresholding methods based on standard deviation; 95th and 90th percentile thresholding methods based on residual percentiles; and thresholding methods based on interquartile range (IQR).

[0126] For each method, provided that the proportion of outliers is limited (e.g., between 1% and 20%) and the number of retained data points is not less than the minimum number required for fitting, data points with residuals exceeding the threshold are marked as outliers and temporarily deleted. The remaining data are then used to re-perform the double exponential function fitting with boundary soft constraints to obtain new parameters and their R² and RMSE indices.

[0127] In this embodiment, representative results of each method include: 3σ method: Remove 3 data points, accounting for 1.00% of the total, and after refitting, R²≈0.950106; 2.5σ method: Remove 12 data points, accounting for 4.00% of the total, and after refitting, R²≈0.967362; 95th percentile method: Remove 15 data points, accounting for 5.00% of the total, and after refitting, R²≈0.970831; 90th percentile method: Remove 30 data points, accounting for 10.00% of the total, and after refitting, R²≈0.987063; IQR method: Remove 43 data points, accounting for 14.33% of the total, and after refitting, R²≈0.993634.

[0128] Under the constraint of outlier ratio, the IQR method yielded the largest improvement in R², increasing from 0.931386 to 0.993634 compared to the result using all data, an improvement of approximately +0.062248 (+6.68%). The RMSE decreased from 0.04075179 to 0.01233394, a reduction of approximately 0.028418 (an improvement of approximately 69.73%). Therefore, the fitting result after removing 43 outliers using the IQR method is taken as the final fitting result in this embodiment.

[0129] In this case, 257 data points were retained, and the final fitting parameters were: a = 0.660543871488996; b = 3119.054418554528183;c=0.339456128511004;d= 11.430550917960563. Corresponding R²=0.99363417, RMSE=0.0123339432. At the boundary, the function values ​​are f(0)=1.000000 and f(1)=0.000004, with deviations from the target values ​​of 1.0 and 0.0 respectively of 0.000000 and 0.000004, thus continuing to meet the boundary soft constraint requirements.

[0130] Step S5: Construction and application of the torsional stiffness distribution equation.

[0131] Substituting the final fitting parameters determined in step S4 into the double exponential function model, we obtain the distribution equation of torsional stiffness as a function of normalized axial position: y kt (x)=0.6605438715exp(-3119.0544185545x)+0.3394561285exp(-11.4305509180x) Where x∈[0,1] is the normalized axial position of the blade, and y kt (x) represents the normalized linear torsional stiffness at the corresponding position.

[0132] In engineering applications, the above distribution equation can be scaled back by combining the actual length of the target blade and the reference torsional stiffness level, and the normalized position can be mapped to the real axis coordinates to obtain the predicted value of torsional stiffness at any axial position. Figure 13 , Figure 14 , Figure 15 and Figure 16 The normalized fitting curve of torsional stiffness distribution with blade length in this embodiment is shown, which can provide a reference for estimating the stiffness distribution of ultra-large wind turbine blades under torsional loads and designing test benches.

[0133] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for constructing blade stiffness and linear density based on finite element analysis, characterized in that, Includes the following steps: Step S1: Based on finite element analysis, perform dimensionality reduction modeling, calculate the linear density of the wind turbine blades, and set a rigid plane to calculate the blade stiffness; Step S2: Normalize and preprocess the calculated blade stiffness and linear density to obtain a discrete point dataset; Step S3: Fit a double exponential function with 15 different initial values, use the Levenberg-Marquardt algorithm to solve for the optimal parameters, and select R... 2 The highest fitting result is taken as the preliminary optimal solution; Step S4: Adaptive outlier detection and removal, comparing the R-squared values ​​of the standard deviation method, interquartile range method, and percentile method. 2 For the degree of improvement, select R. 2 The method that maximizes improvement is taken as the final result; Step S5: Use the double exponential function model for subsequent prediction of the linear density and stiffness of ultra-long blades.

2. The method for constructing blade stiffness and linear density based on finite element analysis according to claim 1, characterized in that, Step S1 includes the following steps: Step S1.1: Based on the multi-source wind turbine blade model data, establish a three-dimensional physical model of the wind turbine blade to be analyzed, and extract the wind turbine blade axis and the geometric dimensions and material parameters of each cross section distributed along the axis; Step S1.2: In the finite element analysis software, the wind turbine blade is equivalent to a variable cross-section beam element model arranged along the axis; the wind turbine blade axis is divided into several short beam elements according to a preset unit length. Step S1.3: Calculate the mass, linear density value and gravitational moment of each short beam unit based on the three-dimensional physical model of the wind turbine blade.

3. The method for constructing blade stiffness and linear density based on finite element analysis according to claim 2, characterized in that, Step S1 further includes the following steps: Step S1.4: Set the short beam element near the leaf root as the fixed end and the end away from the leaf root as the cantilever end. Apply rigid planar constraints to the cantilever end to form simplified cantilever beam boundary conditions. Step S1.5: Calculate the swing stiffness, sway stiffness and torsional stiffness of each short beam element.

4. The method for constructing blade stiffness and linear density based on finite element analysis according to claim 3, characterized in that, In step S1.5, unit rotation angles are applied to the cantilever end about the swinging and oscillation directions, respectively, and the corresponding fixed end bending moments are solved by finite element software. M x , M y Based on the pre-set short beam length and cross-sectional moment of inertia I x I y and elastic modulus E The swing stiffness and oscillation stiffness of each short beam element were calculated.

5. The method for constructing blade stiffness and linear density based on finite element analysis according to claim 3, characterized in that, In step S1.5, a unit rotation angle about the torsional direction is applied to the cantilever end, and the torque at the fixed end is solved by finite element software. M z Based on the short beam length and the polar moment of inertia of the cross section I p and shear modulus G The torsional stiffness of each short beam element was calculated.

6. The method for constructing blade stiffness and linear density based on finite element analysis according to any one of claims 1-5, characterized in that, Step S2 includes the following steps: Step S2.1: Normalize the axial position of each short beam unit obtained in step S1 according to the full length of the wind turbine blade, and map the axial length of the wind turbine blade to the interval [0,1] to obtain the corresponding dimensionless position coordinates; Step S2.2: Normalize the amplitude of the linear density and linear stiffness in each direction of each short beam element to obtain dimensionless data that is easy to fit numerically; export the normalized axial position, linear density and linear stiffness as a CSV file in tabular form. Step S2.3: Read the CSV file and related original parameters, and automatically remove invalid values ​​to ensure that the number of valid data points is not less than the preset threshold, so as to increase the degree of freedom for subsequent fitting. Step S2.4: Establish a correspondence between the normalized axial position and the normalized linear density and linear stiffness to form a discrete point dataset for fitting.

7. The method for constructing blade stiffness and linear density based on finite element analysis according to any one of claims 1-5, characterized in that, Step S3 includes the following steps: Step S3.1: Select the double exponential function as the fitting model for the distribution of wind turbine blade stiffness or linear density along the axial length; Step S3.2: Construct multiple sets of initial parameter sets; Step S3.3: For each set of initial parameters, the Levenberg-Marquardt nonlinear least squares algorithm is used to fit and solve the double exponential function model. Under the premise of satisfying the preset convergence conditions and the maximum number of iterations constraint, the corresponding fitting parameters are obtained.

8. The method for constructing blade stiffness and linear density based on finite element analysis according to claim 7, characterized in that, Step S3 further includes step S3.4: calculate the coefficient of determination R² and root mean square error RMSE for each group of fitting results, and select the group of parameters with the highest R² value as the preliminary optimal solution of the double exponential function model, using R² as the main evaluation index.

9. The method for constructing blade stiffness and linear density based on finite element analysis according to any one of claims 1-5, characterized in that, Step S5 includes the following steps: Step S5.1: Substitute the final fitting parameters obtained in step S4 into the double exponential function model to obtain the continuous distribution equation of the wind turbine blade linear density or stiffness as a function of normalized axial position. Step S5.2: Based on the actual length of the target wind turbine blade and the reference linear density and stiffness level, scale the distribution equation to restore or scale it, and map the normalized position to the true axis coordinates to obtain the predicted values ​​of linear density and stiffness of the target wind turbine blade at any axial position. Step S5.3: Use the prediction results for preliminary estimation of the structural strength and driving capability of the test bench, selection of design parameters for ultra-large blade structures, and related performance evaluation.

10. The method for constructing blade stiffness and linear density based on finite element analysis according to any one of claims 1-5, characterized in that, The parameters of the wind turbine blades are simplified to linear density, flapping stiffness, oscillation stiffness, and torsional stiffness.

Images