Flexible blade structure dynamics system calculation method

Abstract

The invention belongs to the field of helicopter rotor dynamics analysis and design, and discloses a flexible blade structure dynamics system calculation method. A helicopter rotor is of a typical flexible cantilever beam structure with a high aspect ratio, the basic motion characteristic of the helicopter rotor is that large-range rigid body rotation and elastic deformation coexist, the helicopter rotor is a typical rigid-flexible coupling system, and the deformation rule is very complex. For such systems, the dynamics behaviors of the systems are generally described by adopting a multi-body system dynamics theory based on Cartesian coordinates. From the form, the kinetic equation of the flexible blade is generally DAEs with a second-order differential form. In order to ensure the calculation precision of system response and structural dynamic characteristics and solve the defect of constraint default existing in a traditional direct integration method, the invention provides a DAEs solving method which adopts a simple and easy-to-implement default correction strategy to accurately correct the generalized displacement and the generalized speed of the system, so that the accuracy of the system response and structural dynamic characteristics is improved. Therefore, the accurate solution of the blade structure dynamics DAEs is realized.

Description

Technical Field

[0001] This invention belongs to the field of helicopter rotor dynamics analysis and design, and discloses a calculation method for the dynamic system of a flexible blade structure. Background Technology

[0002] A fundamental characteristic of helicopter rotor motion is the coexistence of large-scale rigid body rotation and elastic deformation, making it a typical rigid-flexible coupled system with highly complex deformation laws. Furthermore, as user demands for helicopter performance increase, the aerodynamic shape of helicopter blades becomes increasingly complex, leading to stronger nonlinear effects and making dynamic behavior analysis more difficult. Employing multibody system-based mechanical analysis methods is one of the most effective means of solving the dynamic calculation and analysis of blades with complex aerodynamic shapes. When modeling the motion of flexible helicopter blades using multibody system dynamics theory based on Cartesian coordinates, the resulting dynamic equations are typically a system of differential-algebraic equations (DAEs) in second-order differential form. Due to the significant nonlinear characteristics of the rigid-flexible coupling, analytical solutions to these equations are very difficult, necessitating numerical methods.

[0003] Traditional direct integration methods solve the system's dynamic differential equations and acceleration constraint equations simultaneously to obtain the generalized acceleration, and then calculate the generalized displacement and velocity using Newmark time integration. While this method is easy to implement, it can lead to the generalized velocity and displacement failing to satisfy the system constraints. As the number of integration steps increases, the numerical solution eventually diverges, and even with constraint violation stabilization techniques (such as the Baumgrate constraint violation stabilization method), the computational stability and accuracy of the system cannot be fully guaranteed. To address this issue, researchers have improved the traditional direct integration method by solving the system's dynamic differential equations and the original displacement constraint equations simultaneously. While this effectively solves the violation problem of generalized velocity and displacement, the cost is that the Jacobian coefficient matrix of the entire dynamic system must be calculated at each time iteration step. When the system has strong nonlinear coupling, its Jacobian coefficient matrix cannot be analytically obtained and must be numerically approximated using the finite difference method, which inevitably increases additional computational and storage costs. Summary of the Invention

[0004] Objective: This invention proposes a computational method for the dynamic system of a flexible blade structure. This method avoids the cumbersome computational process of numerically solving the system's Jacobian coefficient matrix, significantly reducing computational and storage costs. Furthermore, compared to the traditional Baumgrate constraint-violated stability method, the new method performs accurate correction calculations for the system's generalized velocity and displacement, effectively solving the problem of velocity and displacement constraint violation.

[0005] The technical solution is as follows: A calculation method for the dynamic system of a flexible blade structure, comprising the following steps: Step 1: Construct the differential-algebraic equation system (DAEs) for the dynamic system of the flexible rotor blade structure using multibody system dynamics theory and helicopter rotor blade structure dynamics theory, and set the time step ∆. t and convergence tolerance ， Given the initial conditions of the system and ; Step two, within a short time interval Inside, according to t The numerical solution of the system at time 1 , and Solve according to the following process t +∆ t Numerical solution at time , and : ①For the first k ( k =0,1,2,...) iteration steps ， Numerical solution and It has been determined by calculation in the previous iteration step; if k =0, then , ; ②Based on the fundamental knowledge of helicopter rotor blade structural dynamics and multibody system dynamics, with and Computational system mass array Generalized force and constrained Jacobian matrix Information such as; ③ with and Given unknown variables, use numerical methods to solve the system of algebraic equations. k Generalized acceleration of +1 iteration step ; ④ Calculate the temporary values ​​of generalized velocity and generalized displacement. , ; ⑤ Based on the fundamental knowledge of helicopter rotor blade structural dynamics and multibody system dynamics, , Calculate the temporary values ​​of the system mass matrix, generalized forces, and constraint Jacobian matrix, respectively. , and ; ⑥ Calculate the generalized velocity using the algorithm proposed in this invention. and generalized displacement ; ⑦ Determine if convergence has occurred. If convergence has occurred, terminate the iteration and obtain the result. t +∆ t The system solution at time t, i.e. , , Otherwise, let k = k +1, return to ①.

[0006] Step 3: Refer to Step 2 to calculate the next time step, until the time expires.

[0007] Furthermore, in step one, the process of constructing the differential-algebraic equation system of the flexible blade structure dynamic system is as follows: First, the blade is divided into several straight beam segments. A floating coordinate system is established at the left endpoint of the region pointing towards the blade root in each straight beam segment. By introducing 6 rigid body degrees of freedom and 14 elastic degrees of freedom, the large deformation motion of the flexible beam structure of the blade is decomposed into a superposition of large-scale rigid body motion in the floating coordinate system and small-to-medium elastic deformation motion relative to the floating coordinate system. Then, relative rotational and relative translational geometric constraints are used to coordinate the connection deformation between adjacent straight beam segments, thereby describing the complex nonlinear deformation behavior of the entire flexible blade. The 6 rigid body degrees of freedom include 3 translational degrees of freedom and 3 rotational degrees of freedom; the 14 elastic degrees of freedom include 3 flapping degrees of freedom, 3 oscillating degrees of freedom, 4 torsional degrees of freedom, and 4 tensile degrees of freedom.

[0008] Furthermore, in step one, the system of equations is as follows:

[0009] in, t Indicates time; Represents the generalized coordinates of a multibody system, where , and These represent the translational, rotational, and elastic degrees of freedom of the system, respectively. ,in , , and These represent the degrees of freedom for stretching, oscillation, swinging, and torsion, respectively. and These are called generalized velocity and generalized acceleration, which represent the generalized coordinates, respectively. First and second derivatives with respect to time; The following text is abbreviated as , represents the generalized mass matrix of a multibody system; The following text is abbreviated as , represents the constraint function vector of a multibody system, which describes the coordination relationship between translation, rotation and elastic deformation among the straight beam segments of the blade; The following text is abbreviated as , representing the constraint function vector For generalized coordinates The Jacobian matrix; The following text is abbreviated as , represents the generalized external force vector; Represents the constrained Lagrange multiplier vector; the superscript "T" indicates matrix transpose.

[0010] Furthermore, in step ③ of step two, the solution is... k Generalized acceleration of +1 iteration step The process is as follows: constraint equations Regarding time t Differentiating, we obtain the generalized velocity constraint equation, i.e.

[0011] here Represents the constraint function vector Regarding time t The partial derivatives; Applying the generalized velocity constraint equation to time t Differentiating again, we obtain the generalized acceleration constraint equation.

[0012] here Representing vectors right Jacobian matrix, express Regarding time t The partial derivatives, Represents the constraint function vector Regarding time t The second-order partial derivatives. For ease of description, the right-hand side of the generalized acceleration constraint equation is represented by... It means, that is .

[0013] The multibody system dynamics differential equations and generalized acceleration constraint equations are combined and expressed as follows:

[0014] Based on the k Numerical solution for each iteration step and ,by and For unknown variables, solving the above system of equations yields the system's... kGeneralized acceleration of +1 iteration step , means as follows

[0015] Furthermore, in step 4 of step two, the generalized acceleration is... Perform Newmark integration to obtain the th k Temporary values ​​of generalized velocity and generalized displacement for +1 iteration step , ,Right now

[0016] here Represents the time step, parameters and It determines the integration accuracy and stability, and its function is to control... The change in acceleration over the time interval is often taken as in actual calculations. .

[0017] Furthermore, in step 6 of step two, the process is as follows: The combined equations of the multibody system dynamics differential equations and the generalized velocity constraint equations are expressed as follows:

[0018] In order to obtain the solution from the above simultaneous equations The Newmark integral formula is transformed as follows:

[0019] Substituting this into the above system of equations and rearranging, we get...

[0020] In this way, and For unknown variables, solving the above system of equations will yield the system's... k Generalized speed of +1 iteration step ; The dynamic differential equations of the multibody system are combined with the displacement constraint equations, expressed as follows:

[0021] Similarly, the Newmark integral formula is first transformed as follows:

[0022] Substituting these equations into the system of simultaneous equations, we can transform it into the following form.

[0023] by and Let be the unknown variable. Solving the above system of equations yields the first... k +1 iteration step generalized displacement .

[0024] Clearly, in solving (or correcting) the generalized velocity... Generalized displacement Simultaneous equations and solving for generalized acceleration The simultaneous equations are identical in form, which simplifies the programming process. It only requires abstracting the solution process of the three equations into a function design, and then inputting specific parameters for different objects. Furthermore, due to the generalized speed... Generalized displacement The solution utilizes velocity constraint equations and displacement constraint equations, thus the numerical solution automatically satisfies the system's displacement and velocity constraints, perfectly resolving the problem of displacement and velocity violations.

[0025] Furthermore, in step 2, step 7, the convergence criterion is as follows: like If the expression converges, then the convergence is complete; otherwise, it does not converge. ε This is the convergence threshold.

[0026] Furthermore, in step ① of step two, the initial iteration step is... k When =0, there is , .

[0027] An electronic device includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of a calculation method for a dynamic system of a flexible blade structure.

[0028] A computer-readable storage medium storing a calculator program thereon, which, when executed by a processor, implements the steps of a calculation method for a dynamic system of a flexible blade structure.

[0029] In summary, the beneficial effects of the present invention are as follows: For the numerical solution of the dynamic system of flexible blade structures, the calculation method proposed in this invention can effectively solve the violation problem of generalized velocity and generalized displacement without sacrificing calculation accuracy. Compared with the currently popular solution methods, this invention is simpler to implement, avoiding the cumbersome calculation process of solving the Jacobian coefficient matrix of the entire system, thus saving computational and storage costs. A comparison of the calculation results and experimental results (reference results) shows that the new method proposed in this invention has high calculation accuracy. The error of the rotational dynamic characteristics of the BO105 blade is within 0.5%, and the error between the calculated and experimental values ​​of the Princeton beam static deformation is within 3%. Even under large deformation conditions, the new method remains stable and reliable, and the calculation results are in consistent agreement with the reference results.

[0030] Unlike popular computational analysis software (such as CAMRAD II, RCAS, etc.) which use finite difference to calculate the unit mass matrix and generalized forces, the new method proposed in this invention can directly derive the analytical expressions of the system mass matrix and generalized forces, further simplifying the calculation process. Attached Figure Description

[0031] Figure 1 This is a schematic diagram comparing the calculated deformation and displacement results of the end of a straight beam under different static loads with the reference results. Detailed Implementation

[0032] Helicopter rotors are typical high-aspect-ratio flexible cantilever beam structures. Their fundamental kinematic characteristic is the coexistence of large-scale rigid body rotation and elastic deformation, making them typical rigid-flexible coupled systems with highly complex deformation laws. For such systems, multibody dynamics theory based on Cartesian coordinates is generally used to describe their dynamic behavior. Formally, the dynamic equations of flexible rotor blades are typically DAEs with second-order differential forms. To ensure the accuracy of system response and structural dynamic characteristic calculations, and to address the constraint violation defects (displacement and velocity) inherent in traditional direct integration methods, this invention proposes a solution method for the dynamic system of flexible rotor blade structures. By employing a simple and easily implemented constraint correction strategy, the generalized displacement and generalized velocity of the system are accurately corrected (calculated), thereby achieving an accurate solution for the rotor blade structure dynamic system. The core steps of this method can be summarized as follows.

[0033] Step 1: Construct the differential-algebraic equations (DAEs) for the dynamic system of the flexible blade structure.

[0034] The complex deformation motion of a high aspect ratio flexible blade is described using the multibody system dynamics theory based on Cartesian coordinates and the helicopter blade structure dynamics theory.

[0035] First, the blade is divided into several straight beam segments. A floating coordinate system is established at the left endpoint of each straight beam segment pointing towards the blade root region. By introducing 6 rigid body degrees of freedom (including 3 translational and 3 rotational degrees of freedom) and 14 elastic degrees of freedom (including 3 flapping, 3 oscillating, 4 torsional, and 4 tensile degrees of freedom), the large deformation motion of the flexible beam structure of the blade is decomposed into a superposition of large-scale rigid body motion in the floating coordinate system and small-to-medium elastic deformation motion relative to the floating coordinate system. Then, relative rotational and relative translational geometric constraints are used to coordinate the connection deformation between adjacent straight beam segments, thereby accurately characterizing the complex nonlinear deformation behavior of the entire flexible blade. The mathematical modeling process is not detailed here; only the final expression of the dynamic equations of the flexible blade structure in the multi-body floating coordinate system is given. (1) in, t Indicates time; Represents the generalized coordinates of a multibody system, where , and These represent the translational, rotational, and elastic degrees of freedom of the system, respectively. ,in , , and These represent the degrees of freedom for stretching, oscillation, swinging, and torsion, respectively. and These are called generalized velocity and generalized acceleration, which represent the generalized coordinates, respectively. First and second derivatives with respect to time; The following text is abbreviated as , represents the generalized mass matrix of a multibody system; The following text is abbreviated as , represents the constraint function vector of a multibody system, which describes the coordination relationship between translation, rotation and elastic deformation among the straight beam segments of the blade; The following text is abbreviated as , representing the constraint function vector For generalized coordinates The Jacobian matrix; The following text is abbreviated as , represents the generalized external force vector; Represents the constrained Lagrange multiplier vector; the superscript "T" indicates matrix transpose.

[0036] Step two, within a short time interval Inside, according to t The numerical solution of the system at time 1 , and Solve according to the following processt +∆ t Numerical solution at time , and : ①For the first k ( k =0,1,2,...) iteration steps ， Numerical solution and It has been determined by calculation in the previous iteration step; if k =0, then , ; ②Based on the fundamental knowledge of helicopter rotor blade structural dynamics and multibody system dynamics, taking the current... k Generalized speed of each iteration step and generalized displacement Calculate the mass matrix of the system Generalized force and constrained Jacobian matrix Information such as; ③ with and Given unknown variables, use numerical methods to solve the system of algebraic equations. k Generalized acceleration of +1 iteration step The specific process is as follows: constraint equations Regarding time t Differentiating, we obtain the generalized velocity constraint equation, i.e.

[0037] here Represents the constraint function vector Regarding time t The partial derivatives; Applying the generalized velocity constraint equation to time t Differentiating again, we obtain the generalized acceleration constraint equation.

[0038] here Representing vectors right Jacobian matrix, express Regarding time t The partial derivatives, Represents the constraint function vector Regarding time t The second-order partial derivatives. For ease of description, the right-hand side of the generalized acceleration constraint equation is represented by... It means, that is .

[0039] The multibody system dynamics differential equations and generalized acceleration constraint equations are combined and expressed as follows:

[0040] Based on the k Numerical solution for each iteration step and ,by and For unknown variables, solving the above system of equations yields the system's... k Generalized acceleration of +1 iteration step ,Right now

[0041] ④ Calculate the temporary values ​​of generalized velocity and generalized displacement. , The specific process is as follows: Generalized acceleration Perform Newmark integration to obtain the th k Temporary values ​​of generalized velocity and generalized displacement for +1 iteration step , ,Right now

[0042] here Represents the time step, parameters and It determines the integration accuracy and stability, and its function is to control... The change in acceleration over the time interval is often taken as in actual calculations. .

[0043] ⑤ Based on the fundamental knowledge of helicopter rotor blade structural dynamics and multibody system dynamics, , Calculate the temporary values ​​of the system mass matrix, generalized forces, and constraint Jacobian matrix, respectively. , and ; ⑥ Calculate generalized velocity and generalized displacement ; Since step ③ solves for DAEs containing acceleration constraint equations, the solution result theoretically only satisfies the acceleration constraint equations. Traditional direct integration methods will... , As The final values ​​of the generalized velocity and generalized displacement of the time-series system are inaccurate because they are only calculated through Newmark time integration. The numerical solution may not satisfy the velocity and displacement constraint equations. As the number of calculation steps increases, the error accumulates and eventually causes the numerical solution to diverge. To solve this problem, this invention employs the following correction strategy to calculate the generalized velocity. and generalized displacement The specific process is as follows: The combined equations of the multibody system dynamics differential equations and the generalized velocity constraint equations are expressed as follows:

[0044] In order to obtain the solution from the above simultaneous equations The Newmark integral formula is transformed as follows:

[0045] Substituting this into the above system of equations and rearranging, we get...

[0046] In this way, and For unknown variables, solving the above system of equations will yield the system's... k Generalized speed of +1 iteration step ; The dynamic differential equations of the multibody system are combined with the displacement constraint equations, expressed as follows:

[0047] Similarly, the Newmark integral formula is first transformed as follows:

[0048] Substituting these equations into the system of simultaneous equations, we can transform it into the following form.

[0049] by and Let be the unknown variable. Solving the above system of equations yields the first... k +1 iteration step generalized displacement .

[0050] Obviously, in solving the generalized velocity Generalized displacement Simultaneous equations and solving for generalized acceleration The simultaneous equations are identical in form, which simplifies the programming process. It only requires abstracting the solution process of the three equations into a function design, and then inputting specific parameters for different objects. Furthermore, due to the generalized speed... Generalized displacement The solution utilizes velocity constraint equations and displacement constraint equations, thus the numerical solution automatically satisfies the system's displacement and velocity constraints, perfectly resolving the problem of displacement and velocity violations.

[0051] ⑦ Determine if convergence has occurred. If the iteration converges, the iteration terminates, and we obtain... t +∆ t The system solution at time t, i.e. , , Otherwise, let k = k +1, return to ①.

[0052] Step 3: Refer to Step 2 to calculate the next time step, until the time expires.

[0053] For dynamics of flexible blade multibody systems (DAEs), the direct integration method, while simple to implement, suffers from violations of generalized velocity and generalized displacement rules. To address this issue without sacrificing computational accuracy, this invention employs the following two key strategies: ① Combine the differential dynamic equations of the system with the velocity constraint equations of the system to calculate the generalized velocity; ② Combine the system differential dynamics equations with the system displacement constraint equations to calculate the generalized displacement; Accurate correction calculations for the system's generalized velocity and generalized displacement were achieved.

[0054] Example 1 To verify the accuracy of the moderately deformable beam model, Dowell and Traybar conducted a series of static and dynamic experiments on uniform cantilever beams under large deformation conditions in 1974, known as the Princeton beam experiments. These experiments later became the standard for evaluating the correctness of nonlinear beam theory. The Princeton beam was made of 7075Al, with the following material parameters: elastic modulus 71 GPa, shear modulus 26.7 GPa, Poisson's ratio 0.33, and density 2800 kg / m³. 3 The beam has a rectangular cross-section, with a length of 0.508m, a width of 0.0127m, and a thickness of 0.003175m. Under different loads, the calculated and experimental values ​​of the static deformation of the Princeton beam have an error of less than 3%, verifying the effectiveness of the calculation method proposed in this invention.

[0055] Table 1 Comparison of calculated and experimental values ​​of static deformation of Princeton beam

[0056] Example 2 The rotor of a BO105-type helicopter has a rotor radius of 4.928m, no pre-twist angle or pre-cone angle, 4 blades, and a rotational speed of 424 rpm. Table 2 below compares the rotor blade dynamic characteristics calculated by this invention with those calculated in the literature. It can be seen that the calculation results of this invention are in high agreement with the calculation results in the published literature, with relative errors all within 0.5%, further verifying the effectiveness of the calculation method proposed in this invention.

[0057] Table 2 Comparison of Calculated Dynamic Characteristics of BO105 Rotor with Reference Results

[0058] Example 3 The beam axis is 1m long, fixed at one end, and the effect of gravity is negligible. Relevant material parameters: linear density 0.2 kg / m³, torsional stiffness 50 Nm. 2 Swing stiffness 50 Nm 2 The oscillation stiffness is 1000 Nm. 2 The tensile stiffness is 106 N. Static loads of 10 N to 150 N were applied to the ends of the beam shaft, and the resulting end deformations are as follows: Figure 1 As shown, the calculation results are in consistent with the reference results (software calculation), demonstrating that the new method remains stable and reliable even under large deformation conditions. Obviously, the embodiments described in the specific implementation details of this application are merely for the purpose of more clearly explaining the technical solutions in the specification, and are only a part of the embodiments of this application, and are not intended to limit this application. All other embodiments obtained by those skilled in the art based on the embodiments in the specific implementation details without creative effort should fall within the protection scope of this application.

Claims

1. A calculation method for the dynamic system of a flexible blade structure, the calculation steps of which are as follows: Step 1: Construct the differential-algebraic equation system (DAEs) for the dynamic system of the flexible rotor blade structure using multibody system dynamics theory and helicopter rotor blade structure dynamics theory, and set the time step ∆. t and convergence tolerance Given the initial conditions of the system and ; Step two, within a short time interval Inside, according to t The numerical solution of the system at time 1 , and Solve according to the following process t +∆ t Numerical solution at time , and : ①For the first k ( k =0,1,2,...) iteration steps ， Numerical solution and It has been determined by calculation in the previous iteration step; if k =0, then , ; ②Based on the fundamental knowledge of helicopter rotor blade structural dynamics and multibody system dynamics, with and Computational system mass array Generalized force and constrained Jacobian matrix ; ③ with and Given unknown variables, use numerical methods to solve the system of algebraic equations. k Generalized acceleration of +1 iteration step ; ④ Calculate the temporary values ​​of generalized velocity and generalized displacement. , ; ⑤ Based on the fundamental knowledge of helicopter rotor blade structural dynamics and multibody system dynamics, , Calculate the temporary values ​​of the system mass matrix, generalized forces, and constraint Jacobian matrix, respectively. , and ; ⑥ Calculate generalized velocity and generalized displacement ; ⑦ Determine if convergence has occurred. If convergence has occurred, terminate the iteration and obtain the result. t +∆ t The system solution at time t. , , Otherwise, let k = k +1, return ①; Step 3: Refer to Step 2 to calculate the next time step, until the time expires.

2. The method according to claim 1, characterized in that: In step one, the process of constructing the differential-algebraic equation system of the flexible blade structure dynamic system is as follows: First, the blade is divided into several straight beam segments. A floating coordinate system is established at the left end of the region pointing towards the blade root of each straight beam segment. By introducing 6 rigid body degrees of freedom and 14 elastic degrees of freedom, the large deformation motion of the flexible beam structure of the blade is decomposed into the superposition of large-scale rigid body motion under the floating coordinate system and small-to-medium elastic deformation motion relative to the floating coordinate system. Then, relative rotation and relative translational geometric constraints are used to coordinate the connection deformation between adjacent straight beam segments, thereby describing the complex nonlinear deformation behavior of the entire flexible blade. The 6 rigid body degrees of freedom include 3 translational degrees of freedom and 3 rotational degrees of freedom; the 14 elastic degrees of freedom include 3 flapping degrees of freedom, 3 oscillating degrees of freedom, 4 torsional degrees of freedom and 4 tensile degrees of freedom.

3. The method according to claim 2, characterized in that: In step one, the system of equations is as follows: in, t Indicates time; Represents the generalized coordinates of a multibody system, where , and These represent the translational, rotational, and elastic degrees of freedom of the system, respectively. ,in , , and These represent the degrees of freedom for stretching, oscillation, swinging, and torsion, respectively. , These are called generalized velocity and generalized acceleration, respectively, representing the generalized coordinates. First and second derivatives with respect to time; The following text is abbreviated as , represents the generalized mass matrix of a multibody system; The following text is abbreviated as , represents the constraint function vector of a multibody system, which describes the coordination relationship between translation, rotation and elastic deformation among the straight beam segments of the blade; The following text is abbreviated as , representing the constraint function vector For generalized coordinates The Jacobian matrix; The following text is abbreviated as , represents the generalized external force vector; Represents the constrained Lagrange multiplier vector; the superscript "T" indicates matrix transpose.

4. The method according to claim 3, characterized in that: In step 2, step ③, solve for the... k Generalized acceleration of +1 iteration step The process is as follows: constraint equations Regarding time t Differentiating, we obtain the generalized velocity constraint equation: here Represents the constraint function vector Regarding time t The partial derivatives; Applying the generalized velocity constraint equation to time t Differentiating again, we obtain the generalized acceleration constraint equation. here Representing vectors right Jacobian matrix, express Regarding time t The partial derivatives, Represents the constraint function vector Regarding time t The second-order partial derivatives; for ease of description, the right-hand side of the generalized acceleration constraint equation is represented by... express, ; The combined equations of dynamics of the multibody system and the generalized acceleration constraint equations are expressed as follows: Based on the k Numerical solution for each iteration step and ,by and For unknown variables, solving the above system of equations yields the system's... k Generalized acceleration of +1 iteration step , means as follows: 。 5. The method according to claim 4, characterized in that: In step 4 of step two, the generalized acceleration is... Perform Newmark integration to obtain the th k Temporary values ​​of generalized velocity and generalized displacement for +1 iteration step , The formula is as follows: here Represents the time step, parameters and It determines the integration accuracy and stability, and its function is to control... The change in acceleration over the time interval is often taken as in actual calculations. .

6. The method according to claim 5, characterized in that: In step 6 of step two, the process is as follows: The combined equations of dynamics of the multibody system and the generalized velocity constraint equations are expressed as follows: In order to obtain the solution from the above simultaneous equations The Newmark integral formula is transformed as follows: Substituting this into the above system of equations and rearranging, we get... by and For unknown variables, solving the above system of equations yields the system's... k Generalized speed of +1 iteration step ; The combined equations of dynamics and displacement constraints of the multibody system are expressed as follows: Similarly, the Newmark integral formula is first transformed as follows: Substituting these equations into the system of simultaneous equations, we can transform it into the following form. by and Let be the unknown variable. Solving the above system of equations yields the first... k +1 iteration step generalized displacement .

7. The method according to claim 6, characterized in that: In step 7 of step two, the convergence criterion is as follows: like If the expression converges, then the convergence is complete; otherwise, it does not converge. ε This is the convergence threshold.

8. The method according to claim 1, characterized in that: In step 2, step ①, the initial iteration step, k When =0, there is , .

9. An electronic device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that: When the processor executes the computer program, it implements the steps of the method according to any one of claims 1-8.

10. A computer-readable storage medium storing a calculator program thereon, characterized in that: When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1-8.

Images