Planetary roller screw dynamics modeling method based on generalized finite element method

Abstract

The invention discloses a planetary roller screw dynamic modeling method based on a generalized finite element method, and belongs to the field of dynamic characteristic analysis of planetary roller screws. The method comprises the following steps: carrying out unit division on a planetary roller screw finite element model; constructing a dynamic model of each type of unit; respectively establishing mass matrixes, stiffness matrixes, damping matrixes and load column vectors of the shaft section unit, the lead screw-roller contact side thread meshing unit, the nut-roller contact side thread meshing unit and the roller tooth-inner gear ring gear meshing unit; and assembling a system dynamics model to complete planet roller screw dynamics modeling. According to the method, the dynamic characteristics of the real physical model can be accurately analyzed, and an effective means is provided for analyzing the vibration response of the system.

Description

Technical Field

[0001] The present invention belongs to the field of dynamic characteristic analysis of planetary roller screws, and particularly relates to a dynamic modeling method of planetary roller screws based on the generalized finite element method. Background Art

[0002] A planetary roller screw mainly consists of a screw, rollers, a nut, an internal gear ring and a cage, and transmits force and motion through the thread meshing between multiple rollers and the screw and the nut. As a new type of actuator, it exhibits the meshing characteristics of multiple pairs (thread pair, gear pair), multiple bodies (multiple rollers are in contact and load-bearing with the screw and the nut simultaneously), and multiple points (the roller teeth and thread teeth of multiple rollers are involved in contact simultaneously), so it has advantages such as high precision, high load-bearing capacity, high rigidity and long service life.

[0003] At present, the electro-mechanical servo technology based on planetary roller screws has been widely applied in many fields such as aircraft rudder surfaces, aircraft brakes, rocket launcher erection, and servo joints of humanoid robots. In order to improve the dynamic comprehensive transmission performance of planetary roller screws, it is very necessary to carry out research on the dynamic contact characteristics of multiple-pair meshing. Establishing a complete and accurate dynamic model of the planetary roller screw system is the basis for its dynamic characteristic and vibration analysis.

[0004] However, the dynamic models in existing research are usually based on certain simplified assumptions, such as using the lumped mass method to model the system. Although this method has high computational efficiency, it is difficult to accurately reflect the complex connection relationships between components, especially for such precision transmission components as planetary roller screws. The lumped mass method discretizes a part into mass points and rigid connections, and cannot characterize the distribution characteristics of shaft segment units, resulting in inherent limitations in the model when analyzing vibration or contact nonlinear behaviors, further affecting the accuracy of dynamic characteristic prediction. Therefore, it is urgent to establish a more refined dynamic model to make up for the above defects. Summary of the Invention

[0005] Technical Problems to be Solved In order to avoid the deficiencies of the prior art, the present invention provides a dynamic modeling method of planetary roller screws based on the generalized finite element method, aiming to establish a dynamic modeling method that can simultaneously consider the synchronous meshing of the thread pair-gear pair and characterize the distribution characteristics of shaft segment units; this method uses the generalized finite element method to divide the system into different types of units (such as shaft segment units, thread meshing units, gear meshing units), and assembles its mass, stiffness, damping matrices and load vectors, and finally forms a system-level dynamic equation. The present invention can accurately analyze the dynamic characteristics of a real physical model and provide an effective means for analyzing the vibration response of the system.

[0006] The technical solution of the present invention is: a dynamic modeling method for planetary roller screws based on the generalized finite element method, and the specific steps are as follows: Perform element division on the finite element model of the planetary roller screw: discretize the screw, nut, rollers, and internal gear ring of the planetary roller screw into nodes, and divide them into shaft segment elements, screw-roller contact side thread meshing elements, nut-roller contact side thread meshing elements, and roller tooth-internal gear ring gear meshing elements; Construct the dynamic models of various types of elements: establish the mass matrix, stiffness matrix, damping matrix, and load column vector of the shaft segment element, screw-roller contact side thread meshing element, nut-roller contact side thread meshing element, and roller tooth-internal gear ring gear meshing element respectively; Assemble the system dynamic model: according to the mapping relationship between the element nodes and the global nodes, assemble each element matrix into the overall mass matrix, overall stiffness matrix, and overall damping matrix at the system level, that is, complete the dynamic modeling of the planetary roller screw.

[0007] A further technical solution of the present invention is: the shaft segment element is modeled using Timoshenko beam elements, and the screw-roller contact side thread meshing element, nut-roller contact side thread meshing element, and roller tooth-internal gear ring gear meshing element are modeled based on the relative displacement relationship of the contact points, and the nodes of each contact element coincide with the nodes of the shaft segment element.

[0008] A further technical solution of the present invention is: the expression of the dynamic model of the shaft segment element is as follows:

[0009] In the formula, q T The generalized coordinates of the two nodes of the shaft segment element; M T Is the consistent mass matrix of the shaft segment element; K T Is the stiffness matrix of the shaft segment element; C T Is the damping matrix of the shaft segment element; The mass matrix M T The expression is as follows:

[0010]

[0011]

[0012]

[0013]

[0014] In the formula, A Is the cross-sectional area of the shaft segment; ρ Is the material density;l is the length of the shaft section; I is the polar moment of inertia; I x is the yz section moment of inertia in the I y coordinate plane, xz section moment of inertia in the The stiffness matrix K T has the following expression:

[0015]

[0016]

[0017]

[0018]

[0019] In the formula, E is the elastic modulus of the material; G is the shear elastic modulus of the material; k is the correction factor introduced to consider that the actual shear strain and shear stress are not uniformly distributed, that is, the shear coefficient of the material; The damping matrix C T has the following expression:

[0020] In the formula, α 0, α 1 are the mass proportionality coefficient and the stiffness proportionality coefficient in Rayleigh damping.

[0021] A further technical solution of the present invention is that the relative displacement relationships of the screw-roller contact side thread meshing unit, the nut-roller contact side thread meshing unit, and the roller tooth-inner gear ring gear meshing unit are as follows: The relative displacement of the screw-roller contact side thread meshing unit is:

[0022]

[0023]

[0024]

[0025] In the formula, are respectively the vibration displacements of the screw in the direction; The torsional displacements of the lead screw in the direction are respectively; The vibration displacements of the rollers in the direction are respectively, The torsional displacements of the rollers in the direction are respectively; is the nominal radius of the lead screw; is the nominal radius of the roller; is the normal force in the plane's projection angle; is the normal force in the plane's projection angle; is the normal force in the plane's projection and the x axis's included angle; is the roller phase angle, j represents the roller's serial number, is the total number of rollers; n is the number of thread starts; represents the angle that the roller part coordinate system rotates around its z axis; The upper parts of the symbols " " and " " mean the lead screw rotates counterclockwise, and the lower parts mean the lead screw rotates clockwise; The relative displacement of the nut-roller contact side thread meshing unit is:

[0026]

[0027] In the formula, The vibration displacements of the nut in the direction are respectively, The torsional displacements of the nut in the direction are respectively; The vibration displacements of the rollers in the direction are respectively, The torsional displacements of the rollers in the direction are respectively; is the nominal radius of the nut; is the nominal radius of the roller; is the normal force in the plane's projection angle; is the normal force in the plane's projection angle; is the normal force in the plane's projection and the x axis's included angle; The relative displacement of the roller tooth - internal gear ring gear meshing unit is as follows:

[0028]

[0029] In the formula, are respectively the vibration displacements of the internal gear ring in the direction, are respectively the torsional displacements of the internal gear ring in the direction; are respectively the vibration displacements of the roller tooth in the direction, are respectively the torsional displacements of the roller tooth in the direction; is the pitch circle radius of the internal gear ring; is the pitch circle radius of the roller tooth; is the helix angle of the gear; is the pressure angle of the gear; is the sum of the gear pressure angle and the roller phase angle.

[0030] A further technical solution of the present invention is that the dynamic model expression of the screw - roller contact side thread meshing unit is as follows:

[0031] In the formula, is the mass matrix of the q th thread meshing unit on the screw - roller contact side; is the displacement matrix of the two nodes of the q th thread meshing unit on the screw - roller contact side; is the stiffness matrix of the q th thread meshing unit on the screw - roller contact side; is the damping matrix of the q th thread meshing unit on the screw - roller contact side; The mass matrix has the following expression:

[0032] In the formula, is the mass of the q th thread tooth of the screw; , and are respectively the moments of inertia of the q th thread tooth of the screw in the x , y , z directions; is the mass of the qQuality of a thread tooth; , and are the moments of inertia of the q th thread tooth of the roller in the x , y , z directions respectively; Displacement matrix has the following expression:

[0033] In the formula, , and are the vibration displacements of the q th thread tooth of the lead screw in the x , y , z directions respectively; , and are the torsional displacements of the q th thread tooth of the lead screw in the x , y , z directions respectively; , and are the vibration displacements of the q th thread tooth of the roller in the x , y , z directions respectively; , and are the torsional displacements of the q th thread tooth of the roller in the x , y , z directions respectively; Stiffness matrix has the following expression:

[0034] In the formula, is the thread meshing stiffness on the lead screw-roller contact side; V S is the projection matrix on the lead screw-roller side, and its expression is as follows:

[0035] Damping matrix has the following expression:

[0036] In the formula, is the thread meshing damping on the lead screw-roller contact side.

[0037] A further technical solution of the present invention is that the dynamic model expression of the nut-roller contact side thread engagement unit is as follows:

[0038] In the formula, is the mass matrix of the q th thread engagement unit on the nut-roller contact side; is the displacement matrix of the two nodes of the q th thread engagement unit on the nut-roller contact side; is the stiffness matrix of the q th thread engagement unit on the nut-roller contact side; is the damping matrix of the q th thread engagement unit on the nut-roller contact side; The mass matrix has the following expression:

[0039] In the formula, is the mass of the q th thread tooth of the nut; , and are the moments of inertia of the q th thread tooth of the nut in the x , y , z directions respectively; is the mass of the q th thread tooth of the roller; , and are the moments of inertia of the q th thread tooth of the roller in the x , y , z directions respectively; The displacement matrix has the following expression:

[0040] In the formula, , , are the vibration displacements of the q th thread tooth of the nut in the x , y , z directions respectively; , and are the vibration displacements of the q th thread tooth of the nut in the x , y ,z The torsional displacement in the 、 and directions are respectively the vibration displacements of the q th thread tooth of the roller in the x 、 y 、 z directions; 、 and directions are respectively the torsional displacements of the q th thread tooth of the roller in the x 、 y 、 z directions; The stiffness matrix has the following expression:

[0041] In the formula, is the thread meshing stiffness on the nut-roller contact side; V N is the projection matrix on the nut-roller side, and the expression is:

[0042] The damping matrix has the following expression:

[0043] In the formula, is the thread meshing damping on the nut-roller contact side.

[0044] A further technical solution of the present invention is that the dynamic model expression of the roller tooth-inner gear ring gear meshing unit is as follows:

[0045] In the formula, is the mass matrix of the q th gear meshing unit of the roller tooth-inner gear ring; is the displacement matrix of the two nodes of the q th gear meshing unit of the roller tooth-inner gear ring; is the stiffness matrix of the q th gear meshing unit of the roller tooth-inner gear ring; is the damping matrix of the q th gear meshing unit of the roller tooth-inner gear ring; is the column vector formed by the components of the normal impact force of the gear pair in each degree of freedom; The mass matrix has the following expression:

[0046] In the formula, is the mass of the q th tooth of the internal gear ring; , and are respectively the moments of inertia of the q th tooth of the internal gear ring in the x , y , z directions; is the mass of the q th roller tooth; , and are respectively the moments of inertia of the q th roller tooth in the x , y , z directions; The displacement matrix is expressed as follows:

[0047] In the formula, , , are respectively the vibration displacements of the q th tooth of the internal gear ring in the x , y , z directions; , , are respectively the torsional displacements of the x th tooth of the internal gear ring in the y , z directions; , , are respectively the vibration displacements of the q th roller tooth on the roller in the x , y , z directions; , , are respectively the torsional displacements of the q th roller tooth of the roller in the x , y , z directions; The stiffness matrix is expressed as follows:

[0048] In the formula, is the meshing stiffness of the gear pair on the contact side of the roller tooth - internal gear ring; V r is the projection matrix on the contact side of the roller tooth - internal gear ring, and the expression is:

[0049] Damping matrix The expression is as follows:

[0050] In the formula, is the meshing damping of the gear pair on the roller tooth-inner gear ring contact side; Load column vector The expression is as follows:

[0051] In the formula, is the normal impact force of the gear pair.

[0052] A further technical solution of the present invention is that when assembling each unit matrix, according to the corresponding relationship between the local number and the global number of the unit nodes, the sub-blocks of each unit stiffness matrix are superimposed on the corresponding positions of the global matrix; The expression for establishing the dynamic motion differential equation of the planetary roller screw is as follows:

[0053] In the formula, is the overall mass matrix of the planetary roller screw system; is the overall damping matrix of the planetary roller screw system; is the overall stiffness matrix of the planetary roller screw system; where,

[0054]

[0055]

[0056] In the formula, , , and , , are respectively the mass matrix, damping matrix, and stiffness matrix of the first shaft segment unit and the q th shaft segment unit; , , and , , are respectively the mass matrix, damping matrix, and stiffness matrix of the first and the q th screw-roller contact side thread meshing unit; , , and , , are respectively the mass matrix, damping matrix, and stiffness matrix of the 1st and the q th roller tooth - internal gear ring contact - side gear meshing unit; , , and , , are respectively the mass matrix, damping matrix, and stiffness matrix of the 1st and the q th nut - roller contact - side thread meshing unit.

[0057] A further technical solution of the present invention is that the dynamic response is solved by the numerical integration method of the dynamic motion differential equation of the planetary roller screw, which specifically includes: Dividing the total load into several load steps and applying them gradually; Updating the contact state after each step of calculation, and dynamically adjusting the stiffness matrix and the adjoint matrix; Judging the convergence by the residual. If not converged, re - iterate.

[0058] A planetary roller screw dynamics modeling system based on the generalized finite element method includes the following modules: Structure discretization module: used to discretize the screw, nut, roller, and internal gear ring of the planetary roller screw into global nodes, and divide them into shaft - segment units, screw - roller contact - side thread meshing units, nut - roller contact - side thread meshing units, and roller tooth - internal gear ring gear meshing units; Unit dynamics modeling module: based on the generalized finite element method, generate the mass matrix, stiffness matrix, damping matrix, and load vector of each unit; System matrix assembly module: according to the global number mapping relationship of the unit nodes, stack the unit matrices to the system - level overall matrix according to the degrees of freedom; Dynamic response solving module: based on the overall mass matrix, overall stiffness matrix, and overall damping matrix, establish the system dynamics differential equation and solve the dynamic response; Model verification module: compare the relative error between the calculated result of the dynamic contact force and the theoretical value to ensure that the model accuracy error is less than the threshold.

[0059] Beneficial effects The beneficial effects of the present invention are as follows: 1. The present invention conducts dynamic modeling on the planetary roller screw system based on the generalized finite element method. Through high-order shape functions and consistent mass matrices, it accurately describes the distribution characteristics of the shaft segment elements between the screw, roller, and nut nodes. The influence of distributed mass is coupled to all relevant nodes through the off-diagonal terms of the stiffness matrix and mass matrix, achieving high-precision prediction of the dynamic behavior of the shaft segment and avoiding model distortion caused by simplification in the traditional lumped mass method; it lays an important foundation for the accuracy of subsequent dynamic characteristic prediction and vibration analysis of the planetary roller screw system.

[0060] 2. In the present invention, the relative errors between the steady-state values of the calculation results of the dynamic model based on the generalized finite element method and the calculation results of the theoretical formula and the finite element simulation results are all within 5%, verifying the accuracy and effectiveness of the established dynamic model of the planetary roller screw, and it can effectively describe the dynamic characteristics of the planetary roller screw.

[0061] 3. For the dynamic model established in the present invention, as long as the basic information parameters of the system are input according to the regulations, various types of element matrices of the system can be generated. This modeling method is also suitable for the dynamic modeling of other types of planetary roller screw drive systems, effectively reducing repeated modeling. Compared with finite element simulation analysis, it can greatly shorten the calculation time and improve the analysis efficiency. Description of the Drawings Figure 1 It is the flow chart for establishing the dynamic model of the planetary roller screw system in the embodiment of the present invention; Figure 2 It is the force analysis diagram of the planetary roller screw drive system in the embodiment of the present invention; Figure 3 It is the finite element model diagram of the planetary roller screw drive system in the embodiment of the present invention; Figure 4 It is the shaft segment unit diagram of the planetary roller screw drive system in the embodiment of the present invention; Figure 5 It is the relative displacement relationship diagram of the contact position of the screw-roller contact side thread meshing unit in the embodiment of the present invention; Figure 6 It is the relative displacement relationship diagram of the contact position of the nut-roller contact side thread meshing unit in the embodiment of the present invention; Figure 7 It is the relative displacement relationship diagram of the contact position of the roller tooth-inner gear ring contact side gear pair meshing unit in the embodiment of the present invention; Figure 8 It is the schematic diagram of the stiffness matrix assembly of the planetary roller screw drive system in the embodiment of the present invention; Figure 9 It is the dynamic contact force curve diagram of the screw-roller contact side in the embodiment of the present invention; Figure 10It is the dynamic contact force curve diagram of the nut-roller contact side in the embodiment of the present invention; Figure 11 It is the comparison diagram of the dynamic contact force of the screw-roller contact side in the embodiment of the present invention; Figure 12 It is the comparison diagram of the dynamic contact force of the nut-roller contact side in the embodiment of the present invention. Detailed implementation manners

[0062] The embodiments described below with reference to the accompanying drawings are exemplary and are intended to explain the present invention and should not be construed as limiting the present invention.

[0063] At present, for the research on the dynamic research methods and dynamic models of planetary roller screws, there are the differential equations of component motion established by Jones M. H., Velinsky S. A., Lasky T. A. Dynamics of the planetary roller screw mechanism. Journal of Mechanisms and Robotics, 2016, 8(1): 014503 according to the Lagrangian method; the bond graph model established by Ma S. J., Zhang T., Liu G., et al. Bond Graph-based dynamic model of planetary roller screw mechanism with consideration of axial clearance and friction. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2018, 232(16): 2899-2911; the dynamic model established by Guo Jianan, He Peng, Huang Hongyan, etc. Analysis and experimental research on the dynamic characteristics of planetary roller screws. Journal of Propulsion Technology, 2018, 39(8): 841-848 using the discrete element method; the lumped mass pure torsion model and the bending-torsion coupling dynamic model established by Wu L. P., Ma S. J., Wan Q., et al. Dynamic modal of planetary roller screw mechanism with considering torsional degree of freedom. The 6th International Conference on Mechatronics and Mechanical Engineering, Japan, Osaka, Sep. 27-30, 2019; and the nonlinear six degrees of freedom dynamic model of planetary roller screw mechanism established by Fu X. J., Liu G., Tong R. T., et al.The six-degree-of-freedom nonlinear rigid body dynamics model established in Mechanism and Machine Theory, 2018, 119:22-36. The dynamics models in existing research are usually based on certain simplified assumptions. For example, the lumped mass method is used to model the system. Although this method has high computational efficiency, it is difficult to accurately reflect the complex connection relationships between components, especially for precision transmission components such as planetary roller screws. The lumped mass method discretizes a part into mass points and rigid connections, and cannot characterize the distribution characteristics of the shaft segment elements, resulting in inherent limitations in the model when analyzing vibration or contact nonlinear behaviors, further affecting the accuracy of dynamic characteristic prediction. Therefore, the present invention provides a dynamics modeling method for planetary roller screws based on the generalized finite element method, and the specific steps are as follows:. Perform element division on the finite element model of the planetary roller screw: Discretize the screw, nut, rollers, and internal gear ring of the planetary roller screw into nodes, and divide them into shaft segment elements, screw-roller contact side thread meshing elements, nut-roller contact side thread meshing elements, and roller tooth-internal gear ring gear meshing elements; Construct the dynamics models of various types of elements: Establish the mass matrix, stiffness matrix, damping matrix, and load column vector of the shaft segment elements, screw-roller contact side thread meshing elements, nut-roller contact side thread meshing elements, and roller tooth-internal gear ring gear meshing elements respectively; Assembly of the system dynamics model: According to the mapping relationship between the element nodes and the global nodes, assemble each element matrix into the overall mass matrix, overall stiffness matrix, and overall damping matrix at the system level, that is, complete the dynamics modeling of the planetary roller screw.

[0064] The present invention provides a dynamics modeling system for planetary roller screws based on the generalized finite element method, including the following modules: Structure discretization module: Used to discretize the screw, nut, rollers, and internal gear ring of the planetary roller screw into global nodes, and divide them into shaft segment elements, screw-roller contact side thread meshing elements, nut-roller contact side thread meshing elements, and roller tooth-internal gear ring gear meshing elements; Element dynamics modeling module: Based on the generalized finite element method, generate the mass matrix, stiffness matrix, damping matrix, and load vector of each element; System matrix assembly module: According to the global number mapping relationship of the element nodes, superimpose each element matrix to the system-level overall matrix according to the degrees of freedom; Dynamic response solution module: Based on the overall mass matrix, overall stiffness matrix, and overall damping matrix, establish the system dynamics differential equation and solve the dynamic response; Model verification module: Compare the relative error between the calculated result of the dynamic contact force and the theoretical value to ensure that the model accuracy error is less than the threshold.

[0065] The above technical solution will be further described below in conjunction with the accompanying drawings: In one embodiment, referring to Figures 1-8 , this embodiment provides a dynamic modeling method for planetary roller screws based on the generalized finite element method, including the following steps: S1: As shown in Figure 3 , according to the structural characteristics, divide the finite element model of the planetary roller screw; discretize the planetary roller screw into a series of nodes and form different types of elements, namely: shaft segment element, screw-roller contact side thread meshing element, nut-roller contact side thread meshing element, roller tooth-inner gear ring gear meshing element. Since the gear pair nodes and the thread pair nodes coincide with the shaft nodes, the shaft nodes in the sub-structures (screw, roller, nut, inner gear ring) are numbered, and through the position and numbering relationship between the sub-structures, the numbers of the shaft nodes and each element in the planetary roller screw are obtained; S2: As shown in Figure 4 , the dynamic modeling of the shaft segment element of the planetary roller screw: Use the Timoshenko beam element to establish the motion equation of the shaft segment element and complete the dynamic modeling of the shaft segment element. The matrix form of the dynamic motion differential equation of the shaft segment element of the planetary roller screw is shown in Equation (1): (1) In the formula, q T The generalized coordinates of the two nodes of the shaft segment, and q T = { x 1, y 1, z 1, θ x1 , θ y1 , θ z1 , x 2, y 2, z 2, θ x2 , θ y2 , θ z2} where: x q , y q , z q ( q = 1, 2) are the displacements of the node q along the local coordinate direction, θ xq , θ yq , θzq ( q = 1, 2) are the rotations of the cross-section at the node about the three coordinate axes; M q is the consistent mass matrix of the shaft element, as shown in Equation (2-6); K T is the stiffness matrix of the shaft element, which is represented here by the element stiffness matrix of the two-node Timoshenko space beam element, as shown in Equation (7-11); T is the damping matrix of the shaft element. The commonly used Rayleigh damping calculation in engineering is shown in Equation (12). C T

[0066] The mass matrix M T has the following expression: (2) (3) (4) (5) (6) In the formula, A is the cross-sectional area of the shaft segment; ρ is the material density; l is the length of the shaft segment; I is the polar moment of inertia; I x and I y are respectively the moments of inertia of the cross-section in the yz coordinate plane and the xz coordinate plane.

[0067] The stiffness matrix K T has the following expression: (7) (8) (9) (10) (11) In the formula, E is the elastic modulus of the material; G is the shear elastic modulus of the material; k is a correction factor introduced to consider that the actual shear strain and shear stress are not uniformly distributed, that is, the shear coefficient of the material. For a circular cross-section, it can be taken as k = 0.9.

[0068] (12)​ In the formula, α 0, α 1 are the mass proportionality coefficient and the stiffness proportionality coefficient in Rayleigh damping.

[0069] S3: As Figure 5 shown, dynamic modeling of the screw - roller contact - side thread meshing unit: According to the force - analysis diagram of the screw - roller contact side, the relative displacement of the screw - roller contact side is deduced to obtain the dynamic motion differential equation of the screw - roller contact - side thread meshing unit.

[0070] The relative displacement of the screw - roller contact side is:

[0071] (14) (15) (16) In the formula, are respectively the vibration displacements of the screw in the direction; are respectively the torsional displacements of the screw in the direction; are respectively the vibration displacements of the roller in the direction, are respectively the torsional displacements of the roller in the direction; is the nominal radius of the screw; is the nominal radius of the roller; is the normal force in the plane projection angle; is the normal force in the plane projection angle; is the projection of the normal force in the plane and the x axis included angle; is the roller phase angle, j represents the roller serial number, is the total number of rollers; n is the number of thread starts; represents the angle that the roller part coordinate system rotates around its z axis; The upper part of the symbols " " and " " means the screw rotates counter - clockwise, and the lower part means the screw rotates clockwise; The dynamic motion differential equation of the screw - roller side thread meshing unit is as shown in Equation (17), and its matrix form is as shown in Equation (18):

[0072] (18) In the formula, is the mass matrix of the q th thread meshing unit on the screw - roller contact side; is the displacement matrix of the two nodes of the q th thread meshing unit on the screw - roller contact side; is the stiffness matrix of the q th thread meshing unit on the screw - roller contact side; is the damping matrix of the q th thread meshing unit on the screw - roller contact side. Among them, the specific expressions of the mass matrix, displacement matrix, stiffness matrix and damping matrix are shown in Equation (19 - 26): (19) In the formula, is the mass of the q th thread tooth of the screw; , and are the moments of inertia of the q th thread tooth of the screw in the x , y , z directions respectively; is the mass of the q th thread tooth of the roller; , and are the moments of inertia of the q th thread tooth of the roller in the x , y , z directions respectively; (20) In the formula, , and are the vibration displacements of the q th thread tooth of the screw in the x , y , z directions respectively; , and are the torsional displacements of the q th thread tooth of the screw in the x , y , z directions respectively; , and are the vibration displacements of the q th thread tooth of the roller in thex , y , z vibrational displacement in the direction; , and are respectively the torsional displacements of the q th thread tooth of the roller in the x , y , z directions; (21) (22) (23)

[0073]

[0074] The stiffness matrix and damping matrix of the screw - roller contact - side thread meshing unit are complex in form. By carefully studying the structural composition of the matrix, a simple calculation form can be obtained: (27) (28) where is the screw - roller contact - side thread meshing stiffness, which can be calculated by Eqs. (29 - 31); is the screw - roller contact - side thread meshing damping, which can be calculated by Eq. (32); V S is the screw - roller side projection matrix, as shown in Eq. (33): (29) (30) (31) (32)

[0075] where is the normal contact force on the screw - roller contact side; is the normal contact deformation on the screw - roller contact side; is the axial contact force on the screw - roller contact side; is the flank angle on the screw - roller contact side; is the helix angle of the roller; k ( e ) is the first - kind elliptic integral, eis the eccentricity; is the equivalent elastic modulus, as shown in Equation (34); is related to the principal curvature function F ( ρ ) and can be obtained by looking up a table. The principal curvature function is calculated through Equations (35 - 37); η is the meshing damping ratio coefficient, taking values from 0.03 to 0.17; and are the masses of the two meshing nodes on the screw - roller contact side respectively.

[0076] (34) (35)

[0077] (37) In the formula, E S and E R are the elastic moduli of the screw and the roller respectively; μ S and μ R are the Poisson's ratios of the screw and the roller respectively; R is the equivalent spherical radius; d S is the pitch diameter of the screw; d R is the pitch diameter of the roller; α is the thread flank angle.

[0078] S4: As shown in Figure 6 , dynamic modeling of the thread meshing unit on the nut - roller contact side: According to the force analysis diagram of the nut - roller contact side, the relative displacement of the nut - roller contact side is deduced, and the dynamic motion differential equation of the thread meshing unit on the nut - roller contact side is obtained. The relative displacement of the nut - roller contact side is:

[0079] (39) In the formula, are the vibration displacements of the nut in the direction respectively, are the torsional displacements of the nut in the direction respectively; are the vibration displacements of the roller in the direction respectively, are the torsional displacements of the roller in the direction respectively; is the nominal radius of the nut; is the nominal radius of the roller; is the normal force at the projection angle on the surface; is the normal force at the projection angle on the surface; is the normal force at the projection on the surface and the x angle between the axis; The dynamic differential equation of motion of the nut-roller side thread engagement unit is shown in Equation (40), and its matrix form is shown in Equation (41):

[0080] (41) In the formula, is the mass matrix of the q th thread engagement unit on the nut-roller contact side; is the displacement matrix of the two nodes of the thread engagement unit on the nut-roller contact side; is the stiffness matrix of the thread engagement unit on the nut-roller contact side; is the damping matrix of the thread engagement unit on the nut-roller contact side. Among them, the specific expressions of the mass matrix, displacement matrix, stiffness matrix, and damping matrix are shown in Equations (42-45):

[0081] In the formula, is the mass of the q th thread tooth of the nut; , and are the moments of inertia of the q th thread tooth of the nut in the x , y , z directions respectively; is the mass of the q th thread tooth of the roller; , and are the moments of inertia of the q th thread tooth of the roller in the x , y , z directions respectively;

[0082] In the formula, , , are the q th thread teeth of the nut in thex , y , z Vibrational displacement in the , and directions respectively represent the torsional displacements of the q th thread tooth of the nut in the x , y , z directions; , and respectively represent the vibrational displacements of the q th thread tooth of the roller in the x , y , z directions; , and respectively represent the torsional displacements of the q th thread tooth of the roller in the x , y , z directions; (44) (45) In the formula, is the thread meshing stiffness on the nut-roller contact side and can be calculated by formula (46); is the thread meshing damping on the nut-roller contact side and can be calculated by formula (47); V N is the projection matrix on the nut-roller side, as shown in formula (48): (46) (47)

[0083] In the formula, and are respectively the normal contact force and normal contact deformation on the nut-roller contact side, and their calculation methods are the same as those on the screw-roller side; is the axial contact force on the nut-roller contact side; α N is the flank angle on the nut-roller contact side; and are respectively the masses of the two meshing nodes on the nut-roller contact side.

[0084] S5: As Figure 7As shown in the figure, the dynamic modeling of the roller tooth-inner gear ring gear meshing unit: According to the force analysis diagram of the roller tooth-inner gear ring, the relative displacement of the roller tooth-inner gear ring is deduced, and the dynamic motion differential equation of the roller tooth-inner gear ring gear pair meshing unit is obtained. The relative displacement of the roller tooth-inner gear ring is as follows:

[0085] (50) In the formula, are respectively the vibration displacements of the inner gear ring in the direction, are respectively the torsional displacements of the inner gear ring in the direction; are respectively the vibration displacements of the roller tooth in the direction, are respectively the torsional displacements of the roller tooth in the direction; is the pitch circle radius of the inner gear ring; is the pitch circle radius of the roller tooth; is the helix angle of the gear; is the pressure angle of the gear; is the sum of the gear pressure angle and the roller phase angle.

[0086] The dynamic motion differential equation of the roller tooth-inner gear ring gear meshing unit is shown in Equation (51), and its matrix form is shown in Equation (52): (52) In the formula, is the mass matrix of the q th gear meshing unit of the roller tooth-inner gear ring; is the displacement matrix of the two nodes of the q th gear meshing unit of the roller tooth-inner gear ring; is the stiffness matrix of the q th gear meshing unit of the roller tooth-inner gear ring; is the damping matrix of the q th gear meshing unit of the roller tooth-inner gear ring; is the column vector formed by the components of the normal impact force of the gear pair in each degree of freedom. Among them, the specific expressions of the mass matrix, displacement matrix, stiffness matrix and damping matrix are shown in Equations (53-57): (53) In the formula, is the mass of the q th tooth of the inner gear ring; , and are the moments of inertia of the q th tooth of the internal gear ring in the x , y , z directions; is the mass of the q th roller tooth; , and are the moments of inertia of the q th roller tooth of the roller in the x , y , z directions; (54) In the formula, , , are the vibration displacements of the q th tooth of the internal gear ring in the x, y, and z directions; , , are the torsional displacements of the q th tooth of the internal gear ring in the x, y, and z directions; , , are the vibration displacements of the q th roller tooth of the roller in the x, y, and z directions; , , are the torsional displacements of the q th roller tooth of the roller in the x, y, and z directions; (55) (56) (57) In the formula, is the average meshing stiffness of the gear pair on the contact side of the roller tooth - internal gear ring, and its calculation method can refer to the literature [Compilation of Parts and Related Standards, Gear and Gear Transmission Volume (Part 1). Beijing: China Standards Press, 2012.]; is the meshing damping of the gear pair on the contact side of the roller tooth - internal gear ring, which can be calculated by formula (58); is the normal impact force of the gear pair; V r is the projection matrix on the contact side of the roller tooth - internal gear ring, as shown in formula (59): (58)

[0087] In the formula, and are the masses of the two meshing nodes on the contact side of the roller teeth and the internal gear ring respectively.

[0088] S6: As Figure 8 shown, the dynamic model of the planetary roller screw system is assembled based on the element numbers with the mass, stiffness, damping matrices and load column vectors of the shaft segment elements, the screw-roller contact side thread pair meshing elements, the nut-roller contact side thread pair meshing elements, and the roller tooth-internal gear ring gear meshing elements. The assembly process of the overall stiffness matrix of the system is the same as that of the overall stiffness matrix in the finite element analysis of structural mechanics, that is, according to the corresponding relationship between the local numbers of the nodes of each element and the overall numbers of the system nodes, the sub-matrices corresponding to each degree of freedom of the element matrix are successively superimposed on the corresponding positions of the overall matrix. The matrix form of the dynamic differential equation of the planetary roller screw system is shown in Equation (60): (60) In the formula, is the overall mass matrix of the planetary roller screw system; is the overall damping matrix of the planetary roller screw system; is the overall stiffness matrix of the planetary roller screw system; its specific expression is shown in Equations (61 - 63).

[0089] (61) (62) (63) In the formula, , , and , , are the mass matrix, damping matrix, and stiffness matrix of the first shaft segment element and the q th shaft segment element respectively; , , and , , are the mass matrix, damping matrix, and stiffness matrix of the first and the q th screw-roller contact side thread meshing elements respectively; , , and , , are the mass matrix, damping matrix, and stiffness matrix of the first and the q th roller tooth-internal gear ring gear meshing elements respectively; , , and , , are the mass matrix, damping matrix, and stiffness matrix of the first and the q th nut-roller contact side thread meshing unit respectively; S7: To verify the accuracy of the planetary roller screw dynamic model based on the generalized finite element method proposed in this paper, two methods are used for comparative verification: First, calculate the static contact force through theoretical formulas. Second, the literature [Dynamic Contact Load Characteristics of Planetary Roller Screw Thread Pair - Gear Pair Synchronous Meshing, Journal of Xi'an Jiaotong University, 2022, 56(10), 11 - 21.] uses finite element simulation to obtain the steady-state value of the dynamic contact force. Finally, the calculation results of these two methods are compared and analyzed with the steady-state value of the dynamic contact force calculated by the dynamic model in this paper.

[0090] According to the structural parameters and load conditions of the planetary roller screw given in Table 1, the stiffness parameters of the planetary roller screw are calculated as shown in Table 2. During the calculation process, the number of rollers is set to 3, and the system is discretized into 18 nodes. Thus, the established dynamic model has 108 degrees of freedom. Based on this model, the time-domain responses of the dynamic contact forces on the screw-roller contact side and nut-roller contact side in the planetary roller screw system are solved, and their variation laws are respectively as Figure 9 and Figure 10 shown.

[0091] On this basis, the static contact forces on the screw-roller contact side and nut-roller contact side are calculated using equations (29) - (31) and equation (46) respectively, and compared with the steady-state value of the dynamic contact force of the thread pair solved by the model in this paper. The results are shown in Table 3. Comparative analysis shows that the relative error between the steady-state value of the dynamic contact force calculated by the model of the present invention and the static contact force does not exceed 2%, verifying the correctness of the proposed dynamic model.

[0092] According to the research of the literature [Dynamic Contact Load Characteristics of Planetary Roller Screw Thread Pair - Gear Pair Synchronous Meshing, Journal of Xi'an Jiaotong University, 2022, 56(10), 11 - 21], under the working conditions of 3 rollers, an axial load of 5000N, and a screw rotation speed of 30 rad / s, the planetary roller screw is simulated and analyzed by the finite element method. On the premise of keeping the same structural parameters, material parameters, and working conditions, the dynamic model established in the present invention numerically solves the dynamic contact forces on the screw-roller contact side and nut-roller contact side, and the result comparison diagram is as Figure 11 , Figure 12As shown. The comparative analysis of the results shown in Table 4 indicates that the relative errors between the steady-state values of the dynamic contact forces of the screw pair calculated by the kinetic model of the present invention and the finite element simulation results in the literature are 0.65% and 0.33% respectively, and the two are highly consistent, further verifying the correctness and reliability of the established kinetic model.

[0093] In summary, the kinetic modeling method of the planetary roller screw proposed based on the generalized finite element idea in the present invention can accurately describe the distribution characteristics of the shaft segment units, accurately reflect the complex connection relationships between components, and make its kinetic model closer to the physical model. At the same time, this modeling method effectively reduces repeated modeling. Compared with the analysis using finite element software, it can greatly shorten the calculation time and improve the analysis efficiency.

[0094] Table 1 Structural parameters of the planetary roller screw

[0095] Table 2 Static contact forces and stiffness parameters of the planetary roller screw

[0096] Table 3 Comparison of the results of static contact forces and steady-state values of dynamic contact forces

[0097] Table 4 Comparison of the finite element simulation results in the literature and the steady-state values of dynamic contact forces

[0098] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limitations of the present invention. Those of ordinary skill in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention without departing from the principles and purposes of the present invention.

Claims

1. A dynamic modeling method for planetary roller screws based on the generalized finite element method, characterized in that The specific steps are as follows: Perform element division on the finite element model of the planetary roller screw: Discretize the screw, nut, rollers, and internal gear ring of the planetary roller screw into nodes, and divide them into shaft segment elements, screw-roller contact side thread meshing elements, nut-roller contact side thread meshing elements, and roller tooth-internal gear ring gear meshing elements; Construct the dynamic models of various types of elements: Establish the mass matrix, stiffness matrix, damping matrix, and load column vector of the shaft segment element, screw-roller contact side thread meshing element, nut-roller contact side thread meshing element, and roller tooth-internal gear ring gear meshing element respectively; System dynamic model assembly: According to the mapping relationship between the element nodes and the global nodes, assemble each element matrix into the overall mass matrix, overall stiffness matrix, and overall damping matrix at the system level, that is, complete the dynamic modeling of the planetary roller screw.

2. The dynamic modeling method of planetary roller screw based on the generalized finite element method according to claim 1, characterized in that: The shaft segment element is modeled using Timoshenko beam elements, and the screw-roller contact side thread meshing element, nut-roller contact side thread meshing element, and roller tooth-internal gear ring gear meshing element are modeled based on the relative displacement relationship of the contact points, and the nodes of each contact element coincide with the nodes of the shaft segment element.

3. The dynamic modeling method of a planetary roller screw based on the generalized finite element method according to claim 2, characterized in that: The expression of the dynamic model of the shaft segment element is as follows: where, q T generalized coordinates of two nodes of the shaft segment element; M T is the consistent mass matrix of the shaft segment element; K T is the stiffness matrix of the shaft segment element; C T is the damping matrix of the shaft segment element; Mass matrix M T The expression is as follows: Wherein, A is the cross-sectional area of the shaft section; ρ is the material density; l is the length of the shaft section; I is the polar moment of inertia; I x is the moment of inertia of the cross-section in the yz coordinate plane, I y is the moment of inertia of the cross-section in the xz coordinate plane; Stiffness matrix K T The expression is as follows: In the formula, E is the elastic modulus of the material; G is the shear elastic modulus of the material; k is a correction factor introduced to account for the fact that the actual shear strain and shear stress are not uniformly distributed, i.e., the shear coefficient of the material; Damping matrix C T The expression is as follows: In the formula, α 0, α 1 are the mass proportionality coefficient and the stiffness proportionality coefficient in Rayleigh damping.

4. A dynamic modeling method for planetary roller screws based on the generalized finite element method according to claim 3, characterized in that: The relative displacement relationships of the screw-roller contact side thread meshing element, nut-roller contact side thread meshing element, and roller tooth-internal gear ring gear meshing element are as follows: The relative displacement of the screw-roller contact side thread meshing element is: Wherein, are respectively the vibration displacements of the lead screw in the direction; are respectively the torsional displacements of the lead screw in the direction; are respectively the vibration displacements of the roller in the direction, are respectively the torsional displacements of the roller in the direction; is the nominal radius of the lead screw; is the nominal radius of the roller; is the normal force in the plane projection angle; is the normal force in the plane projection angle; is the angle between the projection of the normal force in the plane and the x axis; is the roller phase angle, j represents the serial number of the roller, is the total number of rollers; n is the number of thread starts; represents the angle by which the roller part coordinate system rotates around its z axis; The upper part of the symbols " " and " " indicates counterclockwise rotation of the lead screw, and the lower part indicates clockwise rotation of the lead screw; The relative displacement of the nut-roller contact side thread meshing element is: Wherein, are respectively the vibration displacements of the nut in the direction, are respectively the torsional displacements of the nut in the direction; are respectively the vibration displacements of the roller in the direction, are respectively the torsional displacements of the roller in the direction; is the nominal radius of the nut; is the nominal radius of the roller; is the normal force in the plane projection angle; is the normal force in the plane projection angle; is the angle between the projection of the normal force in the plane and the x axis; The relative displacement of the roller tooth-internal gear ring gear meshing element is: Wherein, are respectively the vibration displacements of the internal gear ring in the direction, are respectively the torsional displacements of the internal gear ring in the direction; are respectively the vibration displacements of the roller teeth in the direction, are respectively the torsional displacements of the roller teeth in the direction; is the pitch circle radius of the internal gear ring; is the pitch circle radius of the roller teeth; is the helix angle of the gear; is the pressure angle of the gear; is the sum of the gear pressure angle and the roller phase angle.

5. A dynamic modeling method for planetary roller screws based on the generalized finite element method according to claim 4, characterized in that: The expression of the dynamic model of the screw-roller contact side thread meshing element is as follows: In the formula, is the mass matrix of the q th thread engagement unit on the screw-roller contact side; is the displacement matrix of the two nodes of the q th thread engagement unit on the screw-roller contact side; is the stiffness matrix of the q th thread engagement unit on the screw-roller contact side; is the damping matrix of the q th thread engagement unit on the screw-roller contact side; Mass matrix The expression is as follows: In the formula, is the mass of the q th thread tooth of the lead screw; , and are respectively the moments of inertia of the q th thread tooth of the lead screw in the x , y , z directions; is the mass of the q th thread tooth of the roller; , and are respectively the moments of inertia of the q th thread tooth of the roller in the x , y , z directions; Displacement matrix The expression is as follows: Wherein, , and are respectively the vibration displacements of the q th thread tooth of the lead screw in the x , y , z directions; , and are respectively the torsional displacements of the q th thread tooth of the lead screw in the x , y , z directions; , and are respectively the vibration displacements of the q th thread tooth of the roller in the x , y , z directions; , and are respectively the torsional displacements of the q th thread tooth of the roller in the x , y , z directions; Stiffness matrix The expression is as follows: In the formula, is the screw-roller contact side thread meshing stiffness; V S is the screw-roller side projection matrix, and the expression is as follows: Damping matrix The expression is as follows: In the formula, is the thread meshing damping on the lead screw-roller contact side.

6. The dynamic modeling method of planetary roller screw based on the generalized finite element method according to claim 5, characterized in that: The expression of the dynamic model of the nut-roller contact side thread meshing element is as follows: In the formula, is the mass matrix of the q th thread engagement unit on the nut-roller contact side; is the displacement matrix of the two nodes of the q th thread engagement unit on the nut-roller contact side; is the stiffness matrix of the q th thread engagement unit on the nut-roller contact side; is the damping matrix of the q th thread engagement unit on the nut-roller contact side; Mass matrix The expression is as follows: In the formula, is the mass of the q th thread tooth of the nut; , and are the moments of inertia of the q th thread tooth of the nut in the x , y , z directions respectively; is the mass of the q th thread tooth of the roller; , and are the moments of inertia of the q th thread tooth of the roller in the x , y , z directions respectively; Displacement matrix The expression is as follows: In the formula, , , are respectively the vibration displacements of the q th thread tooth of the nut in the x , y , z directions; , and are respectively the torsional displacements of the q th thread tooth of the nut in the x , y , z directions; , and are respectively the vibration displacements of the q th thread tooth of the roller in the x , y , z directions; , and are respectively the torsional displacements of the q th thread tooth of the roller in the x , y , z directions; Stiffness matrix The expression is as follows: In the formula, is the thread meshing stiffness on the nut-roller contact side; V N is the projection matrix on the nut-roller side, and the expression is: Damping matrix The expression is as follows: In the formula, is the thread meshing damping on the nut-roller contact side.

7. A dynamic modeling method for planetary roller screws based on the generalized finite element method according to claim 6, characterized in that: The expression of the dynamic model of the roller tooth-internal gear ring gear meshing element is as follows: In the formula, is the mass matrix of the q th gear meshing unit of the roller gear - internal gear ring; is the displacement matrix of the two nodes of the q th gear meshing unit of the roller gear - internal gear ring; is the stiffness matrix of the q th gear meshing unit of the roller gear - internal gear ring; is the damping matrix of the q th gear meshing unit of the roller gear - internal gear ring; is the column vector formed by the components of the normal impact force of the gear pair in each degree of freedom; Mass matrix The expression is as follows: Wherein, is the mass of the q th tooth of the internal gear ring; , and are respectively the moments of inertia of the q th tooth of the internal gear ring in the x , y , z directions; is the mass of the q th roller tooth; , and are respectively the moments of inertia of the q th roller tooth in the x , y , z directions; Displacement matrix The expression is as follows: Wherein, 、 、 are respectively the vibration displacements of the q th tooth of the internal gear ring in the x 、 y 、 z directions; 、 、 are respectively the torsional displacements of the q th tooth of the internal gear ring in the x 、 y 、 z directions; 、 、 are respectively the vibration displacements of the q th roller tooth in the x 、 y 、 z directions; 、 、 are respectively the torsional displacements of the q th roller tooth in the x 、 y 、 z directions; Stiffness matrix The expression is as follows: In the formula, is the meshing stiffness of the roller tooth-inner gear ring contact side gear pair; V r is the projection matrix of the roller tooth-inner gear ring contact side, and the expression is: Damping matrix The expression is as follows: In the formula, is the meshing damping of the roller gear-inner gear ring contact side gear pair; Load column vector The expression is as follows: In the formula, is the normal impact force of the gear pair.

8. A dynamic modeling method for planetary roller screws based on the generalized finite element method according to claim 7, characterized in that: When assembling each element matrix, according to the corresponding relationship between the local number and the global number of the element nodes, superimpose the sub-blocks of each element stiffness matrix to the corresponding positions of the global matrix; The expression for establishing the dynamic motion differential equation of the planetary roller screw is as follows: In the formula, is the overall mass matrix of the planetary roller screw system; is the overall damping matrix of the planetary roller screw system; is the overall stiffness matrix of the planetary roller screw system; where wherein and and and and and are respectively the mass matrix, damping matrix, and stiffness matrix of the first shaft segment unit and the q th shaft segment unit; and and and and and are respectively the mass matrix, damping matrix, and stiffness matrix of the first and the q th screw - roller contact side thread meshing unit; and and and and and are respectively the mass matrix, damping matrix, and stiffness matrix of the first and the q th roller tooth - internal gear ring gear meshing unit; and and and and and are respectively the mass matrix, damping matrix, and stiffness matrix of the first and the q th nut - roller contact side thread meshing unit.

9. A dynamic modeling method for planetary roller screws based on the generalized finite element method according to claim 8, characterized in that: The numerical integration method of the dynamic motion differential equation of the planetary roller screw is used to solve the dynamic response, specifically including: Divide the total load into several load steps and apply them gradually; Update the contact state after each step of calculation, and dynamically adjust the stiffness matrix and the adjoint matrix; Judge the convergence through the residual. If it does not converge, re-iterate.

10. A planetary roller screw dynamic modeling system based on the generalized finite element method, which is used to implement the planetary roller screw dynamic modeling method based on the generalized finite element method according to any one of claims 1-9; characterized in that, It includes the following modules: Structure discretization module: Used to discretize the screw, nut, rollers, and internal gear ring of the planetary roller screw into global nodes, and divide them into shaft segment elements, screw-roller contact side thread meshing elements, nut-roller contact side thread meshing elements, and roller tooth-internal gear ring gear meshing elements; Element dynamic modeling module: Based on the generalized finite element method, generate the mass matrix, stiffness matrix, damping matrix, and load vector of each element; System matrix assembly module: According to the global number mapping relationship of the element nodes, superimpose each element matrix to the system-level overall matrix according to the degrees of freedom; Dynamic response solution module: Based on the global mass matrix, global stiffness matrix and global damping matrix, establish the system dynamics differential equation and solve the dynamic response; Model verification module: Compare the relative error between the calculated result of the dynamic contact force and the theoretical value to ensure that the model accuracy error is less than the threshold.

Images