Multi-dimensional aeroelastic calculation method and device based on nonlinear Euler-Bernoulli beam and storage medium

Abstract

The invention discloses a multi-dimensional aeroelastic calculation method and device based on a nonlinear Euler-Bernoulli beam and a storage medium, and belongs to the technical field of wind turbine aerodynamics, the method comprises the following steps: inputting the stiffness of each section of a blade into a nonlinear Euler-Bernoulli beam equation to obtain a modal-displacement transformation matrix; the rigidity, the weight and the external force of each section of the blade are converted to obtain a generalized rigidity matrix, a generalized mass matrix and a generalized force; starting circulation, substituting into a forced coupling vibration equation, solving modal displacement of each order at the current moment, and inputting into a modal-displacement inverse transformation matrix to obtain blade displacement; the generalized force is substituted into a blade element momentum theoretical formula to solve aerodynamic force at the next moment, the generalized force is updated to obtain updated generalized force, and aerodynamic force at each moment is obtained; and calculating to obtain the multi-dimensional aerobomb. According to the method, the calculation of the blade deformation and the load considers the change of the torsional degree of freedom, and more accurate and more efficient load and deformation calculation can be carried out.

Description

Technical Field

[0001] The invention relates to a multi-dimensional aeroelastic calculation method, a device and a storage medium based on a nonlinear Euler Bernoulli beam, and belongs to the technical field of wind turbine aerodynamics. Background Art

[0002] Nowadays, with the need to reduce green gas emissions from fossil energy, wind energy is developing very rapidly around the world. In order to reduce energy costs, wind turbines will become larger and larger in the near future. In the design process of wind turbines, aeroelastic response is an important characteristic, which is related to the stability and safety of wind turbines. Therefore, choosing a suitable aeroelastic model is of great significance to the accuracy and speed of the model. Nowadays, there are many ways to solve aeroelastic problems. The most common ones are the coupling of computational fluid dynamics and structural mechanics, and the coupling of blade element momentum theory and linear Euler Bernoulli beam. Traditional computational mechanics methods require modeling and calculation of blade units, resulting in high computational costs and low computational efficiency. However, the traditional blade element momentum theory and linear Euler Bernoulli beam model do not consider the torsional degree of freedom, and the aerodynamic coefficient is sensitive to the torsional degree of freedom. This shows that the traditional model has gradually become unsuitable for large wind turbine blades, and it is impossible to perform more accurate and efficient load and deformation calculations. Summary of the invention

[0003] The technical problem to be solved by the present invention is to overcome the defects of the prior art and provide a multi-dimensional aeroelastic calculation method, device and storage medium based on the nonlinear Euler Bernoulli beam. The calculation of blade deformation and load takes into account the change of torsional freedom, and can perform more accurate and efficient load and deformation calculations.

[0004] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0005] In a first aspect, the present invention provides a multi-dimensional aeroelastic calculation method based on a nonlinear Euler Bernoulli beam, comprising the following steps:

[0006] Obtaining a nonlinear Euler-Bernoulli beam equation, wherein the nonlinear Euler-Bernoulli beam equation is constructed based on a Newton-Euler method;

[0007] Based on the modal superposition method, the stiffness of each section of the blade is input into the nonlinear Euler-Bernoulli beam equation. By solving the unit flexibility matrix, the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode is obtained. Then, the corresponding modal-displacement inverse transformation matrix is ​​obtained according to the modal-displacement transformation matrix.

[0008] The stiffness, weight and external force of each blade section are converted through the forced coupled vibration equation and the mode-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force;

[0009] Start the cycle, substitute the generalized stiffness matrix, generalized mass matrix and generalized force at the current moment into the forced coupled vibration equation, transform to obtain the forced decoupled vibration equation, solve the modal displacements of each order at the current moment, and then input the modal displacements of each order into the modal-displacement inverse transformation matrix to obtain the blade displacement;

[0010] Substitute the blade displacement into the blade element momentum theory formula to solve the aerodynamic force at the next moment, and obtain the blade external force based on the aerodynamic force, blade gravity and blade centrifugal force at the next moment. Then, the blade external force is updated by the mode-displacement inverse transformation matrix to obtain the updated generalized force. Finally, the updated generalized force is used as the generalized force at the current moment to repeat the above cycle to obtain the aerodynamic force at each moment.

[0011] The multi-dimensional air bomb is obtained according to the aerodynamic force calculation at each moment, and the multi-dimensional air bomb includes the wind wheel thrust, torque and power.

[0012] The construction of the nonlinear Euler-Bernoulli beam equation includes the following steps:

[0013] According to Newton-Euler principle, we can get the force balance and torque balance equations:

[0014] (1),

[0015] in and are the external force and the external torque, and are the structural internal force and the moment around the elastic axis, x is the displacement in the x-axis direction, is the unit vector in the x-axis direction after deformation, and the following specific moment balance equation is obtained:

[0016] (2),

[0017] Among them, M x , M y and M z yes The components in the x, y, and z directions, V x , V y and V z yes The components in the x, y, and z directions, q x ,q y and q z yes The components in the x, y and z directions, v and w represent the displacement in the y and z directions;

[0018] There is also a force balancing equation:

[0019] (3),

[0020] Among them, p x , p y and p z yes For the components in the x, y, and z directions, take the partial differential of the moment balance equation with respect to the x direction and substitute the force balance equation into the equation to obtain:

[0021] (4),

[0022] Combining equation (4) with the geometric relationship of the deformed beam, and ignoring the second-order and higher-order terms in equation (4), the following equation is obtained:

[0023] (5),

[0024] Formula (5) is the moment conservation equation in the deformed coordinate system, where: is the moment about x', y', z', is the airfoil twist angle, is the torsion angle, u is the displacement in the x direction, and They are The speed and acceleration of

[0025] in, Can be made by M y and M z Transformed from:

[0026] (6),

[0027] Because the blade beam structure has a flapping curvature in the x direction and shimmy curvature Geometric relationship with v'', w'':

[0028] (7),

[0029] And according to the relationship between torque and beam curvature:

[0030] (8),

[0031] In the formula, EI1, EI2 and GI p They are the flapping stiffness, shimmy stiffness and torsional stiffness. By discarding the external force in the equation and the internal force along the x direction, the beam equilibrium equation without external force is obtained:

[0032] (9),

[0033] Formula (9) is the constructed nonlinear Euler-Bernoulli beam equation.

[0034] Based on the modal superposition method, the stiffness of each section of the blade is input into the nonlinear Euler-Bernoulli beam equation, and the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode is obtained by solving the unit flexibility matrix, including:

[0035] First, assume that the blade displacement X is described as a linear combination of modes:

[0036] (10),

[0037] In the formula, are the displacement components of each modal order, is the displacement of each mode, [U] 3×n It is a matrix with 3 rows and n columns. The number of rows represents the displacement in three directions, n represents the number of modes, and Y n It represents the displacement vectors of each modal order with n rows and 1 column;

[0038] Then, according to the nonlinear Euler-Bernoulli beam equation, we can obtain For each component in the formula, the mode-displacement transformation matrix is ​​obtained.

[0039] The solution process for each component in the formula includes:

[0040] Treat the blade as a cantilever beam and decompose it along the x direction, i.e. the blade radius direction, into The displacement of each blade segment is described as:

[0041] (11),

[0042] In the vibration displacement equation, flexibility is defined as the deformation of an elastic element under a unit force. Its physical meaning and dimension are exactly opposite to stiffness. When the external force is zero and only the internal force of the system is considered, the vibration equation is transformed into:

[0043] (12),

[0044] In formula (12), the matrix F is the flexibility matrix of the system. For the deformation of the blade, the displacement has the form of simple harmonic motion, and the acceleration is proportional to the displacement:

[0045] (13),

[0046] Substituting the relationship between acceleration and displacement, velocity and displacement into the vibration equation variant, we get:

[0047] (14),

[0048] Finally, the modal-displacement transformation matrix is ​​obtained by solving the flexibility matrix. The modal-displacement transformation matrix obtained by solving the flexibility matrix includes:

[0049] Taking the z direction as an example, a unit force in the z direction is applied to the i-th blade segment, and no force is applied to the other airfoil segments. The z-direction displacement of all blade segments under this condition is obtained by the nonlinear Euler-Bernoulli beam equation. The displacement of each small block is obtained by integrating the curvature and angular deformation in the simple beam theory along the x direction. At the same time, the external force load on the blade is integrated along the x direction to solve the torque and , this part is solved by Taylor integration;

[0050] (2) Because there is Each airfoil segment needs to be subjected to a unit force, and each time it is applied, we get displacement, so only the unit force in the z direction, under the displacement condition in the z direction, there is At the same time, there are nine situations: the blade is deformed in the y direction and around the x direction under a unit force in the z direction; the blade is deformed in the z, y and around the x direction under a unit force in the y direction; the blade is displaced in the z direction, the y direction and the x direction under a unit force around the x direction. These data are stored in a Next, these data need to be stored in the flexibility matrix in the following order:

[0051] (15)

[0052] Among them, fzz in the block matrix represents the effect of the z-direction force in the z-direction, frxrx is the effect of the moment around x in the x-direction, and the lower right corner of each block matrix represents: F1 is the unit force on the first block, FN is the unit force on the Nth block, d1 refers to the deformation on the first block, d2 refers to the deformation on the second block, and dN refers to the deformation on the Nth block;

[0053] Solve the eigenvalue and eigenvector of the determinant of the flexibility matrix and obtain The characteristic frequencies of each mode of a blade segment are:

[0054] (16)

[0055] in, is the eigenvalue, fre is the eigenfrequency, and the solved eigenvector matrix is ​​a The matrix is ​​stored in the following order:

[0056] (17)

[0057] In formula (17), the block matrices BB1, BB2 and BB 3n Represent the displacements of the first, second and 3nth modes respectively, fre1 and fre3N refer to the first and 3Nth modes respectively. The above matrices are normalized to obtain the unit mode component matrix, which is the mode-displacement transformation matrix:

[0058] (18)

[0059] In formula (18), L1, L2 and L 3N are the sum of the maximum modal components under the corresponding modes, That is The mode-displacement transformation matrix is , , is the deformation at a certain characteristic frequency.

[0060] The stiffness, weight and external force of each section of the blade are converted through the forced coupled vibration equation and the modal-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force including:

[0061] According to the blade forced coupling vibration equation, it satisfies:

[0062] (19)

[0063] In the formula, [M] is the mass matrix of each blade section, [C] is the damping matrix of each blade section, [K] is the stiffness matrix of each blade section, {F} is the external force of each blade section, and X changes to , and multiply the left side of the equation by , we get the following relationship:

[0064] (20)

[0065] The above formula is the characteristic mode decoupled vibration equation, and it can be obtained by solving it:

[0066] (twenty one)

[0067] in, is the generalized mass matrix, is the generalized damping matrix, is the generalized stiffness matrix, and the matrix size is , are all diagonal matrices;

[0068] The generalized force matrix is:

[0069] (twenty two).

[0070] The blade displacement is substituted into the blade element momentum theory formula to solve the aerodynamic force at the next moment, and the blade external force is obtained based on the aerodynamic force at the next moment, the blade gravity and the blade centrifugal force. The blade external force is then updated by the mode-displacement inverse transformation matrix to obtain the updated generalized force. Finally, the updated generalized force is used as the generalized force at the current moment to repeat the above cycle, and the aerodynamic force at each moment satisfies the following relationship:

[0071] (twenty three)

[0072] In formula (23) is the relative speed, is the incoming wind speed, is the wake loss wind speed, is the linear velocity of rotation, is the blade structural speed;

[0073] (twenty four)

[0074] In formula (24), is the relative velocity in the z direction, is the relative velocity in the y direction,

[0075] is the inflow angle, which is the angle between the incoming flow and the rotating plane. is the angle of attack, which is the angle between the incoming flow and the blade chord. is the blade twist angle, It is the sum of the blade twist angle and the pitch angle. After obtaining the angle of attack, look up the lift coefficient in the lift and drag coefficient table to obtain the lift coefficient. , drag coefficient and moment coefficient , AC is the aerodynamic center, so the aerodynamic torque is obtained , lift L, drag D, normal force and tangential force :

[0076] (25)

[0077] (26)

[0078] In the formula, is the air density, c is the chord length of each blade segment, the elastic center is closer to the center of the blade, and the moment around the elastic axis and its moment coefficient Written as:

[0079] (27)

[0080] In the formula, is the rotational moment around the elastic axis, e is the distance from the elastic axis to the aerodynamic center, and the conversion from the aerodynamic center torque to the elastic center torque can be obtained by the following formula: the blade lift is positive upward, its center is 1 / 4 chord length from the leading edge, and e is the distance from the aerodynamic center to the elastic axis.

[0081] In the second aspect, the present invention provides a multi-dimensional aeroelastic computing device for a nonlinear Euler Bernoulli beam, comprising:

[0082] A nonlinear Euler-Bernoulli beam equation acquisition module, used for acquiring a nonlinear Euler-Bernoulli beam equation, wherein the nonlinear Euler-Bernoulli beam equation is constructed based on the Newton-Euler method;

[0083] The modal-displacement transformation matrix solving module is used to input the stiffness of each section of the blade into the nonlinear Euler-Bernoulli beam equation based on the modal superposition method, and obtain the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode by solving the unit flexibility matrix, and then obtain the corresponding modal-displacement inverse transformation matrix according to the modal-displacement transformation matrix;

[0084] A conversion module is used to convert the stiffness, weight and external force of each section of the blade through the forced coupled vibration equation and the mode-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force;

[0085] The blade displacement calculation module is used to start the cycle, substitute the generalized stiffness matrix, the generalized mass matrix and the generalized force at the current moment into the forced coupled vibration equation, transform to obtain the forced decoupled vibration equation, solve the modal displacements of each order at the current moment, and then input the modal displacements of each order into the modal-displacement inverse transformation matrix to obtain the blade displacement;

[0086] The generalized force cycle module is used to substitute the blade displacement into the blade element momentum theory formula to solve the aerodynamic force at the next moment, obtain the blade external force based on the aerodynamic force at the next moment, blade gravity and blade centrifugal force, and then update the generalized force of the blade external force through the mode-displacement inverse transformation matrix to obtain the updated generalized force. Finally, the updated generalized force is used as the generalized force at the current moment to repeat the above cycle to obtain the aerodynamic force at each moment;

[0087] The multi-dimensional aeroelastic calculation module is used to obtain the multi-dimensional aeroelastic according to the aerodynamic force calculation at each moment. The multi-dimensional aeroelastic includes the wind wheel thrust, torque and power.

[0088] A third aspect: The present invention provides a computer-readable storage medium having a computer program / instruction stored thereon, which, when executed by a processor, implements the steps of the multi-dimensional aeroelastic calculation method based on a nonlinear Euler Bernoulli beam.

[0089] Beneficial effects of the present invention: The present invention provides a multi-dimensional aeroelastic calculation method, device and storage medium based on a nonlinear Euler Bernoulli beam. Based on the nonlinear Euler Bernoulli beam model, the relationship between blade displacement and characteristic modal displacement is obtained with the modal superposition method, and the blade vibration equation is solved to iteratively obtain a more accurate displacement and load size under aeroelastic coupling, so that the aerodynamic load calculation accuracy of large wind turbines is higher, which is more in line with the needs of large-scale development of wind turbines. At the same time, compared with the traditional method of coupling computational fluid dynamics and computational structural mechanics, the present invention provides similar accuracy and greatly reduces the simulation time, which is suitable for wind turbine design. BRIEF DESCRIPTION OF THE DRAWINGS

[0090] Figure 1 It is a flow chart of the multi-dimensional aeroelastic calculation method based on the nonlinear Euler Bernoulli beam of the present invention;

[0091] Figure 2 It is a coordinate geometric relationship diagram of the blade before and after deformation in the present invention;

[0092] Figure 3 is the incoming flow angular attack velocity triangle of a certain airfoil section in the present invention;

[0093] Figure 4 It is a comparison diagram of blade deformation under NREL5MW, DTU10MW, and IEA15MW rated conditions and different beam model and aerodynamic model combination methods in the present invention;

[0094] Figure 5 It is a blade thrust comparison diagram of different beam model and aerodynamic model combination methods under NREL5MW, DTU10MW, and IEA15MW rated conditions in the present invention. DETAILED DESCRIPTION

[0095] The present invention will be further described below in conjunction with the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and are not intended to limit the protection scope of the present invention.

[0096] Example 1

[0097] like Figure 1 As shown, the present invention provides a flowchart of a multi-dimensional aeroelastic calculation method based on a nonlinear Euler Bernoulli beam, comprising the following steps:

[0098] Step 1: Obtain a nonlinear Euler-Bernoulli beam equation, wherein the nonlinear Euler-Bernoulli beam equation is constructed based on the Newton-Euler method. The construction of the nonlinear Euler-Bernoulli beam equation includes the following steps:

[0099] According to Newton-Euler principle, we can get the force balance and torque balance equations:

[0100] (1),

[0101] in and are the external force and the external torque, and are the structural internal force and the moment about the elastic axis, Non-deformed The unit vector in the coordinate system, After deformation The unit vector in the coordinate system corresponds specifically to Figure 2 The coordinate system diagram of , thus obtaining the following specific moment equilibrium equation:

[0102] (2),

[0103] Among them, M x , M y and M z yes The components in the x, y, and z directions, V x , V y and V z yes The components in the x, y, and z directions, q x ,q y and q z yes The components in the x, y and z directions, v and w represent the displacement in the y and z directions;

[0104] There is also a force balancing equation:

[0105] (3),

[0106] Among them, p x , p y and p z yes For the components in the x, y, and z directions, take the partial differential of the moment balance equation with respect to the x direction and substitute the force balance equation into the equation to obtain:

[0107] (4),

[0108] Combining equation (4) with the geometric relationship of the deformed beam, and ignoring the second-order and higher-order terms in equation (4), the following equation is obtained:

[0109] (5),

[0110] Formula (5) is the moment conservation equation in the deformed coordinate system, where: is the moment about x', y', z', is the airfoil twist angle, is the torsion angle, u is the displacement in the x direction, and They are The speed and acceleration of

[0111] in, Can be made by M y and M z Transformed from:

[0112] (6),

[0113] Because the blade beam structure has a flapping curvature in the x direction and shimmy curvature Geometric relationship with v'', w'':

[0114] (7),

[0115] And according to the relationship between torque and beam curvature:

[0116] (8),

[0117] In the formula, EI1, EI2 and GI p They are the flapping stiffness, shimmy stiffness and torsional stiffness. By discarding the external force in the equation and the internal force along the x direction, the beam equilibrium equation without external force is obtained:

[0118] (9),

[0119] Formula (9) is the constructed nonlinear Euler-Bernoulli beam equation.

[0120] Step 2: Based on the modal superposition method, the stiffness of each section of the blade is input into the nonlinear Euler-Bernoulli beam equation. By solving the unit flexibility matrix, the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode is obtained, and then the corresponding modal-displacement inverse transformation matrix is ​​obtained according to the modal-displacement transformation matrix.

[0121] First, assume that the blade displacement X is described as a linear combination of modes:

[0122] (10),

[0123] In the formula, are the displacement components of each modal order, is the displacement of each mode, [U] 3×n It is a matrix with 3 rows and n columns. The number of rows represents the displacement in three directions, n represents the number of modes, and Y n It represents the displacement vectors of each modal order with n rows and 1 column;

[0124] Then, according to the nonlinear Euler-Bernoulli beam equation, we can obtain For each component in the formula, the mode-displacement transformation matrix is ​​obtained.

[0125] The solution process for each component in the formula includes:

[0126] Treat the blade as a cantilever beam and decompose it along the x direction, i.e. the blade radius direction, into The displacement of each blade segment is described as:

[0127] (11),

[0128] In the vibration displacement equation, flexibility is defined as the deformation of an elastic element under a unit force. Its physical meaning and dimension are exactly opposite to stiffness. When the external force is zero and only the internal force of the system is considered, the vibration equation is transformed into:

[0129] (12),

[0130] In formula (12), the matrix F is the flexibility matrix of the system. For the deformation of the blade, the displacement has the form of simple harmonic motion, and the acceleration is proportional to the displacement:

[0131] (13),

[0132] Substituting the relationship between acceleration and displacement, velocity and displacement into the vibration equation variant, we get:

[0133] (14),

[0134] Finally, the mode-displacement transformation matrix is ​​obtained by solving the flexibility matrix.

[0135] Among them, solving the flexibility matrix to obtain the mode-displacement transformation matrix includes:

[0136] Taking the z direction as an example, a unit force in the z direction is applied to the i-th blade segment, and no force is applied to the other airfoil segments. The z-direction displacement of all blade segments under this condition is obtained by the nonlinear Euler-Bernoulli beam equation. The displacement of each small block is obtained by integrating the curvature and angular deformation in the simple beam theory along the x direction. At the same time, the external force load on the blade is integrated along the x direction to solve the torque and , which is solved by Taylor integration.

[0137] Because there is Each airfoil segment needs to be subjected to a unit force, and each time it is applied, we get displacement, so only the unit force in the z direction, under the displacement condition in the z direction, there is The displacement produced under the displacement condition in the z direction The displacement data, as well as the deformation of the blade in the y direction and around the x direction under the unit force in the z direction; the deformation in the z, y and around the x direction under the unit force in the y direction; the displacement in the z direction, the displacement in the y direction and the deformation around the x direction under the unit force around the x direction are stored in a In the three-dimensional flexibility matrix, the storage order is:

[0138] (15)

[0139] Among them, fzz in the block matrix represents the effect of the z-direction force in the z-direction, frxrx is the effect of the moment around x in the x-direction, and the lower right corner of each block matrix represents: F1 is the unit force on the first block, FN is the unit force on the Nth block, d1 refers to the deformation on the first block, d2 refers to the deformation on the second block, and dN refers to the deformation on the Nth block;

[0140] Solve the eigenvalue and eigenvector of the determinant of the flexibility matrix and obtain The characteristic frequencies of each mode of a blade segment are:

[0141] (16)

[0142] in, is the eigenvalue, fre is the eigenfrequency, and the solved eigenvector matrix is ​​a The matrix is ​​stored in the following order:

[0143] (17)

[0144] In formula (17), the block matrices BB1, BB2 and BB 3n Represent the displacements of the first, second and 3nth modes respectively, fre1 and fre3N refer to the first and 3Nth modes respectively. The above matrices are normalized to obtain the unit mode component matrix, which is the mode-displacement transformation matrix:

[0145] (18)

[0146] In formula (18), L1, L2 and L 3N are the sum of the maximum modal components under the corresponding modes, That is The mode-displacement transformation matrix is , , is the deformation at a certain characteristic frequency.

[0147] Step three, transform the stiffness, weight and external force of each section of the blade through the forced coupled vibration equation and the modal-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force.

[0148] According to the blade forced coupling vibration equation, it satisfies:

[0149] (19)

[0150] In the formula, [M] is the mass matrix of each blade section, [C] is the damping matrix of each blade section, [K] is the stiffness matrix of each blade section, {F} is the external force of each blade section, and X changes to , and multiply the left side of the equation by , we get the following relationship:

[0151] (20)

[0152] The above formula is the characteristic mode decoupled vibration equation, and it can be obtained by solving it:

[0153] (twenty one)

[0154] in, is the generalized mass matrix, is the generalized damping matrix, is the generalized stiffness matrix, and the matrix size is , are all diagonal matrices. From a physical point of view, it can be considered that the modes of different orders are not coupled with each other, and each equation is independent.

[0155] The generalized force matrix is:

[0156] (twenty two).

[0157] By solving the above equations with the fourth-order Runge-Kutta method, we can obtain the displacements of each modal order. Usually, the first few modes are selected for solution. After obtaining the modal displacement, the blade displacement is obtained through the inverse transformation matrix of the modal-displacement.

[0158] Step 4, start the loop, substitute the generalized stiffness matrix, generalized mass matrix and generalized force at the current moment into the forced coupled vibration equation, transform to obtain the forced decoupled vibration equation, solve the modal displacements of each order at the current moment, and then input the modal displacements of each order into the modal-displacement inverse transformation matrix to obtain the blade displacement.

[0159] Step five, substitute the blade displacement into the blade element momentum theory formula to solve for the aerodynamic force at the next moment, obtain the blade external force based on the aerodynamic force at the next moment, blade gravity and blade centrifugal force, and then update the generalized force of the blade external force through the mode-displacement inverse transformation matrix to obtain the updated generalized force, and finally repeat the above cycle using the updated generalized force as the generalized force at the current moment to obtain the aerodynamic force at each moment.

[0160] Satisfies the following relationship:

[0161] (twenty three)

[0162] In formula (23) is the relative speed, is the incoming wind speed, is the wake loss wind speed, is the linear velocity of rotation, is the blade structural speed;

[0163] (twenty four)

[0164] In formula (24), is the relative velocity in the z direction, is the relative velocity in the y direction,

[0165] is the inflow angle, which is the angle between the incoming flow and the rotating plane. is the angle of attack, which is the angle between the incoming flow and the blade chord. is the blade twist angle, It is the sum of the blade twist angle and the pitch angle. After obtaining the angle of attack, look up the lift coefficient in the lift and drag coefficient table to obtain the lift coefficient. , drag coefficient and moment coefficient , AC is the aerodynamic center, so the aerodynamic torque is obtained , lift L, drag D, normal force and tangential force :

[0166] (25)

[0167] (26)

[0168] In the formula, is the air density, c is the chord length of each blade segment, the elastic center is closer to the center of the blade, and the moment around the elastic axis and its moment coefficient Written as:

[0169] (27)

[0170] In the formula, is the rotational moment around the elastic axis, and e is the distance from the elastic axis to the aerodynamic center, such as Figure 3 The conversion from aerodynamic center torque to elastic center torque can be obtained by the following formula, where the blade lift is positive upward and its center is 1 / 4 chord length from the leading edge.

[0171] Step six, obtain the multi-dimensional air bomb according to the aerodynamic force calculation at each moment, and the multi-dimensional air bomb includes the wind wheel thrust, torque and power.

[0172] Example 2

[0173] Based on Example 1, this example provides blade displacement diagrams of NREL5MW, DTU10MW, and IEA15MW wind turbines along the span direction under rated conditions using the method of the present invention. Figure 4 Figure a is a comparison of blade deformation obtained by NREL5MW using the existing model and the computational fluid dynamics plus finite element method. NEBB refers to the nonlinear Euler Bernoulli beam model. It can be seen from the figure that the method proposed in this paper can obtain similar results to the computational fluid dynamics method and save a lot of time. Figure 4 In the figure b, DTU10MW uses the existing model, which is similar to the result of CFD+Shell+ROM. Shell is the shell element method and ROM is the reduced-order model method. These are the calculation methods with the highest assurance at present. TB is the Timoshenko beam model. This model can obtain similar results with both high-fidelity methods and save time. Figure 4 In the figure, c is the IEA15MW model, H2-PTNT and H2-FPM are two beam models of HAWC2 software, the former does not consider torsion, while the latter considers torsion, ElastoDyn and Beamdyn are two beam models of OpenFAST software, the former does not consider torsion, while the latter does, GEBT is a geometrically accurate beam model. As can be seen from the figure, the result obtained by considering the torsional degree of freedom shows that the flapping displacement of the blade is smaller than that of the traditional model, and the conclusion of the traditional model has the problem of exaggerating the flapping displacement of the blade.

[0174] Figure 5 Figures a, b and c are thrust comparison diagrams of NREL5MW, DTU10MW and IEA15MW with the change of incoming wind speed. Compared with the results without considering the torsional degree of freedom, there is a conclusion that the maximum thrust is reduced, and there is a good match in the method of coupling computational fluid dynamics and computational structural mechanics, which verifies the effectiveness of the aeroelastic calculation method of wind turbine blades provided by the present invention.

[0175] Example 3

[0176] This embodiment discloses a device for judging failure of contact between the large end face of a hub bearing roller and a rib, comprising:

[0177] A nonlinear Euler-Bernoulli beam equation acquisition module, used for acquiring a nonlinear Euler-Bernoulli beam equation, wherein the nonlinear Euler-Bernoulli beam equation is constructed based on the Newton-Euler method;

[0178] The modal-displacement transformation matrix solving module is used to input the stiffness of each section of the blade into the nonlinear Euler-Bernoulli beam equation based on the modal superposition method, and obtain the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode by solving the unit flexibility matrix, and then obtain the corresponding modal-displacement inverse transformation matrix according to the modal-displacement transformation matrix;

[0179] A conversion module is used to convert the stiffness, weight and external force of each section of the blade through the forced coupled vibration equation and the mode-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force;

[0180] The blade displacement calculation module is used to start the cycle, substitute the generalized stiffness matrix, the generalized mass matrix and the generalized force at the current moment into the forced coupled vibration equation, transform to obtain the forced decoupled vibration equation, solve the modal displacements of each order at the current moment, and then input the modal displacements of each order into the modal-displacement inverse transformation matrix to obtain the blade displacement;

[0181] The generalized force cycle module is used to substitute the blade displacement into the blade element momentum theory formula to solve the aerodynamic force at the next moment, obtain the blade external force based on the aerodynamic force at the next moment, blade gravity and blade centrifugal force, and then update the generalized force of the blade external force through the mode-displacement inverse transformation matrix to obtain the updated generalized force. Finally, the updated generalized force is used as the generalized force at the current moment to repeat the above cycle to obtain the aerodynamic force at each moment;

[0182] The multi-dimensional aeroelastic calculation module is used to obtain the multi-dimensional aeroelastic according to the aerodynamic force calculation at each moment. The multi-dimensional aeroelastic includes the wind wheel thrust, torque and power.

[0183] Example 4

[0184] This embodiment discloses a computer-readable storage medium on which a computer program is stored. When the program is executed by a processor, the steps of the method described in Embodiment 1 or Embodiment 2 are executed.

[0185] It will be appreciated by those skilled in the art that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment 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 disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0186] 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 flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, 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 flowchart and / or block diagram. 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.

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

[0188] 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 instructions for implementing the process in the computer or other programmable device. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0189] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A multi-dimensional aeroelastic calculation method based on nonlinear Euler Bernoulli beam, characterized by: The following steps are involved: Obtaining a nonlinear Euler-Bernoulli beam equation, wherein the nonlinear Euler-Bernoulli beam equation is constructed based on a Newton-Euler method; Based on the modal superposition method, the stiffness of each section of the blade is input into the nonlinear Euler-Bernoulli beam equation. By solving the unit flexibility matrix, the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode is obtained. Then, the corresponding modal-displacement inverse transformation matrix is ​​obtained according to the modal-displacement transformation matrix. The stiffness, weight and external force of each blade section are converted through the forced coupled vibration equation and the mode-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force; Start the cycle, substitute the generalized stiffness matrix, generalized mass matrix and generalized force at the current moment into the forced coupled vibration equation, transform to obtain the forced decoupled vibration equation, solve the modal displacements of each order at the current moment, and then input the modal displacements of each order into the modal-displacement inverse transformation matrix to obtain the blade displacement; Substitute the blade displacement into the blade element momentum theory formula to solve the aerodynamic force at the next moment, and obtain the blade external force based on the aerodynamic force, blade gravity and blade centrifugal force at the next moment. Then, the blade external force is updated by the mode-displacement inverse transformation matrix to obtain the updated generalized force. Finally, the updated generalized force is used as the generalized force at the current moment to repeat the above cycle to obtain the aerodynamic force at each moment. The multi-dimensional air bomb is obtained according to the aerodynamic force calculation at each moment, and the multi-dimensional air bomb includes the wind wheel thrust, torque and power.

2. The multi-dimensional aeroelastic calculation method based on the nonlinear Euler Bernoulli beam according to claim 1, characterized in that: The construction of the nonlinear Euler-Bernoulli beam equation includes the following steps: According to Newton-Euler principle, we can get the force balance and torque balance equations: （1）， in and are the external force and the external torque, and are the structural internal force and the moment around the elastic axis, x is the displacement in the x-axis direction, is the unit vector in the x-axis direction after deformation, and the following specific moment balance equation is obtained: （2）， Among them, M x , M y and M z yes The components in the x, y, and z directions, V x , V y and V z yes The components in the x, y, and z directions, q x ,q y and q z yes The components in the x, y and z directions, v and w represent the displacement in the y and z directions; There is also a force balancing equation: （3）， Among them, p x , p y and p z yes For the components in the x, y, and z directions, take the partial differential of the moment balance equation with respect to the x direction and substitute the force balance equation into the equation to obtain: （4）， Combining equation (4) with the geometric relationship of the deformed beam, and ignoring the second-order and higher-order terms in equation (4), the following equation is obtained: （5）， Formula (5) is the moment conservation equation in the deformed coordinate system, where: is the moment about x', y', z', is the airfoil twist angle, is the torsion angle, u is the displacement in the x direction, and They are The speed and acceleration of in, Can be made by M y and M z Transformed from: （6）， Because the blade beam structure has a flapping curvature in the x direction and shimmy curvature Geometric relationship with v'', w'': （7）， And according to the relationship between torque and beam curvature: （8）， In the formula, EI1, EI2 and GI p They are the flapping stiffness, shimmy stiffness and torsional stiffness. By discarding the external force in the equation and the internal force along the x direction, the beam equilibrium equation without external force is obtained: （9）， Formula (9) is the constructed nonlinear Euler-Bernoulli beam equation.

3. The multi-dimensional aeroelastic calculation method based on nonlinear Euler Bernoulli beam according to claim 2, characterized in that: Based on the modal superposition method, the stiffness of each section of the blade is input into the nonlinear Euler-Bernoulli beam equation, and the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode is obtained by solving the unit flexibility matrix, including: First, assume that the blade displacement X is described as a linear combination of modes: （10）， In the formula, is the first-order modal displacement component, is the second-order modal displacement component, is the third-order modal displacement component, is the displacement of each mode, [U] 3×n It is a matrix with 3 rows and n columns. The number of rows represents the displacement in three directions, n represents the number of modes, and Y n It represents the displacement vectors of each modal order with n rows and 1 column; Then, according to the nonlinear Euler-Bernoulli beam equation, we can obtain For each component in the formula, the mode-displacement transformation matrix is ​​obtained.

4. The multi-dimensional aeroelastic calculation method based on nonlinear Euler Bernoulli beam according to claim 3, characterized in that: The solution process for each component in the formula includes: Treat the blade as a cantilever beam and decompose it along the x direction, i.e. the blade radius direction, into The displacement of each blade segment is described as: （11）， In the vibration displacement equation, flexibility is defined as the deformation of an elastic element under a unit force. Its physical meaning and dimension are exactly opposite to stiffness. When the external force is zero and only the internal force of the system is considered, the vibration equation is transformed into: （12）， In formula (12), M is the mass matrix, and matrix F is the flexibility matrix of the system. For the deformation of the blade, the displacement has the form of simple harmonic motion, and the acceleration is proportional to the displacement: （13）， Substituting the relationship between acceleration and displacement, velocity and displacement into the vibration equation variant, we get: （14）， Finally, the mode-displacement transformation matrix is ​​obtained by solving the flexibility matrix.

5. The multi-dimensional aeroelastic calculation method based on nonlinear Euler Bernoulli beam according to claim 4, characterized in that: Solving the flexibility matrix to obtain the mode-displacement transformation matrix includes: A unit force under a certain degree of freedom is applied to the i-th blade segment, and the other airfoil segments remain in a state without external force applied. The external force load on the blade is integrated along the x direction to solve the torque under all degrees of freedom; The torque is introduced into the nonlinear Euler Bernoulli beam equation, and combined with the relationship between curvature and angular deformation in simple beam theory, the displacement of all blade segments in all degrees of freedom under this condition is determined. displacement data, Ndof is the number of degrees of freedom, and all degrees of freedom are generated The displacement data, Ndof takes the value of 3, is stored in a In the three-dimensional flexibility matrix, the storage order is: （15）， Among them, fzz in the block matrix represents the effect of the z-direction force in the z-direction, frxrx is the effect of the moment around x in the x-direction, and the lower right corner of each block matrix represents: F1 is the unit force on the first block, FN is the unit force on the Nth block, d1 refers to the deformation on the first block, d2 refers to the deformation on the second block, and dN refers to the deformation on the Nth block; Solve the eigenvalue and eigenvector of the determinant of the flexibility matrix and obtain The characteristic frequencies of each mode of a blade segment are: （16）， in, is the eigenvalue, fre is the eigenfrequency, and the solved eigenvector matrix is ​​a The matrix is ​​stored in the following order: （17）， In formula (17), the block matrices BB1, BB2 and BB 3n Represent the displacements of the first, second and 3nth modes respectively, fre1 and fre3N refer to the first and 3Nth modes respectively. The above matrices are normalized to obtain the unit mode component matrix, which is the mode-displacement transformation matrix: （18）， In formula (18), L1, L2 and L 3N are the sum of the maximum modal components under the corresponding modes, That is The mode-displacement transformation matrix is , , is the deformation at a certain characteristic frequency.

6. The multi-dimensional aeroelastic calculation method based on nonlinear Euler Bernoulli beam according to claim 5, characterized in that: The stiffness, weight and external force of each section of the blade are converted through the forced coupling vibration equation and the modal-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force, including: According to the blade forced coupling vibration equation, it satisfies: （19）， In the formula, [M] is the mass matrix of each blade section, [C] is the damping matrix of each blade section, [K] is the stiffness matrix of each blade section, {F} is the external force of each blade section, and X changes to , and multiply the left side of the equation by , we get the following relationship: （20）， The above formula is the characteristic mode decoupled vibration equation, and it can be obtained by solving it: （21）， in, is the generalized mass matrix, is the generalized damping matrix, is the generalized stiffness matrix, and the matrix size is , are all diagonal matrices; The generalized force matrix is: （22）。 7. The multi-dimensional aeroelastic calculation method based on nonlinear Euler Bernoulli beam according to claim 6, characterized in that: Substituting the blade displacement into the blade element momentum theory formula to solve the aerodynamic force at the next moment, the blade external force is obtained based on the aerodynamic force at the next moment, the blade gravity and the blade centrifugal force, and then the blade external force is updated by the mode-displacement inverse transformation matrix to obtain the updated generalized force, and finally the updated generalized force is used as the generalized force at the current moment to repeat the above cycle, and the aerodynamic force at each moment satisfies the following relationship: （23）， In formula (23) is the relative speed, is the incoming wind speed, is the wake loss wind speed, is the linear velocity of rotation, is the blade structural speed; （24）， In formula (24), is the relative velocity in the z direction, is the relative velocity in the y direction, is the inflow angle, which is the angle between the incoming flow and the rotating plane. is the angle of attack, which is the angle between the incoming flow and the blade chord. is the blade twist angle, It is the sum of the blade twist angle and the pitch angle. After obtaining the angle of attack, look up the lift coefficient in the lift and drag coefficient table to obtain the lift coefficient. , drag coefficient and moment coefficient , AC is the aerodynamic center, so the aerodynamic torque is obtained , lift L, drag D, normal force and tangential force : （25）， （26）， In the formula, is the air density, c is the chord length of each blade segment, the elastic center is closer to the center of the blade, and the moment around the elastic axis and its moment coefficient Written as: （27）， In the formula, is the rotational moment around the elastic axis, e is the distance from the elastic axis to the aerodynamic center, and the conversion from the aerodynamic center torque to the elastic center torque can be obtained by the following formula: the blade lift is positive upward, its center is 1 / 4 chord length from the leading edge, and e is the distance from the aerodynamic center to the elastic axis.

8. A multi-dimensional aeroelastic computing device for nonlinear Euler Bernoulli beams, characterized by: include; A nonlinear Euler-Bernoulli beam equation acquisition module, used for acquiring a nonlinear Euler-Bernoulli beam equation, wherein the nonlinear Euler-Bernoulli beam equation is constructed based on the Newton-Euler method; The modal-displacement transformation matrix solving module is used to input the stiffness of each section of the blade into the nonlinear Euler-Bernoulli beam equation based on the modal superposition method, and obtain the modal-displacement transformation matrix between the blade flapping displacement, swing displacement and torsional displacement and each order mode by solving the unit flexibility matrix, and then obtain the corresponding modal-displacement inverse transformation matrix according to the modal-displacement transformation matrix; A conversion module is used to convert the stiffness, weight and external force of each section of the blade through the forced coupled vibration equation and the mode-displacement transformation matrix to obtain the generalized stiffness matrix, generalized mass matrix and generalized force; The blade displacement calculation module is used to start the cycle, substitute the generalized stiffness matrix, the generalized mass matrix and the generalized force at the current moment into the forced coupled vibration equation, transform to obtain the forced decoupled vibration equation, solve the modal displacements of each order at the current moment, and then input the modal displacements of each order into the modal-displacement inverse transformation matrix to obtain the blade displacement; The generalized force cycle module is used to substitute the blade displacement into the blade element momentum theory formula to solve the aerodynamic force at the next moment, obtain the blade external force based on the aerodynamic force at the next moment, blade gravity and blade centrifugal force, and then update the generalized force of the blade external force through the mode-displacement inverse transformation matrix to obtain the updated generalized force. Finally, the updated generalized force is used as the generalized force at the current moment to repeat the above cycle to obtain the aerodynamic force at each moment; The multi-dimensional aeroelastic calculation module is used to obtain the multi-dimensional aeroelastic according to the aerodynamic force calculation at each moment. The multi-dimensional aeroelastic includes the wind wheel thrust, torque and power.

9. A computer-readable storage medium having a computer program / instruction stored thereon, characterized in that: When the computer program / instruction is executed by a processor, the steps of the multi-dimensional aeroelastic calculation method based on the nonlinear Euler Bernoulli beam described in any one of claims 1 to 7 are implemented.

Images