A flutter state judgment method for a two-dimensional airfoil blade structure

Abstract

The present invention provides a method for judging the flutter state of a binary airfoil blade structure, which belongs to the technical field of wind power generation based on computer data processing. First, a binary airfoil flutter model of the blade is established, and the flutter equation is derived through the Lagrange equation. Then, the frequency domain unsteady aerodynamic force is converted into the time domain using the Jones approximation, and the Theodorsen function is regarded as a filter transfer function to construct a state space model to achieve the time domain expression of the aerodynamic force. Aerodynamic state variables are introduced, and the aerodynamic force is decomposed into non-circular and circular parts, and a time domain aeroelastic equation containing displacement, velocity and aerodynamic state variables is established. By calculating the time response of generalized coordinates under different wind conditions, the aeroelastic stability of the system is judged by the attenuation, equal amplitude oscillation or divergence of the response curve, and the critical state of flutter is determined. The present invention has high accuracy, wide application range, and strong intuitiveness, and can be quickly applied to the flutter analysis of blades of different models.

Description

Technical Field

[0001] The present invention belongs to the technical field of wind power generation based on computer data processing, and in particular relates to a method for judging the flutter state of a binary airfoil blade structure. Background Art

[0002] Offshore wind turbine installed capacity is increasing annually. To better capture wind energy, offshore wind turbine blades are increasingly longer, increasing blade flexibility and aeroelastic coupling, leading to an increased risk of blade flutter. Flutter is a self-excited vibration phenomenon caused by aerodynamic forces, which can lead to blade damage or even fracture in severe cases. Existing flutter analysis techniques based on quasi-steady-state aerodynamic models assume that the airflow is in a steady state at every moment, ignoring dynamic changes and hysteresis effects. This makes them unable to cope with unsteady disturbances such as sudden offshore wind speed changes and dynamic stall, resulting in significant errors in the calculation of wind turbine blade flutter velocity. Traditional CFD / CSD coupling methods, however, can take over 24 hours per simulation, making them inadequate for online wind farm monitoring. In summary, existing technologies struggle to meet the demand for high-precision, real-time flutter warnings for offshore wind turbine blades under unsteady wind conditions. A time-domain, real-time flutter assessment method that integrates unsteady aerodynamic forces is needed to address the challenge of warning offshore blades against aeroelastic instability over a wide wind speed range.

[0003] The invention, entitled "A Method and Apparatus for Calculating the Critical Wind Speed ​​of Wind Turbine Blade Airfoil Flutter," with application number CN108491644B, discloses establishing a blade aeroelastic equation based on the energy balance principle, thereby developing a power flow balance equation. Finally, an analytical formula for the critical flutter speed (including amplitude, damping, and phase parameters) is derived, which is solved by combining a time-domain averaging method with numerical iteration. However, this method relies on a steady aerodynamic model, ignores unsteady aerodynamic effects, and does not involve adaptability verification for complex operating conditions. The invention, entitled "A Method for Predicting Wind Turbine Blade Airfoil Flutter," with application number CN103810341B, discloses using an eigenvalue method to coarsely screen the flutter interval, combining it with a time-domain analysis method to refine it in the critical interval. Finally, the Runge-Kutta method is used to solve the response curve and, combined with tip speed ratio adjustment, to determine the severity of the flutter. However, this method still uses the assumption of a steady aerodynamic model, resulting in low accuracy in calculating flutter wind speed in actual wind fields. The invention, entitled "A Wind Turbine Blade Flutter Protection Method, System, Medium, and Equipment," with application number CN118686735A, discloses determining flutter risk by monitoring the current and torque characteristic values ​​of the variable pitch motor and combining wind speed and blade angle thresholds. However, this method relies on motor signals and cannot determine the critical point, thus representing post-facto protection. None of the above methods incorporates a high-precision unsteady aerodynamic model, nor does it cover complex operating conditions such as turbulence and dynamic stall. Furthermore, they ignore aerodynamic lag effects, resulting in insufficient accuracy in critical state calculations. Summary of the Invention

[0004] To address the above problems, the present invention adopts the time domain analysis method, Theodorsen theory combined with the state space method, and directly solves the aeroelastic coupling equation through unsteady aerodynamic modeling and Laplace domain transformation. It can accurately simulate complex working conditions such as dynamic stall, and directly determine the stability through the response curve (divergence / convergence). It is more comprehensive in adaptability and has higher calculation accuracy.

[0005] The present invention provides a method for determining the flutter state of a binary airfoil blade structure, comprising the following steps:

[0006] S1, obtaining the blade physical parameters in the corresponding scenario, establishing a blade flutter physical model based on the selected binary airfoil model, calculating the system kinetic energy and system potential energy based on the blade flutter physical model, substituting the system kinetic energy and system potential energy into the Lagrange equation to obtain the flutter equation including the mass matrix, stiffness matrix, damping matrix and aerodynamic force vector; wherein the aerodynamic force vector is calculated using unsteady aerodynamic theory, including unsteady aerodynamic forces and unsteady aerodynamic moments using the Theodorsen function;

[0007] S2, using the Jones approximation method (rational function approximation method) to convert the aerodynamic force vector of the S1 flutter equation into a transfer function of a second-order system, and converting the transfer function into a state equation form through Laplace transform to construct a state space model of the second-order system; the rational function approximation algorithm (i.e., the Jones approximation method) used in the present invention is a method based on aeroelasticity theory that performs model order reduction and rapid calculation of complex aerodynamic problems under specific assumptions. Specifically, it ignores high-order nonlinear terms and approximates the coupling relationship between aerodynamic force and structural deformation as a linear relationship, thereby simplifying the unsteady aerodynamic problem into a superposition of a series of quasi-steady states, and is suitable for the field of aeroelastic stability analysis;

[0008] S3, decomposes the transformed aerodynamic force vector in the state-space model into a non-circulating part and a circulating part, introduces aerodynamic state variables including displacement, velocity and filter state, and constructs the time-domain state-space equation based on the flutter equation of S1;

[0009] S4 inputs the real-time wind speed value into the time-domain state-space equation, solves the time-domain state-space equation using the fourth-order Runge-Kutta method, outputs the response curves of the heave displacement and pitch angle, and determines the flutter state based on the oscillation shape of the curve.

[0010] Preferably, the blade physical parameters include cross-sectional mass m , swing stiffness EI , torsional stiffness GJ , Section distance to blade root radius r , center of mass offset coefficient x a , chord length of airfoil sectionc and half chord length b=c / 2 ;

[0011] The elastic axis is established at the shear center of the blade, and the center of mass is located behind the shear center. x a b where x a is the mass center offset coefficient, the aerodynamic force and aerodynamic moment act on the aerodynamic center, and the aerodynamic center is at a distance from the leading edge of the blade. c / 4 The translation displacement of the blade reflecting the up and down vibration h , the downward direction is positive, reflecting the angular displacement of the blade torsional vibration is , the head rising into the wind is positive; the required parameters are obtained through the blade structure design data.

[0012] Preferably, the system kinetic energy in S1 is the sum of translational kinetic energy and rotational kinetic energy, and the translational kinetic energy is given by the cross-sectional mass m , translational speed and the static mass moment of the airfoil per unit length about the elastic axis S EA Determine; the rotational kinetic energy is determined by the moment of inertia around the elastic axis J EA and angular velocity Determine; the system potential energy is the sum of the swing potential energy and the torsional potential energy. The swing potential energy is determined by the swing stiffness k h and translational displacement h Determine; torsional potential energy is determined by torsional stiffness and angular displacement Sure;

[0013] Inputting the system kinetic energy and system potential energy expressions into the Lagrange equation, the flutter equation matrix form is derived:

[0014] ;

[0015] in, is the acceleration in the sinking and floating direction, is the pitch angular acceleration, g h , is the damping coefficient, - L is the unsteady aerodynamic force, is the unsteady aerodynamic moment, - L and Together they form the aerodynamic force vector, which includes the air density , incoming flow velocity V and Theodorsen function C ( k ).

[0016] Preferably, the specific process of S2 includes:

[0017] The Jones approximation method is used to convert the aerodynamic force vector of the S1 flutter equation into a transfer function of a second-order system. The Jones approximation method uses a simple rational function to approximate the complex aerodynamic force function. The numerator describes the response of the aerodynamic force to the motion speed; the denominator describes the aerodynamic hysteresis effect. For a two-dimensional airfoil, the Jones approximation of the Theodorsen function is:

[0018] ;

[0019] Where, P 1 is the coefficient of the numerator quadratic term, and its Jones approximate universal constant is 0.006825. P 2 is the coefficient of the first-order term in the numerator, and its Jones approximate universal constant is 0.1080075. Z 1 is the coefficient of the quadratic term in the denominator, and its Jones approximate universal constant is 0.01365. Z 2 is the coefficient of the first-order term in the denominator, and its Jones approximate universal constant is 0.3455. V is the wind speed, b is the half chord length, s is the Laplace variable;

[0020] After the Theodorsen function is expressed as a transfer function, the unsteady term of the aerodynamic force vector is described by the transfer function and a linear combination of displacement and velocity. The input of the transfer function is defined as the motion combination of the blade, that is, the linear combination of displacement and velocity, and the output is the unsteady term of the aerodynamic force component.

[0021] According to the dynamic characteristics of the transfer function, two internal state variables are automatically generated and recorded as x f1 and x f2 , which is used to characterize the dynamic evolution process of the aerodynamic hysteresis effect. Through these two state variables, the state space model of the second-order filter is constructed. The output of the state space model is the time domain expression of the unsteady term of the aerodynamic force vector.

[0022] Preferably, the specific process of S3 includes:

[0023] Decompose the transformed aerodynamic force vector in the state space model into non-circular components F c and circulation part F nc , the non-circulating part is the steady term in the aerodynamic component, which is composed of the linear combination of acceleration term and velocity term; the circulating part is the unsteady term in the aerodynamic component, which is composed of the linear combination of filter state variables and displacement and velocity term;

[0024] The physical quantities used to constitute the non-circular part and the circular part are combined into a state vector, which contains the displacement ( h , ),speed( , ) and filter state variables ( x f1 , x f2 ), based on the flutter equation of S1 and the state space model of S2, the time domain state space equation of the entire system is constructed:

[0025] ;

[0026] in, A ( V ) is a 6-row 6-column matrix:

[0027] ;

[0028] matrix M is the sum of the mass matrix and the mass matrix derived from the non-circular aerodynamics, the matrix D is the product of the damping matrix and the damping matrix derived from the non-circular aerodynamics and the wind speed. K is the sum of the stiffness matrix and the product of the stiffness matrix derived from the non-circular aerodynamic force and the square of the wind speed, K a is the additional stiffness matrix caused by the circulation aerodynamic force, D a is the additional damping matrix caused by the circulation aerodynamic force, A a is the filter state matrix caused by the circulation aerodynamic force.

[0029] Preferably, the specific process of S4 is:

[0030] Enter real-time wind speed V In the time-domain state-space equation, the fourth-order Runge-Kutta method is used to iteratively solve the time series of the state vector, and the time response curves of the heave displacement and pitch angle are extracted from the state vector solution:

[0031] ;

[0032] ;

[0033] ;

[0034] ;

[0035] ;

[0036] in, is the time step, is the time point of the nth iteration, is the state vector of the nth iteration, k 1 to k 4 is the four incremental steps of the fourth-order Runge-Kutta method calculation process. Each step depends on the previous state and the right-hand side function of the state equation. The state vector is updated by weighted average incremental steps, and the time series data of the heave displacement and pitch angle are extracted from the state vector.

[0037] Preferably, the specific method of determining the flutter state in S4 is:

[0038] If the response curve shows attenuated oscillation, the system is stable; if the response curve shows constant amplitude oscillation, it is judged to be in a critical flutter state; if the response curve shows exponential divergence, the system is unstable.

[0039] Preferably, the wind speed values ​​are scanned incrementally at a certain wind speed interval within a preset wind speed range; a response curve is solved for each wind speed value, and the wind speed corresponding to the first occurrence of equal-amplitude oscillation is identified as the flutter critical wind speed.

[0040] Compared with the prior art, the present invention has the following beneficial effects:

[0041] Existing flutter assessment methods mostly rely on quasi-steady theory or complex hardware devices (such as adding flaps and piezoelectric drive units). The method of the present invention does not require modification of the blade structure, reducing manufacturing and maintenance costs.

[0042] Compared with methods based on eigenvalue method or threshold judgment that require complex calculations, this method directly judges stability through the shape of the time domain response curve. The results are intuitive and computationally efficient, avoiding tedious parameter calibration and threshold setting.

[0043] Some existing methods rely on computational fluid dynamics simulation and have high requirements on computing resources. This method builds a model based on classical theory, has strong versatility, and can be quickly applied to flutter analysis of blades of different models, providing a convenient technical means for engineering practice.

[0044] Specific advantages include:

[0045] (1) Based on the classic Theodorsen theory, the aeroelastic properties of the entire blade can be predicted by taking the two-dimensional airfoil as the research object, without relying on specific hardware modification, which is convenient for the promotion and application of different types of offshore wind turbine blades;

[0046] (2) No complex eigenvalue calculations or energy equations are required; stability can be directly judged by the shape of the generalized coordinate time response curve, and the results are intuitive and easy to understand;

[0047] (3) The complex aeroelastic coupling problem is converted into a standard time-domain state equation for solution, which is suitable for real-time analysis under complex and variable operating conditions of offshore wind power;

[0048] (4) The calculation is simple and fast, with higher efficiency than the traditional finite element method and higher calculation accuracy than the traditional quasi-steady theory. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] In order to more clearly illustrate the technical solutions of the present invention or the prior art, a brief introduction will be given below to the drawings required for use in the embodiments or the description of the prior art. Obviously, what is described below is only one embodiment of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0050] Figure 1 It is a flow chart of the overall process of the present invention.

[0051] Figure 2 It is a schematic diagram of constructing a two-dimensional airfoil flutter model of a blade according to the present invention.

[0052] Figure 3 1 is a response curve diagram of the pitching motion of the system under different wind speeds in an embodiment of the present invention.

[0053] Figure 4 2 is a response curve diagram of the system's sinking and floating motion under different wind speeds in an embodiment of the present invention. DETAILED DESCRIPTION

[0054] This embodiment takes offshore wind turbine blades as an example to further illustrate the specific implementation of the present invention.

[0055] The development of offshore wind turbine blades towards ultra-long lengths has led to an increased risk of flutter. The existing flutter assessment method based on quasi-steady aerodynamic models ignores the unsteady aerodynamic lag effect. Under complex working conditions such as turbulence and dynamic stall, the calculation error of the critical flutter wind speed is large, and the flutter judgment accuracy is low, which makes it difficult to meet the requirements for accurate assessment of the critical flutter state in actual wind fields.

[0056] This study uses an NREL 5 MW blade as the research object. Based on Theodorsen's unsteady aerodynamic theory, the Jones approximation (rational function approximation) is used to transform the frequency-domain unsteady aerodynamic forces into the time domain. This approach, combined with the state-space method, efficiently solves the aeroelastic response in the time domain. Aerodynamic state variables are introduced to construct the time-domain aeroelastic equations. By calculating the time response of the blade's generalized coordinates under different wind conditions, the aeroelastic stability of the system is determined based on the oscillation form of the response curve (attenuation, constant amplitude oscillation, and divergence), thereby enabling analysis and evaluation of blade flutter.

[0057] The overall process is as follows Figure 1 As shown:

[0058] A two-dimensional airfoil flutter model of the blade is established, considering the positional relationship between the elastic axis, center of mass, and aerodynamic center as well as the flapping and torsional dual-degree-of-freedom vibrations, and the flutter equation is derived through the Lagrange equation.

[0059] The Jones approximation method is used to convert the frequency domain unsteady aerodynamic force vector into the time domain. The Theodorsen function is regarded as the filter transfer function, and a state space model is constructed to realize the time domain expression of the aerodynamic force.

[0060] Aerodynamic state variables are introduced, and the transformed aerodynamic force vector in the state space model is decomposed into a non-circular part and a circulating part. The time-domain aeroelastic equation including displacement, velocity and aerodynamic state variables is established.

[0061] By calculating the time response of generalized coordinates under different wind speeds, the aeroelastic stability of the system is judged by the attenuation, constant amplitude oscillation or divergence of the response curve, and the critical flutter state is determined.

[0062] To achieve the above objectives, the embodiments of the present invention provide the following solutions:

[0063] The flutter model is established by taking the typical two-dimensional airfoil section of the blade as the research object. The chord length of the airfoil section is c , half chord length of blade b=c / 2 , the elastic axis is established at the shear center C S The torsional stiffness of the elastic shaft is , the fixed stiffness in the displacement direction is k h , leaf centroid C G Behind the shear center x a b Department, C A The point is the aerodynamic center of the blade, which is located at one quarter of the chord length of the leading edge of the blade. b / 2 At, aerodynamic L and aerodynamic torque Acting on the aerodynamic center, such as Figure 2 The blade has two degrees of freedom of vibration in two directions, and the translational displacement of the blade that reflects the up and down vibration is h , the angular displacement of the blade torsional vibration is .

[0064] The structural parameters selected in this implementation are shown in Table 1:

[0065] Table 1 Blade structure parameters

[0066]

[0067] The vibration amplitude of the blade is small near the equilibrium position, and it can be approximately considered , the airfoil flutter equation is derived through the Lagrange equation. The coordinate origin is selected at the center of mass, then the displacement of any point under the airfoil is: .

[0068] The kinetic energy of the system can be expressed as:

[0069] ;

[0070] Where, m is the mass of the airfoil per unit length, S EA Static mass moment of unit length airfoil about the rigid center; J EA is the moment of inertia about the elastic axis.

[0071] The potential energy of the system is: ;

[0072] Where, k h is the equivalent flapping stiffness of the cantilever blade section, ; is the equivalent torsional stiffness of the section, ; EI Provide blade swing stiffness; GJ is the blade torsional stiffness; r is the radius from the cross section to the blade root.

[0073] Substituting the system kinetic energy and potential energy into the Lagrange equations, we obtain the flutter equations containing the mass matrix, stiffness matrix, damping matrix, and aerodynamic force vector. The aerodynamic force vector is calculated using unsteady aerodynamic theory and includes unsteady aerodynamic forces and unsteady aerodynamic moments using the Theodorsen function.

[0074] Substitute the expressions for kinetic energy and potential energy into the Lagrange equations:

[0075] ;

[0076] The flutter equation of the blade is:

[0077] ;

[0078] Where, is the acceleration in the sinking and floating direction, is the pitch angular acceleration, g h , is the damping coefficient.

[0079] Traditional aeroelastic analysis is typically performed in the frequency domain. With the advancement of control theory, the need to address aeroelastic problems in the time domain is increasing, requiring the establishment of a time-domain aeroelastic state equation for the system. Accurate nonlinear aeroelastic analysis relies on a time-domain unsteady aerodynamic model that describes the arbitrary motion of the system. However, for aeroelastic control applications, model reduction methods based on frequency-domain data, such as the Roger approximation and the minimum state approximation, can be employed. These methods utilize unsteady aerodynamic data corresponding to a finite number of discrete frequency points, extend them to the Laplace domain through analytical extension, and express them as rational functions. Subsequently, the unknown coefficients in the rational functions are determined through function fitting techniques. Auxiliary aerodynamic hysteresis state variables are introduced to obtain the time-domain differential equation, ultimately expressing the dynamics of the entire aeroelastic system in state-space form.

[0080] According to the flutter equation derived above, the Theodorsen function is calculated using the Jones approximation method C ( k ), converted into a transfer function of a second-order system C s , the expression is:

[0081] ;

[0082] Where, P 1 is the coefficient of the numerator quadratic term, and its Jones approximate universal constant is 0.006825. P 2 is the coefficient of the first-order term in the numerator, and its Jones approximate universal constant is 0.1080075. Z 1 is the coefficient of the quadratic term in the denominator, and its Jones approximate universal constant is 0.01365. Z 2 is the coefficient of the first-order term in the denominator, and its Jones approximate universal constant is 0.3455. s is the Laplace variable, V is the wind speed, b Theodorsen function can be regarded as a second-order transfer function of a filter. The coefficients in front of Laplace are replaced by a 0, a 1, b 0, b 1 instead, then a 0= P 1 V 2 / b 2 , a 1= P 2 V / b , b 0= Z 1V 2 / b 2 , b 1= Z 2 V / b , the unsteady part of the aerodynamic force vector containing the Theodorsen function is expressed by a transfer function, whose input is: , the output is: , where and They are and Laplace transform.

[0083] According to the dynamic characteristics of the transfer function, two internal state variables are automatically generated and recorded as x f1 and x f2 , which is used to characterize the dynamic evolution of the aerodynamic hysteresis effect. Through these two state variables, the state space model of the second-order filter is constructed. The input of the state space model is:

[0084] ;

[0085] in, Represents state variables x f1 The rate of change over time, Represents state variables x f2 The rate of change over time reflects the dynamic response of the system under the aerodynamic hysteresis effect.

[0086] The input to the state-space model is: ;

[0087] According to the constructed state space model, the time domain expressions of the aerodynamic force vector are:

[0088] The aerodynamic force expression is: ;

[0089] The aerodynamic torque expression is: ;

[0090] in, is the dimensionless distance from the midpoint to the elastic axis, s p For the exhibition length.

[0091] Decompose the transformed aerodynamic force vector in the state space model into non-circular components F nc and circulation part Fc , the non-circulating part is the steady term in the aerodynamic component, which is represented by the acceleration term and the velocity term:

[0092] ;

[0093] in:

[0094] ;

[0095] The circulation part is the unsteady term in the aerodynamic component, which is represented by the filter state variables, velocity terms and displacement terms:

[0096] ;

[0097] Where:

[0098] .

[0099] Multiply the part in the brackets with the part outside the brackets in the circulation part expression and write it in matrix form to construct the time domain state space equation. The circulation part is obtained as follows:

[0100] ;

[0101] in:

[0102] .

[0103] The circulation part and non-circulation part of the aerodynamic force vector transformed in the state space model are brought into the flutter equation of S1 to form a complete state space equation:

[0104] ;

[0105] in: , , .

[0106] The filtering state variables are also written in the form of state equations: ;

[0107] in: .

[0108] The physical quantities used to constitute the non-circular part and the circular part are combined into a state vector, which contains the displacement ( h , ),speed( , ) and filter state variables ( x f1 , x f2 ), then the state vector is: , at this time the time domain state space equation of the entire system is: ;

[0109] in, ,

[0110] .

[0111] The fourth-order Runge-Kutta method is used to solve the time domain state space equation, and the time response curves of the heave displacement and pitch angle are extracted from the state vector solution:

[0112] ;

[0113] ;

[0114] ;

[0115] ;

[0116] ;

[0117] in, is the time step, is the time point of the nth iteration, is the state vector of the nth iteration, k 1 to k 4 is the four incremental steps of the fourth-order Runge-Kutta method calculation process. Each step depends on the previous state and the right-hand side function of the state equation. The state vector is updated by weighted average incremental steps, and the time series data of the heave displacement and pitch angle are extracted from the state vector.

[0118] The time response of generalized coordinates is calculated in the time domain to determine the aeroelastic stability of the system. Three different inflow velocities of 100m / s, 103m / s, and 106m / s are taken to calculate the aeroelastic response of the blade tip. The velocities are respectively less than the flutter critical speed, the flutter critical speed, and greater than the flutter critical speed. When the system is in a critical flutter state, the motion form of the system is a constant amplitude oscillation; when it is higher than the critical speed, the motion form of the system is a divergent form; when it is lower than the critical speed, the motion form of the system is an attenuated form; Figure 3 and Figure 4 shown.

[0119] As can be seen from the figure, the generalized coordinate response during critical flutter is a standard simple harmonic motion. At 106 m / s, the coordinates diverge, indicating that the system has experienced flutter. At 103 m / s, the coordinates converge, tending to the equilibrium position after the disturbance disappears. The time-domain aeroelastic response verifies the accuracy of the aeroelastic model calculation results.

[0120] The proposed time-domain flutter calculation method achieves efficient analysis and assessment of offshore wind turbine blade flutter by converting frequency-domain aerodynamic forces to the time domain and combining it with state-space methods. This method uses intuitive response curves to determine stability, boasting a solid theoretical foundation and a streamlined calculation process. It provides a reliable solution for flutter risk assessment during blade design and condition monitoring during operation, and has significant engineering application value for improving the safety and reliability of offshore wind turbine equipment.

[0121] The foregoing description is merely a preferred embodiment of the present application and is not intended to limit the present application. Various modifications and variations are readily apparent to those skilled in the art. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present application shall be included within the scope of protection of the present application.

[0122] Although the above describes the specific implementation methods of the present invention, it does not limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art on the basis of the technical solution of the present invention without creative work are still within the scope of protection of the present invention.

Claims

1. A method for determining the flutter state of a binary airfoil blade structure, characterized in that: The following processes are included: S1, obtaining the blade physical parameters in the corresponding scenario, establishing a blade flutter physical model based on the selected binary airfoil model, calculating the system kinetic energy and system potential energy based on the blade flutter physical model, substituting the system kinetic energy and system potential energy into the Lagrange equation to obtain the flutter equation including the mass matrix, stiffness matrix, damping matrix and aerodynamic force vector; wherein the aerodynamic force vector is calculated using unsteady aerodynamic theory, including unsteady aerodynamic forces and unsteady aerodynamic moments using the Theodorsen function; S2, the aerodynamic force vector of the S1 flutter equation is converted into a transfer function of a second-order system using the rational function approximation method, and the transfer function is converted into the state equation form through Laplace transform to construct the state space model of the second-order system; S3, decomposes the transformed aerodynamic force vector in the state-space model into a non-circulating part and a circulating part, introduces aerodynamic state variables including displacement, velocity and filter state, and constructs the time-domain state-space equation based on the flutter equation of S1; S4 inputs the real-time wind speed value into the time-domain state-space equation, solves the time-domain state-space equation using the fourth-order Runge-Kutta method, outputs the response curves of the heave displacement and pitch angle, and determines the flutter state based on the oscillation shape of the curve.

2. The method for determining the flutter state of a dual-element airfoil blade structure according to claim 1, wherein: The blade physical parameters include cross-sectional mass m , swing stiffness EI , torsional stiffness GJ , Section distance to blade root radius r , center of mass offset coefficient x a , chord length of airfoil section c and half chord length b=c / 2 ; in The elastic axis is established at the shear center of the blade, and the center of mass is located behind the shear center. x a b where x a is the mass center offset coefficient, the aerodynamic force and aerodynamic moment act on the aerodynamic center, and the aerodynamic center is at a distance from the leading edge of the blade. c / 4 The translation displacement of the blade reflecting the up and down vibration h , the downward direction is positive, reflecting the angular displacement of the blade torsional vibration is α , looking up into the wind is positive; The required parameters are obtained through blade structure design data.

3. The method for determining the flutter state of a dual-element airfoil blade structure according to claim 1, wherein: The system kinetic energy in S1 is the sum of translational kinetic energy and rotational kinetic energy. The translational kinetic energy is given by the cross-sectional mass. m , translational speed and the static mass moment of the airfoil per unit length about the elastic axis S EA Determine; the rotational kinetic energy is determined by the moment of inertia around the elastic axis J EA and angular velocity Determine; the system potential energy is the sum of the swing potential energy and the torsional potential energy. The swing potential energy is determined by the swing stiffness k h and translational displacement h Determine; torsional potential energy is determined by torsional stiffness and angular displacement Sure; Inputting the system kinetic energy and system potential energy expressions into the Lagrange equation, the flutter equation matrix form is derived: ； in, is the acceleration in the sinking and floating direction, is the pitch angular acceleration, g h , is the damping coefficient, - L is the unsteady aerodynamic force, is the unsteady aerodynamic moment, - L and Together they form the aerodynamic force vector, which includes the air density , incoming flow velocity V and Theodorsen function C ( k ).

4. The method for determining the flutter state of a dual-element airfoil blade structure according to claim 1, wherein: The specific process of S2 includes: The rational function approximation method is used to convert the aerodynamic force vector of the S1 flutter equation into a transfer function of a second-order system. The rational function approximation method uses a simple rational function to approximate the complex aerodynamic force function. The numerator describes the response of the aerodynamic force to the motion speed; the denominator describes the aerodynamic hysteresis effect. For a two-dimensional airfoil, the rational function approximation of the Theodorsen function is: ； Where, P 1 is the coefficient of the quadratic term in the numerator, and its rational function approximation universal constant is 0.006825. P 2 is the coefficient of the first-order term in the numerator, and its rational function approximation universal constant is 0.1080075. Z 1 is the coefficient of the quadratic term in the denominator, and its rational function approximation universal constant is 0.01365. Z 2 is the coefficient of the first-order term in the denominator, and its rational function approximation universal constant is 0.3455. V is the wind speed, b is the half chord length, s is the Laplace variable; After the Theodorsen function is expressed as a transfer function, the unsteady term of the aerodynamic force vector is described by the transfer function and a linear combination of displacement and velocity. The input of the transfer function is defined as the motion combination of the blade, that is, the linear combination of displacement and velocity, and the output is the unsteady term of the aerodynamic force component. According to the dynamic characteristics of the transfer function, two internal state variables are automatically generated and recorded as x f1 and x f2 , which is used to characterize the dynamic evolution process of the aerodynamic hysteresis effect. Through these two state variables, the state space model of the second-order filter is constructed. The output of the state space model is the time domain expression of the unsteady term of the aerodynamic force vector.

5. The method for determining the flutter state of a dual-element airfoil blade structure according to claim 1, wherein: The specific process of S3 includes: Decompose the transformed aerodynamic force vector in the state space model into non-circular components F c and circulation part F nc , the non-circulating part is the steady term in the aerodynamic component, which is composed of the linear combination of acceleration term and velocity term; the circulating part is the unsteady term in the aerodynamic component, which is composed of the linear combination of filter state variables and displacement and velocity term; The physical quantities used to constitute the non-circulating part and the circulating part are combined into a state vector. The state vector contains displacement, velocity and filter state variables. According to the flutter equation of S1 and the state space model of S2, the time domain state space equation of the entire system is constructed: ； in, A ( V ) is a matrix with 6 rows and 6 columns, ,matrix M is the sum of the mass matrix and the mass matrix derived from the non-circular aerodynamics, the matrix D is the product of the damping matrix and the damping matrix derived from the non-circular aerodynamics and the wind speed. K is the sum of the stiffness matrix and the product of the stiffness matrix derived from the non-circular aerodynamic force and the square of the wind speed, K a is the additional stiffness matrix caused by the circulation aerodynamic force, D a is the additional damping matrix caused by the circulation aerodynamic force, A a is the filter state matrix caused by the circulation aerodynamic force.

6. The method for determining the flutter state of a dual-element airfoil blade structure according to claim 1, wherein: The specific process of S4 is: Enter real-time wind speed V In the time-domain state-space equation, the fourth-order Runge-Kutta method is used to iteratively solve the time series of the state vector, and the time response curves of the heave displacement and pitch angle are extracted from the state vector solution: ； ； ； ； ； in, is the time step, is the time point of the nth iteration, is the state vector of the nth iteration, k 1 to k 4 is the four incremental steps of the fourth-order Runge-Kutta method calculation process. Each step depends on the previous state and the right-hand side function of the state equation. The state vector is updated by weighted average incremental steps, and the time series data of the heave displacement and pitch angle are extracted from the state vector.

7. The method for determining the flutter state of a dual-element airfoil blade structure according to claim 6, wherein: The specific method of determining the vibration state in S4 is: If the response curve shows attenuated oscillation, the system is stable; if the response curve shows constant amplitude oscillation, it is judged to be in a critical flutter state; if the response curve shows exponential divergence, the system is unstable.

8. The method for determining the flutter state of a dual-element airfoil blade structure according to claim 7, wherein: The wind speed value is scanned incrementally with a set wind speed interval step size within a preset wind speed range; the response curve is solved for each wind speed value, and the wind speed corresponding to the first occurrence of equal-amplitude oscillation is identified as the flutter critical wind speed.

Images