Helicopter rotor system dynamics simulation efficient solving method based on three-dimensional entity unit

Abstract

The invention discloses a helicopter rotor system dynamics simulation efficient solving method based on a three-dimensional entity unit, and belongs to the technical field of helicopter rotor system dynamics simulation analysis. Comprising the following steps: dividing a helicopter rotor geometric model into three-dimensional entity finite element models, and endowing unit attributes; and according to the finite element model, dividing co-rotation analysis components, and realizing effective butt joint among the components. Establishing a dynamic differential algebraic equation of the rotor system, and discretizing the dynamic differential algebraic equation; two key items in the efficient solving process are calculated, calculation of iteration increments is completed, and the dynamic time-history response of the rotor system is obtained. According to the thought of a component-level co-rotation coordinate method, extracting the elastic deformation of the rotor component, completing the stress tensor calculation of the rotor structure, and outputting a stress nephogram. According to the method, each iteration step linearized equation set is solved by adopting a Schel complement method, so that the solving efficiency of the large-scale linear equation set is improved; on the premise of keeping high precision, the dynamic simulation analysis of the finite element model of the large-scale rotor system can be efficiently completed.

Description

Technical Field

[0001] The invention belongs to the technical field of helicopter rotor system dynamics simulation analysis, and relates to an efficient solution method for helicopter rotor system dynamics simulation based on a three-dimensional solid finite element component-level co-rotating coordinate method (C-CR). Background Art

[0002] As an aircraft capable of vertical take-off and landing and hovering, helicopters are widely used in military, rescue, transportation, exploration and other fields. One of the core components of a helicopter is the rotor system, which is a key component that provides lift and thrust. The performance of the rotor directly affects the flight stability, maneuverability and safety of the helicopter. Therefore, it is of great significance to conduct in-depth research on the structural dynamics of the rotor system.

[0003] The rotor system is mainly composed of rotor blades, rotor hub, rotor propeller and rotor drive system. During flight, the rotor system is subjected to complex dynamic loads, such as centrifugal force, inertial force and elastic deformation. The interaction of these loads makes the dynamic characteristics of the rotor system very complex. In order to ensure the safety and reliability of the helicopter, the structural dynamics of the main rotor system must be accurately simulated and analyzed. Traditional analysis methods mainly include the following aspects:

[0004] (1) Multibody dynamics analysis: Simulate the relative motion and interaction of the components in the rotor system and analyze its overall dynamic characteristics.

[0005] (2) Finite element analysis: Detailed modeling of the main rotor system is performed to analyze its stress-strain distribution and fatigue life under complex loads.

[0006] (3) Experimental testing: The accuracy and reliability of the simulation model are verified through physical experiments, but this method is usually time-consuming and costly.

[0007] With the development of computer technology and simulation technology, the use of computer simulation technology for structural dynamics analysis of the main rotor system has become a key research direction. Simulation technology can predict and optimize the dynamic characteristics of the main rotor system during the design stage, reduce the number and cost of physical tests, and improve design efficiency and product performance. Although the current simulation analysis method has solved the problem of rotor system dynamics analysis to a certain extent, it still faces many challenges in dealing with complex nonlinear dynamic responses and high-precision numerical simulations. For example, the rotor blades will generate complex centrifugal forces and inertial forces during rotation, resulting in nonlinear responses of the structure, which puts higher requirements on simulation analysis.

[0008] In the field of helicopter rotor system simulation and analysis, beam elements are usually used to simulate blade structures. However, with the continuous updating and iteration of helicopter configurations, more blades with advanced configurations and complex geometric shapes have emerged. Beam elements are difficult to meet the needs of accurate analysis when facing these complex helicopter rotor structures, especially in the modeling and analysis of auxiliary components such as hubs and bearings. Therefore, the development of an efficient and accurate rotor system structural dynamics simulation and analysis method has important theoretical significance and engineering application value.

[0009] In order to achieve high-fidelity dynamic analysis of helicopter rotor systems, the co-rotating coordinate method (Full-CR) based on three-dimensional solid elements is the best way to solve this problem. Three-dimensional solid elements can meet the full-scale modeling analysis of any complex rotor geometry, and its complete post-processing analysis technology can obtain the stress and strain distribution state of the helicopter rotor structure. This is of great significance to the design and optimization of advanced airfoils. However, the main problems faced by using three-dimensional solid elements for calculations are:

[0010] (1) A large number of units: In order to accurately describe complex geometric shapes and capture subtle structural features, the model needs to be divided into a large number of three-dimensional solid units, which leads to a sharp increase in the size of the model.

[0011] (2) High consumption of computing resources: Due to the large number of units, the storage and memory requirements required for calculations increase significantly, placing high demands on computing resources.

[0012] (3) Long calculation time: The huge calculation scale leads to a significant increase in calculation time, especially in nonlinear dynamic analysis, where the iterative solution process is complex and time-consuming.

[0013] (4) Numerical stability problems: Numerical stability problems are prone to occur during the solution process, especially when dealing with nonlinear and complex contact problems. Summary of the invention

[0014] To solve the above problems, the present invention proposes an efficient simulation analysis algorithm for rotor system dynamics based on the three-dimensional solid unit component-level co-rotating coordinate method (C-CR). The method uses three-dimensional solid units to perform finite element meshing on the rotor system geometry model, and divides the finite element model into multiple co-rotating analysis components. The units in each component share the same co-rotating coordinate system, and combine the multi-body concept to achieve effective docking between the various components of the system through constraint equations of multiple constraints. In the nonlinear iterative solution process, the linearized equation group of each iteration step is solved by the Schur complement method, and an efficient calculation method for key items in the solution process is derived to improve the solution efficiency of large-scale linear equation groups. Through this method, the dynamic simulation analysis of the finite element model of a large-scale rotor system can be efficiently completed under the premise of ensuring high precision. The method of the present invention can not only accurately simulate the complex geometric shape and dynamic response of the rotor system, but also effectively solve the shortcomings of the traditional method in terms of calculation scale and efficiency, and provide strong support for the design and optimization of complex rotor systems. Thereby further improving the design level and operating performance of the main rotor system, and promoting the innovation and development of helicopter rotor systems.

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

[0016] A method for efficiently solving the dynamics simulation of a helicopter rotor system based on three-dimensional solid units, wherein the method first discretizes the helicopter rotor geometry model into a three-dimensional solid finite element model through finite element meshing. Secondly, the co-rotation analysis components (hereinafter referred to as components) of the three-dimensional solid finite element model are divided. Thirdly, time-invariant constraints or time-varying constraints are used to achieve effective connection of the various components. Fourthly, the dynamics equations of the rotor system are established. Fifthly, the dynamics equations of the helicopter rotor system are efficiently solved. Sixthly, two key items in the solution process of the linearized equation group are calculated. The innovation points of the present invention focus on the second and sixth steps. The second step proposes a new component analysis method based on three-dimensional solid finite elements. The sixth step proposes an efficient solution calculation method. Specifically, the following steps are included:

[0017] The first step is to discretize the helicopter rotor geometry model into a three-dimensional solid finite element model through finite element meshing.

[0018] The helicopter rotor system usually includes structures such as blades, hub arms, and rotating shafts. First, the helicopter rotor system geometric model is meshed using three-dimensional solid units (such as eight-node hexahedron units, twenty-node hexahedron units, etc.) to obtain a finite element model of the rotor system. Secondly, a model file such as the model node coordinates and unit node numbers of the finite element model is generated. Finally, the unit material properties are assigned to the corresponding units.

[0019] Step 2: Divide the co-rotation analysis components of the 3D solid finite element model

[0020] In the first step, the helicopter rotor system geometric model is discretized into a three-dimensional solid finite element model through finite element meshing, and then the three-dimensional solid finite element model is divided into several components, as shown in the attached figure. Figure 1 As shown. In the figure, the units in each component share the same coordinate system, which is called the component co-rotation coordinate system. The component co-rotation coordinate system is generated by sequentially finding three nodes that are far apart and not on the same straight line on each component, which are called component representative nodes, and quickly constructing the corresponding component co-rotation coordinate system through the edge parallel method.

[0021] The construction method of the component co-rotation coordinate system basis vector is as follows:

[0022]

[0023] Wherein, e1 represents the unit vector formed by the second representative node and the first representative node on the component; e3 represents the unit vector formed by the third and first representative nodes on the component and the cross product of the e1 vector; e2 represents the unit vector obtained by the cross product of the unit vector e3 and e1; x1 represents the coordinate column vector of the first representative node; x2 represents the coordinate column vector of the second representative node; x3 represents the coordinate column vector of the third representative node; ||…|| represents the modulus of the calculated vector;

[0024] The rotation matrix of the component can be expressed as The establishment of the component co-rotational coordinate system eliminates the solution process of each unit co-rotational coordinate system in the traditional unit-level co-rotational coordinate (Full CR) method, which greatly reduces the program calculation and solution time.

[0025] Step 3: Use time-invariant constraints or time-varying constraints to achieve effective connection of each component

[0026] In the second step, after the helicopter rotor system is divided into several components, the multi-body system dynamics concept is adopted to achieve the connection between the components through constraints. The constraints include time-invariant constraints and time-varying constraints.

[0027] 3.1 Constraint equations and constraint Jacobian matrix of time-invariant constraints

[0028] The time-invariant constraints can achieve effective docking between components. The finite element nodes between two adjacent components are coupled by coupling node degrees of freedom, as shown in the attached figure. Figure 2As shown. The method of coupling nodal degrees of freedom is achieved by limiting the relative coordinates of each pair of coincident nodes between two components in three dimensions of space. The constraint equation generated by the method of coupling nodal degrees of freedom is the same as the constraint equation of the spherical joint constraint in multibody dynamics, and the constraint equation can be expressed as Φ0=0. The subscript 0 represents the term related to the time-invariant constraint. To further obtain the constraint Jacobian matrix of the constraint equation, let the node coordinates of the component Generalized coordinates as components The constraint Jacobian matrix Φ of the spherical joint constraint can be calculated by differentiating the constraint equations with respect to generalized coordinates. q_0 Φ q_0 It is a sparse constant matrix that presents a value of 1 or -1 only at the generalized coordinate position related to a certain component, and all other positions are 0. Therefore, the constraint Jacobian matrices of all time-invariant constraints can be pre-calculated and stored, and then these matrices can be directly back-substituted in each iterative solution process of the simulation to improve the computational efficiency.

[0029] 3.2 Constraint Equations and Constraint Jacobian Matrix of Time-Varying Constraints

[0030] The time-varying constraints include angle-driven constraints, vertical constraints, and distance constraints. The combination and superposition of these constraints can achieve effective docking between the helicopter hub structure and components. The angle-driven constraint is used to control the rotational motion of the three-dimensional solid finite element model. The vertical constraint and distance constraint are used to constrain the vertical angle and relative distance of the components. In order to implement the three time-varying constraints, external auxiliary nodes need to be introduced into the three-dimensional solid finite element model. The purpose of introducing external auxiliary nodes is to combine the auxiliary nodes with the finite element nodes in the three-dimensional solid finite element model to jointly complete the writing of the time-varying constraint equation Φ1. As shown in the attached figure Figure 4 As shown, the generalized coordinates of the auxiliary node are denoted by q aux . By q aux and the generalized coordinates of the components described in Section 3.1 Vectors A, B, C, and D are generated between the auxiliary nodes and the nodes of the three-dimensional solid finite element model of the component. Vector A represents the arm in the hub structure, vector C is the helicopter shaft, and vectors B and D are two fixed vectors on the helicopter body, which are used to fix the entire rotor system. These vectors are used to write the constraint equation Φ1, where the angle drive constraint is achieved by controlling the change of the angle θ between vectors A and B. The constraint equation Φ1 of the time-varying constraint is expressed as follows:

[0031] Φ1 A =A·B-|A||B|cosθ (2)

[0032] Φ1 V =C·D (3)

[0033] Φ1 D =|A|-d (4)

[0034] Among them, the subscript 1 represents the time-varying constraint; Φ1 A represents the constraint equation of the angle drive constraint; θ represents the drive angle; Φ1 V represents the vertical constraint equation; Φ1 D Represents the distance constraint equation; d represents the control distance of the distance constraint.

[0035] The constraint equation Φ1 of the time-varying constraint is transformed into the generalized coordinate q aux and The derivative can obtain the constraint Jacobian matrix Φ of the time-varying constraint q_1 The constraint Jacobian matrix of the time-varying constraint Φ q_1 It is highly sparse.

[0036] Step 4: Establish the dynamic equation of the rotor system

[0037] After the above three steps, the generalized coordinates q of the rotor system are the combination of the generalized coordinates of all three-dimensional solid finite element nodes and the generalized coordinates of the introduced auxiliary nodes:

[0038]

[0039] in, represents the generalized coordinate column vector of the 3D solid finite element nodes in the first component; represents the generalized coordinate column vector of the 3D solid finite element nodes in the second component; Indicates the Nth b The generalized coordinate column vector of the 3D solid finite element nodes within each component; Represents the column vector of generalized coordinates of auxiliary nodes.

[0040] Using the first kind of Lagrangian equations and considering all constraint equations (2-4), the dynamic control equations of the helicopter rotor system can be written as the following index 3 differential algebraic equations (DAEs):

[0041]

[0042] in, represents the generalized acceleration corresponding to all generalized coordinates in the helicopter rotor system; Φ(q,t) represents the total constraint equation of the helicopter rotor system, which includes the time-invariant constraint equation Φ0 in Section 3.1 and the time-varying constraint equation Φ1 in Section 3.2; M is the mass matrix of the helicopter rotor system; Φ q =[Φ q_0 ; Φ q_1 ] is the system constraint Jacobian matrix; λ is the Lagrange multiplier vector reflecting the constraint force physically; fINT is the internal force column vector of all nodes, which is calculated from the elastic internal force of the component; EXT It is the external force of the helicopter rotor system, including gravity and aerodynamic force.

[0043] Furthermore, the internal force vector f INT And the mass matrix M of the helicopter rotor system can be expressed as:

[0044]

[0045] in, represents the mass matrix of the first component; represents the mass matrix of the second component; Indicates the Nth b The mass matrix of a component; M aux Represents the auxiliary node mass matrix.

[0046] Step 5: Efficiently solve the helicopter rotor system dynamics equations

[0047] 5.1 Calculation of the Jacobian Matrix of the Helicopter Rotor System

[0048] The rotor system dynamics equations (6) in the fourth step describe a set of differential algebraic equations with highly nonlinear characteristics and typical rigid properties. The typical implicit Newmark direct integration method is used to discretize the system dynamics equations into a set of nonlinear algebraic equations at each integral / time step. If the solution at the nth time step is obtained, the differential algebraic equations are discretized and the Newmark hypothesis is introduced in the n+1th step to produce the following nonlinear equations:

[0049]

[0050] Where h is the integration step; α and δ represent two parameters in the Newmark algorithm; t n+1 represents the n+1th moment after time discretization; represents the acceleration column vector at time t=n+1; q n+1 Represents the generalized coordinate column vector at time t=n+1; Represents the velocity column vector at time t=n+1.

[0051] After derivation, we get is a nonlinear equation system with unknown variables. The linear equation system in the iteration step obtained by the Newton-Raphson method is:

[0052] JΔX=G RES (10)

[0053] Where J is the Jacobian matrix; is the iteration increment; G RES =[F T ,Ψ T ] T is the residual of the nonlinear equation;

[0054] Furthermore, the Jacobian matrix J is:

[0055]

[0056] The Jacobian matrix can be abbreviated as:

[0057]

[0058] 5.2 Fast Solution of Linear Equations

[0059] In Section 5.1, we obtained the Jacobian matrix J of the linear equation system. By observing formula (12), we can find that the C matrix is ​​a large-scale block diagonal matrix, expressed as

[0060]

[0061] in,

[0062]

[0063] in is the rotation matrix of the component described in Section 2.

[0064] To achieve a fast solution to the linear equation system (10), the Schur complement method can be used. First, the linear equation system (10) is divided into blocks. After the blocks are divided, the equation system (10) can be rewritten as:

[0065]

[0066] The first equation of formula (15) can be solved as:

[0067]

[0068] Combining formula (16) with the second equation of formula (15), we can solve the Lagrangian multiplier increment of the nonlinear algebraic equation:

[0069]

[0070] Finally, substituting formula (17) back into formula (16) yields the acceleration increment: The final result. This term is the Schur complement of the submatrix C in the matrix J. Therefore, the problem of solving a large-scale linearized system of equations is transformed into solving in each iteration step by the Schur complement method. and C-1 F these two key items.

[0071] Step 6: Calculate the two key terms in the process of solving the linearized equations

[0072] 6.1 Efficient calculation of the first key term C -1 F

[0073] From formula (13), we can see that the large-scale matrix C has a block diagonal structure. Directly using the parallel mode to solve it can obtain the following expression:

[0074]

[0075] In the formula It is C (i) The lower triangular matrix obtained after the Choleshki decomposition of the matrix. Formula (18) involves two types of matrix calculations. One type is the multiplication of a large-scale block diagonal matrix formed by the rotation matrix as a submatrix and a column vector, such as the first and fourth steps in expression (18). This step of calculation can be uniformly expressed as diag(R)V. The other type is the multiplication of the inverse matrix of a large-scale sparse lower triangular matrix and a column vector, such as the second and third steps in expression (18), which can be recorded as The efficient calculation of these two types of matrix multiplication is the key to improving the algorithm's calculation efficiency. The specific operation process is shown in the attached Figure 6 First, when dealing with the first matrix multiplication problem, V 3m×1 The vector is converted into V using the Reshape function 3×m The subscript m represents the total number of finite element nodes. The rotation matrix R is then used to perform the multiplication operation.

[0076] Finally, the obtained matrix is ​​reshaped into a 3m×1 column vector using the Reshape function.

[0077] In calculation When you multiply matrices of the form It can be pre-calculated and stored, and only needs to be back-calculated during each iterative calculation process, without the need for reassembly and decomposition.

[0078] 6.2 Efficient calculation of the second key term

[0079] In key items In this example, the inverse matrix of the C matrix can be expressed as

[0080]

[0081] The D matrix is ​​the rotation matrix R of each component. (i) The large-scale rotation matrix formed. and D are both block diagonal matrices. Substituting expression (19) into the expression of key term 2, we can get the following expression:

[0082]

[0083] It can be seen from expression (20) that as long as the term Then perform matrix multiplication on the transposed matrix and itself to obtain the result of key item 2.

[0084] As attached Figure 7 As shown, we first consider the time-invariant constraints. For each component, there is the following expression

[0085]

[0086] In the formula Represents a matrix whose elements with indices m and n are 1 and all other elements are 0.

[0087] From Section 3.1, we can see that in formula (21) It is a completely constant matrix, which can be pre-calculated and stored before the iterative operation. During the iteration process, the rotation matrix of a component at the current moment is calculated. Then, substituting into formula (21), we can get the Go through all the parts one by one and finally assemble them The above algorithm avoids the time consumption caused by repeated large-scale sparse matrix calculation assembly and multiplication operations.

[0088] Next, consider the time-varying constraint part in formula (20): Calculation of . From Section 3.2, we know that this term is about D and Φ q_1 The time-varying constraint part terms may not necessarily lead to higher computational efficiency in program implementation if they are pre-calculated and back-substituted. In fact, the calculation containing double time-varying terms does not need to be pre-processed by formula (21). Usually, it is a very sparse matrix with a small number of columns. Putting this item directly into the iterative calculation for repeated calculation will not have a significant impact on the calculation time. Therefore, the calculation of this item can be completed in each iteration step.

[0089] After calculating all the above steps, the dynamic equation (6) of the helicopter rotor system can be solved. Then the dynamic time history response is obtained. The dynamic time history response includes global displacement, elastic displacement, velocity and acceleration, etc. So far, the efficient simulation solution method of the helicopter rotor system dynamics based on the C-CR method proposed by the present invention is completed.

[0090] Beneficial and positive effects of the present invention:

[0091] (1) The helicopter rotor system dynamics modeling method based on three-dimensional solid finite elements proposed in the present invention can accurately simulate the helicopter rotor system blade model with any complex shape, and has good modeling adaptability to any hinge component in the rotor system.

[0092] (2) The efficient solution method for the rotor system dynamics based on the component-level co-rotating coordinate method (C-CR) proposed in the present invention can significantly improve the computational efficiency of nonlinear iterative solutions while ensuring the computational accuracy, compared with the traditional unit-level co-rotating coordinate (Full-CR) formulation, full Lagrangian (TL) formulation and updated Lagrangian (UL) formulation.

[0093] (3) The present invention uses a component-level stress solution method to calculate the stress-strain distribution state of the rotor system blade structure based on the concept of component-level co-rotating coordinates. Compared with the traditional beam unit modeling and analysis method, the method of the present invention has higher analysis accuracy and calculation efficiency, and can help engineers to efficiently and high-quality complete the dynamic design and analysis of the rotor system of a new configuration helicopter. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 It is a flow chart of the implementation of the present invention.

[0095] Figure 2 It is the division of components of the finite element model of the rotor system.

[0096] Figure 3 It is a schematic diagram of the method for docking time-invariant constraint components.

[0097] Figure 4 It is a schematic diagram of the docking method of time-varying constraint components.

[0098] Figure 5 The Schur complement method is used to calculate the linear equations.

[0099] Figure 6 It is a fast calculation process for key item one.

[0100] Figure 7 This is the fast calculation process for key item two.

[0101] Figure 8 It is a three-dimensional solid finite element model of the tiltrotor.

[0102] Fig. 9 It is a simulated cloud map of the rotor tilting process.

[0103] Fig.10 It is a comparison of the global displacement time history curves of the blade tip node during the rotor tilt process.

[0104] Fig.11 It is a comparison of the global velocity time history curves of the blade tip node during the rotor tilt process.

[0105] Fig.12 It is a comparison of the global acceleration time history curves of the blade tip node during the rotor tilt process.

[0106] Fig.13 It is a comparison of the time history curves of the local x-direction elastic displacement of the blade tip node during the rotor tilt process.

[0107] Fig.14 It is a comparison of the time history curves of the local y-direction elastic displacement of the blade tip node during the rotor tilt process.

[0108] Fig.15 It is a comparison of the time history curves of the local z-direction elastic displacement of the blade tip node during the rotor tilt process.

[0109] Fig.16 This is a comparison of simulation calculation time using different methods. DETAILED DESCRIPTION

[0110] The present invention is described below in conjunction with specific embodiments.

[0111] Combined with the helicopter rotor dynamics analysis example, the specific implementation method is described and the calculation accuracy and efficiency of the method are compared and analyzed. The specific implementation method of the present invention is as follows:

[0112] The first step is to discretize the helicopter rotor system geometry model into a three-dimensional solid finite element model through finite element meshing.

[0113] In this example, the Figure 8 The geometric model and three-dimensional solid finite element model of the helicopter tilt-rotor system are given. The geometric model contains three blades in total. The length of a single blade is 4.5m, made of aluminum, and has a pre-twist angle of -12° along the span. The rotor system hub is a hingeless rigid hub. The three-dimensional solid finite element model of the rotor system contains 3450 units and 6090 nodes, with a total of 18270 degrees of freedom. The above finite element data is used to generate a model file with model node coordinates and unit node numbers, and the unit material properties are assigned to the units.

[0114] Step 2: Divide the co-rotation analysis components of the 3D solid finite element model

[0115] After the helicopter rotor system geometric model is discretized into a three-dimensional solid finite element model through finite element meshing in the first step, the three-dimensional solid finite element model is divided into several components. The example model is divided into 6 co-rotation analysis components, and the units in each component share the same coordinate system, which is called the component co-rotation coordinate system. The component co-rotation coordinate system is generated by first finding three nodes on each component that are far apart and not on the same straight line, which are called component representative nodes. The corresponding component co-rotation coordinate system is quickly constructed through the edge parallel method of the three component representative nodes.

[0116] The construction method of the component co-rotation coordinate system basis vector is as follows:

[0117]

[0118] Wherein, e1 represents the unit vector formed by the second representative node and the first representative node on the component; e3 represents the unit vector formed by the third and first representative nodes on the component and the cross product of the e1 vector; e2 represents the unit vector obtained by the cross product of the unit vector e3 and e1; x1 represents the coordinate column vector of the first representative node; x2 represents the coordinate column vector of the second representative node; x3 represents the coordinate column vector of the third node; ||…|| represents the modulus of the calculated vector;

[0119] The rotation matrix of the component can be expressed as The establishment of the component co-rotational coordinate system eliminates the solution process of each unit co-rotational coordinate system in the traditional unit-level co-rotational coordinate method, which greatly reduces the program calculation and solution time.

[0120] The third step is to use time-invariant constraints or time-varying constraints to achieve effective connection of various components.

[0121] In the second step, after the three-dimensional solid finite element model of the helicopter rotor system is divided into several components, the multi-body system dynamics concept is adopted to achieve the docking between components through constraints. The constraints include time-invariant constraints and time-varying constraints. In this example, two components on the same blade are docked again through time-invariant constraints. Components 1, 3, and 5 represent the component numbers at the root of each blade. These components are docked through time-varying constraints to form an equivalent rigid hub structure.

[0122] The time-varying constraints and the time-invariant constraints are specifically as follows:

[0123] 3.1 Constraint equations of time-invariant constraints and calculation of constraint Jacobian matrix

[0124] The time-invariant constraint can achieve effective docking between the components. The finite element nodes between two adjacent components are coupled by the method of coupling node degrees of freedom. The method of coupling node degrees of freedom is achieved by limiting the relative coordinates of each pair of coincident nodes between the two components in three dimensions of space. The constraint equation generated by the method of coupling node degrees of freedom is the same as the constraint equation of the ball joint constraint in multibody dynamics. The constraint equation can be expressed as Φ0=0, where the subscript 0 represents the term related to the time-invariant constraint. In this equation, the node coordinates of the component are Generalized coordinates as components After obtaining the constraint Jacobian matrix Φ0 of the spherical joint constraint, the constraint Jacobian matrix Φ of the time-invariant constraint can be calculated by differentiating the constraint equation with respect to the generalized coordinates of the component. q_0 Φ q_0 It is a sparse constant matrix that presents a value of 1 or -1 only at the generalized coordinate position related to a certain component, and all other positions are 0. Therefore, the Jacobian matrices of all spherical joint constraints can be pre-calculated and stored, and these matrices can be directly substituted in each iterative solution process of the simulation, thereby improving the computational efficiency.

[0125] 3.2 Constraint equations for time-varying constraints and calculation of constraint Jacobian matrix

[0126] The time-varying constraints include angle-driven constraints, vertical constraints, and distance constraints. The combination and superposition of these constraints can achieve effective docking between the helicopter hub structure and components. The angle-driven constraint is used to control the rotational motion of the three-dimensional solid finite element model. The vertical constraint and distance constraint are used to constrain the vertical angle and relative distance of the three-dimensional solid finite element model. In order to realize the three time-varying constraints, it is necessary to introduce external auxiliary nodes in the three-dimensional solid finite element model. The purpose of introducing external auxiliary nodes is to combine the auxiliary nodes with the finite element nodes in the three-dimensional solid finite element model to jointly complete the writing of the time-varying constraint equation Φ1. The generalized coordinates of the auxiliary nodes are denoted by q aux . By q aux and the generalized coordinates of the components described in Section 3.1 Vectors A, B, C, and D are generated between the auxiliary nodes and the nodes of the three-dimensional solid finite element model of the component. Vector A represents the arm in the hub structure, vector C is the rotating shaft of the helicopter rotor system, and vectors B and D are two fixed vectors on the helicopter body, which are used to fix the entire rotor system. These vectors are used to write the constraint equation Φ1. The angle drive constraint is achieved by controlling the change of the angle θ between vectors A and B. The constraint equation Φ1 of the time-varying constraint is expressed as follows:

[0127] Φ1 A =A·B-|A||B|cosθ (23)

[0128] Φ1 V =C·D (24)

[0129] Φ1 D =|A|-d (25)

[0130] Among them, the subscript 1 represents the time-varying constraint; Φ1 A represents the constraint equation of the angle drive constraint; θ represents the drive angle; Φ1 V Indicates vertical constraint; Φ1 D Represents distance constraint; d represents the control distance of the distance constraint.

[0131] The constraint equation Φ1 of the time-varying constraint is transformed into the generalized coordinate q aux and The derivative can obtain the constraint Jacobian matrix Φ of the time-varying constraint q_1 The constraint Jacobian matrix Φ of the time-varying constraint q_1 It is highly sparse.

[0132] Step 4: Establish the dynamic equation of the rotor system.

[0133] After the above three steps, the generalized coordinates q of the rotor system are the combination of the generalized coordinates of all three-dimensional solid finite element nodes and the generalized coordinates of the introduced auxiliary nodes:

[0134]

[0135] in, represents the generalized coordinate column vector of the 3D solid finite element nodes in the first component; represents the generalized coordinate column vector of the 3D solid finite element nodes in the second component; Indicates the Nth b The generalized coordinate column vector of the 3D solid finite element nodes within each component; Represents the column vector of generalized coordinates of auxiliary nodes.

[0136] Using the first kind of Lagrangian equations and considering all constraint equations (23-25), the dynamic control equations of the helicopter rotor system can be written as the following index 3 differential algebraic equations (DAEs):

[0137]

[0138] in, represents the generalized acceleration corresponding to all generalized coordinates in the helicopter rotor system; Φ(q,t) represents the total constraint equation of the helicopter rotor system, which includes the time-invariant constraint equation Φ0 in Section 3.1 and the time-varying constraint equation Φ1 in Section 3.2; M is the mass matrix of the helicopter rotor system; Φ q =[Φq_0 ; Φ q_1 ] is the system constraint Jacobian matrix; λ is the Lagrange multiplier vector reflecting the constraint force physically. INT is the internal force column vector of all nodes, which is calculated from the elastic internal force of the component; f EXT It is the external force of the helicopter rotor system, including gravity and aerodynamic force.

[0139] Furthermore, the internal force vector f INT And the mass matrix M of the helicopter rotor system can be expressed as:

[0140]

[0141] in, represents the mass matrix of the first component; represents the mass matrix of the second component; Indicates the Nth b The mass matrix of a component; M aux Represents the auxiliary node mass matrix.

[0142] Step 5: Efficiently solve the dynamic equations of the helicopter rotor system and obtain the structural dynamic response.

[0143] 5.1 Calculation of the Jacobian Matrix of the Helicopter Rotor System

[0144] In order to calculate the Jacobian matrix of the helicopter rotor system, this example sets the calculation parameters. The physical time of the simulation process is 10s, which is divided into 10,000 time integration steps, with a step length of 0.001s. The convergence tolerance of the Newton-Raphson iteration is set to 5×10 -3 The Newmark integral parameters are 0.56 and 0.25 respectively. The maximum rotor speed is 20 rad / s, and it accelerates uniformly from rest to the maximum speed within the first second. At t=1.2s, the rotor begins to tilt, and the rotor tilt angular velocity is 0.157 rad / s.

[0145] The typical implicit Newmark direct integration method is further used to discretize the system dynamics equation into a set of nonlinear algebraic equations at each integral / time step. If the solution of the first n-th time step is obtained, the differential algebraic equation is discretized and the Newmark hypothesis is introduced in the n+1th step to produce the following nonlinear equations:

[0146]

[0147] Where h is the integration step; α and δ represent two parameters in the Newmark algorithm; t n+1 represents the n+1th moment after time discretization; represents the acceleration column vector at time t=n+1; q n+1 Represents the generalized coordinate column vector at time t=n+1; Represents the velocity column vector at time t=n+1.

[0148] After derivation, we get is a nonlinear equation system with unknown variables. The linear equation system in the iteration step obtained by the Newton-Raphson method is:

[0149] JΔX=G RES (31)

[0150] Where J is the Jacobian matrix; is the iteration increment; G RES =[F T ,Ψ T ] T is the residual of the nonlinear equation;

[0151] Furthermore, the Jacobian matrix J is:

[0152]

[0153] The Jacobian matrix can be abbreviated as:

[0154]

[0155] 5.2 Fast Solution of Linear Equations

[0156] In formula (33), the C matrix is ​​a large-scale block diagonal matrix, expressed as

[0157]

[0158] in,

[0159]

[0160] in is the rotation matrix of the component described in Section 2.

[0161] In order to quickly solve the large-scale linear equation system (31), the Schur complement method can be used. First, the linear equation system (31) is divided into blocks. After the block division, the equation system (31) can be rewritten as:

[0162]

[0163] The first equation of formula (36) can be solved as:

[0164]

[0165] Combining formula (37) with the second equation of formula (36), we can solve the Lagrangian multiplier increment of the nonlinear algebraic equation:

[0166]

[0167] Finally, substituting formula (38) back into formula (37) yields the acceleration increment: The final result. This term is the Schur complement of the submatrix C in the matrix J. Therefore, the problem of solving a large-scale linearized system of equations is transformed into solving in each iteration step by the Schur complement method. and C -1 F these two key items.

[0168] Step 6: Calculate the two key terms in the process of solving the linearized equations

[0169] 6.1 Efficient calculation of the first key term C -1 F

[0170] From formula (34), we can see that the large-scale matrix C has a block diagonal structure. Directly using the parallel mode to solve it can obtain the following expression:

[0171]

[0172] In the formula It is C (i) The lower triangular matrix obtained after the Choleshki decomposition of the matrix. Formula (39) involves two types of matrix calculations. One type is the multiplication of a large-scale block diagonal matrix formed by the rotation matrix as a submatrix and a column vector, such as the first and fourth steps in expression (18). This step of calculation can be uniformly expressed as diag(R)V. The other type is the multiplication of the inverse matrix of a large-scale sparse lower triangular matrix and a column vector, such as the second and third steps in expression (39), which can be recorded as The efficient calculation of these two types of matrix multiplication is the key to improving the algorithm's calculation efficiency. The specific operation process is shown in the attached Figure 6 First, when dealing with the first matrix multiplication problem, V 3m×1 The vector is converted into V using the Reshape function 3×m The subscript m represents the total number of finite element nodes. Then the rotation matrix R is used to perform the multiplication operation. Finally, the obtained matrix is ​​converted into a 3m×1 column vector using the Reshape function again.

[0173] In calculation When you multiply matrices of the form It can be pre-calculated and stored, and only needs to be back-calculated during each iterative calculation process, without the need for reassembly and decomposition.

[0174] 6.2 Efficient calculation of the second key term

[0175] In key items In this example, the inverse matrix of the C matrix can be expressed as

[0176]

[0177] in and D are both block diagonal matrices, from which we can get the following expression

[0178]

[0179] It can be seen from expression (41) that as long as the term Then perform matrix multiplication on the transposed matrix and itself to obtain the result of key item 2.

[0180] First, calculate the time-invariant part of this term. Figure 7 As shown, for each component, there is the following expression

[0181]

[0182] In the formula represents a matrix in which the element value at the position with index m,n is 1 and the element value at other positions is 0. As can be seen from Section 3.1, the value in formula (42) It is a completely constant matrix, which can be pre-calculated and stored before the iterative operation. When the rotation matrix of a component at the current moment is calculated during the iteration process Then, substituting into formula (21), we can get the Go through all the parts one by one and finally assemble them The above algorithm avoids the time consumption caused by repeated large-scale sparse matrix calculation assembly and multiplication operations.

[0183] Next, consider the time-varying part of formula (41): Calculation of . From Section 3.2, we know that this term is about D and Φ q_1 The time-varying constraint part terms may not necessarily lead to higher computational efficiency in program implementation if they are pre-calculated and back-substituted. In fact, the calculation containing double time-varying terms does not need to be pre-processed by formula (21). Usually, it is a very sparse matrix with a small number of columns. Putting this item directly into the iterative calculation for repeated calculation will not have a significant impact on the calculation time. Therefore, the calculation of this item can be completed in each iteration step.

[0184] After calculating all the above steps, the dynamic differential algebraic equation (27) of the helicopter rotor system can be solved. Then the dynamic time history response is obtained. The dynamic time history response includes global displacement, elastic displacement, velocity and acceleration, etc. So far, the helicopter rotor system dynamics efficient simulation solution method based on the C-CR method proposed by the present invention is completed.

[0185] Step 7: Analysis of calculation results

[0186] In order to compare and analyze the results of the example, we established the same model in the commercial finite element software ANSYS. The rigid hub structure was modeled using the MPC-184 rigid beam element, while the rotor rotary drive hinge and tilt hinge were modeled using the MPC-184 pin connection element. At the same time, the unit-level co-rotating coordinate (Full CR) method was also used to verify the results.

[0187] Attached Fig. 9 The dynamic simulation effect diagram of the tilt-rotor system calculated by the component-level co-rotating coordinate method is shown. The C-CR algorithm can calculate the motion trajectory of the rotor tip node during the entire tilting process. We project the spatial motion trajectory, motion velocity and acceleration of the node to the three coordinate axes respectively, and compare them with the calculation results of ANSYS and the calculation results of the unit-level co-rotating coordinate (Full CR) method. The relevant curves are shown in the attached figure. Fig.10 , Attachment Fig.11 and attached Fig.12 The calculation results show that the displacement, velocity and acceleration of the blade tip node obtained by the three different methods are basically consistent, and the errors are within the allowable range. Due to the influence of the elastic effect, the acceleration curve experiences high-frequency oscillation and large overall movement. Fig.13 , Attachment Fig.14 and attached Fig.15 The time history curves of the local elastic displacement of the blade tip node in the x, y, and z directions obtained by the three calculation methods are shown respectively. The vibration trends of these curves are basically consistent, and the results verify the correctness of the C-CR calculation method.

[0188] In terms of computational efficiency, ANSYS takes 32030 seconds, the unit-level co-rotating coordinate method (Full-CR) takes 24066 seconds, and the C-CR method takes 4358 seconds, which is about 1 / 7 of the time taken by ANSYS and 1 / 5 of the time taken by the Full-CR algorithm. It can be seen that the C-CR method has a higher computational efficiency when dealing with helicopter rotor system dynamics problems.

[0189] The above-described embodiments merely express the implementation methods of the present invention, but they cannot be understood as limiting the scope of the patent 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. An efficient solution method for helicopter rotor system dynamics simulation based on three-dimensional solid elements, characterized in that: The efficient solution method for helicopter rotor system dynamics simulation firstly discretizes the helicopter rotor geometry model into a three-dimensional solid finite element model through finite element meshing; secondly, divides the co-rotation analysis components of the three-dimensional solid finite element model; thirdly, uses time-invariant constraints or time-varying constraints to achieve effective connection of each co-rotation analysis component; Fourth, establish the dynamic equations of the rotor system. Fifth, efficiently solve the dynamic equations of the helicopter rotor system. Sixth, calculate the two key terms in the process of solving the linearized equations.

2. The efficient solution method for helicopter rotor system dynamics simulation based on three-dimensional solid elements according to claim 1 is characterized in that: The specific steps include: The first step is to discretize the helicopter rotor geometry model into a three-dimensional solid finite element model through finite element meshing, specifically: The helicopter rotor system geometric model is meshed by finite element using three-dimensional solid units to obtain the finite element model of the rotor system; secondly, the model node coordinates and unit node number model file of the finite element model are generated; the unit material properties are assigned to the corresponding units; The second step is to divide the co-rotation analysis components of the three-dimensional solid finite element model, specifically: The three-dimensional solid finite element model is divided into several co-rotation analysis components. The units in each co-rotation analysis component share the same coordinate system, which is called the component co-rotation coordinate system. The component co-rotation coordinate system is generated by finding three nodes that are far apart and not on the same straight line on each component in turn, which are called component representative nodes. The corresponding component co-rotation coordinate system is quickly constructed by the edge parallel method. The third step is to use time-invariant constraints or time-varying constraints to achieve effective connection of various components, specifically: Adopting the multi-body system dynamics concept, the docking between the co-rotation analysis components is realized through constraints; the constraints include time-invariant constraints and time-varying constraints, and the time-invariant constraint equation Φ0 and the time-varying constraint equation Φ1 are obtained; Step 4: Establish the dynamic equation of the rotor system After the first to third steps, the generalized coordinates q of the helicopter rotor system are the combination of the generalized coordinates of all three-dimensional solid finite element nodes and the generalized coordinates of the introduced auxiliary nodes: in, represents the generalized coordinate column vector of the 3D solid finite element nodes in the first component; represents the generalized coordinate column vector of the 3D solid finite element nodes in the second component; Indicates the Nth b The generalized coordinate column vector of the 3D solid finite element nodes within each component; represents the generalized coordinate column vector of the auxiliary nodes; Using the first kind of Lagrangian equations and considering all constraint equations, the dynamic control equations of the helicopter rotor system are written as the following index 3 differential algebraic equations DAEs: in, represents the generalized acceleration corresponding to all generalized coordinates in the helicopter rotor system; Φ(q,t) represents the total constraint equation of the helicopter rotor system, including the time-invariant constraint equation Φ0 and the time-varying constraint equation Φ1; M is the mass matrix of the helicopter rotor system; Φ q =[Φ q_0 ; Φ q_1 ] is the system constraint Jacobian matrix; λ is the Lagrange multiplier vector reflecting the constraint force physically; f INT is the internal force column vector of all nodes, which is calculated from the elastic internal force of the component; EXT It is the external force of the helicopter rotor system, including gravity and aerodynamic force; Step 5: Efficiently solve the helicopter rotor system dynamics equations 5.1 Calculation of the Jacobian Matrix of the Helicopter Rotor System At each integral / time step, the rotor system dynamics equations as shown in formula (6) are discretized into a set of nonlinear algebraic equations; if a solution is obtained at the nth time step, the differential algebraic equations are discretized, and the Newmark hypothesis is introduced at the n+1th step to generate nonlinear equations; Then get is a nonlinear system of unknowns; where q n+1 represents the generalized coordinate column vector at time t=n+1, then represents the acceleration column vector at time t=n+1, λ n+1 represents the Lagrangian multiplier column vector at time t=n+1; The linear equations in the iteration step obtained by the Newton-Raphson method are: JΔX=G RES (10) Where J is the Jacobian matrix; is the iteration increment; G RES =[F T ,Ψ T ] T is the residual of the nonlinear equation; The Jacobian matrix J is: The Jacobian matrix is ​​abbreviated as: 5.2 Fast Solution of Linear Equations The C matrix is ​​represented as: in, in, The rotation matrix of the component for the corotation analysis; The Schur complement method is used to solve the linear equation system (10); the linear equation system (10) is processed in blocks and written as: Solving the first equation of formula (15) yields: Combine formula (16) with the second equation of formula (15) to solve the Lagrangian multiplier increment of the nonlinear algebraic equation: Substituting formula (17) back into formula (16) yields the acceleration increment: The final result; and C -1 F is the two key terms for solving large-scale linearized equations; Step 6: Calculate the two key terms in the process of solving the linearized equations 6.1 Efficient calculation of the first key term C -1 F From formula (13), we can see that the large-scale matrix C has a block diagonal structure. Using the parallel mode to solve it, we can get the following expression: In the formula, It is C (i) The lower triangular matrix obtained after the Choleshsky decomposition of the matrix; Formula (18) involves two types of matrix calculations. One type is the multiplication of a large-scale block diagonal matrix formed by a rotation matrix as a submatrix and a column vector, such as the first and fourth steps in expression (18), which are uniformly expressed as diag(R)V; the other type is the multiplication of the inverse matrix of a large-scale sparse lower triangular matrix and a column vector, such as the second and third steps in expression (18), which are expressed as The efficient calculation process of these two types of matrix multiplication is as follows: First, when dealing with the first type of matrix multiplication problem, V 3m×1 The vector is converted into V using the Reshape function 3×m Matrix, the subscript m represents the total number of finite element nodes; secondly, the rotation matrix R is used to complete the multiplication operation; finally, the obtained matrix is ​​again formed into a 3m×1 column vector using the Reshape function; In calculation When you multiply matrices of the form Can be pre-calculated and stored, only back-calculation is required during each iterative calculation; 6.2 Efficient calculation of the second key term In key items In this example, the inverse matrix of the C matrix is ​​expressed as Among them, the D matrix is ​​composed of the rotation matrix R of each component (i) The large-scale rotation matrix formed; and D are both block diagonal matrices; Substituting expression (19) into the expression of key term 2, we get the following expression First, consider the time-invariant constraints; for each component, there is the following expression In the formula, Represents a matrix whose element values ​​are 1 and other elements are 0. The rotation matrix of a component at the current moment is calculated during the iteration process Then, substitute into formula (21) to obtain the value of component i Go through all the parts one by one and finally assemble them Secondly, in each iteration step, the time-varying constraint part of formula (20) is calculated Finally, the dynamic equation (6) of the helicopter rotor system is solved and the dynamic time history response is obtained.

3. The efficient solution method for helicopter rotor system dynamics simulation based on three-dimensional solid elements according to claim 2 is characterized in that: In the second step, the construction method of the component co-rotation coordinate system basis vector is as follows: Among them, e1 represents the unit vector formed by the second representative node and the first representative node on the component; e3 represents the unit vector formed by the third and first representative nodes on the component and the cross product of the e1 vector; e2 represents the unit vector obtained by the cross product of the unit vector e3 and e1; x1 represents the coordinate column vector of the first representative node; x2 represents the coordinate column vector of the second representative node; x3 represents the coordinate column vector of the third representative node; ||…|| represents the modulus of the calculated vector; the rotation matrix of the component is expressed as 4. The efficient solution method for helicopter rotor system dynamics simulation based on three-dimensional solid elements according to claim 2 is characterized in that: In the third step, the time-invariant constraints and the time-varying constraints are specifically: 3.1 Constraint equations and constraint Jacobian matrix of time-invariant constraints The time-invariant constraint achieves effective docking between the co-rotation analysis components; the finite element nodes between two adjacent co-rotation analysis components are coupled by using the method of coupling node degrees of freedom; The constraint equation generated by the coupled nodal degree of freedom method is the same as the constraint equation of the spherical joint constraint in multibody dynamics. The constraint equation is expressed as Φ0=0, where the subscript 0 represents the term related to the time-invariant constraint; let the node coordinates of the component Generalized coordinates as components The constraint equations are differentiated with respect to generalized coordinates to calculate the constraint Jacobian matrix Φ of the spherical joint constraint. q_0 ; 3.2 Constraint Equations and Constraint Jacobian Matrix of Time-Varying Constraints The time-varying constraints include angle drive constraints, vertical constraints and distance constraints; the effective docking between the helicopter hub structure and the components is achieved through the combination and superposition of these constraints; external auxiliary nodes are introduced into the three-dimensional solid finite element model; The auxiliary node is combined with the finite element node in the three-dimensional solid finite element model to jointly complete the writing of the time-varying constraint equation Φ1; the generalized coordinate of the auxiliary node is denoted by q aux , through q aux and the generalized coordinates of the components Vectors A, B, C, and D are generated between the auxiliary nodes and the nodes of the three-dimensional solid finite element model of the component; vector A represents the arm in the hub structure, vector C is the helicopter shaft, and vectors B and D are two fixed vectors on the helicopter body, which are used to fix the entire rotor system; the above vectors are used to write the constraint equation Φ1, in which the angle drive constraint is realized by controlling the change of the angle θ between the vectors A and B; the constraint equation Φ1 of the time-varying constraint is expressed as follows: Φ1 A =A·B-|A||B|cosθ (2) Φ1 V =C·D (3) Φ1 D =|A|-d (4) Among them, the subscript 1 represents the time-varying constraint; Φ1 A The constraint equation representing the angle-driven constraint; θ represents the driving angle; Φ1 V represents the vertical constraint equation; Φ1 D represents the distance constraint equation; d represents the control distance of the distance constraint; The constraint equation Φ1 of the time-varying constraint is transformed into the generalized coordinate q aux and The derivative can obtain the constraint Jacobian matrix Φ of the time-varying constraint q_1 ; The constraint Jacobian matrix of the time-varying constraint Φ q_1 It is highly sparse.

5. The efficient solution method for helicopter rotor system dynamics simulation based on three-dimensional solid elements according to claim 4 is characterized in that: The angle drive constraint is used to control the rotational movement of the three-dimensional solid finite element model; the vertical constraint and distance constraint are used to constrain the vertical angle and relative distance of the components.

6. The efficient solution method for helicopter rotor system dynamics simulation based on three-dimensional solid elements according to claim 2 is characterized in that: In the fourth step, the internal force vector f INT And the mass matrix M of the helicopter rotor system can be expressed as: in, represents the mass matrix of the first component; represents the mass matrix of the second component; Indicates the Nth b The mass matrix of each component; M aux Represents the auxiliary node mass matrix.

Images