Wind turbine blade flutter characteristic calculation method based on joint simulation

Abstract

The invention discloses a wind turbine blade flutter characteristic calculation method based on joint simulation, and relates to the technical field of wind power generation. The method comprises the following steps: S1, reading structural body data of a wind turbine blade, and extracting parameters such as a spanwise position, an airfoil profile coordinate and aerodynamic characteristics of the blade; s2, full-automatic parametric modeling of the blade is achieved through a parametric modeling tool, and a related input file is generated; s3, based on an unsteady blade element momentum theory, aerodynamic force and torque responses of the blade at different spanwise positions are calculated step by step; s4, standardized calculation steps are packaged through a programmed method, personal errors are reduced, and calculation efficiency is improved; s5, designing a cabin interface, a tower interface and a floating platform interface, and realizing flexible data interaction and expansion among modules; and S6, dynamic response analysis is combined, and flutter characteristics of the wind turbine in a complex load environment are simulated. According to the method, the calculation process of wind turbine blade flutter characteristic analysis is accelerated, and the analysis precision and the calculation efficiency are improved. Through programmed and modular design, the flutter simulation of the wind turbine blade can be quickly and accurately completed, and a scientific basis is provided for structural optimization of the blade and safety design of a wind turbine system.

Description

Technical Field

[0001] The present invention relates to the technical field of wind power generation, and in particular to a method for calculating flutter characteristics of wind turbine blades based on joint simulation. Background Art

[0002] With the continuous increase in the capacity of wind turbines and their widespread application in offshore floating environments, the structure and dynamic characteristics of wind turbine systems have become increasingly complex. The increase in blade length and flexibility makes blades more prone to flutter under complex aerodynamic loads, leading to energy loss, fatigue damage, and even blade failure, which directly affects the safety and reliability of wind turbines. Therefore, accurate analysis of blade flutter characteristics has become a key link in wind turbine design.

[0003] Traditional blade flutter analysis methods are mainly based on aeroelastic theory and simplified rigid body models, which makes it difficult to fully consider the dynamic characteristics of blades under unsteady aerodynamic loads. To this end, existing technologies gradually introduce joint simulation methods to integrate the aerodynamic, structural and dynamic characteristics of wind turbine blades to improve the accuracy of flutter characteristic calculations. Through the unsteady blade element momentum theory, this method gradually calculates the aerodynamic force and torque response of the blade at different spanwise positions to achieve a more realistic flutter simulation, providing a basis for blade structure optimization and safety design.

[0004] In order to solve the above problems, the present invention is based on the wind turbine blade flutter characteristic calculation method of joint simulation, which programs the calculation process, uses MATLAB and other tools to realize automatic operation, and encapsulates various standardized calculation steps into functions or script modules to avoid manual repetitive operations, improve efficiency and reduce human errors. At the same time, the method designs the nacelle interface, tower interface and floating platform interface to facilitate model expansion and iteration. Summary of the invention

[0005] The purpose of the present invention is to propose a method for calculating the flutter characteristics of wind turbine blades based on joint simulation to solve the problems raised in the background technology. The present invention adopts the theory of unsteady blade element momentum, combined with the aerodynamic characteristic parameters of the blade along the span direction, and gradually calculates the aerodynamic force and torque response of the blade at different positions to accurately simulate the flutter behavior of the blade under complex environmental loads. By programming the calculation process and using tools such as MATLAB to realize automated operation, various standardized calculation steps are encapsulated as functions or script modules to avoid manual repetitive operations, improve efficiency, reduce human errors, and can be used for flutter parameter sensitivity detection. At the same time, the method designs the cabin interface, tower interface and floating platform interface, so that data can be flexibly exchanged between modules, which is convenient for the expansion and iteration of the model, thereby providing a scientific basis for the optimization design of blade structure.

[0006] To achieve the above object, the present invention provides the following solutions:

[0007] A method for calculating the flutter characteristics of a wind turbine blade based on joint simulation comprises the following steps:

[0008] S1. Read parameterized data from wind turbine benchmark models published by IEA, NREL, DUT, etc. or other wind turbine files that comply with the IEA WT Ontology standard, and create wind turbine structure data and NuMAD blade object data in MATLAB.

[0009] S2, using the NuMAD blade object data created in S1, and using the improved NuMAD software (denoted as NuMAD++) to realize fully automatic parametric modeling of wind turbine blades, and at the same time automatically generating a FAST input file that meets the requirements, and an APDL script file for creating a wind turbine blade shell model;

[0010] S3, using the ANSYS APDL script described in S2, calling ANSYS software through MATLAB programming, and creating a blade plate shell finite element model for subsequent simulation analysis;

[0011] S4. Using the structural data read in S1, MATLAB programming is used to automatically generate an ANSYS APDL script file for outputting the MNF modal neutral file. This script can create multiple reference airfoils on the wind turbine blade plate shell model in S3, and rigidly connect the nodes at the quarter chord length position (aerodynamic center) of each reference airfoil to its surrounding structure, perform modal analysis on the blade using the LANB method, and output a modal neutral file that conforms to the ADMAS reading format;

[0012] S5, using the ANSYS APDL script for outputting the MNF modal neutral file generated in S4, calling the ANSYS software through MATLAB programming, and generating a blade modal neutral file that meets the ADMAS reading format;

[0013] S6. Based on the structural data read in S1, an ADAMS / View command file (cmd file) for creating a blade dynamics model can be automatically generated through MATLAB programming. The cmd file can call the modal neutral file described in S5 to build a blade flexible body model in ADAMS. The model design includes setting a drive at the blade root to rotate it around the hub center and applying a concentrated force limit load at the aerodynamic center of each reference airfoil of the blade.

[0014] At the same time, the file also sets the input and output units: the displacement and velocity in the flapping direction and the swing direction are used as the data of the output unit, and the aerodynamic force calculated by the external function is used as the data of the input unit. In the blade root coordinate system, the concentrated force load of the wind turbine blade can be specifically expressed as:

[0015]

[0016] In the formula, F f represents the load in the swinging direction, F e Indicates the load in the swing direction, F z represents the gravity load and centrifugal load generated by the node and its surrounding mass, θ is the deflection angle of the flapping-shimmy load direction relative to the blade root coordinate system;

[0017] S7. Obtain the airfoil model and two-dimensional coordinates of each blade reference section from the structural data in S1, calculate the aerodynamic performance of the airfoil using XFOIL software, and then calculate the aerodynamic force at each section using the unsteady blade element momentum theory (UBEM) through MATLAB programming. The load F in the flapping direction described in S6 f and the load F in the swing direction e It can be expressed as:

[0018]

[0019] Where T is the thrust on the reference airfoil of the wind turbine, acting in the radial direction of the blade rotation plane; Q is the torque on the reference airfoil, acting in the tangential direction of the blade rotation plane; β is the offset angle between the flapping-swing array coordinate system and the blade coordinate system, and β is generally taken as 0;

[0020] S8. Using MATLAB Simulink software, the data of the aerodynamic part (aerodynamic load) in S7 and the structural part (blade displacement, velocity) in S6 can be exchanged in real time to preliminarily establish the aeroelastic coupling model of the flexible blade. In Simulink, the core settings include the MATLAB function module, which is based on the unsteady aerodynamic model, takes the wind speed time series and the blade spanwise aerodynamic characteristics (such as lift coefficient and drag coefficient) as input, and uses the Gauss-Legendre numerical integration method to gradually calculate the aerodynamic force of the blade. The calculated aerodynamic data is transmitted to ADAMS in real time through the interface to drive the structural response calculation of the blade. In ADAMS, the dynamic response of the blade (including displacement, velocity and stress distribution) is calculated according to the geometric parameters and material properties of the blade, and the response data is fed back to Simulink in real time through the interface module to update the next aerodynamic calculation. Simulink and ADAMS achieve synchronous interaction through a modular interface to ensure the high accuracy and consistency of the aerodynamic and structural coupling model;

[0021] S9. In S8, the aeroelastic coupling model of the flexible blade also reserves the nacelle interface, tower interface and floating platform interface to support model expansion and modular design.

[0022] Preferably, S1 is mainly used to automatically read wind turbine benchmark model files or industry standard model files published by IEA, NREL and DUT, etc. The specific operation steps are as follows:

[0023] S101 uses a file format that complies with the IEA WT Ontology standard. This format is widely used in the wind energy industry. It has a simple structure and is easy to read by scripts. It contains all the geometry, aerodynamics, and material information required to define the wind turbine structure.

[0024] S102, using MATLAB programming, read the file content described in S101, and create structural data of different parts according to the classification of blades, nacelles, towers and floating platforms. Secondly, convert the structural data of the blades into object data required by NuMAD, an open source software of the Sandia National Key Laboratory in the United States.

[0025] S103, through MATLAB programming, the data in S102 is automatically written into an Excel table, and the Excel table can be read by NuMAD software to create blade or blade object data. This step is optional, and is intended to facilitate reading the material, shape and ply information of the blade and quickly modify the blade parameters. At the same time, it completes the automatic reading and writing process of the object for NuMAD software, which is convenient for iterative design of blades.

[0026] Preferably, S2 mainly realizes fully automatic parametric modeling of wind turbine blades through improved NuMAD software (NuMAD++), and the specific operation steps are as follows:

[0027] S201. The improved NuMAD++ software adds parameterization of ply angles based on the original NuMAD, and can design and model the ply angles of blades.

[0028] S202, the angle parameterization method of the plies includes, first, giving the angles of different plies at all spanwise positions in the IEA WT Ontology file, and second, filling in the laminate notation (e.g. [-45 / 0 / 45]) under the fiber angle label in the Excel table described in S103. Next, NuMAD++ is used to read the above file, obtain the blade parameterized model and output the FAST input file and the ANSYS APDL script for the blade finite element model.

[0029] Preferably, the process of using the unsteady blade element momentum theory (UBEM) calculation in S7 is:

[0030] S701, data initialization:

[0031] Extract the spanwise position of the blade and the two-dimensional coordinates of the airfoil of the reference section from the structural data read in S1. Associate the airfoil model at the reference section with the corresponding aerodynamic performance data, and initialize the aerodynamic force, pitch angle, twist angle, etc. of each reference section.

[0032] S702, numbering and position initialization:

[0033] All reference sections are numbered in order from blade root to blade tip, and the starting position of the calculation is set from the blade root, processing each reference section in turn.

[0034] S703. Calculate the inflow angle and the angle of attack:

[0035] For the current reference section, according to the incoming wind speed U ∞ and blade rotation speed Ω, calculate the inflow angle Φ and angle of attack α:

[0036]

[0037] Where U ∞ is the upstream undisturbed wind speed, Ω is the blade rotation speed, α is the angle of attack, Φ is the inflow angle, θ is the sum of the pitch angle and the twist angle, a is the axial induction factor, a' is the tangential induction factor, V flap represents the disturbance velocity in the swinging direction, V edge Indicates the disturbance velocity in the swing direction;

[0038] S704, renewal inducing factor:

[0039] According to the aerodynamic characteristics of the current section (including the lift coefficient C l and the drag coefficient C d ), iteratively update the axial induction factor a and the circumferential induction factor a':

[0040]

[0041] Where F is the Prandtl tip correction factor, which is used to correct the error caused by selecting a limited number of sections to solve the aerodynamic force on the complete blade, B represents the number of blades, and c(r) represents the chord length of the airfoil at a distance r from the blade root;

[0042] S705, convergence judgment:

[0043] Determine whether the induction factors a and a' meet the convergence condition. If they have converged, proceed to the next step; otherwise, return to S703 and S704 to continue iterative updating;

[0044] S706, Calculate aerodynamic force:

[0045] Based on the convergence of the induction factor, the aerodynamic force at the current reference section is calculated.

[0046]

[0047] Where dT represents the thrust acting on the blade element, dQ represents the bending moment acting on the blade, and W is the relative inflow velocity, which is obtained from the following formula:

[0048]

[0049] S707, output result:

[0050] The aerodynamic results calculated for the current section are stored and updated to the output variable. Check whether all reference sections have been processed; if not, proceed to the calculation of the next reference section and repeat S703 to S707; the calculation formula for aerodynamics is as follows:

[0051]

[0052] In the formula, R hub Indicates the hub radius, r i and r i+1 Indicates the spanwise distance of the blade between the previous and current reference sections;

[0053] Preferably, the cabin interface, tower interface and floating platform interface in S9 specifically include the following contents:

[0054] S901, Cabin interface:

[0055] In the Simulink model, relevant entrances are reserved for the cabin interface to realize data transmission. The cabin is simplified as a rigid body with concentrated mass and rotational inertia by default, acting on the top of the tower and allowing it to rotate around the central axis of the tower. The wind rotor is connected to the cabin through a revolute pair. If the dynamic effect of the cabin needs to be introduced, the cabin structure can be imported into the multi-body dynamics model of the wind rotor through the ADAMS model extension module to form a dynamic coupling relationship between the cabin and the wind rotor. At the same time, a data parameter interface is reserved in the MATLAB function module. When the interface is not enabled, the influence of the cabin on the unsteady aerodynamic load is not considered; when the interface is enabled, the externally transmitted cabin dynamic data (such as vibration and rotation characteristics) will be used as input, affecting the calculation results of the aerodynamic load.

[0056] S902, tower interface:

[0057] An entry is reserved in the Simulink model for receiving dynamic tower data from the outside. By default, the effect of the tower on the aerodynamic load is ignored, but when necessary, a multi-body dynamic model including the tower can be constructed through the ADAMS extension module. The bottom of the tower is connected to the floating platform through a hinge or flexible connection (small movement is allowed), the top is fixedly constrained to the cabin, and the cabin is regarded as a flexible body, so as to achieve dynamic coupling between the tower and other modules. Tower-related parameters are also reserved in the MATLAB function module. When the interface is not enabled, the tower has no effect on the calculation of aerodynamic loads; after the interface is enabled, the dynamic response data of the tower (such as tower gravity, vibration mode and wind load effect) will be input through the Simulink model and used to update the calculation formula of the aerodynamic load. The specific parameter relationship can be described by the following formula:

[0058]

[0059] In the formula, C p is the power coefficient, D is the rotor diameter, ω is the rotor speed, A is the windward area, P is the output power, and ρ is the air density.

[0060] S903, floating platform interface:

[0061] The floating platform interface is used to transmit the dynamic response data of the floating platform to achieve interaction with the tower module. The dynamic response data includes:

[0062] Translational motion: displacement, velocity and acceleration of the floating platform in the X, Y and Z directions;

[0063] Rotational motion: rotation angle and angular velocity around the X, Y, and Z axes.

[0064] Furthermore, the dynamic response data of the floating platform can be processed in the following two ways:

[0065] S9031. The first method is to use the coordinate transformation matrix between the floating platform and the tower to map the motion response data of the floating platform to the coordinate system at the bottom of the tower.

[0066] S9032. The second method is to regard the movement of the floating platform as the disturbance speed to the wind field, and simulate the influence of the dynamic movement of the platform on the wind field distribution, thereby indirectly affecting the force and dynamic response of the tower and blades.

[0067] In the floating platform coordinate system, the velocity disturbance U of the floating platform on the flow field is platform , the specific formula used is:

[0068]

[0069] It can also be expressed as:

[0070]

[0071] Where U surge represents the velocity in the longitudinal direction, U sway represents the sway velocity, U heave represents the heave direction velocity, represents the angular velocity in the rolling direction, represents the angular velocity of the heading direction, represents the angular velocity in the pitch direction, Indicates the position of a point in space in the floating platform coordinate system.

[0072] Furthermore, the velocity, acceleration, force and other vector information and coordinate point position vectors between the blades (wind rotors), nacelles, towers and floating platforms mentioned above are affected by coordinate transformation to calculate aerodynamic loads. The coordinate system includes 9 coordinate systems, namely floating platform coordinate system 0, tower bottom coordinate system 1, tower top coordinate system 2, yaw coordinate system 3, hub coordinate system 4, wind rotor coordinate system 5, blade root coordinate system 6, wing surface coordinate system 7 and wind direction coordinate system 8. The origin of floating platform coordinate system 0 is located at the global center of gravity of the floating platform, the coordinate axis of tower bottom coordinate system 1 coincides with floating platform coordinate system 0, and the origin is located at the bottom of the tower, the coordinate axis of tower top coordinate system 2 coincides with the tower bottom coordinate system, and the origin is located at the top of the tower at the same position as the nacelle foundation, the origin of yaw coordinate system 3 coincides with the tower top coordinate system, but differs from the tower top coordinate system by a yaw angle θ in the xy plane. Nyaw The origin of the hub coordinate system 4 coincides with the center of the wind rotor, the y-axis coincides with the y-axis of the yaw coordinate system, the x-axis coincides with the main axis of the wind rotor and points downstream in the wind direction, and the origin of the wind rotor coordinate system 5 coincides with the hub coordinate system and differs from it by θ in the yz plane. wing , the origin of the blade root coordinate system 6 is located at the center of the blade root, the y axis coincides with the y axis of the wind rotor coordinate system, and the cone angle θ differs in the zx plane cone The origin of the wing surface coordinate system 7 is located at the intersection of the pitch axis and the wing surface, the z axis coincides with the z axis of the blade root coordinate system, and the difference in the xy plane is a pitch angle θ twist The origin of wind direction coordinate system 8 coincides with coordinate system 0, the x-axis is consistent with the wind direction, and the wind direction angle θ differs from the floating platform coordinate system wind . Its coordinate system and transformation matrix are as follows:

[0073] S9033. For vectors such as speed and force, the following formula can be used for conversion:

[0074] X B =a BA X A

[0075] In the formula, X A =(x A ,y A ,zA ), X B =(x B ,y B ,z B ), which represents the value of the vector in different coordinate systems, a BA It is expressed as a transformation matrix from coordinate system A to coordinate system B.

[0076] S9034. For the coordinate point position vector, the following formula can be used for conversion:

[0077] M B =a BA M A +r BA ×F A

[0078] Where M A , M B They represent the radius vector from the origin of coordinate system A to the origin of coordinate system B in coordinate system A respectively.

[0079] S9035, the transformation matrix between the floating platform and the tower bottom coordinate system. Here we only analyze the semi-submersible floating platform, and the wind turbine is located in the center of the floating platform:

[0080]

[0081] Where d is the distance from the origin of the floating platform coordinate system to the origin of the tower bottom coordinate system.

[0082] S9036, the transformation matrix between the tower bottom coordinate system and the tower top coordinate system is,

[0083]

[0084] In the formula, θ 1x and θ 1y They respectively represent the inclination angles of the tower due to bending or external forces around the x-axis and y-axis of the tower base coordinate system.

[0085] S9037, the transformation matrix between the tower top coordinate system and the yaw coordinate system is,

[0086]

[0087] In the formula, θ Nyaw is the yaw angle.

[0088] S9038. The transformation matrix between the yaw coordinate system and the hub coordinate system is:

[0089]

[0090] In the formula, θ tilt is the pitch angle.

[0091] S9039, the transformation matrix between the hub coordinate system and the wind rotor coordinate system is:

[0092]

[0093] In the formula, θ wing is the azimuth.

[0094] S9040, the transformation matrix between the wind rotor coordinate system and the blade root coordinate system is,

[0095]

[0096] In the formula, θ cone is the cone angle.

[0097] S9041. The transformation matrix between the blade root coordinate system and the wing surface coordinate system is:

[0098]

[0099] In the formula, θ twist is the pitch angle.

[0100] S9042. The transformation matrix between the wind direction coordinate system and the tower bottom coordinate system is:

[0101]

[0102] In the formula, θ wind is the wind direction angle.

[0103] Furthermore, the present invention also protects the application of the above method in sensitivity analysis of wind turbine flutter characteristics and parameters thereof, specifically including the following contents:

[0104] S1001. Rapidly generate or adjust multiple design variables of wind turbine blades through NuMAD++, including the geometric and structural parameters of key components such as leading edge, leading edge panel, leading edge reinforcement area, trailing edge, trailing edge panel, trailing edge reinforcement area, main beam and web.

[0105] S1002. Rapidly generate relevant scripts and files through methods from S1 to S9 for subsequent joint simulation.

[0106] S1003. Based on S1001, the programmed joint simulation method is used to calculate the flutter time domain response of the wind turbine blades, and the flutter response of the wind turbine blades under different wind speeds, airflow environments and blade dynamic behaviors is analyzed, including characteristics such as flutter frequency, amplitude and stability.

[0107] S1004. Based on the flutter time domain characteristics calculated in S1003, sensitivity analysis is performed to evaluate the impact of various design variables on the flutter characteristics of wind turbine blades. By gradually changing the leading edge, trailing edge, main beam and other parameters, the impact of different variables on the flutter amplitude, frequency and stability is analyzed.

[0108] S1005. Based on the results of the sensitivity analysis, further calculate the flutter critical speed of the wind turbine blades.

[0109] S1006. The sensitivity analysis and critical speed data obtained in S1004 and S1005 provide guidance for the optimal design of wind turbine blades. Based on the analysis results, the geometric shape and structural design of the blades are optimized, especially the adjustment of key parameters, and the aerodynamic performance and structural stability of the blades are optimized. Ultimately, it is ensured that the wind turbine blades can maintain stable flutter characteristics under various working conditions, thereby improving the operating efficiency, reliability and safety of the wind turbine.

[0110] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0111] The wind turbine blade flutter characteristic calculation method based on joint simulation provided by the present invention greatly improves the efficiency and accuracy of wind turbine blade flutter characteristic analysis by using the unsteady blade element momentum theory combined with programmed calculation and multi-module interface design. The present invention realizes flexible interaction and expansion between modules by reserving the nacelle interface, tower interface and floating platform interface, and supports iterative optimization of the model. The present invention helps to quickly and accurately simulate the flutter behavior of blades under complex environmental loads, and provides a scientific basis for the optimal design of blades and the improvement of wind turbine safety. At the same time, through the programmed calculation method, the standardized calculation steps are encapsulated as functions or script modules, which significantly reduces manual repetitive operations, improves calculation efficiency and reduces human errors. The present invention effectively solves the problems of low efficiency and poor model scalability in the prior art of wind turbine flutter characteristic analysis, and further promotes the design optimization and cost reduction of wind power generation. BRIEF DESCRIPTION OF THE DRAWINGS

[0112] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings involved in the embodiments are briefly introduced. Obviously, the drawings in the following description are only schematic illustrations of some embodiments of the present invention. For those skilled in the art, other forms of drawings can also be constructed based on these drawings without creative work.

[0113] Figure 1 It is a flow chart of a method for calculating the flutter characteristics of wind turbine blades based on joint simulation mentioned in an embodiment of the present invention;

[0114] Figure 2It is a flow chart of the joint simulation aerodynamic calculation mentioned in the embodiment of the present invention;

[0115] Figure 3 This is a wind turbine coordinate system diagram mentioned in an embodiment of the present invention. DETAILED DESCRIPTION

[0116] The following is a clear and comprehensive description of the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings in the embodiments of the present invention. It should be noted that the embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0117] The purpose of the present invention is to provide a programmed calculation method for the flutter characteristics of wind turbine blades based on joint simulation. Through programmed and modular calculation methods, the flutter characteristics of wind turbine blades under complex loads can be analyzed efficiently and accurately, while achieving model expansion and iteration.

[0118] In order to make the above-mentioned objects, features and advantages of the present invention more clear, the technical solution of the present invention is described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0119] like Figure 1 As shown, the programmed calculation method of wind turbine blade flutter characteristics based on joint simulation provided by an embodiment of the present invention includes the following steps:

[0120] S1: Using MATLAB programming, the parameterized data is read from the IEA-15-240-RWT offshore reference wind turbine model published in IEA Wind Task 37 and converted into the open source software NuMAD object of Sandia National Key Laboratory in the United States.

[0121] S2: Using the converted data from S1, NuMAD++ is used to realize the fully automatic parametric modeling of wind turbine blades, and the ply angles at the leading and trailing edges are set to [±45] n / 0 /

[45] n , the ply angle at the main beam is [±45 / -45] m , where n and m represent the number of repetitions, depending on the number of layers. Through this method, an APDL script file for creating a wind turbine blade shell model is automatically generated.

[0122] S3: Based on the APDL script file for creating the wind turbine blade shell model generated in S2, the wind turbine blade shell model is automatically generated in ANSYS using MATLAB programming, and the correctness of the model parameters such as shape, ply, and material are preliminarily verified.

[0123] S4: Automatically generate an APDL script that outputs a modal neutral file that meets the ADMAS reading format through MATLAB programming. This script can achieve: automatically create multiple reference airfoils on the wind turbine blade shell model in S3 in ANSYS, and rigidly connect the nodes at the quarter chord length position (aerodynamic center) of each reference airfoil to its surrounding structure, and finally output a modal neutral file that meets the ADMAS reading format. The parameters of the reference section are shown in the following table:

[0124] Reference cross section r(m) Airfoil Model Relative thickness (%) Chord length(m) Torsion angle (°) 0 0 Circle 1 5.2 - 1 17.55 SNL-FFA-W3-500 0.5 5.65 11.03 2 28.68 FFA-W3-360 0.36 5.70 7.18 3 38.47 FFA-W3-330blend 0.33 5.14 4.73 4 51.38 FFA-W3-301 0.301 4.482 2.56 5 62.91 FFA-W3-270blend 0.27 3.96 1.22 6 74.67 FFA-W3-241 0.241 3.501 0.19 7 90.29 FFA-W3-211 0.211 2.90 -1.58 8 102.38 FFA-W3-211 0.211 2.39 -2.16 9 114.08 FFA-W3-211 0.211 1.848 -1.57 10 117 FFA-W3-211 0.211 0.500 -1.24

[0125] S5: Automatically generate ADAMS / View command cmd file through MATLAB programming. This script can realize: automatically read the modal neutral file described in S4 in ADAMS and create a flexible blade model, and set concentrated forces at each reference section of the flexible blade. Secondly, set the input and output units:

[0126] S501: Input unit includes the aerodynamic force of the blade at each reference section, including the aerodynamic force F in the flapping direction f and the aerodynamic force F in the swing direction e .

[0127] S502: The output unit includes the speed of the blade in the flapping direction and the swing direction at each section, by obtaining the speed and angular velocity of the aerodynamic center in the X direction (chord direction) in the airfoil coordinate system at each section, and calculating:

[0128]

[0129] In the formula, v x (d) represents the velocity in the chord direction at the section with a spanwise distance d in the blade coordinate system, which is obtained by the ADAMS function VX(), ω x (d) represents the angular velocity in the chord direction of the section with a spanwise distance of d in the blade coordinate system, which is obtained by the ADAMS function WX(), and ω represents the angular velocity of the wind wheel.

[0130] S6: The parameterized data read in S1 include the aerodynamic parameters of the cross section at each spanwise position of the blade. The aerodynamic loads are calculated based on these aerodynamic parameters and using the unsteady blade element momentum theory.

[0131] Furthermore, the process of using unsteady blade element momentum theory calculation is described in S6.

[0132] S601: Data initialization:

[0133] Extract the spanwise position of the blade and the two-dimensional coordinates of the airfoil of the reference section from the structural data read in S1. Use XFOIL software to calculate the aerodynamic performance of the airfoil at the reference section, associate the airfoil model at the reference section with the corresponding aerodynamic characteristics database, and initialize the aerodynamic parameters of each reference section. Set the incoming wind speed U ∞ =10.59m / s, blade rotation speed Ω=7.56rpm.

[0134] S602: Numbering and position initialization:

[0135] All reference sections are numbered in order from blade root to blade tip, and the starting position of the calculation is set from the blade root, processing each reference section in turn.

[0136] S603: Calculate the inflow angle and angle of attack:

[0137] For the current reference section, the inflow angle Φ and the angle of attack α are calculated according to the incoming wind speed and the blade rotation speed:

[0138]

[0139] S604: Renewal inducing factor:

[0140] According to the aerodynamic characteristics of the current section (including the normal force coefficient C n and the tangential force coefficient C t ), iteratively update the axial induction factor a and the tangential induction factor a':

[0141]

[0142] S605: Convergence judgment:

[0143] It is determined whether the induction factors a and a' meet the convergence condition. If they have converged, proceed to the next step; otherwise, return to S603 and S604 to continue iterative updating.

[0144] S606: Calculate aerodynamic forces:

[0145] Based on the convergence of the induction factor, the aerodynamic force at the current reference section is calculated, and the load F in the flapping direction is obtained. f and the load F in the swing direction e .

[0146]

[0147] In the formula, R hub =3.97m.

[0148]

[0149] In the formula, β=0.

[0150] S607: Output result:

[0151] The aerodynamic results calculated for the current section are stored and updated to the output variables. Check whether all reference sections have been processed; if not, proceed to the calculation of the next reference section and repeat S603 to S607.

[0152] S7: Using MATLAB Simulink software, real-time data exchange between the aerodynamic part (aerodynamic load) in S6 and the structural part (blade displacement, speed) in S5 can be achieved, and an aeroelastic coupling model of the flexible blade is preliminarily established.

[0153] S8: Considering the offshore monopile wind turbine in the IEA Wind 15-MW reference wind turbine model, that is, the bottom is fixed in the seabed, the nacelle interface and the tower interface are added.

[0154] Furthermore, the specific operations of the cabin interface and the tower interface in S8 are as follows.

[0155] S801: For the nacelle interface, the nacelle is regarded as a rigid body, simplified as concentrated mass and moment of inertia acting on the top of the tower, and allowed to rotate around the central axis of the tower. The wind rotor is connected to the nacelle through a revolute pair.

[0156] S802: Using the multi-body dynamics method in ADAMS software, the tower is simplified into a thin-walled conical cylindrical flexible structure with a fixed bottom and a free top, and fixed constraints are imposed on the cabin and the tower top.

[0157] Furthermore, the specific operations of using the multi-body dynamics method in ADAMS software in S801 and S802 are as follows.

[0158] S8011: In S801, the nacelle is regarded as a rigid body. Through the parametric data of the wind turbine read in S1, a simplified nacelle model of the corresponding size is created in ADAMS. The model is generally set as a cuboid to facilitate modeling and simplify calculations. The front end of the nacelle is connected to the wind rotor (hub center) and a revolute pair is set.

[0159] S8021: In S802, the tower is simplified into a thin-walled conical cylindrical flexible structure with a fixed bottom and a free top. The wind turbine parametric data read in S1 is used to construct a thin-walled conical cylinder tower model using the top section dimensions and bottom section dimensions of the tower. The tower model is then imported into ANSYS to generate a modal neutral file. Finally, the flexible model of the tower is created in ADAMS using the modal neutral file. The center of the top section of the tower is connected to the bottom of the nacelle, and a revolute pair is set.

[0160] S9: Based on the above, a parameter sensitivity test can also be performed. Taking the leading edge thickness as an example, the leading edge thickness can be set to 2, 3, and 4 times the original value, respectively, to compare the flutter characteristics of the blades and obtain the flutter critical speed.

[0161] The remaining technical features in this embodiment can be flexibly selected by those skilled in the art according to actual conditions to meet different specific practical needs. However, it is obvious to those of ordinary skill in the art that it is not necessary to adopt these specific details to implement the present invention. In other examples, in order to avoid confusing the present invention, the composition, structure or components of the formula are not specifically described, and they are all within the technical protection scope defined by the technical solution claimed for protection in the claims of the present invention.

[0162] Modifications and changes made by those skilled in the art that do not depart from the spirit and scope of the present invention should be within the scope of protection of the claims attached to the present invention. In the above description, a large number of specific details are set forth in order to provide a thorough understanding of the present invention. However, it is obvious to those of ordinary skill in the art that these specific details are not necessary to practice the present invention. In other examples, in order to avoid confusing the present invention, well-known technologies, such as specific construction details, operating conditions and other technical conditions, are not specifically described.

[0163] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for calculating the flutter characteristics of wind turbine blades based on joint simulation, characterized in that: The specific steps include: S1. Use MATLAB programming to automatically read and preprocess data, automatically read parameterized data from wind turbine benchmark models or other wind turbine files that meet the IEA WT Ontology standard, and create wind turbine structure data in MATLAB; S2. Use MATLAB programming to convert the structural data in S1 into object data, and use the improved open source software NuMAD to achieve fully automatic parametric modeling of wind turbine blades, automatically generate FAST input files that meet the requirements, and APDL script files for creating wind turbine blade shell models, providing standardized input data support for simulation analysis and structural optimization; S3, based on the structure data read in S1, automatically generate APDL script files through MATLAB programming; Use the generated script file to create a reference airfoil on the wind turbine blade shell model obtained in S2, and rigidly connect the nodes of the reference airfoil at a position one quarter of the chord length from the leading edge to its surrounding structure, and output a modal neutral file that meets the ADMAS reading format; S4. Using MATLAB programming to automatically generate an ADAMS / View command cmd file; using the command cmd file to read the blade modal neutral file in S3 in ADAMS and create a flexible blade model, and load a concentrated force limit load at a reference section of the flexible blade; In the blade root coordinate system, the calculation formula of the concentrated force load of the wind turbine blade is expressed as: In the formula, F f represents the load in the swinging direction, F e Indicates the load in the swing direction, F z represents the gravity load and centrifugal load generated by the node and its surrounding mass, θ is the deflection angle of the flapping-shimmy load direction relative to the blade root coordinate system; S5. Obtain the airfoil model and two-dimensional coordinates of each blade reference section from the structural data in S1, calculate the aerodynamic performance of the airfoil using XFOIL software, and then calculate the aerodynamic force at each reference section using unsteady blade element momentum theory through MATLAB programming, and obtain the load F in the flapping direction. f and the load F in the swing direction e , the specific calculation formula is: Where T is the thrust on the reference airfoil of the wind turbine, acting in the radial direction of the blade rotation plane; Q is the torque on the reference airfoil, acting in the tangential direction of the blade rotation plane; β is the offset angle between the flapping-swing array coordinate system and the blade coordinate system; S6. Using the Simulink function module, based on the aerodynamic calculation method described in S5, the aerodynamic force acting on the flexible blade is calculated according to the wind speed time series and the aerodynamic characteristics along the blade span direction; the calculated aerodynamic data is transmitted to the ADAMS software in real time through the ADAMS / Simulink coupling interface, and in ADAMS, the dynamic response of the blade is calculated by combining the flexible blade model, blade geometric parameters and material properties; the dynamic response data is fed back to Simulink in real time through the modular interface for updating the calculation of the aerodynamic force; Synchronous data interaction between Simulink and ADAMS to achieve high-precision real-time coupling between the pneumatic and structural parts; S7, construct a wind turbine model, reserve a nacelle interface, a tower interface and a floating platform interface on the wind turbine model, and further construct a wind turbine model coordinate system to realize the conversion of various physical quantities such as load, displacement, velocity and acceleration between various components and interfaces; the wind turbine model coordinate system includes 9 coordinate systems, namely: Floating platform coordinate system, ① tower bottom coordinate system, ② tower top coordinate system, ③ yaw coordinate system, ④ hub coordinate system, ⑤ rotor coordinate system, ⑥ blade root coordinate system, ⑦ wing surface coordinate system and ⑧ wind direction coordinate system.

2. The wind turbine blade flutter characteristic calculation method based on joint simulation according to claim 1, characterized in that: The specific improvements of the improved open source software NuMAD described in S2 compared with the conventional NuMAD are as follows: S201, writing the object data in S2 into an Excel table: the generated table content and format meet the requirements of NuMAD software file reading, so that users can view and modify parameter information, and with the help of the function of reading Excel table in NuMAD software, read and write Excel, so as to be used for iterative design and optimization of blades; S202. Adding ply angle parameterization based on NuMAD: reading ply angle data from an industry standard model or a wind turbine reference model, wherein the angle data includes ply angle configurations of different regions of the leading edge, trailing edge, and main beam; S203, the user manually defines the ply angles of different areas according to specific design requirements; S204, accurately allocating the ply angles to the corresponding blade regions according to the reference data or user-defined data read; specifically, automatically matching the angle configurations of the leading edge, trailing edge, and main beam regions to the segmented model of the blade to ensure the accuracy and consistency of ply distribution; S205, verifying the input ply angle to ensure that the data meets the material mechanics constraints and design requirements; S206. Expand the input and output functions to generate the aerodynamic, structural and dynamic input files of the blades with one click for subsequent joint simulation analysis.

3. The method for calculating wind turbine blade flutter characteristics based on joint simulation according to claim 1, characterized in that: The calculation of the aerodynamic force at each reference section using the unsteady blade element momentum theory described in S5 specifically includes the following contents: S501, using the object data in S2, obtaining the airfoil model at each blade reference section and the airfoil coordinates in the airfoil coordinate system; S502, numbering the reference airfoils in S501 in order from blade root to blade tip; S503, initializing the position of the reference section and starting calculation from the blade root; S504. Determine whether the number of the current reference section exceeds the number of the reference section at the blade tip. If so, it indicates that all reference sections have been calculated, and then output the thrust T and torque Q on the blade at each reference section position; otherwise, continue to use the unsteady blade element momentum theory to calculate the aerodynamic force of the current reference section.

4. The method for calculating wind turbine blade flutter characteristics based on joint simulation according to claim 3, characterized in that: The calculation of the aerodynamic force of the current reference section using the unsteady blade element momentum theory described in S504 specifically includes the following contents: S5041, initializing the axial induction factor a and the circumferential induction factor a'; S5042, update the inflow angle Φ and the angle of attack α, the calculation formulas are as follows: Where U ∞ is the upstream undisturbed wind speed, Ω is the blade rotation speed, α is the angle of attack, Φ is the inflow angle, θ is the sum of the pitch angle and the twist angle, a is the axial induction factor, and a' is the tangential induction factor; S5043. Update the axial induction factor a and the circumferential induction factor a'. The calculation formula is as follows: In the formula, C l is the lift coefficient, C d is the drag coefficient, F is the Prandtl tip correction factor, which is used to correct the error caused by selecting a finite number of sections to solve the aerodynamic force on the complete blade, B is the number of blades, and c(r) is the chord length of the airfoil at a distance r from the blade root; S5044. Determine whether a and a' in S5043 converge. If so, re-execute S504. Otherwise, repeat S5042 and S5043.

5. The method for calculating the flutter characteristics of wind turbine blades based on joint simulation according to claim 4, characterized in that: The calculation of the inflow angle Φ described in S5042 takes into account the influence of the blade bending velocity and introduces the disturbance velocity V in the blade flapping direction on the axial induced velocity. flap , the disturbance velocity V in the blade swing direction is introduced into the tangential induced velocity term edge .

6. The method for calculating the flutter characteristics of wind turbine blades based on joint simulation according to claim 1, characterized in that: The S6 specifically includes the following contents: S601, using a customized MATLAB function module, based on the method of using the unsteady blade element momentum theory to calculate the aerodynamic force at each reference cross section described in S5, and using the Gauss-Legendre numerical method to calculate the aerodynamic force of the flexible blade according to the wind speed time series and the blade spanwise aerodynamic characteristics; S602, transmitting the aerodynamic force calculated in S601 to ADAMS in real time through a modular interface, and calculating the dynamic response of the blade in ADAMS by combining the flexible blade model, blade geometric parameters and material properties, wherein the dynamic response includes displacement, velocity and stress distribution; S603. Feedback the blade dynamic response to Simulink in real time through the modular interface to update the aerodynamic calculation, and realize high-precision real-time coupling between the aerodynamic part and the structural part through synchronous data interaction between Simulink and ADAMS. The ode45 algorithm is used to solve the aeroelastic coupling to ensure the calculation accuracy of unsteady aerodynamic force and structural dynamic response.

7. The method for calculating wind turbine blade flutter characteristics based on joint simulation according to claim 1, characterized in that: S7 constructs a wind turbine model coordinate system to achieve the conversion of various physical quantities such as load, displacement, velocity and acceleration between various components and interfaces, specifically including the following contents: S701. For vectors such as velocity and force, use the following formula for conversion: X B =a BA X A Where, X A =(x A ,y A ,z A ), X B =(x B ,y B ,z B ), which represents the value of the vector in different coordinate systems, a BA It is expressed as the transformation matrix from coordinate system A to coordinate system B; S702: For the coordinate point position vector, use the following formula to convert: M B =a BA M A +r BA ×F A Where M A , M B They represent the radius vector from the origin of coordinate system A to the origin of coordinate system B in coordinate system A respectively; S703, construct the transformation matrix between the floating platform and the tower bottom coordinate system, only analyze the semi-submersible floating platform, and the wind turbine is located at the center of the floating platform. The specific function is expressed as: Where d is the distance from the origin of the floating platform coordinate system to the origin of the tower bottom coordinate system; S704, constructing a transformation matrix between the tower bottom coordinate system and the tower top coordinate system, the specific function is expressed as: In the formula, θ 1x and θ 1y They respectively represent the inclination angles of the tower frame around the x-axis and y-axis of the tower base coordinate system due to bending or external forces; S705, constructing a transformation matrix between the tower top coordinate system and the yaw coordinate system. The specific function is expressed as: In the formula, θ Nyaw is the yaw angle; S706, constructing a transformation matrix between the yaw coordinate system and the hub coordinate system. The specific function is expressed as: In the formula, θ tilt is the pitch angle; S707, construct the transformation matrix between the hub coordinate system and the wind rotor coordinate system. The specific function is expressed as: In the formula, θ wing is the azimuth; S708, constructing a transformation matrix between the wind rotor coordinate system and the blade root coordinate system. The specific function is expressed as: In the formula, θ cone is the cone angle; S709, construct a transformation matrix between the blade root coordinate system and the wing surface coordinate system. The specific function is expressed as: In the formula, θ twist is the pitch angle; S710, construct a transformation matrix between the wind direction coordinate system and the tower bottom coordinate system, and the specific function is expressed as: In the formula, θ wind is the wind direction angle.

8. The method for calculating wind turbine blade flutter characteristics based on joint simulation according to claim 1, characterized in that: In S7, the nacelle interface, tower interface and floating platform interface are reserved on the wind turbine model for model expansion and modular design, including the following: S7011. In the Simulink model, an entry is reserved for the cabin interface to realize data transmission; the cabin is simplified as a rigid body with concentrated mass and rotational inertia by default, acting on the top of the tower and rotating around the central axis of the tower; the wind rotor and the cabin are connected through a revolute pair; if the dynamic effect of the cabin is introduced, the cabin structure is imported into the multi-body dynamics model of the wind rotor through the ADAMS model extension module to form a dynamic coupling relationship between the cabin and the wind rotor; a data parameter interface is reserved in the MATLAB function module. When the interface is not enabled, the influence of the cabin on the unsteady aerodynamic load is not considered; when the interface is enabled, the cabin dynamic data transmitted externally is used as input to consider the influence of the aerodynamic load on the calculation results; S7012. In the Simulink model, an entry is reserved for the tower interface to receive dynamic tower data from the outside. By default, the influence of the tower on the aerodynamic load is ignored. For situations where the influence cannot be ignored, a multi-body dynamic model including the tower is constructed through the ADAMS extension module. The bottom of the tower is connected to the floating platform through a hinge or flexible connection, and the top is fixedly constrained to the cabin to achieve dynamic coupling between the tower and other modules. Tower-related parameters are reserved in the MATLAB function module. When the interface is not enabled, the tower has no effect on the calculation of aerodynamic loads. After the interface is enabled, the dynamic response data of the tower is input through the Simulink model and used to update the calculation formula of the aerodynamic load. S7013. A floating platform interface is reserved in the Simulink model for transmitting the dynamic response data of the floating platform. By default, the influence of the floating platform on the aerodynamic load is ignored. When the interaction of the floating platform with the aerodynamic load needs to be introduced, its motion response data is received through the Simulink interface, including dynamic information of six degrees of freedom: translation in three directions and rotation around three axes, specifically referring to the displacement, velocity and acceleration in the X, Y and Z directions in the floating platform coordinate system, as well as the rotation angle and angular velocity around the X, Y and Z axes. The above data are transmitted to Simulink through external simulation software for updating the calculation of the aerodynamic load. A parameter entry for the floating platform is reserved in the MATLAB function module. When it is not enabled, the influence of the floating platform on the aerodynamic load is ignored by default. After enabling, the externally input floating platform data dynamically participates in the correction calculation of the aerodynamic load, forming a two-way coupling between the aerodynamic load and the dynamic response of the floating platform.

9. Application of the method according to claims 1 to 8 to the analysis of wind turbine flutter characteristics and parameter sensitivity thereof, characterized in that: Specifically, it includes the following contents: by adjusting or generating test variables, calculating the flutter time-domain characteristics of the wind turbine, determining the sensitivity of key parameters, and further obtaining the flutter critical speed, providing guidance for the optimal design of wind turbine blades.

Images