A modeling method for helicopter blades with complex configuration based on multi-body dynamics

Abstract

The present invention provides a method for modeling helicopter blades with complex configurations based on multi-body dynamics, which belongs to the field of helicopter blade modeling. First, according to the blade configuration characteristics, it is divided into several sub-modules, and the dynamic models of each sub-module are independent of each other. Secondly, according to the motion and deformation characteristics of the sub-modules, the dynamic models of each sub-module are established. Third, according to the multi-body dynamics idea, the displacement coordination boundary conditions between the sub-modules are introduced, and the constraint equations are established to enforce the continuity of the boundary node displacement. Fourth, based on the first-class Lagrangian equations, the dynamic equations and constraints of each sub-module are combined to assemble a blade dynamic model of a modular group. Finally, the generalized α method is used to numerically discretize the blade dynamic model of the modular group, and the Newton iteration method is used to solve the discretized linear equation group. The present invention utilizes the modular multi-body dynamics idea, which is of great significance for helicopter blade structure design and configuration analysis.

Description

Technical Field

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

[0002] Modern helicopter blades have evolved into complex configurations such as upward reversal, downward reversal, forward projection, and swept back, such as the Blue Edge blade (forward projection and swept back type) and the BERP blade (swept back type). The blade tips of these blades are all 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 blades of different types. Studies have shown that the swept back 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 result in larger errors. The present invention conducts relevant research on the problem of efficient and accurate unified modeling of blades of different configurations.

[0003] Among the current approaches for modeling helicopter blades of varying configurations, three-dimensional solid finite element modeling is one of the most accurate. This method can reveal more structural details of the blades based on the blade configuration through finite element meshing, allowing for modeling of blades of varying configurations. For example, Chinese invention patent [CN111339607B] provides an automated three-dimensional modeling method for D-shaped beam blades. This method utilizes ActiveX automation technology to rapidly model blades by invoking CATIACOM components via the VBWindows application platform. However, while this modeling method provides a highly accurate model, computational efficiency is significantly reduced. This is because establishing a precise three-dimensional solid model results in a large number of degrees of freedom in the dynamic model, which, when combined with aerodynamic model simulation, results in a significantly low solution rate. To address this computational efficiency issue, Chinese invention patent [CN116305589B] proposes a blade structural reduction analysis method based on a three-dimensional unit cell finite element mesh model. This method establishes an efficient, simplified, low-degree-of-freedom model of a helicopter blade, suitable for efficient simulation and stress analysis of the blade's structural response to stress and deformation. Although this method can reduce the degrees of freedom of the blade model to a certain extent and improve the simulation efficiency to a certain extent, it requires the re-establishment of the finite element mesh of the three-dimensional unit cell for blades of different configurations, which is a relatively complicated task. Moreover, this method is 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 property matrix, which limits its scope of application. Another type of helicopter blade modeling method is to approximate the blade as a slender beam and model the blade based on geometric nonlinear beam theory. This 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] The modeling of blades with different configurations, especially those with complex configurations, has various difficulties, 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 response to 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 to make the displacements between different sub-modules 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. The coefficient matrices between adjacent sub-modules are then 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 multibody dynamics includes the following steps:

[0008] 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. 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 point 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 motions 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 established at the rightmost end, which is recorded as the right end node A; a cross-sectional coordinate system B is established at any point on the internal axis of the sub-module, with the origin located on the sub-module axis, and the origin is recorded as point P, which is used to describe the position of any point D on the cross section.

[0010] Step 2: Based on 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] Where, is the origin of the module coordinate system E of the i-th submodule Position coordinates expressed in the inertial coordinate system I; is the origin of the module coordinate system E of the i-th submodule Euler angle coordinates expressed in the 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] Where, 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 in the module coordinate system E. The axis length; The origin of the cross-section 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] Where, 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-section 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 ; 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] Where, 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 formula (3) Replace 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 known as 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] Where, 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; are vectors 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] Where, 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] Where, is the Euler angle coordinate in formula (1) The first two elements in .

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

[0036] (9)

[0037] Where, 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 generalized inertial force expression (5) to obtain the generalized inertial force expression of the i-th submodule, and then organize it to obtain:

[0041] (11)

[0042] Where, is the beam material density; are the generalized mass matrix and nonlinear term of the i-th submodule respectively, which are expanded as follows:

[0043] (12)

[0044] In the formula, the variable superscript T 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 i-th submodule is obtained by taking the partial derivative of the submodule's strain energy with respect to the generalized coordinates. Therefore, the strain energy of the i-th submodule must be derived. According to the definition of strain energy, strain energy is calculated from stress and strain, so the stress and strain expressions for the submodule must be calculated. The following derivations are performed sequentially.

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

[0048] (13)

[0049] Where, is the axial strain; is the shear strain; The origin of the cross-section 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, superscript That is, it represents the first-order arc-length derivative; Indicates deformation The second-order arc-length derivative of ; Indicates deformation The second-order arc-length derivative of ; is the coordinate of any point D on the i-th 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 cross-section coordinate system B relative to the module coordinate system E in the x-axis direction; Indicates the elastic torsion angle of the cross-section coordinate system B relative to the module coordinate system E in the x-axis direction The 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] Where, 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 organized into the form of section load:

[0054] (15)

[0055] Where, are the generalized force strain and generalized moment strain of the submodule respectively; are the 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. Calculate the strain energy For generalized coordinates The first-order derivative of is used to obtain the generalized elastic force of the i-th submodule and write it in matrix form:

[0057] (16)

[0058] Where, is the stiffness matrix of the i-th 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 the Kane 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:

[0061] (17)

[0062] Where, is the rth partial 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 Expressed as:

[0064] (18)

[0065] Where, represents the r-th generalized coordinate of the i-th submodule.

[0066] Substituting Equation (18) into Equation (17), we obtain 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] Where, 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 generalized dynamic forces for a single submodule. However, the key to modular multibody modeling for various blade configurations is establishing constraint equations between the submodules, dynamically connecting them and ultimately obtaining complete dynamic equations. For blades, these constraint equations fall into the following four categories.

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

[0074] The first type of constraint is introduced by the connection 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 include: bearingless rotor, hingeless rotor, and articulated rotor. The three types of constraint equations are expressed as:

[0075] (twenty one)

[0076] Where, The constraint equations for three types of rotors are respectively expressed as 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.

[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 or rotation between the submodules. According to multi-body dynamics theory, the connection needs to be established through constraint equations. Therefore, the following constraint equation exists between the i-th submodule and the i+1-th submodule:

[0079] (twenty two)

[0080] Where, 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; Represents 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] Where, represents the second type of total constraint equation, whose elements are the subscript i in equation (22) rewritten as (1,2,…,N-1); the superscript T represents the transpose 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. Let the rotor rotation angular velocity be , time t, the driving constraint is written as:

[0086] (twenty four)

[0087] Where, 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, 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 degrees, and this angle is introduced as a constraint. 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] Where, It represents the second element of the Euler angle coordinates of the i+1th submodule in equation (1).

[0092] Steps 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 Lagrangian equation, introduce the Lagrangian multiplier, combine the dynamic equations and constraint equations of each submodule, and form a modular blade dynamic model. Specifically:

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

[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 submodules.

[0099] According to the Kane equation, the blade balance equation of the 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] Where, 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 the constraint equation (26). Therefore, the coefficient matrix of the blade dynamics model formula (29) is obtained: and vector Different submodules are decoupled, so adding or deleting submodules will not affect the coefficients of existing submodules, 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 discretization;

[0107] The generalized -α method is used to numerically discretize the blade dynamic model shown in formula (29). Assuming that the blade dynamic response within the 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] Where, Indicates the nth time step, and the subscript n of other variables also indicates 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 at 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] Where, Represents auxiliary variables The initial value of is the generalized acceleration The initial value of Are algorithm parameters that satisfy:

[0115] (33)

[0116] Where, 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 of 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 k-th iteration of the system of equations at the n-th time step.

[0122] also, The specific expression is:

[0123] (35)

[0124] In finding 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 Lagrangian 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 the time period T have been solved, that is, determine 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 , obtain the discrete formula of the blade dynamic 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 the 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 the modeling of blades with different complex configurations by adjusting the constraint parameters.

[0133] (2) The present invention uses the implicit integration algorithm 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 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 protruding and swept blade tip 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 diagram of the submodule.

[0140] Figure 7 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 by the present invention.

[0142] Figure 9 Schematic diagram of the assembly of different sub-module elements in the nonlinear term of the blade dynamics model system established for the present invention.

[0143] Figure 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 is subjected to a concentrated external force F in Example 1 of the present invention.

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

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

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

[0147] Figure 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 is subjected to a concentrated external force F in Example 1 of the present invention.

[0148] Figure 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 is subjected to a concentrated external force F in Example 1 of the present invention.

[0149] Figure 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] Figure 17 This is a structural side view of the blade with a swept-tip configuration in Example 2 of the present invention.

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

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

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

[0154] The following embodiments of the present invention are described in further detail with reference to 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. Figures 2 to 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 submodule structure, and different modules are connected using Figure 7 The submodules 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 solutions and advantages of the present invention more clear, the following two specific embodiments are combined with the attached Figures 5-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 in FIG, this embodiment is a rectangular straight blade subjected to a dynamic load F at the blade tip, where 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 with no turning section along the blade axis. It is divided into 4 submodules, such as Figure 5 At the same time, Figure 6 , establish an inertial coordinate system I, and all the motions of each submodule are converted 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 established at the rightmost end; a cross-sectional coordinate system B is established 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: Based on 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 i-th submodule in the module coordinate system E before and after deformation, as shown in formula (2); finally, based on formulas (1) and (2), define the position coordinates of any point D on the i-th 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 i-th 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 i-th submodule, then calculate the strain energy, and finally obtain the generalized inertial force of the i-th 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 i-th 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 the Kane method, the generalized active force term of the i-th submodule is derived Use Equation (18) to calculate the eccentric velocity of point P, and substitute it into Equation (17) to obtain the generalized active force of the i-th 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 multibody modeling is establishing constraint equations between submodules, achieving dynamic connections between them and ultimately obtaining complete dynamic equations. For propeller blades, constraint equations fall into 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 blades 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 equation (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 in the figure, the connection method between different submodules is given. According to the multi-body dynamics theory, the connection needs to be established through constraint equations. Therefore, the constraint equations between the i-th submodule and the i+1-th submodule are established using formula (22), and the final set is in the form of formula (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 blade axis. Since this embodiment uses a straight blade, there is no axial turning of the blade, so the constraints in this type of constraint are .

[0181] Steps 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 Lagrangian equation, introduce the Lagrangian multiplier, combine the dynamic equations and constraint equations of each submodule, and form a modular blade dynamic model. Specifically:

[0184] From Equations (11), (16), and (20), we consider dividing the blade into N submodules of generalized inertial force, generalized elastic force, and generalized active force, and group them to obtain the generalized force of Equation (27); then, according to the Kane equation, we can obtain the blade balance equation (28) of the N submodules; further, considering the constraint equation (26), we obtain a set of differential algebraic equations, and organize them to obtain the blade dynamics model of the N submodules, as shown in Equation (29).

[0185] According to formula (27), each submodule is completely independent and is only connected by the constraint equation (26). Therefore, the coefficient matrix of the blade dynamics model formula (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 Figure 9 As shown, it is more convenient to model blades with 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 discretization;

[0188] The generalized-α method is used to numerically discretize the blade dynamic model shown in formula (29). In this embodiment, the simulation time 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 the numerical discretization of the method, a set of linear equations is obtained. In each time step, the Jacobian matrix of formula (35) is first calculated using the difference method, and then the Newton method of formula (34) is used to solve the equations. Finally, it is judged 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 Lagrangian multiplier at that moment are obtained. At this point, 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 , obtain the discrete formula of the blade dynamic 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 Figures 10 to 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 Ω = 70 rad / s, and the external force F = 50sin (20t) acts on the blade tip and changes along the z-axis. Figures 10 to 12 is the linear displacement of the blade tip in the axial, shimmy and flapping directions, Figures 13 to 15 The angular displacements of the blade tip in the axial, shimmy, and flapping directions are shown in the figure. As can be seen from the figure, the calculation results of the modeling method of the present invention are basically consistent with those of the multi-body dynamics software DYMORE. Under high-speed rotation, the changes in the axial deformation curve and the torsion angle curve show a clear bending-stretching-torsion coupling phenomenon, which is also the typical dynamic characteristic of helicopter blades in operation. Therefore, the ability to accurately describe this phenomenon demonstrates 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 complex configuration swept blades and 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 with embodiment 1. The first is the number of sub-modules divided according to the blade configuration in step 1. For the tip-swept blade, 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 as 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 adjusting the value of , the dynamic models 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 Figure 16 and Figure 17 The following are the top and side views of the blade with swept tip used in this embodiment. The length of the straight blade portion is L1 = 0.8636m, the length of the swept portion is L2 = 0.1524m, the width is b1 = 0.0254m, and the height is h1 = 0.0016m. In order to better compare the effectiveness of the modeling method of the present invention, the sweep angle is Two working conditions, 15° and 45°, were selected for verification.

[0198] like Figure 18 The figure shows the comparison between the displacement of the free end of the blade and the calculation result of Ansys software. The calculation condition is the sweep angle. , the speed gradually increases from 0 to 1 second. 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. Figure 19 and Figure 20As shown in the figure, the first six natural frequency curves of the two sweep angles of 15° and 45° at 500r / min and 750r / min speeds are respectively given, 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 all increased to a certain extent. Because the increase in speed, the centrifugal force of the blade will also increase, thereby increasing the natural frequency. 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-back configuration blades.

[0199] In summary, by performing effective numerical simulations on the above two typical embodiments and comparing them with the results of Ansys software or experimental results, the feasibility and development necessity of the method for modeling helicopter blades with complex configurations based on multi-body dynamics proposed by the present invention can be well illustrated. 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 to a blade dynamics model of any connection method. Therefore, the method for modeling helicopter blades with complex configurations 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 intended to express specific implementation methods of the present invention, but should not be understood as limiting the patent scope of the present invention. It should be pointed out that for those skilled in the art, several variations and improvements can be made without departing from the concept of the present invention, and these all fall within the scope of protection of the present invention.

Claims

1. A method for modeling helicopter blades with complex configurations based on multi-body dynamics, characterized in that: The following steps are involved: Step 1: Divide the blade into several submodules based on 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: Based on 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, introduce the displacement coordination boundary conditions between sub-modules, establish constraint equations, and enforce the continuity of boundary node displacements; specifically: 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: Where, Φ 11 ,Φ 12 ,Φ 13 The constraint equations for three types of rotors are respectively expressed as 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: Where, 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; Represents the left end node O of the i+1th submodule respectively E 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: Where Φ2 represents the second type of total constraint equation, whose elements are obtained by rewriting the subscript i of Equation (22) 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 of the rotor. The rotor rotates at a steady speed, so the rotor speed is introduced into the blade dynamics model through the drive constraint. Assuming the rotor rotation angular velocity Ω and time t, the drive constraint is written as: Where, Express the third Euler angle coordinate of the first submodule given in equation (1), 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 transferred 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 last submodule after the turning point is the i+1th submodule, it is at a downward angle α0 with the i-th submodule. The constraint equation is written as: Where, The second element of the Euler angle coordinates of the i+1th submodule in equation (1); Steps 3-5: Total constraint equation; Combining equations (21), (23), (24), and (25), we can obtain the total constraint equation of the dynamic model: Step 4: Based on the first-kind Lagrangian equation, Lagrangian multipliers are introduced to combine the dynamic equations and constraint equations of each submodule to form a modular blade dynamic model. Step 5: Use the generalized-α method to numerically discretize the blade dynamics model of the modular group and use Newton's method to solve the discretized linear equations; specifically: Step 5-1: numerical discretization using the generalized-α method; The generalized -α method is used to numerically discretize the blade dynamic model shown in formula (29). Assuming that the blade dynamic response within the 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 t n =nh The blade dynamics model is expressed as the residual g n In the form of: Where, t n Indicates the nth time step, and the subscript n of other variables also indicates the numerical results at the nth time step; p n , Respectively represent the generalized coordinates, generalized velocity, and generalized acceleration of each submodule at the nth time step; λ n represents the Lagrange multiplier of the nth time step; (F non ) n represents the nonlinear term at the nth time step; Φ(p n ,t n ) represents the constraint equation of the nth time step; According to the discretization format of the generalized-α method, p n , Discrete into: Where h is the integration time step, that is, t n =t n-1 +h; 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; a n is an auxiliary variable of the generalized-α method, and a n-1 ,a n Satisfaction relationship: In the formula, a0 represents the auxiliary variable a n-1 The initial value of is the generalized acceleration The initial value of α m ,α f ,∈,γ is the algorithm parameter, satisfying: Where μ∈[0,1] 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; After numerical discretization using the generalized-α 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 the problem. The iterative format is as follows: Where, the superscript k represents the kth iteration of Newton's method; 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: In finding the increment value of the current iteration step After that, the unknown variables of the blade dynamics model are updated as follows: Where, the subscript n represents the variable of the nth time step; the superscript k+1 represents the variable of 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), g n 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 Lagrangian multiplier at that moment are obtained. At this point, complete t n 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 the time period T have been solved, that is, determine Is it true? If so, 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 t n+1 =t n +h, and obtain the discrete formula of the blade dynamic 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.

2. The method for modeling helicopter blades with complex configurations based on multi-body dynamics according to claim 1, 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 motions of each submodule need to be converted to this inertial coordinate system for description; a rigid body motion node is established 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 O E ; An elastic motion node is established at the rightmost end, which is recorded as the right end node A; a cross-sectional coordinate system B is established at any point on the internal axis of the sub-module, with the origin located on the sub-module axis, and the origin is recorded 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, 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 p i , represented by the position coordinates and Euler angle coordinates of the left end node and the deformation of the right end node: Where, is the origin O of the module coordinate system E of the i-th submodule E Position coordinates expressed in the inertial coordinate system I; is the origin O of the module coordinate system E of the i-th submodule E Euler angle coordinates expressed in the inertial coordinate system I; q 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; Then, define the position coordinates of any point D on the i-th submodule in the module coordinate system E, expressed as: Where, 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 between point D and the coordinate origin O in the module coordinate system E. E The axis length; (u, v, w) is the deformation of the origin of the cross-sectional coordinate system B in the module coordinate system E along the axis (x, y, z), u represents the deformation along the x-axis, v represents the deformation along the y-axis, and w represents the deformation along the z-axis; (0, η, ζ) 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, and ζ represents the z-axis coordinate of the cross-sectional coordinate system B; T EB is the transformation matrix from the cross-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: Where, is the position coordinate of point D in the inertial coordinate system I, T IE 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: 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, are the velocity and acceleration expressed in the inertial coordinate system I; ρ is the material density of the submodule; is the first-order time derivative of the generalized coordinates of the i-th submodule, is the generalized rate; dA represents the differential of the cross-sectional area; dx 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 Where, The origin O of the module coordinate system E of the i-th submodule is E The first and second time derivatives of the position coordinates expressed in the inertial coordinate system I; ω, They are the module coordinate system E origin O E Angular velocity and angular acceleration; are vectors ω, 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: Where, They are The first and second time derivatives of G, They are the projection matrix and its first-order time derivative respectively. The projection matrix G is specifically written as: Where, is the Euler angle coordinate in formula (1) The first two elements in ; From equations (1), (6) and (7), we can derive the velocity D of any point on the i-th submodule: For generalized rate The first derivative of : Where I3 is the three-dimensional unit matrix; H q 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: Finally, substitute equations (9) and (10) into the generalized inertial force expression (5) to obtain the generalized inertial force expression of the i-th submodule, and then organize it to obtain: Where ρ is the density of the beam material; M i , are the generalized mass matrix and nonlinear term of the i-th submodule respectively, which are expanded as follows: In the formula, the variable superscript T is the transposition symbol, and the 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 i-th 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 obtained by calculating the 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 the Kane method, the generalized active force term F of the i-th submodule is derived f,i ; Discuss any point P on the axis of the submodule, the active force acting is f, then the rth component of the generalized active force is It can be expressed as: Where, v (r) is the rth partial velocity of point P in the inertial coordinate system I; a 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 v of point P is (r) Expressed as: Where, represents the rth generalized coordinate of the i-th submodule; Substituting Equation (18) into Equation (17), we obtain the discrete formula of the rth component of the generalized active force: Therefore, if all a generalized force components are written in vector form, the generalized active force of the i-th submodule can be expressed as: Where, F f,i represents the generalized active force column vector; represents the radius vector of any point P on the axis of the i-th submodule; p i Represents the column vector of generalized coordinates of the i-th submodule.

4. The method for modeling helicopter blades with complex configurations based on multi-body dynamics according to claim 3, characterized in that: In step 2-1, T EB The specific expression is as follows: Where, Represents the Euler angle from the deformed section coordinate system B to the module coordinate system E. Specifically: θ1 is the total rotation angle of the section coordinate system B relative to the module coordinate system E in the x-axis direction, including the pitch angle θ of the blade C , pre-twist angle θ tw and elastic torsion angle φ, i.e. θ1=θ C +θ tw +φ; β 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 helicopter blades with complex configurations based on multi-body dynamics according to claim 3, 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: Where, ε 11 is the axial strain; ε 12 ,ε 13 is the shear strain; (u ′ ,v′,w′) are the first-order arc-length derivatives of the deformation of the origin of the cross-sectional coordinate system B along the axis (x, y, z) in the module coordinate system E with respect to the axial coordinate, and the superscript ′ indicates the first-order arc-length derivative; v″ indicates the second-order arc-length derivative of the deformation v; w″ indicates the second-order arc-length derivative of the deformation w; (0,η,ζ) are the coordinates of any point D on the i-th submodule in the cross-sectional coordinate system B, η indicates the y-axis coordinate of the cross-sectional coordinate system B, ζ indicates the z-axis coordinate of the cross-sectional coordinate system B; θ1 indicates the total rotation angle of the cross-sectional coordinate system B relative to the module coordinate system E in the x-axis direction; φ′ represents the first-order arc length derivative of the elastic torsion angle φ in the x-axis direction of the cross-section coordinate system B relative to the module coordinate system E 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: Where σ 11 is the axial stress; σ 12 ,σ 13 is the shear stress; Q 11 is Young's modulus; Q 55 ,Q 66 is the shear modulus; Furthermore, the strain energy of the submodule is calculated; According to the definition of strain energy, the strain energy U of the i-th submodule is obtained using formulas (13) and (14): i The expression is organized into the form of section load: Where, γ,κ are the generalized force strain and generalized moment strain of the submodule respectively; F,M are the section force and section moment respectively; the matrix S is the 6-dimensional section stiffness matrix; l is the length of the submodule; dA represents the differential of the cross-sectional area; dx represents the differential of the submodule axis length; Finally, the generalized elastic force of the submodule is calculated; the strain energy U is calculated i For the generalized coordinate p i The first-order derivative of is used to obtain the generalized elastic force of the i-th submodule and write it in matrix form: Where K i is the stiffness matrix of the ith submodule, and the numerical solution is calculated using the difference method; F e,i is a nonlinear term.

6. The method for modeling helicopter blades with complex configurations based on multi-body dynamics according to claim 1, characterized in that: The step 4 is specifically as follows: According to equations (11), (16) and (20), the blade is divided into N submodules of generalized inertial force, generalized elastic force and generalized active force, and the groups are: Among them, F T Represents the total generalized inertial force of N submodules, F e is the total generalized elastic force of N submodules, F f is the total generalized active force of N submodules; According to the Kane equation, the blade balance equation of the N submodule groups is obtained: F T +F e +F f =0 (28) Substitute equation (27) into equation (28), and consider the total constraint equation shown in equation (26). Based on the first kind of Lagrangian equation, we introduce the Lagrangian multiplier to obtain a set of differential algebraic equations to extract the generalized acceleration And the linear term coefficient of the generalized coordinate p, and sort out the blade dynamics model of N submodules: Where M, K, F non are the total mass matrix, stiffness matrix and nonlinear terms of the N submodules respectively; Φ(p,t) represents the constraint equation and is a function of the generalized coordinate p and time t; is the transpose of the Jacobian matrix of the constraint equation with respect to the generalized coordinate p, that is, λ is the Lagrange multiplier.

Images