Helicopter complex configuration blade modeling method based on multi-body dynamics

Abstract

The invention provides a helicopter complex configuration blade modeling method based on multi-body dynamics, and belongs to the field of helicopter blade modeling. Firstly, according to the configuration characteristics of the blade, the blade is divided into a plurality of sub-modules, and kinetic models of the sub-modules are mutually independent. Secondly, according to the motion and deformation characteristics of the sub-modules, a kinetic model of each sub-module is established; thirdly, according to the multi-body dynamics thought, displacement coordination boundary conditions between the submodules are introduced, a constraint equation is established, and the continuity of boundary node displacement is forced; fourthly, on the basis of the first-class Lagrange equation, combining the kinetic equations and constraint conditions of all the sub-modules, and assembling to form a modular assembly blade kinetic model. And finally, a generalized-alpha method is used for numerical discretization of the modular assembled blade kinetic model, and a Newton iteration method is used for solving a discretized linear equation set. According to the method, the modular multi-body dynamics thought is utilized, and the method is of great significance to helicopter blade structure design, configuration analysis and the like.

Description

Technical Field

[0001] The invention belongs to the field of helicopter blade modeling and relates to a modeling method of a helicopter blade with complex configuration based on multi-body dynamics. Background Art

[0002] Modern helicopter blades have evolved into complex configurations such as upward, downward, forward, and swept structures, such as the Blue Edge blade (forward and swept type), BERP blade (swept type), etc. The tip parts of these blades are designed to be non-straight. Although this has improved the overall performance of the helicopter to some extent, it has also brought certain difficulties to the unified modeling of different types of blades. Studies have shown that the swept configuration has a greater impact on the structural loads at the root and turning point of the blade. If straight blades are still used for modeling, it will cause large errors. The present invention conducts relevant research on the problem of efficient and accurate unified modeling of blades with different configurations.

[0003] Nowadays, among the modeling methods for helicopter blades of different configurations, the three-dimensional solid finite element modeling method is one of the most accurate methods. This method can present more structural details of the blades according to the blade configuration through finite element mesh division, and can achieve modeling for blades of different configurations. For example, the Chinese invention patent [CN111339607B] provides a three-dimensional model automation modeling method for D-shaped beam blades, which calls CATIACOM components through the VBwindows application platform and uses ActiveX automation technology to achieve rapid modeling of blades. However, while this modeling method brings high-precision models, the calculation efficiency will be greatly reduced, because the establishment of an accurate three-dimensional solid model will lead to a large degree of freedom of the dynamic model, and when combined with the aerodynamic model simulation, the solution rate will be particularly low. In order to solve the problem of calculation efficiency, the Chinese invention patent [CN116305589B] proposes a blade structure reduction analysis method based on a finite element mesh model of a three-dimensional unit cell, and establishes an efficient low-degree-of-freedom simplified model of helicopter blades, which can be applied to the efficient simulation and stress solution of structural response when the blade is subjected to force and deformation. Although this method can reduce the freedom of the blade model to a certain extent and improve the simulation efficiency to a certain extent, it is necessary to re-establish the finite element mesh of the three-dimensional unit cell for blades of different configurations, which is a relatively complicated task. In addition, this method is also difficult to couple with the aerodynamic model. After the three-dimensional unit cell is established, it is necessary to apply external loads of different working conditions to the unit cell and calculate the equivalent attribute matrix, which limits its scope of application. Another type of helicopter blade modeling method is to approximate the blade as a slender beam for modeling. The blade modeling and analysis based on geometric nonlinear beam theory greatly improves the calculation efficiency and does not result in large calculation errors. It is also the most commonly used blade modeling method. However, this method is based on traditional finite element assembly. For different blade configurations, the traditional finite element modeling method needs to re-derive some formulas, which brings a large workload. Therefore, it is not suitable for unified modeling of different blade configurations.

[0004] For the modeling of blades with different configurations, especially those with complex configurations, there are various difficult problems, such as the low efficiency of automatic three-dimensional modeling with large degrees of freedom, the need to re-establish the three-dimensional unit cell finite element mesh for blades with different configurations using the reduced-order modeling method, and the complexity of the nonlinear beam theory based on traditional finite elements for different configurations. Summary of the invention

[0005] In view of the problems existing in the prior art, the present invention utilizes the modular modeling concept and proposes a modeling method for helicopter blades with complex configurations based on multi-body dynamics. The present invention divides the blades into multiple sub-modules and introduces displacement coordination boundary conditions between the sub-modules, so that the displacements between different sub-modules are continuous. In particular, for blades with complex configurations, multiple sub-modules can be divided, and the structural forms of blades with different configurations can be achieved only by adjusting the boundary condition parameters, thereby being adaptable to the modeling of blades with complex configurations. The blade dynamics model based on multi-body dynamics can expand the scope of analysis, use the Kane method to derive the generalized forces of each sub-module, and establish the dynamic connection between adjacent sub-modules by introducing constraint conditions. Then, the coefficient matrices between adjacent sub-modules are decoupled, and adding or deleting sub-modules does not require re-deriving the coefficient matrices of other sub-modules, and the mechanical properties of the sub-modules can be analyzed independently.

[0006] In order to achieve the above object, the technical solution adopted by the present invention is:

[0007] A method for modeling a helicopter blade with complex configuration based on multi-body dynamics comprises the following steps:

[0008] Step 1: Divide the blade into several sub-modules according to its configuration characteristics, and establish a coordinate system to describe different motions. The dynamic models of each sub-module are independent of each other. Specifically:

[0009] The blade is divided into several submodules according to different configurations. The division criterion is: the blade is divided into different submodules before and after the turning position along the axis, that is, the axis of each submodule is guaranteed to be straight. At the same time, an inertial coordinate system I is established, and all the movements of each submodule need to be converted to this inertial coordinate system; a rigid body motion node is set up at the leftmost node of each submodule, and a module coordinate system E is established at this rigid body motion node. This rigid body motion node is used as the origin of the module coordinate system E and is recorded as ; An elastic motion node is set up at the rightmost end, denoted as the right end node A; a cross-sectional coordinate system B is set up at any point on the internal axis of the submodule, with the origin located on the axis of the submodule, and the origin is denoted as point P, which is used to describe the position of any point D on the cross section.

[0010] Step 2: According to the motion and deformation characteristics of the submodules, the generalized forces of the submodules are derived based on the Kane method, and the dynamic models of each submodule are established. Specifically:

[0011] Step 2-1: Definition of coordinates and motion description of each submodule;

[0012] First, define the coordinate variables of each node of the submodule, and set the generalized coordinates of the i-th submodule to be , represented by the position coordinates and Euler angle coordinates of the left end node and the deformation of the right end node:

[0013] (1)

[0014] In the formula, is the origin of the module coordinate system E of the i-th submodule Position coordinates expressed in inertial coordinate system I; is the origin of the module coordinate system E of the i-th submodule Euler angle coordinates expressed in inertial coordinate system I; is the elastic deformation degree of freedom of the right end node of the i-th submodule in the module coordinate system E; the superscript T represents the transpose operation of the vector.

[0015] Then, define the position coordinates of any point D on the i-th submodule in the module coordinate system E, expressed as:

[0016] (2)

[0017] In the formula, is the position coordinate of point D in the module coordinate system E. The superscript E indicates that the position coordinate is measured in the module coordinate system E. The subscripts i and D indicate point D on the i-th submodule. x is the distance from point D to the origin of the coordinate system in the module coordinate system E. The length of the axis; The origin of the cross-sectional coordinate system B is along the axis in the module coordinate system E. The deformation in the direction, represents the deformation along the x-axis, represents the deformation along the y-axis, represents the deformation along the z-axis; is the coordinate of point D in the cross-sectional coordinate system B, represents the y-axis coordinate of the cross-sectional coordinate system B, represents the z-axis coordinate of the cross-sectional coordinate system B; is the transformation matrix from the section coordinate system B to the module coordinate system E after the submodule is deformed. The specific expression is as follows:

[0018] (3)

[0019] In the formula, It represents the Euler angle from the deformed section coordinate system B to the module coordinate system E. Specifically: It is the total rotation angle of the cross-sectional coordinate system B relative to the module coordinate system E in the x-axis direction, including the pitch angle of the blade , pre-twist angle and elastic torsion angle ,Right now ; It is the rotation angle of the cross-section coordinate system B relative to the module coordinate system E in the y-axis direction; It is the rotation angle of the cross-section coordinate system B relative to the module coordinate system E in the z-axis direction.

[0020] According to equations (1) and (2), the position coordinates of any point D on the i-th submodule after deformation in the inertial coordinate system I are:

[0021] (4)

[0022] In the formula, is the position coordinate of point D in the inertial coordinate system I, is the transformation matrix from the module coordinate system E to the inertial coordinate system I. The specific expression is to transform the Euler angle in equation (3) Replace it with .

[0023] Step 2-2: derive the generalized inertia force of the i-th submodule;

[0024] According to the Kane method, the generalized inertia force of the i-th submodule is expressed as:

[0025] (5)

[0026] in, are the first-order and second-order time derivatives of the position coordinates of any point D on the i-th submodule after deformation in the inertial coordinate system I, that is, the velocity and acceleration expressed in the inertial coordinate system I; is the submodule material density; is the first-order time derivative of the generalized coordinates of the i-th submodule, also called the generalized rate; represents the differential of the cross-sectional area; Represents the differential of the submodule axis length.

[0027] Therefore, the velocity of any point D on the i-th submodule in the inertial coordinate system I is and acceleration :

[0028] (6)

[0029] In the formula, The origin of the module coordinate system E of the i-th submodule The first and second time derivatives of the position coordinates expressed in the inertial coordinate system I; They are the origin of the module coordinate system E Angular velocity and angular acceleration; The vectors are The antisymmetric matrix of .

[0030] According to the projection relationship between the time derivative of the Euler angle and the angular velocity in the multi-body dynamics theory, we have:

[0031] (7)

[0032] In the formula, They are The first and second time derivatives of ; are the projection matrix and its first-order time derivative, the projection matrix Specifically, it can be written as:

[0033] (8)

[0034] In the formula, is the Euler angle coordinate in formula (1) The first two elements in .

[0035] From equations (1), (6) and (7), we can deduce the velocity D of any point on the i-th submodule: For generalized rate The first derivative of :

[0036] (9)

[0037] In the formula, is the three-dimensional identity matrix; is the interpolation function matrix.

[0038] According to equations (6) and (7), the acceleration of any point D on the i-th submodule in the inertial coordinate system I can be Arranged into the following form:

[0039] (10)

[0040] Finally, substitute equations (9) and (10) into the expression (5) of the generalized inertial force to obtain the expression of the generalized inertial force of the i-th submodule, and then organize it to obtain:

[0041] (11)

[0042] In the formula, is the beam material density; are the generalized mass matrix and nonlinear term of the ith submodule, respectively, which are expanded as follows:

[0043] (12)

[0044] In the formula, the variable T in the upper right corner is the transposition symbol, and other variables have been defined in the above formulas.

[0045] Step 2-3: Derive the generalized elastic force of the i-th submodule;

[0046] The generalized elastic force of the ith submodule is obtained by taking the partial derivative of the strain energy of the submodule with respect to the generalized coordinates, so it is necessary to derive the strain energy of the ith submodule. According to the definition of strain energy, strain energy is calculated from stress and strain, so it is necessary to calculate the stress-strain expression of the submodule, which is deduced in turn below.

[0047] First, the strain of the submodule is derived. For the blade, the strain is small. Using the small strain assumption, it is considered that the Green strain is approximately the engineering strain, so the strain-displacement relationship is:

[0048] (13)

[0049] In the formula, is the axial strain; is the shear strain; The origin of the cross-sectional coordinate system B is along the axis in the module coordinate system E. The first-order arc-length derivative of the deformation in the direction with respect to the axial coordinate, with the superscript That is, it represents the first-order arc-length derivative; Indicates deformation The second arc-length derivative of ; Indicates deformation The second arc-length derivative of ; is the coordinate of any point D on the ith submodule in the cross-sectional coordinate system B, represents the y-axis coordinate of the cross-sectional coordinate system B, represents the z-axis coordinate of the cross-sectional coordinate system B; It represents the total rotation angle of the section coordinate system B relative to the module coordinate system E in the x-axis direction; It represents the elastic torsion angle of the section coordinate system B relative to the module coordinate system E in the x-axis direction. First arc-length derivative with respect to the axial coordinate.

[0050] Then, the stress of the submodule is derived. The blade is a slender structure, and the uniaxial stress assumption is used, that is, the stress other than the axial stress is ignored, that is, , so the stress-strain relationship is expressed as:

[0051] (14)

[0052] In the formula, is the axial stress; is the shear stress; is Young's modulus; is the shear modulus.

[0053] Furthermore, the strain energy of the submodule is calculated. According to the definition of strain energy, the strain energy of the ith submodule can be obtained using equations (13) and (14): The expression is transformed into the form of cross-sectional load:

[0054] (15)

[0055] In the formula, They are the generalized force strain and generalized moment strain of the submodule respectively; are section force and section moment respectively; the matrix is the 6-dimensional section stiffness matrix; is the length of the submodule; represents the differential of the cross-sectional area; Represents the differential of the submodule axis length.

[0056] Finally, the generalized elastic force of the submodule is calculated. Calculation of strain energy For generalized coordinates The first-order derivative of is used to obtain the generalized elastic force of the i-th submodule and is written in matrix form:

[0057] (16)

[0058] In the formula, is the stiffness matrix of the ith submodule, and the numerical solution is calculated using the difference method; is a nonlinear term.

[0059] Step 2-4: Derive the generalized active force of the i-th submodule;

[0060] According to the discrete format of the active force in Kane's method, the generalized active force term of the i-th submodule is derived: . For any point P on the axis of the submodule, the main force acting is , then the rth component of the generalized active force is It can be expressed as:

[0061] (17)

[0062] In the formula, is the rth eccentric velocity of point P in the inertial coordinate system I; is the number of generalized coordinates of the submodule.

[0063] remember is the radius vector of point P in the inertial coordinate system I, then the eccentric velocity of point P is It is expressed as:

[0064] (18)

[0065] In the formula, represents the r-th generalized coordinate of the ith submodule.

[0066] Substituting equation (18) into equation (17), we get the discrete formula of the rth component of the generalized active force:

[0067] (19)

[0068] Therefore, all The generalized force components are written in vector form, and the generalized active force of the i-th submodule is expressed as:

[0069] (20)

[0070] In the formula, represents the generalized active force column vector; represents the radius vector of any point P on the axis of the i-th submodule; Represents the column vector of generalized coordinates of the i-th submodule.

[0071] Step 3: Based on the theory of multi-body dynamics, introduce the displacement coordination boundary conditions between sub-modules, establish constraint equations, and enforce the continuity of boundary node displacements. Specifically:

[0072] Step 2 describes the derivation of the generalized dynamic forces of a single submodule. The key to modular multibody modeling applicable to blades of various configurations is to establish constraint equations between submodules, realize the dynamic connection between different submodules, and thus obtain complete dynamic equations. For blades, constraint equations include the following four categories.

[0073] Step 3-1: First type constraint equation;

[0074] The first type of constraint is introduced by the connection condition between the blade and the hub, that is, the connection between the first submodule of the blade and the hub. For helicopter blades, the connection between the blade and the hub is related to the rotor type. The most common rotor types currently include: bearingless rotor, hingeless rotor, and articulated rotor. The three types of constraint equations are expressed as:

[0075] (twenty one)

[0076] In the formula, The constraint equations of three types of rotors are respectively represented: bearingless rotor, hingeless rotor, and articulated rotor; In expression (1) The element of , represented by the i-th submodule in formula (1), and the first type of constraint equation is the connection between the first submodule of the blade and the hub, so the subscript i is written as 1.

[0077] Step 3-2: Second type constraint equation;

[0078] The second type of constraint is introduced by the connection between submodules. In step 1, the blade is divided into multiple submodules. There is no relative displacement and rotation between the submodules. According to the multi-body dynamics theory, it is necessary to establish a connection through constraint equations. Therefore, there is the following constraint equation between the i-th submodule and the i+1-th submodule:

[0079] (twenty two)

[0080] In the formula, They represent the position coordinates and Euler angle coordinates of the right end node A of the i-th submodule in the inertial coordinate system I respectively; Respectively represent the left end node of the i+1th submodule Position coordinates and Euler angle coordinates in inertial coordinate system I.

[0081] Assuming that the blade is divided into N submodules, there are N-1 constraint equations such as equation (22), which can be written as:

[0082] (twenty three)

[0083] In the formula, represents the second type of total constraint equation, whose elements are the subscript i of equation (22) rewritten as (1, 2, …, N-1); the superscript T represents the transposition symbol.

[0084] Step 3-3: The third type of constraint equation;

[0085] The third type of constraint is introduced by the motion of the rotor. The rotor works at a constant speed, so the rotor speed is introduced into the blade dynamics model through the drive constraint. Suppose the rotor rotation angular velocity , time t, the driving constraint is written as:

[0086] (twenty four)

[0087] In the formula, The third Euler angle coordinate of the first submodule given in equation (1) is expressed as follows. Since the drive constraint only needs to apply the rotational angular velocity to the first submodule of the blade, the rotational angular velocity can be transferred to all subsequent submodules through the second type constraint equation.

[0088] Step 3-4: The fourth type of constraint equation;

[0089] The fourth type of constraint is the constraint between the front and rear submodules of blades of different configurations at the turning point of the blade axis. In step 1, two different submodules are divided before and after the turning point of the blade axis. Due to the turning point, the angle between the front and rear submodules is less than 180°, and this angle is introduced in the form of constraints. Assume that the submodule after the turning point is the i+1th submodule, which is at a negative angle to the i-th submodule. , the constraint equation is written as:

[0090] (25)

[0091] In the formula, It represents the second element of the Euler angle coordinates of the i+1th submodule in equation (1).

[0092] Step 3-5: Total constraint equation;

[0093] Combining equations (21), (23), (24) and (25), we can obtain the total constraint equation of the dynamic model:

[0094] (26)

[0095] Step 4: Based on the first kind of Lagrangian equation, Lagrangian multipliers are introduced to combine the dynamic equations and constraint equations of each submodule to form a blade dynamic model of modular group. Specifically:

[0096] According to equations (11), (16) and (20), the blade is divided into N sub-modules of generalized inertial force, generalized elastic force and generalized active force, and the groups are as follows:

[0097] (27)

[0098] in, represents the total generalized inertia force of N submodules, is the total generalized elastic force of N submodules, is the total generalized main driving force of N sub-modules.

[0099] According to the Kane equation, the blade balance equation of N submodule groups is obtained:

[0100] (28)

[0101] Substituting equation (27) into equation (28), while considering the total constraint equation shown in equation (26), and introducing Lagrange multipliers based on the first-kind Lagrange equation, we obtain a set of differential algebraic equations to extract the generalized acceleration and generalized coordinates The linear term coefficients are calculated and the blade dynamics model of N submodules is obtained:

[0102] (29)

[0103] In the formula, They are the total mass matrix, stiffness matrix and nonlinear terms of N submodules respectively; represents the constraint equation and is the generalized coordinate and a function of time t; is the constraint equation for generalized coordinates The transpose of the Jacobian matrix is ; is the Lagrange multiplier.

[0104] According to formula (27), each submodule is completely independent and is only connected by constraint equation (26). Therefore, the coefficient matrix of the blade dynamics model (29) is obtained: and vector Different sub-modules are decoupled, so adding or deleting sub-modules will not affect the coefficients of existing sub-modules, making it easier to model blades with different configurations.

[0105] Step 5: Use generalized- Methods The blade dynamics model of the modular assembly was numerically discretized and the discretized linear equations were solved using Newton's method. Specifically:

[0106] Step 5-1: Use generalized- Method numerical discrete;

[0107] The generalized-α method is used to numerically discretize the blade dynamics model shown in formula (29). Assuming that the blade dynamic response within a time period T needs to be solved, the integral time step is h, that is, the time T is divided into steps of h. time steps, then the nth time step The blade dynamics model is expressed as the residual In the form of:

[0108] (30)

[0109] In the formula, represents the nth time step, and the subscript n of other variables also represents the numerical results at the nth time step; Respectively represent the generalized coordinates, generalized velocity, and generalized acceleration of each submodule at the nth time step; represents the Lagrange multiplier of the nth time step; represents the nonlinear term at the nth time step; Represents the constraint equation for the nth time step.

[0110] According to the broad- The discrete form of the method, Discrete into:

[0111] (31)

[0112] Where h is the integration time step, that is, ; The subscript n-1 or n of each variable represents the value of the variable at the n-1th or nth time step; are algorithm parameters; is an auxiliary variable of the generalized-α method, and , Satisfaction relationship:

[0113] (32)

[0114] In the formula, Represents auxiliary variables The initial value of is the generalized acceleration The initial value of is the algorithm parameter, satisfying:

[0115] (33)

[0116] In the formula, is the spectral radius of the generalized-α method.

[0117] Substituting equations (31) and (32) into equation (30), we can obtain the complete discrete formula of the blade dynamics model.

[0118] Step 5-2: Use Newton's method to solve the linear equations discretized by the generalized-α method;

[0119] Use generalized - After the numerical discretization of the method, a set of linear equations is obtained, and the variables in the nth time step are recorded as , in each time step, the Newton method is used to iteratively solve, and the iterative format is as follows:

[0120] (34)

[0121] In the formula, the superscript Represents Newton's method Iterations; is the residual of the kth iteration of the equation system at the nth time step; is the increment of the kth iteration; is the Jacobian matrix of the kth iteration of the system of equations at the nth time step.

[0122] also, The specific expression is:

[0123] (35)

[0124] Find the increment value of the current iteration step After that, the unknown variables of the blade dynamics model are updated as follows:

[0125] (36)

[0126] In the formula, the subscript Indicates the variable at the nth time step; superscript The variable representing the k+1th iteration of the current time step; and are the process quantities to be solved, respectively: .

[0127] Substitute the variables in formula (36) into formula (34), satisfy , that is, the residual The 2-norm of is less than or equal to the set error limit , the iteration step is completed, and the increments of the generalized coordinates, generalized velocity, generalized acceleration and Lagrange multiplier at that moment can be obtained. .

[0128] So far, completed The task of solving the blade dynamics model at the moment.

[0129] Step 5-3: Determine whether the solution of all discrete time steps is completed;

[0130] Finally, determine whether all time steps in time period T have been solved. Is it true? If it is true, then end the calculation. If not, let the variable value of the nth time step (Equation (36)) be the initial value of the n+1th time step, and substitute it into equations (30), (31) and (32) in step 5-1, and update the time step at the same time, that is, let , and obtain the discrete formula of the blade dynamics model at the n+1th time step; then use the Newton method iterative formulas (34) and (35) in step 5-2 to solve the variable value of the kth iteration step of the n+1th step; finally, execute step 5-3 to determine whether the solution of all time steps is completed, and repeat this cycle until the entire solution task is completed.

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

[0132] (1) Compared with the traditional finite element modeling method, the present invention does not need to re-derive the dynamic equations for different blade configurations. It only needs to adjust the constraint equation parameters between sub-modules to achieve modeling of blades with different complex configurations. The present invention divides the complex configuration blade into several sub-modules, especially at the turning point of the complex configuration blade, establishes the constraint equations between the sub-modules based on the multi-body dynamics theory, and achieves modeling of blades with different complex configurations by adjusting the constraint parameters.

[0133] (2) The present invention uses an implicit integration algorithm, the generalized-α method, to solve the established differential algebraic equations. The traditional explicit integration algorithm needs to use methods such as extracting independent degrees of freedom to eliminate algebraic equations for differential algebraic equations, and the processing process is relatively complicated. The generalized-α method provided by the present invention can introduce Lagrange multiplier terms and directly use a discrete format to process differential algebraic equations, discretize them into a system of linear equations, and then use Newton's method to iteratively solve them. Combined with the module-level difference algorithm, the efficiency of equation solving is further improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0134] Figure 1 It is a calculation flow chart of the present invention.

[0135] Figure 2 This is an example of module division of the rectangular blade configuration listed in the present invention.

[0136] Figure 3 This is an example of module division of the blade configuration with a protruding tip and swept back listed in the present invention.

[0137] Figure 4 This is an example of module division of the reverse blade configuration under the blade tip listed in the present invention.

[0138] Figure 5 Schematic diagram of the rotating straight blade structure under dynamic load in Example 1 of the present invention.

[0139] Figure 6 The inertial coordinate system I, module coordinate system E and cross-sectional coordinate system B established for the present invention, the defined left and right end nodes and internal points of the submodule, and the overall deformation schematic diagram of the submodule.

[0140] Figure 7 It is a schematic diagram of the constraint types between adjacent submodules of the present invention.

[0141] Figure 8 Schematic diagram of the assembly of different submodule elements in the coefficient matrix of the blade dynamics model established for the present invention.

[0142] Fig. 9 A schematic diagram of the assembly of different submodule elements in the nonlinear terms of the blade dynamics model system established for the present invention.

[0143] Fig.10 This is a graph showing how the axial displacement of the blade tip changes with time when the blade tip of a rotating straight blade under dynamic load in Example 1 of the present invention is subjected to a concentrated external force F.

[0144] Fig.11 This is a graph showing the variation of the oscillation displacement of the blade tip with time when the blade tip of the rotating straight blade under dynamic load in Example 1 of the present invention is subjected to a concentrated external force F.

[0145] Fig.12 This is a graph showing the variation of the flapping displacement of the blade tip with time when the blade tip of the rotating straight blade under dynamic load in Example 1 of the present invention is subjected to a concentrated external force F.

[0146] Fig.13 This is a graph showing how the torsion angle of the blade tip changes with time when the blade tip of a rotating straight blade under dynamic load in Example 1 of the present invention is subjected to a concentrated external force F.

[0147] Fig.14 This is a graph showing the variation of the flap angle of the blade tip with time when the blade tip of the rotating straight blade under dynamic load in Example 1 of the present invention is subjected to a concentrated external force F.

[0148] Fig.15 This is a graph showing the variation of the flapping angle of the blade tip with time when the blade tip of the rotating straight blade under dynamic load in Example 1 of the present invention is subjected to a concentrated external force F.

[0149] Fig.16 This is a top view of the structure of the blade with a swept-tip configuration in Example 2 of the present invention.

[0150] Fig.17 This is a structural side view of a blade with a swept-tip configuration in Example 2 of the present invention.

[0151] Fig.18 is the sweep angle of the blade with swept tip configuration in Example 2 of the present invention Comparison of the X and Y displacements of the blade tip when the angle is 45° and the calculation results of Ansys software.

[0152] Fig.19 is the sweep angle of the blade with swept tip configuration in Example 2 of the present invention When the angle is 15°, the comparison chart between the first to sixth order natural frequencies and the test results at speeds of 500r / min and 750r / min.

[0153] Fig. 20 is the sweep angle of the blade with swept tip configuration in Example 2 of the present invention When the angle is 45°, the comparison chart between the first to sixth order natural frequencies and the test results at speeds of 500r / min and 750r / min. DETAILED DESCRIPTION

[0154] The following is a further detailed description of the embodiments of the present invention in conjunction with the accompanying drawings and examples. The following examples are used to illustrate the present invention, but are not intended to limit the scope of the present invention.

[0155] The calculation flow chart of the present invention is as follows: Figure 1 shown. Figure 2~Figure 4 Three common configurations of helicopter blades are listed, and the corresponding blade structure submodule division method is given. For rectangular straight blades, the number of submodules can be divided arbitrarily, while for complex configuration blades, such as Figure 3 and Figure 4 When dividing the submodules, it is necessary to divide them at the turning position along the axial direction of the blade to ensure that each submodule is straight. For any blade configuration, module division can be carried out according to this method; then according to Figure 1 The implementation process of any blade configuration is to model the sub-module structure, and different modules are connected using Figure 7 The sub-modules are connected in such a way that they can be assembled into a complete blade dynamics model.

[0156] In order to make the purpose, technical solution and advantages of the present invention more clear, the following two specific embodiments are combined with the attached Figures 5 to 20 The advantages, accuracy and effectiveness of the present invention are further described in detail.

[0157] (1) Example 1: Dynamic response analysis of a rotating straight blade with a dynamic load on the blade tip;

[0158] like Figure 5 As shown, this embodiment is a rectangular straight blade subjected to a dynamic load F at the blade tip, wherein F=50sin(20t) is a simple harmonic force that varies with time, to verify the applicability and effectiveness of the modeling method of the present invention under the straight blade configuration, wherein the blade structural parameters are shown in Table 1.

[0159] Table 1. Structural parameters of blades subjected to dynamic loads

[0160]

[0161] Step 1: According to the blade configuration characteristics, it is a straight blade configuration, and there is no turning section along the blade axis. It is divided into 4 sub-modules, such as Figure 5 At the same time, Figure 6 , establish an inertial coordinate system I, and convert all movements of each submodule to this coordinate system; establish a module coordinate system E at the leftmost node of each submodule, and use this node as the origin of the module coordinate system E ; An elastic motion node A is set up at the rightmost end; a cross-sectional coordinate system B is set up at any point P on the internal axis of the submodule, with the origin P located on the axis of the submodule, describing the position of any point D on the cross section.

[0162] Step 2: According to the motion and deformation characteristics of the submodules, the generalized forces of the submodules are derived based on the Kane method, and the dynamic models of each submodule are established. Specifically:

[0163] Step 2-1: Definition of coordinates and motion description of each submodule;

[0164] First, define the coordinate variables of each node of the submodule, and set the generalized coordinates of the i-th submodule to be , as shown in formula (1); then, define the position coordinates of any point D on the ith submodule in the module coordinate system E before and after deformation, as shown in formula (2); finally, according to formulas (1) and (2), define the position coordinates of any point D on the ith submodule in the inertial coordinate system I after deformation, as shown in formula (4).

[0165] Step 2-2: derive the generalized inertia force of the i-th submodule;

[0166] According to the Kane method, the generalized inertia force of the i-th submodule is derived, as shown in formula (5). Then, by substituting formulas (9) and (10) into formula (5), the generalized inertia force expression of the i-th submodule can be obtained and organized into the form of formula (11).

[0167] Step 2-3: Derive the generalized elastic force of the i-th submodule;

[0168] The generalized elastic force of the ith submodule is obtained by taking the partial derivative of the strain energy of the submodule with respect to the generalized coordinates. Therefore, it is necessary to first derive the stress and strain of the ith submodule, and then calculate the strain energy, and finally obtain the generalized inertia force of the ith submodule. Use equations (13), (14), and (15) in sequence to calculate the strain, stress, and strain energy of the submodule, and finally use equation (16) to calculate the generalized elastic force of the ith submodule.

[0169] Step 2-4: Derive the generalized active force of the i-th submodule;

[0170] According to the discrete format of the active force in Kane's method, the generalized active force term of the i-th submodule is derived: The eccentric velocity of point P is calculated using equation (18), and then substituted into equation (17) to obtain the generalized active force of the ith submodule as equation (20).

[0171] Step 3: Based on the theory of multi-body dynamics, introduce the displacement coordination boundary conditions between sub-modules, establish constraint equations, and enforce the continuity of boundary node displacements. Specifically:

[0172] The key to modular multi-body modeling is to establish constraint equations between sub-modules, realize the dynamic connection between different sub-modules, and thus obtain complete dynamic equations. For blades, constraint equations include the following four categories.

[0173] Step 3-1: First type constraint equation;

[0174] The first type of constraint is introduced by the connection condition between the blade and the hub, that is, the connection between the first submodule of the blade and the hub. Figure 3 As shown, the blade and the hub in this embodiment are fixedly constrained, that is, the connection mode of the bearingless rotor. Therefore, the first type of constraint equation in this embodiment uses the formula (21) .

[0175] Step 3-2: Second type constraint equation;

[0176] The second type of constraints is introduced by the connections between submodules. Figure 7 As shown, the connection method between different sub-modules is given. According to the multi-body dynamics theory, it is necessary to establish the connection through constraint equations. Therefore, equation (22) is used to establish the constraint equations between the i-th sub-module and the i+1-th sub-module, and the final set is in the form of equation (23).

[0177] Step 3-3: The third type of constraint equation;

[0178] The third type of constraint is introduced by the motion of the rotor. The rotor is rotating at a constant speed, so the rotor speed is introduced into the blade dynamics model through the drive constraint. In this embodiment, the blade rotation angular velocity , the simulation time is t=1s, and the constraint is given by formula (24).

[0179] Step 3-4: The fourth type of constraint equation;

[0180] The fourth type of constraint is the constraint between the front and rear submodules of blades of different configurations at the turning point of the blade axis. Since this embodiment uses a straight blade, there is no turning point in the blade axis, so the constraints in this type are .

[0181] Step 3-5: Total constraint equation;

[0182] Finally, by combining equations (21), (23), and (24), we can obtain the total constraint equation of the blade dynamics model, as shown in equation (26).

[0183] Step 4: Based on the first kind of Lagrangian equation, Lagrangian multipliers are introduced to combine the dynamic equations and constraint equations of each submodule to form a blade dynamic model of modular group. Specifically:

[0184] According to equations (11), (16) and (20), the blade is divided into N sub-modules of generalized inertial force, generalized elastic force and generalized active force, and the generalized force of equation (27) is obtained by grouping them. Then, according to the Kane equation, the blade balance equation (28) of the N sub-modules can be obtained. Furthermore, considering the constraint equation (26), a set of differential algebraic equations is obtained, and the blade dynamics model of the N sub-modules is obtained, as shown in equation (29).

[0185] According to formula (27), each submodule is completely independent and is only connected by constraint equation (26). Therefore, the coefficient matrix of the blade dynamics model (29) is obtained: and vector , different submodules are decoupled, so adding or deleting submodules will not affect the coefficients of existing submodules, such as Figure 8 and Fig. 9 As shown, it is more convenient to model blades of different configurations.

[0186] Step 5: Use generalized- Methods The blade dynamics model of the modular assembly was numerically discretized and the discretized linear equations were solved using Newton's method. Specifically:

[0187] Step 5-1: Use generalized- Method numerical discrete;

[0188] The generalized-α method is used to numerically discretize the blade dynamics model shown in formula (29). The simulation time of this embodiment is T = 1s, and the blade dynamic response is calculated. The integration time step is h = 1×10 -3 s, that is, time T is divided into 1000 time steps with a step size of h. Use equations (31) and (32) to discretize and substitute them into equation (30) to calculate the residual , the spectral radius of the generalized-α method is .

[0189] Step 5-2: Use Newton's method to solve the linear equations discretized by the generalized-α method;

[0190] Use generalized - After numerical discretization, a set of linear equations is obtained. In each time step, the Jacobian matrix of formula (35) is calculated by the difference method, and then the Newton method of formula (34) is used to solve the equations. Finally, it is determined whether the residual meets the error limit, and the variables in formula (36) are substituted into formula (34). satisfy , that is, the residual The 2-norm of is less than or equal to the set error limit In this embodiment, the error limit is =1×10 -6 , the iteration step is completed, and the increments of the generalized coordinates, generalized velocity, generalized acceleration and Lagrange multiplier at that moment can be obtained. So far, the task of solving the blade dynamics model at that moment is completed.

[0191] Step 5-3: Determine whether the solution of all discrete time steps is completed;

[0192] Finally, determine whether all time steps within the time period T = 1s have been solved. Is it true? If it is true, then end the calculation. If not, let the variable value of the nth time step (Equation (36)) be the initial value of the n+1th time step, and substitute it into equations (30), (31) and (32) in step 5-1, and update the time step at the same time, that is, let , and obtain the discrete formula of the blade dynamics model at the n+1th time step; then use the Newton method iterative formulas (34) and (35) in step 5-2 to solve the variable value of the kth iteration step of the n+1th step; finally, execute step 5-3 to determine whether the solution of all time steps is completed, and repeat this cycle until the entire solution task is completed.

[0193] like Figure 10~Figure 15 As shown in the figure, the dynamic simulation results of the rotating straight blade under the action of external force are compared with the calculation results of the multi-body dynamics software DYMORE. The blade rotates at a speed of Ω=70rad / s, and the external force F=50sin(20t) acts on the blade tip and changes along the z-axis direction. Figure 10~Figure 12 is the linear displacement of the blade tip in the axial, shimmy and flapping directions, Figure 13~Figure 15 It is the angular displacement of the blade tip in the axial, swing and flapping directions. As can be seen from the figure, the calculation results of the modeling method of the present invention are basically consistent with the calculation results of the multi-body dynamics software DYMORE. Under high-speed rotation, the changes in the axial deformation curve and the torsion angle curve reflect obvious bending-stretching-torsion coupling phenomenon, which is also the typical dynamic characteristic of the helicopter blade in the working state. Therefore, the ability to accurately describe this phenomenon shows the effectiveness of the modeling method of the present invention.

[0194] (2) Example 2: Dynamic response and modal analysis of a complex swept blade;

[0195] Similarly, steps 1 to 5 are performed to carry out structural dynamic modeling of swept blades with complex configurations, and to solve the dynamic response. Since the blade tip in this embodiment is a swept configuration, there are two slight differences in executing steps 1 to 5 compared to embodiment one. The first is the number of sub-modules divided according to the blade configuration in step 1. For tip-swept blades, the straight blade section is divided into three sub-modules, and the swept blade section is divided into two sub-modules. It should be noted that different sub-modules must be divided at the turning point between the straight blade section and the swept section to ensure the subsequent application of constraints; the second difference is the establishment of the fourth type of constraints in steps 3-4. For this embodiment, the sweep angle is taken and , to calculate different working conditions.

[0196] The advantage of the present invention is that when modeling different blade configurations, it is only necessary to divide the submodules according to the configuration and modify the constraint parameters. By using the value of , the dynamic model of blades with different configurations can be obtained without any other changes, thus reflecting the advantages of the modeling method of the present invention.

[0197] like Fig.16 and Fig.17 The following are the top view and side view of the structure of the swept tip blade used in this embodiment. The length of the straight blade part is L 1 =0.8636m, the length of the swept part is L 2 =0.1524m, width b 1 =0.0254m, height h 1 =0.0016m, in this embodiment, in order to better compare the effectiveness of the modeling method of the present invention, the sweep angle Two working conditions, 15° and 45°, were selected for verification.

[0198] like Fig.18 The figure shows the displacement of the swept free end of the blade compared with the calculation results of Ansys software. The calculation condition is the swept angle. , the speed gradually increases from 0 to 1 second. After that, it keeps rotating at a constant speed. As can be seen from the figure, the calculation results of the modeling method of the present invention completely coincide with the Ansys result curve, which shows the accuracy of the modeling method of the present invention in calculating the large-scale rigid body displacement of the swept blade. Fig.19 and Fig. 20As shown, the first six natural frequency curves of the two sweep angles of 15° and 45° at 500r / min and 750r / min speeds are given respectively, of which the first five are flapping modes and the sixth is a torsional mode. First of all, it can be seen intuitively that compared with the speed of 500r / min, at the speed of 750r / min, the natural frequencies of each order at the same sweep angle are increased to a certain extent. Because the centrifugal force of the blade increases with the increase in speed, the natural frequency is increased. The modeling method proposed in the present invention also calculates this change trend that is consistent with the experimental results. Therefore, the comparison of modal calculation results also proves the accuracy and effectiveness of the modeling method of the present invention for modeling swept configuration blades.

[0199] In summary, by performing effective numerical simulation on the above two typical embodiments and comparing with the results of Ansys software or experimental results, the feasibility and development necessity of the modeling method of helicopter complex configuration blades based on multi-body dynamics proposed by the present invention can be well explained. The precise comparison of numerical results also illustrates the effectiveness and accuracy of the modeling method proposed by the present invention. Although the two embodiments given in the present invention are numerical simulations of bearingless blades, according to the constraint equation of formula (21) in step 3-1, it is easy to expand it into a blade dynamics model of any connection method. Therefore, the modeling method of helicopter complex configuration blades based on multi-body dynamics proposed by the present invention is an efficient modeling algorithm with great development potential.

[0200] The above-described embodiments are only for expressing specific implementation methods of the present invention, but they cannot be understood as limiting the patent scope of the present invention. It should be pointed out that for those skilled in the art, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A method for modeling helicopter blades with complex configuration based on multi-body dynamics, characterized in that: The following steps are involved: Step 1: Divide the blade into several submodules according to its configuration characteristics, and establish a coordinate system to describe different motions. The dynamic models of each submodule are independent of each other. Step 2: According to the motion and deformation characteristics of the submodules, the generalized forces of the submodules are derived based on the Kane method, and the dynamic models of each submodule are established; Step 3: Based on the multi-body dynamics theory, the displacement coordination boundary conditions between sub-modules are introduced, and the constraint equations are established to enforce the continuity of boundary node displacements; Step 4: Based on the first kind of Lagrangian equation, Lagrangian multipliers are introduced to combine the dynamic equations and constraint equations of each submodule to form a blade dynamic model of modular group set; Step 5: Use generalized- Methods The blade dynamics model of modular groups is numerically discretized and the discretized linear equations are solved by Newton's method.

2. The method for modeling a helicopter blade with complex configuration based on multi-body dynamics according to claim 1 is characterized in that: The step 1 is specifically as follows: The blade is divided into several submodules according to different configurations. The division criteria are as follows: the blade is divided into different submodules before and after the turning position along the axis to ensure that the axis of each submodule is straight; at the same time, an inertial coordinate system I is established, and all the movements of each submodule need to be converted to the inertial coordinate system for description; a rigid body motion node is set up at the leftmost node of each submodule, and a module coordinate system E is established at the rigid body motion node, and the rigid body motion node is used as the origin of the module coordinate system E and recorded as ; An elastic motion node is set up at the rightmost end, denoted as the right end node A; a cross-sectional coordinate system B is set up at any point on the internal axis of the submodule, with the origin located on the axis of the submodule, and the origin is denoted as point P, which is used to describe the position of any point D on the cross section.

3. The method for modeling helicopter blades with complex configurations based on multi-body dynamics according to claim 2 is characterized in that: The step 2 is specifically as follows: Step 2-1: Definition of coordinates and motion description of each submodule; First, define the coordinate variables of each node of the submodule, and set the generalized coordinates of the i-th submodule to be , represented by the position coordinates and Euler angle coordinates of the left end node and the deformation of the right end node: (1) In the formula, is the origin of the module coordinate system E of the i-th submodule Position coordinates expressed in inertial coordinate system I; is the origin of the module coordinate system E of the i-th submodule Euler angle coordinates expressed in inertial coordinate system I; is the elastic deformation freedom of the right end node of the ith submodule in the module coordinate system E; the superscript T represents the transposition operation of the vector; Then, define the position coordinates of any point D on the i-th submodule in the module coordinate system E, expressed as: (2) In the formula, is the position coordinate of point D in the module coordinate system E. The superscript E indicates that the position coordinate is measured in the module coordinate system E. The subscripts i and D indicate point D on the i-th submodule. x is the distance from point D to the origin of the coordinate system in the module coordinate system E. The length of the axis; The origin of the cross-sectional coordinate system B is along the axis in the module coordinate system E. The deformation in the direction, represents the deformation along the x-axis, represents the deformation along the y-axis, represents the deformation along the z-axis; is the coordinate of point D in the cross-sectional coordinate system B, represents the y-axis coordinate of the cross-sectional coordinate system B, represents the z-axis coordinate of the cross-sectional coordinate system B; is the transformation matrix from the section coordinate system B to the module coordinate system E after the submodule is deformed; According to equations (1) and (2), the position coordinates of any point D on the i-th submodule after deformation in the inertial coordinate system I are: (4) In the formula, is the position coordinate of point D in the inertial coordinate system I, is the transformation matrix from the module coordinate system E to the inertial coordinate system I; Step 2-2: derive the generalized inertia force of the i-th submodule; According to the Kane method, the generalized inertia force of the i-th submodule is expressed as: (5) in, are the first-order and second-order time derivatives of the position coordinates of any point D on the i-th submodule after deformation in the inertial coordinate system I, and are the velocity and acceleration expressed in the inertial coordinate system I; is the submodule material density; is the first-order time derivative of the generalized coordinates of the ith submodule, is the generalized rate; represents the differential of the cross-sectional area; represents the differential of the submodule axis length; Therefore, the velocity of any point D on the i-th submodule in the inertial coordinate system I is and acceleration : (6) In the formula, The origin of the module coordinate system E of the i-th submodule The first and second time derivatives of the position coordinates expressed in the inertial coordinate system I; They are the origin of the module coordinate system E Angular velocity and angular acceleration; The vectors are The antisymmetric matrix of ; According to the projection relationship between the time derivative of the Euler angle and the angular velocity in the multi-body dynamics theory, we have: (7) In the formula, They are The first and second time derivatives of ; are the projection matrix and its first-order time derivative, the projection matrix Specifically written as: (8) In the formula, is the Euler angle coordinate in formula (1) The first two elements in ; From equations (1), (6) and (7), we can deduce the velocity D of any point on the i-th submodule: For generalized rate The first derivative of : (9) In the formula, is the three-dimensional identity matrix; is the interpolation function matrix; According to equations (6) and (7), the acceleration of any point D on the i-th submodule in the inertial coordinate system I is Arranged into the following form: (10) Finally, substitute equations (9) and (10) into the expression (5) of the generalized inertial force to obtain the expression of the generalized inertial force of the i-th submodule, and then organize it to obtain: (11) In the formula, is the beam material density; are the generalized mass matrix and nonlinear term of the ith submodule, respectively, which are expanded as follows: (12) In the formula, the variable T in the upper right corner is the transposition symbol, and other variables have been defined in the above formulas; Step 2-3: Derive the generalized elastic force of the i-th submodule; The generalized elastic force of the ith submodule is obtained by taking the partial derivative of the strain energy of the submodule with respect to the generalized coordinates, and the strain energy is calculated by stress and strain; Step 2-4: Derive the generalized active force of the i-th submodule; According to the discrete format of the active force in Kane's method, the generalized active force term of the i-th submodule is derived: ; Discuss any point P on the axis of the submodule, the main force acting is , then the rth component of the generalized active force is It can be expressed as: (17) In the formula, is the rth eccentric velocity of point P in the inertial coordinate system I; is the number of generalized coordinates of the submodule; remember is the radius vector of point P in the inertial coordinate system I, then the eccentric velocity of point P is It is expressed as: (18) In the formula, represents the rth generalized coordinate of the i-th submodule; Substituting equation (18) into equation (17), we get the discrete formula of the rth component of the generalized active force: (19) Therefore, all The generalized force components are written in vector form, and the generalized active force of the i-th submodule is expressed as: (20) In the formula, represents the generalized active force column vector; represents the radius vector of any point P on the axis of the i-th submodule; Represents the column vector of generalized coordinates of the i-th submodule.

4. The method for modeling a helicopter blade with complex configuration based on multi-body dynamics according to claim 3 is characterized in that: In the step 2-1, The specific expression is as follows: (3) In the formula, It represents the Euler angle from the deformed section coordinate system B to the module coordinate system E. Specifically: It is the total rotation angle of the cross-sectional coordinate system B relative to the module coordinate system E in the x-axis direction, including the pitch angle of the blade , pre-twist angle and elastic torsion angle ,Right now ; It is the rotation angle of the cross-section coordinate system B relative to the module coordinate system E in the y-axis direction; It is the rotation angle of the cross-section coordinate system B relative to the module coordinate system E in the z-axis direction.

5. The method for modeling a helicopter blade with complex configuration based on multi-body dynamics according to claim 3 is characterized in that: The steps 2-3 are specifically as follows: First, the strain of the submodule is derived; the strain-displacement relationship of the blade is: (13) In the formula, is the axial strain; is the shear strain; The origin of the cross-sectional coordinate system B is along the axis in the module coordinate system E. The first-order arc-length derivative of the deformation in the direction with respect to the axial coordinate, with the superscript That is, it represents the first-order arc-length derivative; Indicates deformation The second arc-length derivative of ; Indicates deformation The second arc-length derivative of ; is the coordinate of any point D on the ith submodule in the cross-sectional coordinate system B, represents the y-axis coordinate of the cross-sectional coordinate system B, represents the z-axis coordinate of the cross-sectional coordinate system B; It represents the total rotation angle of the section coordinate system B relative to the module coordinate system E in the x-axis direction; It represents the elastic torsion angle of the section coordinate system B relative to the module coordinate system E in the x-axis direction. The first-order arc-length derivative with respect to the axial coordinate; Then, the stress of the submodule is derived; The blade is a slender structure, and the uniaxial stress assumption is used. The stress-strain relationship is expressed as: (14) In the formula, is the axial stress; is the shear stress; is Young's modulus; is the shear modulus; Furthermore, the strain energy of the submodule is calculated; According to the definition of strain energy, the strain energy of the i-th submodule is obtained using formulas (13) and (14): The expression is transformed into the form of cross-sectional load: (15) In the formula, They are the generalized force strain and generalized moment strain of the submodule respectively; are section force and section moment respectively; the matrix is the 6-dimensional section stiffness matrix; is the length of the submodule; represents the differential of the cross-sectional area; represents the differential of the submodule axis length; Finally, the generalized elastic force of the submodule is calculated; the strain energy is calculated For generalized coordinates The first-order derivative of is used to obtain the generalized elastic force of the i-th submodule and is written in matrix form: (16) In the formula, is the stiffness matrix of the ith submodule, and the numerical solution is calculated using the difference method; is a nonlinear term.

6. The method for modeling a helicopter blade with complex configuration based on multi-body dynamics according to claim 5 is characterized in that: The step 3 is specifically as follows: Step 3-1: First type constraint equation; The first type of constraint is introduced by the connection condition between the blade and the hub, that is, the connection between the first submodule of the blade and the hub; for helicopter blades, the connection between the blade and the hub is related to the rotor type. Currently, common rotor types include: bearingless rotor, hingeless rotor, and articulated rotor. The three types of constraint equations are expressed as: (21) In the formula, The constraint equations of three types of rotors are respectively represented: bearingless rotor, hingeless rotor, and articulated rotor; In expression (1) The element of, the i-th submodule represented in formula (1), and the first type of constraint equation is the connection between the first submodule of the blade and the hub, so the subscript i is written as 1; Step 3-2: Second type constraint equation; The second type of constraint is introduced by the connection between submodules. In step 1, the blade is divided into multiple submodules. There is no relative displacement or rotation between the submodules. The following constraint equation exists between the i-th submodule and the i+1-th submodule: (22) In the formula, They represent the position coordinates and Euler angle coordinates of the right end node A of the i-th submodule in the inertial coordinate system I respectively; Respectively represent the left end node of the i+1th submodule Position coordinates and Euler angle coordinates in inertial coordinate system I; Assuming that the blade is divided into N submodules, there are N-1 constraint equations such as equation (22), which can be written as: (23) In the formula, represents the second type of total constraint equation, whose elements are the subscript i of equation (22) rewritten as (1, 2, ..., N-1); the superscript T represents the transposition symbol; Step 3-3: The third type of constraint equation; The third type of constraint is introduced by the motion mode of the rotor. The rotor works at a constant speed, so the rotor speed is introduced into the blade dynamics model through the drive constraint. Suppose the rotor rotation angular velocity , time t, the driving constraint is written as: (24) In the formula, The third Euler angle coordinate of the first submodule given in equation (1) is expressed as follows, because the drive constraint only needs to apply the rotational angular velocity to the first submodule of the blade, that is, the rotational angular velocity is transmitted to all subsequent submodules through the second type constraint equation; Step 3-4: The fourth type of constraint equation; The fourth type of constraint is the constraint between the front and rear submodules of blades of different configurations at the turning point of the blade axis; assuming that the latter submodule that turns is the i+1th submodule, it is at a downward angle to the i-th submodule , the constraint equation is written as: (25) In the formula, The second element of the Euler angle coordinates of the i+1th submodule in equation (1); Step 3-5: Total constraint equation; Combining equations (21), (23), (24) and (25), we get the total constraint equation of the dynamic model: (26)。 7. The method for modeling a helicopter blade with complex configuration based on multi-body dynamics according to claim 6 is characterized in that: The step 4 is specifically as follows: According to equations (11), (16) and (20), the blade is divided into N sub-modules of generalized inertial force, generalized elastic force and generalized active force, and the groups are as follows: (27) in, represents the total generalized inertia force of N submodules, is the total generalized elastic force of N submodules, is the total generalized active force of N submodules; According to the Kane equation, the blade balance equation of N submodule groups is obtained: (28) Substituting equation (27) into equation (28), while considering the total constraint equation shown in equation (26), introducing Lagrange multipliers based on the first-kind Lagrange equation, we obtain a set of differential algebraic equations to extract the generalized acceleration and generalized coordinates The linear term coefficients are calculated and the blade dynamics model of N submodules is obtained: (29) In the formula, They are the total mass matrix, stiffness matrix and nonlinear terms of N submodules respectively; represents the constraint equation and is the generalized coordinate and a function of time t; is the constraint equation for generalized coordinates The transpose of the Jacobian matrix is ; is the Lagrange multiplier.

8. The method for modeling helicopter blades with complex configurations based on multi-body dynamics according to claim 7 is characterized in that: The step 5 is specifically as follows: Step 5-1: Use generalized- Method numerical discrete; The generalized-α method is used to numerically discretize the blade dynamics model shown in formula (29). Assuming that the blade dynamic response within a time period T needs to be solved, the integral time step is h, that is, the time T is divided into steps of h. time steps, then the nth time step The blade dynamics model is expressed as the residual In the form of: (30) In the formula, represents the nth time step, and the subscript n of other variables also represents the numerical results at the nth time step; Respectively represent the generalized coordinates, generalized velocity, and generalized acceleration of each submodule at the nth time step; represents the Lagrange multiplier of the nth time step; represents the nonlinear term at the nth time step; represents the constraint equation for the nth time step; According to the broad- The discrete form of the method, Discrete into: (31) Where h is the integration time step, that is, ; The subscript n-1 or n of each variable represents the value of the variable at the n-1th or nth time step; are algorithm parameters; is an auxiliary variable of the generalized-α method, and , Satisfaction relationship: (32) In the formula, Represents auxiliary variables The initial value of is the generalized acceleration The initial value of is the algorithm parameter, satisfying: (33) In the formula, is the spectral radius of the generalized-α method; Substituting equations (31) and (32) into equation (30), we can obtain the complete discrete formula of the blade dynamics model; Step 5-2: Use Newton's method to solve the linear equations discretized by the generalized-α method; Use generalized - After the numerical discretization of the method, a set of linear equations is obtained, and the variables in the nth time step are recorded as , in each time step, the Newton method is used to iteratively solve, and the iterative format is as follows: (34) In the formula, the superscript Represents Newton's method Iterations; is the residual of the kth iteration of the equation system at the nth time step; is the increment of the kth iteration; is the Jacobian matrix of the kth iteration of the system of equations at the nth time step; also, The specific expression is: (35) Find the increment value of the current iteration step After that, the unknown variables of the blade dynamics model are updated as follows: (36) In the formula, the subscript Indicates the variable at the nth time step; superscript The variable representing the k+1th iteration of the current time step; and are the process quantities to be solved, respectively: ; Substitute the variables in formula (36) into formula (34), satisfy , that is, the residual The 2-norm of is less than or equal to the set error limit , the iteration step is completed, and the increments of the generalized coordinates, generalized velocity, generalized acceleration and Lagrange multiplier at that moment are obtained ; So far, completed The task of solving the blade dynamics model at the moment; Step 5-3: Determine whether the solution of all discrete time steps is completed; Finally, determine whether all time steps in time period T have been solved. Is it true? If it is true, end the calculation; if not, let the variable value of the nth time step (Equation (36)) be the initial value of the n+1th time step, substitute it into equations (30), (31) and (32) in step 5-1, and update the time step at the same time, that is, let , and obtain the discrete formula of the blade dynamics model at the n+1th time step; then use the Newton method iterative formulas (34) and (35) in step 5-2 to solve the variable value of the kth iteration step of the n+1th step; finally, execute step 5-3 to determine whether the solution of all time steps is completed, and repeat this cycle until the entire solution task is completed.

Images