Geometric precise beam dynamics numerical modeling and calculating method

Abstract

The invention belongs to the technical field of beam structure systems, and particularly relates to a geometric precise beam dynamics numerical modeling and calculating method. According to the method, the elastic motion of a beam is described through geometric nonlinearity, a Gaussian integral and trapezoidal integral method is adopted to calculate a high-order term and a derivative thereof, and an elastic motion equation and a rigid body motion equation of the beam are established. In specific implementation, variables are discretized, a shape function and a derivative thereof are calculated, and then an elastic variable and a derivative thereof are solved. The method is suitable for medium and large deformation beams with deformation angles smaller than 90 degrees, and has good numerical stability and calculation efficiency. By means of the method, the coupling effect of elastic motion and rigid body motion of the blade and the tensile torsion coupling effect of elastic deformation can be accurately described, and powerful support is provided for structural design of helicopter blades, wind driven generator blades and the like.

Description

Technical Field

[0001] The invention belongs to the technical field of beam structure systems, and in particular relates to a geometrically accurate beam dynamics numerical modeling and calculation method. Background Art

[0002] Due to the increasing application of advanced configuration rotors, hingeless rotors and bearingless rotors are simpler in structure and have lower maintenance costs than articulated rotors. Hingeless rotors and bearingless rotors replace the hinges on traditional rotors with flexible components at the root of the blades. The flapping and shimmying motions of the blades are achieved through the elastic deformation of the root of the blades, which complicates the rotor dynamics problem. In order to reduce R&D costs and shorten the design cycle, during the rotor design stage, the dynamic characteristics and stability analysis level of the rotor depends on the theoretical prediction ability of the rotor aeroelasticity. A high-precision aeroelastic dynamics model of the rotor structure is the key to better developing advanced configurations, and a high-precision blade model is the key to the aeroelastic dynamics model of the rotor structure.

[0003] Conventional rigid blade models and linear elastic blade models are difficult to describe complex nonlinear problems in rotor dynamics, such as high-order terms in the dynamic equations, which cannot accurately describe the coupling effect between the elastic motion of the blade and the rigid body motion, as well as the tensile-torsion coupling effect of elastic deformation. With the widespread use of blades made of composite materials, the large deformation and various characteristics of the blades are more obvious. The insufficient accuracy of the blade model leads to the accuracy of the aeroelastic analysis of the rotor structure not meeting the requirements. In order to describe the nonlinear dynamic characteristics of blades with large deformation, a more accurate beam model is needed for numerical simulation calculations. Summary of the invention

[0004] The purpose of the present invention is to provide a geometrically accurate beam dynamics numerical modeling and calculation method to solve the nonlinear dynamics problem of large deformation beams in structures such as advanced configuration rotors proposed in the background technology.

[0005] The present invention achieves the above-mentioned purpose through the following technical solutions:

[0006] The present invention proposes a geometrically accurate beam dynamics numerical modeling and calculation method, the method comprising:

[0007] S1, obtain the shape function h(ξ) vector at the current Gaussian integration point and its related derivatives with respect to the beam axis coordinate x;

[0008] S2, calculating the elastic variables at the current Gaussian integration point of the beam, and its first-order derivative and second-order derivative with respect to time or beam axis coordinate x;

[0009] S3, using the elastic variable and its derivative at the current Gaussian integration point, calculating the torsional deformation high-order term and axial deformation high-order term generated by bending, and the time derivative of the high-order term through trapezoidal integration, and discretizing the coefficient matrix after the generalized coordinates;

[0010] S4, according to the Green-Lagrange strain tensor, using the calculation results of steps S1, S2, and S3 to calculate the strain measure of the cross section at the corresponding Gaussian integral point, obtain the coefficient matrix after discretizing the generalized coordinates, and discretize the equivalent stiffness matrix after the generalized coordinates;

[0011] S5, using the Gaussian integral method, integrating the calculation result of step S4 along the axial direction of the beam, calculating the strain energy term obtained by discretizing the variational generalized coordinates after the variation, and the equivalent stiffness matrix after discretizing the generalized coordinates;

[0012] S6. According to the D'Alembert principle, the acceleration and virtual displacement at the corresponding Gaussian integral point are calculated using the calculation results of steps S1, S2, and S3, and an equivalent mass matrix is ​​discretized;

[0013] S7, using the Gaussian integration method, integrating the product of the virtual displacement, acceleration and gravity obtained in step S6 along the beam axis, calculating the variation, and discretizing the kinetic energy term and equivalent mass matrix obtained by the variational generalized coordinates;

[0014] S8, calculating the work done by all external loads acting on the elastic axis according to the calculation results of steps S1, S2, and S3, and discretizing the external load terms of variational generalized coordinates after variation;

[0015] S9, according to Hamilton's principle, combined with the calculation results of steps S5, S7, and S8, establish the elastic motion equation and rigid body motion equation of the beam;

[0016] S10. Newton's iteration method is used, with generalized coordinates as the degrees of freedom. The generalized velocity and generalized acceleration are calculated using the Newmark-β method. The motion equations in step S9 are solved iteratively until the generalized coordinates of the current and next two iteration steps converge, and then the calculation is stopped.

[0017] Further, step S1 includes: transforming the elastic motion of the beam by the variable u e , w, v, φ are described by geometric nonlinearity, and the variables are discretized using generalized coordinates q(t) and shape function h(ξ), and the shape function vector at the current Gaussian integration point and its related derivatives with respect to the beam axis coordinate x are calculated.

[0018] Furthermore, in step S1, the shape function is specifically a polynomial under the length of the beam unit, the maximum shape function of bending is 3rd order, and the maximum shape function of axial and torsional deformation is 4th order.

[0019] Further, in step S2, the torsion angle around the axis is represented by a derivative related to elastic bending.

[0020] Furthermore, in step S3, the derivative of the high-order term with respect to time is calculated by using trapezoidal integral numerical calculation until the integral converges.

[0021] Furthermore, in the method, the beam is specifically a deformable beam with a deformation angle less than 90 degrees.

[0022] The beneficial effects of the present invention are:

[0023] The geometrically accurate beam dynamics numerical modeling and calculation method proposed in the present invention significantly improves the accuracy and efficiency of solving medium and large deformation beam dynamics problems. Specifically, through geometric nonlinear description and the discretization of generalized coordinates and shape functions, it can accurately capture the coupling effect of elastic motion and rigid body motion of the beam, as well as the tensile torsion coupling effect of elastic deformation. In addition, the method uses numerical calculation methods such as trapezoidal integral and Gaussian integral to effectively deal with the high-order terms of torsional deformation and axial deformation caused by bending, ensuring the stability and accuracy of the calculation. It is particularly suitable for the design and analysis of deformable beams with a deformation angle less than 90 degrees, such as helicopter blades, wind turbine blades and other structures. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 A schematic diagram of a flow chart of a geometrically accurate beam dynamics numerical modeling and calculation method provided in an embodiment of the present application;

[0025] Figure 2 A schematic diagram of another flow chart of the geometrically accurate beam dynamics numerical modeling and calculation method provided in the embodiment of the present application;

[0026] Figure 3 Schematic diagram of elastic deformation of the beam model in the application embodiment. DETAILED DESCRIPTION

[0027] The present application is further described in detail below in conjunction with the accompanying drawings. It is necessary to point out here that the following specific implementation methods are only used to further illustrate the present application and cannot be understood as limiting the scope of protection of the present application. Technical personnel in this field can make some non-essential improvements and adjustments to the present application based on the above application content.

[0028] Example 1

[0029] like Figure 1-3 As shown, this embodiment proposes a geometrically accurate beam dynamics numerical modeling and calculation method, the method comprising the following steps:

[0030] S1. Obtain the shape function h(ξ) vector at the current Gaussian integration point and its related derivatives with respect to the beam axis coordinate x; specifically, the elastic motion of the beam is transformed from the variable u e , w, v, φ are described by geometric nonlinearity, and the variables are discretized using generalized coordinates q(t) and shape function h(ξ), and the shape function vector at the current Gaussian integration point and its related derivatives with respect to the beam axis coordinate x are calculated;

[0031] S2, calculating the elastic variables at the current Gaussian integration point of the beam, and its first-order derivative and second-order derivative with respect to time or beam axis coordinate x;

[0032] S3, using the elastic variables and their derivatives at the current Gaussian integration point, calculate the high-order terms of torsional deformation and axial deformation caused by bending through trapezoidal integration, and the derivatives of the high-order terms with respect to time, and discretize the coefficient matrix after the generalized coordinates are obtained;

[0033] S4, according to the Green-Lagrange strain tensor, using the calculation results of steps S1, S2, and S3 to calculate the strain measure of the cross section at the corresponding Gaussian integral point, obtain the coefficient matrix after discretizing the generalized coordinates, and discretize the equivalent stiffness matrix after the generalized coordinates;

[0034] S5, using the Gaussian integral method, integrating the calculation result of step S4 along the axial direction of the beam, calculating the strain energy term obtained by discretizing the variational generalized coordinates after the variation, and the equivalent stiffness matrix after discretizing the generalized coordinates;

[0035] S6. According to the D'Alembert principle, the acceleration and virtual displacement at the corresponding Gaussian integral point are calculated using the calculation results of steps S1, S2, and S3, and an equivalent mass matrix is ​​discretized;

[0036] S7, using the Gaussian integration method, integrating the product of the virtual displacement, acceleration and gravity obtained in step S6 along the beam axis, calculating the variation, and discretizing the kinetic energy term and equivalent mass matrix obtained by the variational generalized coordinates;

[0037] S8, calculating the work done by all external loads acting on the elastic axis according to the calculation results of steps S1, S2, and S3, and discretizing the external load terms of variational generalized coordinates after variation;

[0038] S9, according to Hamilton's principle, combined with the calculation results of steps S5, S7, and S8, establish the elastic motion equation and rigid body motion equation of the beam;

[0039] S10. Newton's iteration method is used, with generalized coordinates as the degrees of freedom. The generalized velocity and generalized acceleration are calculated using the Newmark-β method. The motion equations in step S9 are solved iteratively until the generalized coordinates of the current and next two iteration steps converge, and then the calculation is stopped.

[0040] Further preferably, in step S1, the shape function is specifically a polynomial under the beam unit length, the maximum shape function of bending is 3rd order, and the maximum shape function of axial and torsional deformation is 4th order.

[0041] It can be understood that by adopting higher-order shape functions (3rd order for bending, 4th order for axial and torsion), the elastic deformation of the beam can be described more accurately, especially when dealing with large deformations and complex nonlinear effects. This helps to capture more subtle deformation features, thereby improving the overall calculation accuracy. In addition, the use of higher-order shape functions helps to reduce the accumulation of numerical errors, especially in long-term simulation or iterative solution processes. This can significantly enhance the stability of the numerical method and ensure the reliability and consistency of the calculation results.

[0042] Further preferably, in step S2, the torsion angle around the axis is represented by the derivative related to elastic bending; in practical applications, through this processing, the torsion deformation of the beam when subjected to external force can be more accurately described, and this deformation is not only related to the magnitude of the external force, but also closely related to the elastic bending state of the beam. Representing the torsion angle as the derivative related to elastic bending can accurately capture this complex deformation relationship, thereby improving the accuracy of the dynamic equation.

[0043] Further preferably, in step S3, the calculation of the derivative of the high-order term with respect to time is specifically performed by trapezoidal integral numerical calculation until the integral converges; in practical applications, this feature ensures the accuracy and stability when calculating the derivative of the high-order term with respect to time. Since beam dynamics problems often involve complex high-order terms, direct solutions may face difficulties. The use of trapezoidal integral numerical calculation not only simplifies the calculation process, but also can obtain a high-precision solution by continuous iteration until the integral converges.

[0044] Further preferably, the beam in the method is specifically a deformable beam with a deformation angle less than 90 degrees.

[0045] Specific combination Figure 2 According to the above embodiment of the present invention, the implementation of the solution includes the following specific calculation and solution steps:

[0046] Step 1: The elastic motion of the beam is determined by the variable u e , w, v, φ are described by geometric nonlinearity, and the variables are discretized using generalized coordinates q(t) and shape functions h(ξ), see equation (1-1). Calculate the shape function h(ξ) vector at the current Gaussian integration point and its related derivatives with respect to the beam axis coordinate x. The shape function and its derivatives are shown in Table 1 and Table 2 for the shape functions and their derivatives for axial tension and torsion.

[0047]

[0048] The shape function is a polynomial under ξ = x / l (l is the length of the beam element). In the implementation, the maximum shape function of the bending is of order 3, and in order to exclude rigid motion, h(0) = h'(0) = 0 must be satisfied.

[0049] Table 1 Bending shape function and its derivatives

[0050]

[0051] The maximum shape function for axial and torsional deformation is of order 4:

[0052] Table 2 Shape functions and their derivatives for axial tension and torsion

[0053]

[0054]

[0055] Step 2: Calculate the elastic variables at the current Gaussian integration point of the beam, as well as their first-order and second-order derivatives with respect to time or the beam axis coordinate x. The torsion angles β and ζ around the y and z axes are expressed by the related derivatives of the elastic bending v and w, as shown in formula (1-2), where the leading S indicates the sine, C indicates the cosine, and the superscript + represents the derivative of the flexure curve r.

[0056]

[0057] Step 3: Through the variables and their derivatives at the current Gaussian integration point in step 2, trapezoidal integration is used to calculate the second-order and higher-order terms of axial stretch and torsion caused by v and w bending, as well as the time derivatives and coefficient vectors of the higher-order terms after the variations are discretized into generalized coordinates.

[0058] The precise calculation of the beam axial tension (U) and torsion (θ) variables is shown in formula (1-3):

[0059]

[0060] The calculation method of the second-order terms U2 and Θ2 can be expressed by the shape function calculated in step 1, as shown in formula (1-4):

[0061]

[0062] The calculation methods of the high-order terms ΔU, ΔΘ and their derivatives and variations with respect to time are shown in equations (1-5) and (1-7):

[0063]

[0064] Where D is the formula:

[0065]

[0066] For high-order terms, trapezoidal integration is used for numerical calculation. The formula is obtained by using M equally spaced trapezoidal integration.

[0067] (1-8):

[0068]

[0069] (1-8) In the formula:

[0070]

[0071]

[0072] Now we have:

[0073]

[0074] (1-11) In the formula:

[0075]

[0076] Then double the integration step:

[0077]

[0078] Starting from M=2, repeat this process until the integral converges, formula (1-14):

[0079] |I 2M -I M |<ε(referece) (1-14)

[0080] Step 4: According to the Green-Lagrange strain tensor, the results of steps 1, 2, and 3 are substituted into the corresponding formula to calculate the strain measure at the Gaussian integral point and the coefficient matrix after the variational discretization of the generalized coordinates.

[0081] The strain measurement and the calculation method of generalized coordinates of the dispersion are shown in formula (1-15):

[0082]

[0083] In formula (1-15), the calculation method of D is formula (1-16):

[0084]

[0085] Perform variational calculation on strain, see formula (1-17):

[0086]

[0087] According to the material constitutive relationship, the relationship between force and strain is obtained. For isotropic materials, the calculation method is shown in formula (1-18):

[0088]

[0089] The calculation method of each cross-sectional parameter in the above formula is shown in formula (1-19):

[0090]

[0091] For various anisotropic materials, the calculation method is shown in formula (1-20):

[0092]

[0093] The right-hand side of equations (1-18) and (1-20) can be decomposed as follows, see equation (1-21):

[0094]

[0095] Step 5: Use the Gaussian integration method to integrate the stress and strain in step 4 along the axial direction of the beam, calculate the variation of the strain energy, discretize the strain energy term of the variational generalized coordinates, and the equivalent stiffness matrix after discretizing the generalized coordinates.

[0096] The Gaussian integral of the entire beam element along the beam axis is calculated as:

[0097]

[0098] Where:

[0099]

[0100] where ξ i is the Gaussian integration point, w i is the weight at the integration point.

[0101] The strain energy term is calculated as:

[0102]

[0103] From formula (1-17), K in formula (1-24) can be written as:

[0104]

[0105] From formula (1-21), formula (1-24) can be written as:

[0106]

[0107] In the formula, S[i] represents the i-th row of the matrix S in equations (1-18) and (1-20), and MatrixK represents the equivalent stiffness matrix.

[0108] Step 6: Substitute the results of steps 1, 2, and 3 into the corresponding formula to calculate the acceleration and virtual displacement at the Gaussian integral point, and at the same time, the equivalent mass matrix can be discretized.

[0109] According to D'Alembert's principle, if the inertial acceleration is regarded as a kind of physical force, the virtual work of inertial force and gravity is:

[0110]

[0111] where g I is the acceleration in the inertial frame, usually gravity.

[0112] Where r I The calculation method of its first and second derivatives with respect to time is:

[0113]

[0114] R I Find the variation:

[0115]

[0116] M is used to represent the inertia term in the beam motion equation. Using the virtual displacement δx on the body axis BI / B and δψ BI / B Substituting the axis of inertia, the rigid body equations about the B axis are obtained. The displacement, velocity, acceleration and virtual displacement need to be calculated at a point on the beam relative to the body axis:

[0117]

[0118] Where η b , b are the coordinates of any point on the beam cross section after deformation.

[0119] The virtual displacement of the beam can also be expressed as:

[0120]

[0121] Therefore, R in formula (1-29) can be written as:

[0122]

[0123] Where R I The virtual displacement at the center of gravity is given.

[0124] Step 7: Integrate the product of the virtual displacement in step 6 and the inertial force and gravity along the beam axis, use the Gaussian integration method to calculate the variation of the kinetic energy, and discretize the variational generalized coordinate kinetic energy term and the equivalent mass matrix.

[0125] Therefore, the virtual work is:

[0126]

[0127] The inertia term in the beam motion equation can be obtained as follows:

[0128]

[0129] The elastic motion term in equation (1-34) can be decomposed into the following form according to (1-30) and (1-32):

[0130]

[0131] Where MatrixM is the equivalent mass matrix, and VectorR does not contain generalized acceleration terms.

[0132] Step 8. Substitute the results of steps 1, 2, and 3 into the corresponding formula, calculate the variation of the work of all external loads acting on the elastic axis, and discretize the variational generalized coordinates.

[0133] Discretize the variational generalized coordinates δq T The posterior external load term is:

[0134]

[0135] where F E 、M E is the force and moment at a certain point of the loaded elastic axis, The calculation method is as follows:

[0136]

[0137] in:

[0138]

[0139] Step 9: According to Hamilton's principle, equation (1-41), the elastic motion equation and rigid body motion equation of the beam are established based on the results calculated in steps 5, 7 and 8;

[0140]

[0141] The equation of motion for elasticity is:

[0142]

[0143] in is the structural damping, and the diagonal damping matrix is ​​calculated from the linear mass matrix and the spring matrix as M and K are shown in equations (1-26) and (1-35), where the structural damping constant g is twice the critical damping ratio. The rigid body motion equation is:

[0144]

[0145] Step 10. The present invention uses generalized coordinates as degrees of freedom, and the generalized velocity and generalized acceleration are calculated using the Newmark-β method, and the Newton iteration method is used for numerical solution. Given the initial value of the generalized coordinates, the generalized velocity and generalized acceleration are calculated using the Newmark-β method, and steps 1 to 9 are repeated to perform Newton iteration. When the generalized coordinates of the first and second iteration steps converge, the calculation is stopped.

[0146] According to the above-mentioned embodiment of the present invention, the method obtains the shape function vector and derivative at the Gaussian integral point, calculates the elastic variable and its time or space derivative, and then uses trapezoidal integration to calculate the high-order deformation terms and derivatives, and constructs the strain energy term and the kinetic energy term. In addition, combining the Green-Lagrange strain tensor and the Hamiltonian principle, the elastic and rigid body motion equations of the beam are established. This method not only has good numerical stability and computational efficiency, but also can accurately describe the nonlinear dynamic characteristics of large deformation structures such as blades. Therefore, the present invention has broad application prospects in the structural design of helicopter blades, wind turbine blades, etc., and can significantly improve the accuracy and efficiency of aeroelastic analysis.

[0147] The above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit it. Although the present application has been described in detail with reference to the aforementioned embodiments, a person skilled in the art should understand that the technical solutions described in the aforementioned embodiments can still be modified, or some of the technical features can be replaced by equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A geometrically accurate beam dynamics numerical modeling and calculation method, characterized in that: The method comprises: S1, obtain the shape function vector at the current Gaussian integration point and its related derivatives with respect to the beam axis coordinate x; S2, calculating the elastic variables at the current Gaussian integration point of the beam, and its first-order derivative and second-order derivative with respect to time or beam axis coordinate x; S3, using the elastic variable and its derivative at the current Gaussian integration point, calculating the torsional deformation high-order term and axial deformation high-order term caused by bending, and the time derivative of the high-order term through trapezoidal integration, and discretizing the coefficient matrix after the generalized coordinates; S4, according to the Green-Lagrange strain tensor, using the calculation results of steps S1, S2, and S3 to calculate the strain measure of the cross section at the corresponding Gaussian integral point, obtain the coefficient matrix after discretizing the generalized coordinates, and discretize the equivalent stiffness matrix after the generalized coordinates; S5, using the Gaussian integral method, integrating the calculation result of step S4 along the axial direction of the beam, calculating the strain energy term obtained by discretizing the variational generalized coordinates after the variation, and the equivalent stiffness matrix after discretizing the generalized coordinates; S6. According to the D'Alembert principle, the acceleration and virtual displacement at the corresponding Gaussian integral point are calculated using the calculation results of steps S1, S2, and S3, and an equivalent mass matrix is ​​discretized; S7, using the Gaussian integration method, integrating the product of the virtual displacement, acceleration and gravity obtained in step S6 along the axial direction of the beam, calculating the variation, and discretizing the kinetic energy term and equivalent mass matrix obtained by the variational generalized coordinates; S8, calculating the work done by all external loads acting on the elastic axis according to the calculation results of steps S1, S2, and S3, and discretizing the external load terms of variational generalized coordinates after variation; S9, according to Hamilton's principle, combined with the calculation results of steps S5, S7, and S8, establish the elastic motion equation and rigid body motion equation of the beam; S10. Newton's iteration method is used, generalized coordinates are used as degrees of freedom, and generalized velocity and generalized acceleration are calculated using the Newmark-β method. The motion equation in step S9 is solved iteratively until the generalized coordinates of the current and next two iteration steps converge, then the calculation is stopped.

2. A geometrically accurate beam dynamics numerical modeling and calculation method according to claim 1, characterized in that: Step S1 includes: transforming the elastic motion of the beam by the variable u e , w, v, φ are described by geometric nonlinearity, and the variables are discretized using generalized coordinates q(t) and shape function h(ξ), and the shape function vector at the current Gaussian integration point and its related derivatives with respect to the beam axis coordinate x are calculated.

3. The geometrically accurate beam dynamics numerical modeling and calculation method according to claim 1 is characterized in that: In step S1, the shape function is specifically a polynomial under the length of the beam unit, the maximum shape function of bending is the third order, and the maximum shape function of axial and torsional deformation is the fourth order.

4. The geometrically accurate beam dynamics numerical modeling and calculation method according to claim 1 is characterized in that: In step S2, the torsion angle around the axis is represented by the relevant derivative of elastic bending.

5. The geometrically accurate beam dynamics numerical modeling and calculation method according to claim 1 is characterized in that: In step S3, the derivative of the high-order term with respect to time is calculated by using trapezoidal integral numerical calculation until the integral converges.

6. The geometrically accurate beam dynamics numerical modeling and calculation method according to claim 2 is characterized in that: In the method, the beam is specifically a deformable beam with a deformation angle less than 90 degrees.

Images