Long flexible blade structure-pneumatic nonlinear dynamics reduced order simulation method and system

Abstract

The invention belongs to the technical field of dynamic simulation, and provides a long flexible blade structure-aerodynamic nonlinear dynamic reduced-order simulation method and system, the technical scheme is that a dynamic model considering blade composite material characteristics and variable cross-section characteristics at the same time is established based on ANCF, and an aerodynamic force calculation method on a cross-section frame is deduced; the method specifically comprises the following steps: calculating a chord length vector on an airfoil section based on position information of a front edge and a rear edge on the airfoil section, projecting a relative velocity of a three-dimensional blade unit to the airfoil section to obtain a two-dimensional relative velocity, and calculating to obtain an attack angle; the aerodynamic force acting on the blade airfoil section is calculated by combining the chord length vector, the two-dimensional relative speed and the attack angle, the differences of the elastic center, the mass center and the aerodynamic center on the composite material blade airfoil section are considered, aerodynamic force translation is corrected, and the analysis precision of a blade dynamic model is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of dynamics simulation, and in particular relates to a method and system for reducing the order of aerodynamic nonlinear dynamics of a long flexible blade structure. Background Art

[0002] The statements in this section merely provide background information related to the present invention and do not necessarily constitute prior art.

[0003] As core components for energy capture, the dynamic characteristics of wind turbine blades directly impact the turbine's power generation efficiency, structural safety, and operational lifespan. However, as composite blades increase in size and experience the strong coupling of aerodynamic, inertial, and elastic forces in complex wind fields, they can exhibit significant instabilities such as large deformation and nonlinear vibration. Establishing efficient, high-fidelity blade dynamics models to accurately predict their dynamic response has become a key challenge in improving wind turbine reliability and design.

[0004] Currently, there are two main approaches to modeling wind turbine blade dynamics: three-dimensional finite element models based on blade solid or shell elements, and beam models based on beam elements. While these three-dimensional finite element models offer high accuracy, their high computational cost and complex modeling process limit their application. In contrast, beam models based on beam elements balance computational speed and accuracy, making them the most widely used dynamic modeling approach. Wind turbine blades are composed of components made of different materials, such as the skin and reinforcing web. These materials are anisotropic along the blade's span, and their geometry undergoes significant variations along the blade's span, leading to nonlinear variations in linear density, bending-torsional stiffness, and aerodynamic coefficients along the span. Traditional beam model simulations based on beam elements, which use linear interpolation to calculate material and aerodynamic coefficients, are unable to account for the nonlinear variations in material and aerodynamic parameters introduced by streamlined blade geometry. Typical beam models include the Euler-Bernoulli beam model, the Timoshenko beam model, and the geometrically accurate beam model. Among them, the Euler-Bernoulli beam model does not consider the torsional degree of freedom, ignores shear deformation, and is not accurate enough, resulting in the calculated deformation being lower than the actual value; the Timoshenko beam model takes shear deformation into account, and the calculated deformation is closer to the actual situation, but it has shortcomings when dealing with large deformations or finite cross-sectional rotations; the geometrically precise beam model assumes that the frame field on the beam unit cross-section is an orthogonal frame field, so the rotation variables need to be interpolated, which makes the calculation more complicated and has singularity problems.

[0005] The Absolute Nodal Coordinate Formulation (ANCF) method uses gradient fields instead of orthogonal frames. It directly characterizes element deformations using node positions and gradient vectors in a global coordinate system, avoiding singularity issues. Its constant mass matrix significantly improves the efficiency of solving dynamic equations. For flexible, slender structures like blades, the ANCF method simultaneously considers bending, axial deformation, and their coupled effects, making it more suitable for dynamic analysis of large deformations and rotations.

[0006] The inventors discovered that some people have used the absolute node coordinate method to analyze the blade dynamics model. One method is to construct a three-dimensional blade dynamics model based on plate units, and the other method is to combine the absolute node coordinate method with the finite element method to construct beam units to establish a blade model. However, the former does not take into account the composite material properties of the blade and has a huge amount of calculation. The latter finite element model has too many nodes and requires a complex numerical integration process after export, which is time-consuming and difficult to ensure accuracy. In addition, the blade dynamics models constructed by these two methods do not take into account the differences between the center of mass, elastic center and aerodynamic center within the airfoil section. Instead, they assume that aerodynamic forces act on the center of mass line, thereby ignoring the torque term attached to the force translation, which further reduces the accuracy of the dynamic model. Summary of the Invention

[0007] In order to solve at least one technical problem existing in the above-mentioned background technology, the present invention provides a long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method and system, which takes into account the composite material characteristics and variable cross-section characteristics of the blade at the same time, and the model accuracy is relatively high.

[0008] In order to achieve the above object, the present invention adopts the following technical solutions:

[0009] A first aspect of the present invention provides a method for reducing the order of aerodynamic nonlinear dynamics of a long flexible blade structure, comprising the following steps:

[0010] A cross-sectional frame is established at the centroid of the blade airfoil section, and the position information of the leading and trailing edges of the blade airfoil section is extracted based on the acquired blade data.

[0011] The fan blade is divided into several units, the kinetic energy of each blade unit is calculated, and the unit mass matrix of the corresponding unit is obtained based on the kinetic energy calculation of each blade unit;

[0012] According to the blade airfoil section frame vector, the quasi-curvature vector related to bending and torsion on the beam element is determined, and the constant coefficient tensor of each blade element is calculated in combination with the curvature vector. The generalized elastic force of each element of the wind turbine blade is obtained by combining the constant coefficient tensor of each blade element with the generalized coordinates;

[0013] The chord length vector on the airfoil section is calculated based on the position information of the leading edge and the trailing edge on the airfoil section. The relative velocity of the blade unit in three dimensions is projected onto the airfoil section to obtain the two-dimensional relative velocity. The angle of attack is calculated. The aerodynamic force acting on the blade airfoil section is calculated by combining the chord length vector, the two-dimensional relative velocity and the angle of attack. The generalized external force acting on the blade unit is obtained according to the absolute nodal coordinate method.

[0014] The mass matrix, generalized elastic force matrix and generalized external force matrix of each unit are assembled to obtain the wind blade dynamic model of the slender beam model in the absolute nodal coordinate system.

[0015] Furthermore, the constant coefficient tensor of the blade unit is calculated by combining the curvature vector, specifically based on the quasi-curvature vector related to bending and torsion on the beam unit, and the constant coefficient tensor of the blade unit is obtained according to the shape function and its partial derivative with respect to the arc length coordinate x and the stiffness parameter function of the blade unit along the arc length coordinate x of the center of mass line.

[0016] Furthermore, the calculation formula for the quasi-curvature vector related to bending and torsion on the beam element is:

[0017]

[0018] Among them, k y With k z Curvature-like vectors related to bending and torsion, r x represents the partial derivative of the centroid function r(x,t) with respect to x, which is related to the beam arc length coordinate x. y With r z are the two direction vectors of the beam section composed of the gradient field, r xy For r x The partial derivative with respect to y, r xz For r x The partial derivative with respect to z.

[0019] Furthermore, the calculation formula of the constant coefficient tensor of each blade unit is:

[0020]

[0021] Among them, C, C W ,C yy ,C zz is the stiffness parameter function of the beam element along the arc length coordinate x of the centroid line, s α , s β ,y α ,y β , z α , z β represents the shape function, α, β=1,2,3...,10, and the subscript x represents the partial derivative of the shape function with respect to x, n-1 is the left node coordinate of unit n, xn is the right node coordinate of element n.

[0022] Furthermore, the calculation formula for the generalized elastic force of each unit of the fan blade is:

[0023]

[0024] in, c is the viscoelastic damping coefficient, δU is the virtual work of the elastic force obtained according to the virtual work theorem, δq α is the virtual displacement, q α ,q β ,q μ ,q ν is a generalized coordinate, and the subscript range is consistent with the range of α, α=1,2,3...,10, is the first derivative of the generalized coordinate with respect to time.

[0025] Furthermore, the calculation formula of the angle of attack is:

[0026]

[0027] c(x)=(y te -y le )r y +(z te -z le )r z ,

[0028]

[0029] Among them, y te 、y le 、z te 、z le is the coordinate of the leading edge and trailing edge of the airfoil section in the section coordinate system, v rel is the relative velocity of the blade unit in three dimensions, is the two-dimensional relative velocity, c is the chord length vector on the airfoil section, v ∞ is the wind speed at infinity, is the relative velocity of the leading edge of the blade airfoil, G=(r x ,r y ,r z ) T , θ is With r x The angle, r x represents the partial derivative of the centroid function r(x,t) with respect to x, which is related to the beam arc length coordinate x. y With r z are the two direction vectors of the beam section composed of the gradient field, α=1,2,3...,10.

[0030] A second aspect of the present invention provides a long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation system, comprising:

[0031] A position information extraction module is used to establish a cross-section frame at the centroid of the blade airfoil section and extract the position information of the leading edge and the trailing edge of the blade airfoil section in combination with the acquired blade data;

[0032] A mass matrix calculation module is used to divide the wind turbine blade into a number of units, calculate the kinetic energy at each blade unit, and obtain the unit mass matrix of the corresponding unit based on the kinetic energy calculation at each blade unit;

[0033] A generalized elastic force calculation module is used to determine the quasi-curvature vector related to bending and torsion on the beam element based on the blade airfoil section frame vector, calculate the constant coefficient tensor of each blade element based on the curvature vector, and obtain the generalized elastic force of each element of the wind turbine blade by combining the constant coefficient tensor of each blade element with the generalized coordinates;

[0034] A generalized external force calculation module is used to calculate the chord length vector on the airfoil section based on the position information of the leading edge and trailing edge on the airfoil section, project the relative velocity of the blade unit in three dimensions onto the airfoil section to obtain a two-dimensional relative velocity, calculate the angle of attack, calculate the aerodynamic force acting on the blade airfoil section by combining the chord length vector, two-dimensional relative velocity, and angle of attack, and obtain the generalized external force acting on the blade unit according to the absolute node coordinate method;

[0035] The dynamic simulation model construction module is used to assemble the mass matrix, generalized elastic force matrix and generalized external force matrix of each unit to obtain the wind turbine blade dynamic model of the slender beam model in the absolute node coordinate system.

[0036] A third aspect of the present invention provides a computer-readable storage medium.

[0037] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method as described above.

[0038] A fourth aspect of the present invention provides a computer device.

[0039] A computer device comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the steps in the above-mentioned long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method are implemented.

[0040] A fifth aspect of the present invention provides a program product.

[0041] A program product is a computer program product, comprising a computer program. When the computer program is executed by a processor, the steps in the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method as described above are implemented.

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

[0043] 1. The present invention provides a method and system for modeling wind turbine blade dynamics based on the absolute nodal coordinate method. By calculating the blade's mass and stiffness matrices and equating the composite blade model's attribute parameters to functions that vary along the arc-length coordinate, the generalized elastic forces of each unit of the wind turbine blade are calculated. The aerodynamic forces acting on the blade airfoil are then translated to the center of mass of the airfoil cross-section, ensuring that the blades are subjected to a common force. The generalized forces acting on the blades at the same point are then substituted into the absolute nodal coordinate method for solution, thereby expressing the blade's nonlinear characteristics in a more robust and efficient manner. Compared to the traditional beam element method, this method considers the composite material properties and structural force characteristics, improving model accuracy.

[0044] 2. The present invention achieves a vector representation of the chord length of the blade airfoil section on the section frame by obtaining the arc-length-varying function of the leading and trailing edges of the airfoil section within the section gradient frame. The material properties of the blade are then used to construct material functions along the arc-length coordinate, including linear density, section inertia, and stiffness functions, thereby solving the blade's mass and stiffness matrices. In the global coordinate system, the blade is discretized into multiple finite elements, and the shape functions are constructed using the position and gradient information of the nodes at each end of each element. The axial deformation, curvature, and strain state of the element are derived, and the generalized elastic force is calculated in conjunction with the material functions. Simultaneously, the angle of attack is calculated by projecting the relative velocity formed by the incoming wind speed at infinity and the blade structural velocity onto the airfoil section. Existing tables of lift, drag, and moment coefficients are then used to interpolate and solve for the aerodynamic forces at Gaussian points. This is combined with gravity to obtain the generalized external forces on the blade. Compared to the finite element method, the present method is computationally faster.

[0045] Advantages of additional aspects of the present invention will be given in part in the following description and in part will be obvious from the following description, or will be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] The accompanying drawings, which constitute a part of the present invention, are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations on the present invention.

[0047] Figure 1 This is a flow chart of a method for reducing the order of aerodynamic nonlinear dynamics of a long flexible blade structure provided by an embodiment of the present invention;

[0048] Figure 2is a schematic diagram of an absolute node coordinate system for a slender beam provided by an embodiment of the present invention;

[0049] Figure 3 is the lift and drag acting on the airfoil cross section provided by the embodiment of the present invention;

[0050] Figure 4 This is the displacement of the blade tail end point in three directions when the pitch angle is 0°, as shown in the calculation example 1 provided in the embodiment of the present invention;

[0051] Figure 5 This is a calculation example 1 provided by an embodiment of the present invention, showing deformation of the blade tail end in the x and y directions at different pitch angles;

[0052] Figure 6 The deformation and error rate of the blade tail end in the y direction at different pitch angles in the calculation example 1 provided by the embodiment of the present invention are as follows;

[0053] Figure 7 This is the displacement of the blade tail end point in three directions when the pitch angle is 90° in the example 2 provided by the embodiment of the present invention;

[0054] Figure 8 This is the deformation of the blade tail end point in the x and y directions under different pitch angles in the second calculation example provided by the embodiment of the present invention;

[0055] Figure 9 This is the deformation and error rate of the blade tail end point in the y direction at different pitch angles in the second calculation example provided by the embodiment of the present invention. DETAILED DESCRIPTION

[0056] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0057] It should be noted that the following detailed descriptions are illustrative and intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present invention belongs.

[0058] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or combinations thereof.

[0059] Modeling wind turbine blades using 3D solid elements consumes enormous computational resources. During the initial iterative design phase, 1D beam elements are often used to simplify the 3D wind turbine blade model and quickly generate computational results. As wind turbines evolve, their rated power increases, and blades become longer to capture more wind energy.

[0060] Existing wind turbine designs use simplified one-dimensional Timoshenko beam elements to calculate wind turbine blades. However, Timoshenko beams are suitable for calculating short, thick beams. However, they are no longer accurate enough for the large blade lengths, displacements, and rotations of current wind turbines. One of the current challenges is establishing a unified dynamic model that can describe these large blade lengths, displacements, and rotations.

[0061] The present invention proposes a long-flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method based on the absolute node coordinate method. The one-dimensional beam element uses displacement and gradient vectors as a description of blade rotation and deformation, relying only on the slender body assumption. With as few approximations as possible, it retains the characteristics of wind turbine blade stretching, shearing, bending and torsion. Mathematical language can be used to describe situations involving geometric nonlinearity, such as large rotation and large deformation that may exist in existing blades, and is more universal.

[0062] Wind turbine blades are composed of composite components made of different materials, such as the outer skin and reinforcing webs. These materials are anisotropic along the blade's extension. Furthermore, the blade's geometric configuration undergoes significant changes along its span, leading to nonlinear variations in linear density, bending-torsional stiffness, and aerodynamic coefficients along the span. Traditional simulation methods that use linear interpolation to calculate material and aerodynamic coefficients are unable to describe the nonlinear variations in material and aerodynamic parameters introduced by streamlined changes in blade geometry. Accurately describing the changes in material aerodynamic properties caused by the blade's complex geometry is a second current challenge.

[0063] The present invention takes into account that although the frame of the absolute nodal coordinate method does not have the orthogonality assumption, for the problem of material anisotropy, the strain field described by the gradient frame can be converted to the material frame through matrix transformation for calculation. Based on the slender body assumption, the strain on each cross section is calculated, and the blade cross section composed of multiple materials is pre-integrated. Compared with the equivalent stiffness given by the Timoshenko beam, the method of the present invention can calculate a more accurate stiffness matrix, improving the calculation accuracy. Compared with linear interpolation, it is more suitable for describing the continuous changes in the material and aerodynamic properties of wind turbine blades.

[0064] The existing absolute node coordinate method introduces a gradient frame to describe the shearing, bending and torsion of the blade. Although it improves the calculation accuracy and avoids the problem of rotational singularity, it increases the number of generalized coordinates and increases the amount of calculation. How to improve the calculation efficiency of elastic force in simulation is the third existing problem.

[0065] In the process of calculating elastic potential energy, the present invention adopts a pre-integration method to calculate the coefficient matrix related to the elastic force in advance, separates the elastic parameters and generalized coordinates, derives the generalized elastic force and the corresponding Jacobian matrix in the form of tensor contraction, avoids Gaussian integral and numerical difference operations that occupy a large amount of computing resources, and greatly improves the calculation efficiency of the elastic force.

[0066] Example 1

[0067] like Figure 1 As shown, this embodiment provides a long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method, including the following steps:

[0068] Step 1: Establish a cross-section frame at the centroid of the blade airfoil section, and extract the position information of the leading edge and trailing edge of the blade airfoil section based on the acquired blade data;

[0069] It should be noted that, in this embodiment, when establishing the cross-section frame at the centroid of the blade airfoil cross-section, the cross-section frame is established at the centroid of the blade airfoil cross-section in combination with the geometric characteristics of the blade and the characteristics of the absolute node coordinate method;

[0070] Specifically, the geometric characteristics of the blade include the presence of geometric positions such as the center of mass, elastic center, aerodynamic center, and shear center on the blade airfoil section;

[0071] Specifically, if Figure 2 As shown, a slender beam of uniform material with a length of L, in the absolute node coordinate method proposed in this embodiment, the position r of any point p on the beam element at time t is p (x, y, z, t), can be described as formula (1):

[0072] r p (x,y,z,t)=r(x,t)+yr y (x,t)+zr z (x,t) (1),

[0073]

[0074] In formula (1), r(x,t) represents the centroid function related to the beam arc length coordinate x, r y (x,t) and r z (x, t) represents the two direction vectors of the beam section, which form a coordinate system. The coordinates of point p in this coordinate system are (y, z). α (x), y α (x) and z α(x) represents the interpolation shape function, and here we follow the Einstein summation convention, that is, the dot product of two vectors can be written as the sum of the corresponding multiplications of the internal elements. According to the Einstein summation convention, the summation sign is omitted for representation, q α is the generalized coordinate of the unit node, which is defined as shown in formula (3):

[0075]

[0076]

[0077] Where, x n-1 is the left node coordinate of unit n, x n is the right node coordinate of element n; is the coordinate of the midpoint of unit n, ξ is the normalized arc length coordinate within the unit, ξ=(xx n-1 ) / l n , l n represents the length of the beam element. In formula (4), only the shape function s is listed. α (x), y α (x), z α (x) is the non-zero term.

[0078] When extracting the position information of the leading edge and the trailing edge on the blade airfoil section, the acquired blade data can be combined with a direct combined search method or parameterized airfoil determination. This is a prior art and will not be described in detail in this embodiment.

[0079] Step 2: Divide the wind turbine blade into several units, calculate the kinetic energy of each blade unit, and obtain the unit mass matrix of the corresponding unit based on the kinetic energy calculation of each blade unit;

[0080] The specific steps include:

[0081] Step 201: Calculate the kinetic energy of each blade unit;

[0082] When the material of the beam is uniform and the centroid passes through the center of each section, the kinetic energy at the nth element is:

[0083]

[0084] Where G is the frame field of the center of gravity of the airfoil section, the symbol ^ represents the reference configuration, is the first derivative of the centroid function r(x,t) with respect to time related to the beam arc length coordinate x, and is the first-order derivative of the two direction vectors of the beam section composed of the gradient field with respect to time, ρ A is the linear density function of the beam element along the arc length coordinate x of the centroid line, J yy,J zz ,J yz is the cross-sectional mass inertia function of the beam element along the arc length coordinate x of the centroid line;

[0085] Step 202: Calculate the unit mass matrix of each blade unit based on the kinetic energy of the blade unit.

[0086] because Substituting into formula (5), we can get Based on this, the unit mass matrix at the unit is obtained I3 is a third-order unit matrix, and the specific calculation formula is shown in formula (7):

[0087]

[0088] Among them, ρ A ,J yy ,J zz ,J yz are the linear density function and cross-sectional mass inertia function of the beam element along the arc length coordinate x of the centroid line, s α , s β ,y α ,y β , z α , z β Represents shape function, α, β=1,2,3...,10.

[0089] Step 3: Based on the blade airfoil section frame vector, determine the quasi-curvature vector related to bending and torsion on the beam element. Combine the curvature vector to calculate the constant coefficient tensor of the blade element. Combine the constant coefficient tensor of the blade element and the generalized coordinates to obtain the generalized elastic force of each element of the wind turbine blade.

[0090] Specifically include:

[0091] Step 301: Determine the quasi-curvature vector related to bending and torsion on the beam element based on the gradient vector of the blade airfoil cross-section frame;

[0092]

[0093] Among them, k y With k z Curvature-like vectors related to bending and torsion, r x represents the partial derivative of the centroid function r(x,t) with respect to x, which is related to the beam arc length coordinate x. y With r z are the two direction vectors of the beam section composed of the gradient field, r xy For r x The partial derivative with respect to y, r xz For r x Partial derivative with respect to z;

[0094] Step 302: Based on the curvature vector related to bending and torsion on the beam element, the constant coefficient tensor of the blade element is obtained according to the shape function and its partial derivative with respect to the arc length coordinate x and the stiffness parameter function of the blade element along the arc length coordinate x of the centroid line.

[0095] In this embodiment, the constant coefficient tensor of the blade unit The calculation formula is:

[0096]

[0097] Where the subscript x represents the partial derivative of the shape function with respect to x, C, C W ,C yy ,C zz is the stiffness parameter function of the beam element along the arc length coordinate x of the centroid line;

[0098] Step 303: Finally, the elastic potential energy U and the unit generalized elastic force of the wind turbine blade at the nth unit are obtained by combining the generalized coordinates.

[0099] In this embodiment, the elastic potential energy U at the nth unit of the fan blade can be expressed by the virtual work of the elastic force, and the calculation formula is:

[0100]

[0101] The generalized elastic force at the nth element The calculation formula is:

[0102]

[0103] in, c is the viscoelastic damping coefficient, δU is the virtual work of the elastic force obtained according to the virtual work theorem, δq α is the virtual displacement, q α ,q β ,q μ ,q ν is a generalized coordinate, and the subscript range is consistent with the range of α, α=1,2,3...,10, is the first derivative of the generalized coordinate with respect to time.

[0104] Step 4: Calculate the chord length vector on the airfoil section based on the position information of the leading edge and trailing edge on the airfoil section. Combine the relative velocity of the blade unit in three dimensions and project it onto the airfoil section to obtain the two-dimensional relative velocity. Calculate the angle of attack. Combine the chord length vector, two-dimensional relative velocity, and angle of attack to calculate the aerodynamic force acting on the blade airfoil section. Calculate the generalized external force acting on the blade unit using the absolute nodal coordinate method.

[0105] The specific steps include:

[0106] Step 401: First, calculate the chord length vector on the airfoil section based on the leading and trailing edge position coordinates and gradient vectors of the airfoil section frame;

[0107] The chord length vector c(x) on the airfoil section points from the leading edge to the trailing edge, and its position in the frame G = (r x ,r y ,r z ) T The vector expression on is:

[0108] c(x)=(y te -y le )r y +(z te -z le )r z (12),

[0109] Among them, y te 、y le 、z te 、z le are the coordinates of the leading edge and trailing edge of the airfoil section in the section coordinate system;

[0110] Step 402: The relative velocity v of the blade unit in three dimensions can be obtained based on the infinite incoming wind speed and the velocity of the blade airfoil leading edge. rel , and project it onto the airfoil section to obtain the two-dimensional relative velocity v rel ;

[0111] Aerodynamic force is the main source of power acting on the fan blades. Assume that the wind speed at infinite distance is v ∞ , then its relative speed to the leading edge of the blade airfoil is:

[0112]

[0113] in, is the relative velocity of the leading edge of the blade airfoil;

[0114] The relative velocity needs to be projected onto the airfoil section to calculate the angle of attack. The frame field of the center of gravity of the airfoil section is G = (r x ,r y ,r z ) T ,like Figure 3 As shown. The cross-section attack angle is in the two-dimensional frame (r y ,r z ) T According to the above calculation, the incoming flow velocity of the two-dimensional airfoil section is:

[0115]

[0116] in

[0117] Step 403: According to the two-dimensional relative velocity Calculate the angle of attack with the chord length vector;

[0118] The angle of attack is the angle between the chord of the airfoil section and the relative velocity, and its magnitude and direction are:

[0119]

[0120] Where θ is With r x If θ is less than π / 2, the angle of attack is positive, otherwise it is negative.

[0121] Step 404: Finally, the aerodynamic force acting on the blade airfoil section is calculated based on the chord vector, relative velocity, air density, and aerodynamic coefficient;

[0122] The drag acting on the airfoil section is shown in formula (16):

[0123]

[0124] Where ρ is the air density, C D (x,α) is the drag coefficient, and the same C L (x,α),C M (x,α) are the lift coefficient and moment coefficient.

[0125] The lift direction on the airfoil is perpendicular to the drag direction, i = r y ×r z / ||r y ×r z || is the defined rotation axis, and the rotation angle is π / 2. is the rotation matrix obtained according to the Rodris formula, where is the corresponding antisymmetric matrix. If the rotation angle is -π / 2, the rotation matrix is

[0126] According to the rotation matrix and the drag formula on the airfoil section, the lift on the airfoil section can be obtained as:

[0127]

[0128] The combined equations (16) and (17) can then be used to obtain the resultant force on the airfoil section, as shown in equation (18):

[0129]

[0130] C DL =(C L (x,α)A+CD (x,α)I3) (19),

[0131] The moment of the resultant force on the airfoil section relative to the leading edge or 1 / 4 of the chord is called the pitching moment. According to the pitching moment coefficient, the pitching moment on the airfoil can be calculated as:

[0132]

[0133] In formula (20), the resultant aerodynamic force acts at a quarter of the chord length, which is not at the same position as the center of mass. f is the translational torque generated by the aerodynamic force translated to the center of mass.

[0134] Step 405: Combine the aerodynamic force and gravity on the airfoil section and the blade element shape function to obtain the generalized external force acting on the blade element according to the absolute node coordinate method;

[0135] The specific calculation process is:

[0136] The generalized external force at the nth unit can be obtained by the principle of virtual work:

[0137]

[0138] By taking the partial derivative of the virtual work δW in formula (21) with respect to the virtual displacement δq, the generalized external force acting on the blade element can be obtained:

[0139]

[0140] where f(x, t) and t(x, t) are the distributed force and moment of the beam element along the arc length coordinate x of the centroid line, where the distributed force is the sum of gravity and aerodynamic force, δr(x, t) is the virtual displacement, and δΠ(x, t) is an infinitesimal frame rotation.

[0141] Step 5: Assemble the mass matrix, generalized elastic force matrix, and generalized external force matrix of each unit to obtain the mass matrix and generalized force matrix of the entire model. Apply the constraint equations to obtain the wind turbine blade dynamics model of the slender beam model in the absolute nodal coordinate system.

[0142] Based on the calculated loads of all blade units and the mass of all units, they are distributed to the nodes of the blade units according to the Boolean matrix. The loads and masses acting on the nodes are substituted into the motion equation of the blade to complete the modeling. The wind turbine blade dynamic model based on the absolute node coordinate method is obtained. The final dynamic equation is shown in the following equation (23):

[0143]

[0144] Where Φ is the constraint vector, Φ q is the Jacobian matrix of Φ with respect to q, λ is the Lagrange multiplier, Represents the acceleration vector of all nodes.

[0145] To further integrate and solve the blade motion equations, the absolute nodal coordinate method increases the number of generalized coordinates, which reduces the efficiency of elastic force calculations. Therefore, in the elastic potential energy calculation process, the present invention uses a pre-integration method to precalculate the coefficient matrix related to the elastic force. This separates the elastic parameters from the generalized coordinates, avoids Gaussian integration and numerical difference operations within the integrator, and significantly improves the efficiency of elastic force calculations. Finally, the displacement, velocity, and acceleration of all blade nodes within the integration time domain are obtained, thereby analyzing the blade's dynamic response.

[0146] Furthermore, in order to verify the accuracy of the blade dynamics model established by the present invention, this embodiment compares the dynamics model simulation results with the results of the software Bladed. Some parameters of the case blade data are shown in Table 1 below:

[0147] Table 1 Parameters of blade segments

[0148]

[0149]

[0150] Example 1: Deadweight horizontal shutdown condition

[0151] No wind load is applied to the wind turbine blades. One of the blades is placed horizontally with a pitch angle of 0 degrees. The blade is only affected by gravity in the shutdown condition. The displacement of its tail end point is as follows: Figure 4 As shown, the displacement in the x direction is the deformation perpendicular to the blade rotation plane, the displacements in the y and z directions are within the rotation plane, and z is along the blade span direction.

[0152] Depend on Figure 5 As can be seen, the deformation of the blade tip in the x and y directions obtained by solving the dynamic model at different pitch angles is compared with the software results. The dashed line is the Bladed simulation result, and the circled line is the simulation result of the blade dynamic model developed based on ANCF. The simulation results of the Bladed deformation of the blade tip at a pitch angle of 0 are shown in Table 2.

[0153] Table 2 Deformation of the Bladed Blade End Node

[0154] Working conditions x direction / m y direction / m z direction / m Deadweight level stop 0.447127 -1.1255 0.024959 Wind load horizontal shutdown 1.01603 -1.17867 0.072387

[0155] Since the torsional stiffness has been corrected, and the deformation in the x and z directions is smaller than that in the y direction under the shutdown condition, the deformation in the y direction is mainly analyzed here. The y-direction deformation obtained by the dynamic model is smaller than the Bladed result. The y-direction deformation of the blade tail node of the former is 1.1147m, and the deformation result of the latter is 1.1255m, with an error of 0.96%, showing a good agreement. The deformation error in the y direction of the two at different pitch angles is as follows Figure 6 It can be seen that the dynamic simulation results are consistent with the software results.

[0156] Example 2: Wind load horizontal shutdown condition

[0157] Apply wind load to the wind turbine, with a uniform wind speed of 10 m / s, set the pitch angle to 90 degrees, i.e. feathering state, and the viscoelastic damping coefficient to 0.02. Other conditions remain unchanged. The displacement of the blade tail end under shutdown conditions is as follows: Figure 7 As shown, the displacement in the x direction is the deformation perpendicular to the blade rotation plane, the displacements in the y and z directions are within the rotation plane, and z is along the blade span direction.

[0158] When the blade is stopped, it is continuously subjected to a steady wind of 10 m / s. The steady wind comes perpendicular to the rotation plane. Under the action of the wind load, the blade produces a greater deformation along the x direction. The deformation of the blade will increase compared to the no-wind condition. Figure 8 It can be seen that the variation along the x-direction obtained by the dynamic model simulation is quite different from the Bladed simulation result. The reason is that the torsional stiffness in the material properties of the dynamic model is corrected by the axis shift, and the resulting error leads to a large difference with the Bladed result. At the same time, the deformation in the y-direction is also slightly increased compared to the no-wind condition, as shown in Figure 2. Figure 9 It can be seen that the simulation results of the blade model developed based on ANCF are very consistent with the Bladed results. However, due to the errors in material properties, the error between the two is larger under wind conditions than under windless conditions, but most of them remain within 10%.

[0159] Based on the absolute nodal coordinates (ANCF), the present invention develops a high-precision simplified engineering model that takes into account the geometric nonlinearity of wind turbine blades, achieving efficient calculation of the three-dimensional nonlinear deformation problem of long and flexible blades based on one-dimensional beam elements. This method comprehensively considers the nonlinear changes in the linear density, stiffness, and aerodynamic center position of the blade along the span direction, constructs the material and structural function equations along the arc length direction based on the actual two-dimensional airfoil cross-section of the blade, and determines the airfoil cross-section chord vector through the cross-sectional coordinates of the leading and trailing edges of the airfoil cross-section to calculate the angle of attack and aerodynamic load. Finally, the axis shift theorem is used to correct the stiffness data at the shear center and elastic center to the center of mass, so that the entire dynamic model is established based on a single principal axis.

[0160] The results show that the wind turbine blade dynamic model and aerodynamic load calculation method based on the absolute nodal coordinate method can accurately describe the blade's nonlinear deformation and aerodynamic characteristics. Under the conditions of no wind load and zero pitch angle, the deformation of the single blade's tail node in the y-direction is 1.1147m, while the software simulation result is 1.1255m, with an error of 0.96%. Under wind load, the deformation in the x-direction has a larger error than the Bladed software result due to the partial error introduced by the torsional stiffness offset correction, and it changes with the pitch angle. However, the deformation error in the y-direction remains within 10%.

[0161] The blade dynamics model developed in the present invention has very good guiding significance and reference value for the initial simple dynamic analysis of the blade, and can help to better design and optimize the blade.

[0162] Example 2

[0163] This embodiment provides a long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation system, including:

[0164] A position information extraction module is used to establish a cross-section frame at the centroid of the blade airfoil section and extract the position information of the leading edge and the trailing edge of the blade airfoil section in combination with the acquired blade data;

[0165] A mass matrix calculation module is used to divide the wind turbine blade into a number of units, calculate the kinetic energy at each blade unit, and obtain the unit mass matrix of the corresponding unit based on the kinetic energy calculation at each blade unit;

[0166] A generalized elastic force calculation module is used to determine the quasi-curvature vector related to bending and torsion on the beam element based on the blade airfoil section frame vector, calculate the constant coefficient tensor of each blade element based on the curvature vector, and obtain the generalized elastic force of each element of the wind turbine blade by combining the constant coefficient tensor of each blade element with the generalized coordinates;

[0167] A generalized external force calculation module is used to calculate the chord length vector on the airfoil section based on the position information of the leading edge and trailing edge on the airfoil section, project the relative velocity of the blade unit in three dimensions onto the airfoil section to obtain a two-dimensional relative velocity, calculate the angle of attack, calculate the aerodynamic force acting on the blade airfoil section by combining the chord length vector, two-dimensional relative velocity, and angle of attack, and obtain the generalized external force acting on the blade unit according to the absolute node coordinate method;

[0168] The dynamic simulation model construction module is used to assemble the mass matrix, generalized elastic force matrix and generalized external force matrix of each unit to obtain the wind turbine blade dynamic model of the slender beam model in the absolute node coordinate system.

[0169] It should be noted that the specific implementation method of the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation system in the embodiment of the present invention is similar to the specific implementation method of the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method in the embodiment of the present invention. Please refer to the description of the method part for details. In order to reduce redundancy, it will not be repeated here.

[0170] Example 3

[0171] This embodiment provides a computer-readable storage medium having a computer program stored thereon. When the program is executed by a processor, the steps in the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method as described above are implemented.

[0172] Example 4

[0173] This embodiment provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, the steps in the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method described above are implemented.

[0174] Example 5

[0175] This embodiment provides a program product, which is a computer program product, including a computer program. When the computer program is executed by a processor, the steps in the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method as described above are implemented.

[0176] Those skilled in the art will appreciate that embodiments of the present invention may be provided as methods, systems, or computer program products. Thus, the present invention may take the form of hardware embodiments, software embodiments, or embodiments combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage and optical storage, etc.) containing computer-usable program code.

[0177] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowcharts and / or block diagrams, as well as combinations of processes and / or blocks in the flowcharts and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowcharts and / or block diagrams. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0178] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0179] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0180] Those skilled in the art will appreciate that all or part of the processes in the above-described method embodiments can be implemented by instructing related hardware through a computer program. The program can be stored in a computer-readable storage medium, and when executed, the program can include the processes in the above-described method embodiments. The storage medium can be a magnetic disk, an optical disk, a read-only memory (ROM), or a random access memory (RAM).

[0181] The foregoing description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Those skilled in the art will readily appreciate that various modifications and variations of the present invention are possible. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention are intended to be within the scope of protection of the present invention.

Claims

1. Long flexible blade structure - aerodynamic nonlinear dynamics reduced order simulation method, characterized by: include: A cross-sectional frame is established at the centroid of the blade airfoil section, and the position information of the leading and trailing edges of the blade airfoil section is extracted based on the acquired blade data. The fan blade is divided into several units, the kinetic energy of each blade unit is calculated, and the unit mass matrix of the corresponding unit is obtained based on the kinetic energy calculation of each blade unit; According to the blade airfoil section frame vector, the quasi-curvature vector related to bending and torsion on the beam element is determined, and the constant coefficient tensor of each blade element is calculated in combination with the curvature vector. The generalized elastic force of each element of the wind turbine blade is obtained by combining the constant coefficient tensor of each blade element with the generalized coordinates; The chord length vector on the airfoil section is calculated based on the position information of the leading edge and the trailing edge on the airfoil section. The relative velocity of the blade unit in three dimensions is projected onto the airfoil section to obtain the two-dimensional relative velocity. The angle of attack is calculated. The aerodynamic force acting on the blade airfoil section is calculated by combining the chord length vector, the two-dimensional relative velocity and the angle of attack. The generalized external force acting on the blade unit is obtained according to the absolute nodal coordinate method. The mass matrix, generalized elastic force matrix and generalized external force matrix of each unit are assembled to obtain the wind blade dynamic model of the slender beam model in the absolute nodal coordinate system.

2. The long flexible blade structure-aerodynamic nonlinear dynamics reduction simulation method according to claim 1, characterized in that: The constant coefficient tensor of the blade unit is calculated by combining the curvature vector, specifically based on the quasi-curvature vector related to bending and torsion on the beam unit, and the constant coefficient tensor of the blade unit is obtained according to the shape function and its partial derivative with respect to the arc length coordinate x and the stiffness parameter function of the blade unit along the arc length coordinate x of the center of mass line.

3. The long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method according to claim 1, characterized in that: The formula for calculating the curvature-like vector related to bending and torsion on a beam element is: Among them, k y With k z Curvature-like vectors related to bending and torsion, r x represents the partial derivative of the centroid function r(x,t) with respect to x, which is related to the beam arc length coordinate x. y With r z are the two direction vectors of the beam section composed of the gradient field, r xy For r x The partial derivative with respect to y, r xz For r x The partial derivative with respect to z.

4. The long flexible blade structure-aerodynamic nonlinear dynamics reduction simulation method according to claim 3, characterized in that: The calculation formula of the constant coefficient tensor of each blade element is: Among them, C, C W ,C yy ,C zz is the stiffness parameter function of the beam element along the arc length coordinate x of the centroid line, s α , s β ,y α ,y β , z α , z β represents the shape function, α, β=1,2,3...,10, and the subscript x represents the partial derivative of the shape function with respect to x, n-1 is the left node coordinate of unit n, x n is the right node coordinate of element n.

5. The long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method according to claim 1, characterized in that: The calculation formula for the generalized elastic force of each unit of the fan blade is: in, c is the viscoelastic damping coefficient, δU is the virtual work of the elastic force obtained according to the virtual work theorem, δq α is the virtual displacement, q α ,q β ,q μ ,q ν is a generalized coordinate, and the subscript range is consistent with the range of α, α=1,2,3...,10, is the first derivative of the generalized coordinate with respect to time.

6. The long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method according to claim 1, characterized in that: The angle of attack is calculated as: Among them, y te 、y le 、z te 、z le is the coordinate of the leading edge and trailing edge of the airfoil section in the section coordinate system, v rel is the relative velocity of the blade unit in three dimensions, is the two-dimensional relative velocity, c is the chord length vector on the airfoil section, v ∞ is the wind speed at infinity, is the relative velocity of the leading edge of the blade airfoil, G=(r x ,r y ,r z ) T , θ is With r x The angle, r x represents the partial derivative of the centroid function r(x,t) with respect to x, which is related to the beam arc length coordinate x. y With r z are the two direction vectors of the beam section composed of the gradient field, α=1,2,3...,10.

7. Long flexible blade structure - aerodynamic nonlinear dynamics reduced order simulation method, characterized by: include: A position information extraction module is used to establish a cross-section frame at the centroid of the blade airfoil section and extract the position information of the leading edge and the trailing edge of the blade airfoil section in combination with the acquired blade data; A mass matrix calculation module is used to divide the wind turbine blade into a number of units, calculate the kinetic energy at each blade unit, and obtain the unit mass matrix of the corresponding unit based on the kinetic energy calculation at each blade unit; A generalized elastic force calculation module is used to determine the quasi-curvature vector related to bending and torsion on the beam element based on the blade airfoil section frame vector, calculate the constant coefficient tensor of each blade element based on the curvature vector, and obtain the generalized elastic force of each element of the wind turbine blade by combining the constant coefficient tensor of each blade element with the generalized coordinates; A generalized external force calculation module is used to calculate the chord length vector on the airfoil section based on the position information of the leading edge and trailing edge on the airfoil section, project the relative velocity of the blade unit in three dimensions onto the airfoil section to obtain a two-dimensional relative velocity, calculate the angle of attack, and calculate the aerodynamic force acting on the blade airfoil section by combining the chord length vector, two-dimensional relative velocity, and angle of attack. The generalized external force acting on the blade unit is obtained using the absolute node coordinate method; The dynamic simulation model construction module is used to assemble the mass matrix, generalized elastic force matrix and generalized external force matrix of each unit to obtain the wind turbine blade dynamic model of the slender beam model in the absolute node coordinate system.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the steps of the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method according to any one of claims 1 to 6 are implemented.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the steps of the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method according to any one of claims 1 to 6 are implemented.

10. A program product, wherein the program product is a computer program product, comprising a computer program, characterized in that: When the computer program is executed by a processor, the steps of the long flexible blade structure-aerodynamic nonlinear dynamics reduced-order simulation method according to any one of claims 1 to 6 are implemented.

Images