A method for solving internal force of articulated constraint based on newton-euler equation

Abstract

The application discloses a hinged constraint internal force solving method based on Newton-Euler equation, relates to the technical field of structural dynamics numerical simulation analysis in engineering, and aims to solve the problems of insufficient constraint precision and unstable constraint calculation in the prior art in finite element analysis. The method comprises the following steps: simplifying a hinged structure, constructing an equivalent constraint equation, solving a rigid body dynamics equation, solving the constraint equation, and calculating constraint internal force. The method effectively improves the calculation precision and significantly enhances the calculation stability. The hinged force calculation method proposed by the application theoretically avoids the introduction of artificial parameters, so that the calculation result is not affected by parameter selection.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of numerical simulation analysis of structural dynamics in engineering, and particularly relates to a hinge constraint internal force solving method based on Newton-Euler equation. BACKGROUND

[0002] In explicit dynamics analysis of mechanical structure impact, vehicle collision, etc., there are mainly two methods for existing technology to process hinged structures: a fine modeling method based on physical shape and a simplified modeling method based on constraint equations. The fine modeling method based on physical shape divides the grid according to the real geometric structure of the hinge, and defines the contact relationship between the connected parts to transfer force and motion. This method has high calculation accuracy, but the calculation cost is significantly increased, and the calculation stability is low. The simplified modeling method based on constraint equations pays less attention to the local deformation details of the hinge structure itself in many dynamic conditions, such as whole vehicle collision analysis. In order to improve the calculation efficiency and the stability of the local position, a simplified equivalent constraint model is often used in modeling. The behavior of the constraint model can be represented by a constraint equation, and the solving method is the key. The commonly used constraint equation solving methods currently mainly include the Lagrange multiplier method and the penalty function method.

[0003] The Lagrange multiplier method solves the motion parameters and internal forces by introducing the Lagrange multiplier to solve the constraint equations related to the hinge element and the system equations. However, in the explicit finite element analysis process, this method destroys the inherent advantage of the explicit algorithm "no need to solve the simultaneous equations", significantly increases the calculation complexity, and needs to occupy more storage space. Therefore, this method is generally not suitable for explicit dynamics analysis. The penalty function method realizes the constraint by introducing a penalty term in the constraint equation deviation. When the constraint condition is not strictly satisfied, a penalty force is applied on the related node to force the constraint to approach satisfaction. This method does not need to solve the simultaneous equations, the calculation logic is simple, and it is easy to implement in parallel. Therefore, the mainstream explicit dynamics analysis software, such as LS-DYNA, ABAQUS, Radioss, generally uses this method to process the hinge constraint.

[0004] At present, when analyzing the dynamic response of the hinged structure in the explicit dynamics program, the penalty function method is the most widely used method. However, this method has obvious technical defects when processing the hinge constraint: first, the penalty function method relies on the penalty factor for constraint correction, resulting in insufficient constraint accuracy of the hinge element; second, in the explicit integration process, the penalty factor is highly sensitive to the value, which is easy to cause solution instability and even program interruption. It is a very difficult task to select a reasonable penalty factor, and different values are often needed for different conditions. The above problems significantly restrict the application effect of the simplified hinge model in high-precision and large-scale whole vehicle collision simulation. Therefore, the insufficient constraint accuracy and the unstable constraint calculation in finite element analysis have become a scientific problem to be solved. SUMMARY

[0005] The application aims to provide a hinge constraint internal force solving method based on Newton-Euler equation to solve the problems of insufficient constraint accuracy and unstable constraint calculation in the prior art in finite element analysis.

[0006] The application provides a hinge constraint internal force solving method based on Newton-Euler equation, comprising:

[0007] Step one, simplifying the hinge structure and constructing equivalent constraint equations; wherein the hinge structure comprises a spherical hinge and a revolute pair, the constraint characteristics of the spherical hinge are that there is a fixed spatial coincidence relationship between two mechanical components, a specific node N1 on component A and a corresponding node N2 on component B always coincide in space, and the equivalent constraint equation of the spherical hinge is expressed by the coordinates of the two nodes:

[0008] (1);

[0009] Wherein x1 is the coordinate of node N1, and x2 is the coordinate of node N2;

[0010] The equivalent constraint equation of the revolute pair is expressed according to the position vectors of two groups of four constraint nodes, and one group of coincident constraint nodes N3 and N4 is added compared with the spherical hinge, the first group of constraint equations is the same as formula (1), and the second group of constraint equations is:

[0011] (2);

[0012] Wherein x3 is the coordinate of node N3, and x4 is the coordinate of node N4;

[0013] Step two, solving the rigid body dynamics equation; the dynamic response of a single rigid body under the action of external force follows Newton equation and Euler equation, which respectively correspond to the translational freedom and rotational freedom of the rigid body, after the translational acceleration and rotational acceleration of the rigid body are calculated according to the mathematical expressions of the translational freedom and rotational freedom, the central difference algorithm is used to update the displacement increment of the translation and rotation of the rigid body in the explicit time integration framework; the transformation matrix is calculated based on the rotational displacement increment of the rigid body and the inertia matrix at the time step , the inertia tensor of the rigid body at the time step is updated, and the coordinates of the dependent nodes on the rigid body are updated;

[0014] Step three, solving the constraint equation and calculating the constraint internal force.

[0015] Further, in step two, the translational and rotational degrees of freedom of the rigid body are mathematically expressed as:

[0016] (3)

[0017] (4)

[0018] where, m represents the mass of the rigid body, a represents the acceleration of the center of mass of the rigid body, F represents the resultant external force; I represents the inertia tensor of the rigid body, w represents the angular velocity of the rigid body, a represents the angular acceleration of the rigid body, M represents the resultant external moment.

[0019] Further, the mass m of the rigid body is the sum of the masses of the nodes in the dependent node set of the rigid body :

[0020] (5)

[0021] The inertia tensor I of the rigid body is given by the following formula:

[0022] (6)

[0023] where, r represents the position of the node.

[0024] Further, in step two, under the explicit time integration framework, the central difference algorithm is used to update the translational and rotational displacement increments of the rigid body. The update formula for the translational coordinates of the center of mass of the rigid body is as follows:

[0025] (7)

[0026] where, r represents the position of the center of mass of the rigid body at time step n, r represents the position of the center of mass of the rigid body at time step n+1, d represents the translational displacement increment.

[0027] Further, in step two, based on the Hughes-Winget algorithm, the rotational displacement increment of the rigid body is used to calculate the transformation matrix :

[0028] (8) ​​

[0029] (9)

[0030] where, is the identity matrix, is the intermediate matrix, is an arbitrary vector.

[0031] Further, in step two, based on the transformation matrix and the inertia matrix of the rigid body at time step , the inertia tensor of the rigid body at time step is updated based on equation (10) and the coordinates of the slave nodes on the rigid body are updated by equation (11):

[0032] (10)

[0033] (11)

[0034] where, and denote the coordinates of node i at time step n and n+1, respectively, denotes the translational displacement increment of the rigid body, denotes the position vector of node i relative to the center of mass of the rigid body.

[0035] Further, in step three, the constraint equation of the spherical joint is equivalently represented by the following equation:

[0036] (12)

[0037] where, denotes the node displacement at time step n+1, denotes the node displacement at time step n, denotes the position vector of relative to ;

[0038] (13)

[0039] The dynamic equation of the rigid body has been given in equations (3) and (4) and is expressed in the following form:

[0040] (14)

[0041] (15)

[0042] where, is the resultant external force, ​​denotes the mass of the rigid body, a is the translational acceleration of the rigid body, is the angular acceleration of the rigid body, denotes the inertia tensor matrix of the rigid body in the global coordinate system, which is updated at each time step according to equation (10), denotes the equivalent external moment, while denotes the resultant moment acting on the rigid body, and the relationship between them is:

[0043] (16).

[0044] Further, in step three, the method further comprises:

[0045] In explicit time integration, using central difference algorithm to discretize in time, equation (14) and equation (15) are discretized as:

[0046] (17)

[0047] (18)

[0048] where, is the displacement of the center of mass of the rigid body, is the rotational displacement of the rigid body, is the velocity of the rigid body, is the angular velocity of the rigid body, represents the sum of external forces other than the hinge internal force, is the rigid body torque caused by the hinge internal force, is the rigid body torque caused by other forces other than the hinge internal force, is the hinge internal force, is the time increment between the nth step and the n+1 step, denotes the time increment between the nth step and the n+1 step; and are the translational displacements of the rigid body at time steps n and n+1, respectively, and are the rotational displacements of the rigid body at time steps n and n+1, respectively, is the translational velocity of the rigid body at time step n+1 / 2, is the rotational velocity of the rigid body at time step n+1 / 2;

[0049] For any point on the rigid body, the displacement of the point and the center of mass displacement have the following relationship:

[0050] (19)

[0051] where, is the position vector of the constraint node relative to the center of mass of the rigid body, , is the rotational displacement of the rigid body.

[0052] Further, in step three, the method further comprises: according to the rigid body kinematics principle, the displacement increment between time step n and time step n+1 is derived from the center of mass translation displacement and the rotational displacement:

[0053] (20)

[0054] where, is the position vector of the constraint node relative to the center of mass of the rigid body, the formula (20) is approximately expressed as:

[0055] (21)

[0056] Hinge internal force refers to the constraint force of rigid body A and rigid body B at the hinge constraint points 1, 2, is the constraint force of rigid body B on rigid body A at node 1 on rigid body A, is the constraint force of rigid body A on rigid body B at node 2 on rigid body B, and is the interaction force, which satisfies:

[0057] (22)

[0058] Total force acting on the rigid body and the constraint force from other hinge units are coupled with each other.

[0059] Further, in step three, the method further comprises:

[0060] is approximately expressed as:

[0061] (23)

[0062] where, is the total force acting on the rigid body at the current time, i=1 represents the total force acting on rigid body A, i=2 represents the total force acting on rigid body B, is the current study of the spherical hinge constraint internal force at time;

[0063] Substituting formula (17), formula (18), formula (21)-(23) into formula (12), we get:

[0064] (24)

[0065] wherein, is the coefficient matrix, denotes the increment of the hinged constraint force between time steps and ; is a three-dimensional vector calculated by acting on the rigid body by the resultant force and the moment as well as the velocity and displacement at the constraint node

[0066]

[0067] is the inverse matrix of the inertia tensor of rigid body A, is the inverse matrix of the inertia tensor of rigid body B.

[0068] The present application has the following beneficial effects: the hinged constraint internal force solving method based on Newton-Euler equation of the present application directly solves the internal force between multiple rigid body constraint nodes based on a simplified hinged structure, combines an explicit time integration algorithm and a rigid body Newton-Euler equation, and realizes the transient dynamic behavior analysis of the hinged structure in the explicit format without introducing artificial parameters. Compared with the traditional penalty function method and the Lagrange multiplier method, the present application fundamentally avoids the uncertainty and calculation error caused by the selection of the penalty factor in the traditional penalty function method, and ensures that the constraint equation is accurately satisfied. Compared with the Lagrange multiplier method, the present application directly constructs the constraint equation, bypasses the solving process of the augmented system equation, and significantly improves the calculation efficiency through the explicit solving scheme, and is naturally compatible with the parallel computing architecture. Since no artificial parameters are introduced, the calculation stability is no longer affected by the parameter value, and the present application is more stable than the traditional penalty function method. BRIEF DESCRIPTION OF DRAWINGS

[0069] In order to more clearly illustrate the technical solutions of the present application, the drawings needed in the embodiments will be briefly introduced below, and obviously, other drawings can also be obtained by those skilled in the art without creative labor on the basis of these drawings.

[0070] Figure 1 is a schematic diagram of a spherical hinge structure, wherein (a) is a schematic diagram of a spherical hinge structure, (b) is a schematic diagram of a cross-sectional view of the spherical hinge structure;

[0071] Figure 2 is a schematic diagram of a revolute pair structure, wherein (a) is a schematic diagram of a revolute pair structure, (b) is a schematic diagram of the positions of nodes N1, N2, N3 and N4. DETAILED DESCRIPTION

[0072] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions will be described below in connection with the specific embodiments of the present application and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without any creative work fall within the protection scope of the present application. The technical solutions provided by the embodiments of the present application will be described in detail below in connection with the drawings.

[0073] Please refer to Figure 1 and Figure 2 , a hinge constraint internal force solving method based on Newton-Euler equation in an embodiment of the present application, comprising the following steps:

[0074] Step one: simplifying the hinge structure and constructing the equivalent constraint equation.

[0075] In mechanical structures, spherical hinge and rotary pair are the most common two kinds of hinge structures, both of which can be simplified and the equivalent constraint equation can be constructed. The basic structure of spherical hinge is shown in Figure 1 . The constraint characteristics of spherical hinge are that there is a fixed spatial coincidence relationship between two mechanical components, that is, a specific node N1 on component A must always coincide with the corresponding node N2 on component B in spatial position. Under this constraint condition, the mechanical structure can rotate around the constraint point, but the translational freedom is completely constrained. The equivalent constraint equation can be expressed by the coordinates of two nodes:

[0076] (1)

[0077] Wherein, x1 is the coordinate of node N1, and x2 is the coordinate of node N2.

[0078] As shown in Figure 2 , the rotary pair constrains three translational freedoms and two rotational freedoms between components, and only one rotational freedom around the central axis is reserved. In essence, it can be equivalent to the combination of two spherical hinges. Therefore, the equivalent constraint equation of the rotary pair can be expressed according to the position vectors of two groups of four constraint nodes. Compared with the spherical hinge, it increases a group of coincident constraint nodes N3 and N4. The first constraint equation is the same as formula (1), and the second constraint equation is

[0079] (2)

[0080] Wherein, x3 is the coordinate of node N3, and x4 is the coordinate of node N4.

[0081] Step two: solving the rigid body dynamics equation.

[0082] The dynamic response of a single rigid body under external forces follows Newton's equation and Euler's equation, corresponding to the translational and rotational degrees of freedom of the rigid body, respectively, and the mathematical expressions are:

[0083] (3)

[0084] (4)

[0085] where, represents the mass of the rigid body, is the position of the center of mass of the rigid body, is the acceleration of the center of mass, is the inertia tensor of the rigid body, is the angular velocity of the rigid body, is the angular acceleration of the rigid body, is the total external force, is the total external moment. The mass of the rigid body is the sum of the masses of the nodes in the rigid body's dependent node set :

[0086] (5)

[0087] The inertia tensor of the rigid body is given by the following formula:

[0088] (6)

[0089] where is the position of the node .

[0090] After calculating the translational and rotational accelerations of the rigid body using equations (3) and (4), the central difference algorithm is used to update the displacement increments of the rigid body's translation and rotation in the explicit time integration framework. The update formula for the translational coordinates of the center of mass of the rigid body is as follows:

[0091] (7)

[0092] where, is the position of the center of mass of the rigid body at time step n, is its position at time step n+1, is the translational displacement increment.

[0093] Based on the Hughes-Winget algorithm, the rotation displacement increment of the rigid body can be used to calculate the transformation matrix :

[0094] (8)

[0095] (9)

[0096] in, It is the identity matrix. It is the intermediate matrix. It is any vector.

[0097] Based on the transformation matrix and time step Inertia matrix at time step Update time step based on equation (10) Inertial tensor of a rigid body at time t And update the coordinates of the subordinate nodes on the rigid body using equation (11):

[0098] (10)

[0099] (11)

[0100] in, and Let represent the coordinates of node i at time step n and time step n+1, respectively. This represents the increment of the translational displacement of a rigid body. This represents the position vector of node i relative to the center of mass of the rigid body.

[0101] Step 3: Solve the constraint equations and calculate the constraint internal forces.

[0102] The constraint equations for the ball joint can be found here (1), and can be equivalently expressed by the following equations:

[0103] (12)

[0104] in, This represents the node displacement at time step n+1, while This represents the nodal displacement at step n. express relatively The position vector, that is:

[0105] (13)

[0106] The dynamic equations of the rigid body have been given in equations (3) and (4). For ease of calculation, they can be expressed in the following form.

[0107] (14)

[0108] (15)

[0109] In this expression, a is the translational acceleration of the rigid body, is the angular acceleration of the rigid body, denotes the inertia tensor matrix of the rigid body in the global coordinate system, which is updated at each time step according to equation (10). denotes the equivalent external moment, while denotes the resultant moment acting on the rigid body, and their relationship is:

[0110] (16)

[0111] In explicit time integration, using the central difference algorithm to discretize in time, equations (14) and (15) can be discretized as

[0112] (17)

[0113] (18)

[0114] where, is the displacement of the center of mass of the rigid body, is the rotational displacement of the rigid body, is the velocity of the rigid body, is the angular velocity of the rigid body, represents the sum of external forces other than the hinge internal force, is the rigid body torque caused by the hinge internal force, is the rigid body torque caused by other forces other than the hinge internal force, is the hinge internal force, is the time increment between the nth step and the n+1 step, while denotes the time increment between the n step and the n+1 step. and are the translational displacements of the rigid body at time steps n and n+1, respectively, and are the rotational displacements of the rigid body at time steps n and n+1, respectively, is the translational velocity of the rigid body at time step n+1 / 2, is the rotational velocity of the rigid body at time step n+1 / 2.

[0115] For any point on the rigid body, the displacement of the point and the center of mass displacement have the relationship

[0116] (19)

[0117] where, ​​represents the position vector of the point relative to the center of mass of the rigid body, i.e. . represents the rotational displacement of the rigid body.

[0118] According to the kinematics of rigid body, the displacement increment between time step n and time step n+1 can be derived from the translational displacement and rotational displacement of the center of mass:

[0119] (20)

[0120] where, represents the position vector of the constraint node relative to the center of mass of the rigid body. Since the time step is very small in the explicit time integration format, the change of the position vector before and after a single time step can be ignored. Therefore, equation (20) can be approximated as:

[0121] (21)

[0122] The hinge internal force refers to the constraint force of rigid body A and rigid body B at the hinge constraint points 1,2, represents the constraint force of rigid body B on rigid body A at node 1 on rigid body A, represents the constraint force of rigid body A on rigid body B at node 2 on rigid body B, and are interaction forces, satisfying:

[0123] (22)

[0124] The total force acting on the rigid body (excluding the hinge internal force currently being studied) mainly comes from the constraint forces of other hinge elements and external forces such as gravity, contact forces, etc. However, since multiple hinge structures may exert constraints on the same rigid body, the total force is essentially coupled with the constraint forces from other hinge elements.

[0125] Considering that the change of the constraint force of each hinge element within a time increment step can be ignored, it can be approximated as:

[0126] (23)

[0127] where, is the total force acting on the rigid body at the current time, i=1 represents the total force acting on rigid body A, and i=2 represents the total force acting on rigid body B.​ denotes The momentary research spherical hinge constraint internal force.

[0128] Substitute formula (17), formula (18), formula (21)-(23) into formula (12), obtain:

[0129] (24)

[0130] Wherein, is a coefficient matrix, Denote the increment of the hinge constraint force between time step And . It is a three-dimensional vector, which can be calculated by acting on the rigid body force And the moment And the velocity at the constraint node And displacement :

[0131]

[0132] Is the inverse matrix of the inertia tensor of rigid body A, Is the inverse matrix of the inertia tensor of rigid body B.

[0133] In summary, the present application is aimed at the problem of insufficient calculation accuracy and poor calculation stability in the transient dynamics analysis of the articulated structure of the automobile, and a kind of articulated structure constraint force calculation method suitable for explicit time integration format is proposed.The method effectively improves the calculation accuracy while significantly enhancing the calculation stability.The articulated force calculation method proposed in the present application theoretically avoids the introduction of artificial parameters, so that the calculation result is not affected by the parameter selection.Therefore, the method has stronger adaptability under different working conditions, and is more suitable for engineering application.

[0134] The above-mentioned embodiments of the present application do not constitute a limitation on the scope of protection of the present application.

Claims

1. A method for solving the internal force of a hinged constraint based on Newton-Euler equations, characterized in that, Comprise: Step one, simplify the hinged structure, build equivalent constraint equation; wherein, the hinged structure includes spherical hinge and rotating pair, the constraint characteristics of spherical hinge is that there is fixed space coincidence relation between two mechanical components, the specific node N1 on component A and the corresponding node N2 on component B always coincide in space position, the equivalent constraint equation of spherical hinge is expressed by the coordinates of two nodes: (1)； Wherein, x1 is the coordinate of node N1, x2 is the coordinate of node N2; The equivalent constraint equation of rotating pair is expressed according to the position vectors of two groups of four constraint nodes, compared with spherical hinge, it increases a group of coincident constraint nodes N3 and N4, the first group of constraint equations is the same as formula (1), the second group of constraint equations is: (2)； Wherein, x3 is the coordinate of node N3, x4 is the coordinate of node N4; Step two: Solve the rigid body dynamics equations. Under external forces, the dynamic response of a single rigid body follows Newton's and Euler's equations, corresponding to the translational and rotational degrees of freedom, respectively. After calculating the translational and rotational accelerations based on their mathematical expressions, the translational and rotational displacement increments are updated using a central difference algorithm within an explicit time integration framework. The rotational displacement increments of the rigid body are then utilized. To calculate the transformation matrix Based on the transformation matrix and time step Inertia matrix at time step Update time step Inertial tensor of a rigid body at time t And update the coordinates of the subordinate nodes on the rigid body; Step three, solve the constraint equation, calculate the constraint internal force.

2. The method of claim 1, wherein the Newton-Euler equation based method of solving for the internal forces of a joint constrained system is characterized by, In step two, the translational and rotational degrees of freedom of rigid body are mathematically expressed as: (3) (4) wherein, denotes the mass of the rigid body, is the acceleration of the center of mass of the rigid body, is the resultant external force; is the inertia tensor of the rigid body, is the angular velocity of the rigid body, is the angular acceleration of the rigid body, is the resultant external moment.

3. The method of claim 2, wherein the Newton-Euler equation-based method of solving for the internal forces of a joint-constrained body is characterized by, mass of the rigid body set of slave nodes for the rigid body intermediate node mass sum of the masses: (5) Inertia tensor of a rigid body is given by the formula: (6) wherein is the position of the node .

4. The method of claim 3, wherein the Newton-Euler equation-based method of solving for the internal forces of a joint-constrained body is characterized by, In step two, under the explicit time integration framework, the central difference algorithm is used to update the displacement increment of the translation and rotation of the rigid body, and the update formula of the translational coordinate of the center of mass of the rigid body is as follows: (7) wherein, is the position of the center of mass of the rigid body at time step n, is the position of the center of mass of the rigid body at time step n+1, is the translational displacement increment.

5. The method of claim 4, wherein the Newton-Euler equation-based method of solving for the internal forces of a joint-constrained body is characterized by, In step two, the transformation matrix is calculated based on the Hughes-Winget algorithm using the incremental rotation displacement of the rigid body : T = R(θ) = I + (sin θ)S + (1 - cos θ)S2 : where I is the identity matrix, S is the skew symmetric matrix of the vector v, and R(θ) is the rotation matrix (8) (9) wherein is the identity matrix, is an intermediate matrix, is an arbitrary vector.

6. The method of claim 5, wherein the Newton-Euler equation-based method of solving for the internal forces of a joint-constrained body is characterized by, In step two, based on the transformation matrix and time step Inertia matrix at time step Update time step based on equation (10) Inertial tensor of a rigid body at time t And update the coordinates of the subordinate nodes on the rigid body using equation (11): (10) (11) wherein, and Xi(n) and Xi(n+1) denote the coordinates of node i at time step n and time step n+1, respectively, denotes the incremental translational displacement of the rigid body, denotes the position vector of node i relative to the center of mass of the rigid body.

7. The method of claim 6, wherein the Newton-Euler equation based method of solving for the internal forces of a joint constrained system is characterized by, In step three, the constraint equation of spherical hinge is equivalently expressed by the following equation: (12) wherein, denotes the node displacement at time step n+1, denotes the node displacement at step n, denotes the position vector of the relative the position vector of the relative (13) The dynamic equation of rigid body has been given in formula (3) and formula (4), and is expressed in the following form: (14) (15) wherein, is the resultant external force, denotes the mass of the rigid body, a is the translational acceleration of the rigid body, is the angular acceleration of the rigid body, denotes the inertia tensor matrix of the rigid body in the global coordinate system, which is updated at each time step according to equation (10), denotes the equivalent external moment, while denotes the resultant moment acting on the rigid body, the relationship between them being: (16)。 8. The method of claim 7, wherein the Newton-Euler equation based method of solving for the internal forces of a joint constrained system is characterized by, In step three, the method further comprises: In explicit time integration, the central difference algorithm is used to discretize in time, and formula (14) and formula (15) are discretized as: (17) (18) where, is the displacement of the center of mass of the rigid body, is the rotational displacement of the rigid body, is the velocity of the rigid body, is the angular velocity of the rigid body, represents the sum of external forces other than the hinge internal force, is the rigid body torque due to the hinge internal force, is the rigid body torque due to other forces other than the hinge internal force, is the hinge internal force, is the time increment between the nth step and the n+1th step, denotes the time increment between the n th step and the n+1 th step; and are the translational displacements of the rigid body at time steps n and n+1, respectively, and are the rotational displacements of the rigid body at time steps n and n+1, respectively, is the translational velocity of the rigid body at time step n+1 / 2, is the rotational velocity of the rigid body at time step n+1 / 2. For any point on the rigid body, its displacement is related to the displacement of the center of mass by (19) wherein, represents the position vector of the constraint node relative to the center of mass of the rigid body, , represents the rotational displacement of the rigid body.

9. The method of claim 8, wherein the Newton-Euler equation based method of solving for internal forces of a joint constrained system is characterized by, In step three, the method further comprises: (20) wherein, denotes the position vector of the constraint node relative to the center of mass of the rigid body, which is approximated by equation (20) as (21) hinged internal forces denote the constraint forces of rigid body A and rigid body B at the hinged constraint points 1,2, denote the constraint forces of rigid body B on rigid body A at node 1 on rigid body A, denote the constraint forces of rigid body A on rigid body B at node 2 on rigid body B, and are interaction forces, satisfying: (22) resulting force acting on the rigid body with the constraint forces from other articulated units are mutually coupled.

10. The method of claim 9, wherein the Newton-Euler equation based method of solving for internal forces of a joint constrained system is characterized by, According to the kinematics principle of rigid body, the displacement increment from time step n to time step n+1 is derived from the translational displacement and rotational displacement of the center of mass: Approximately expressed as: (23) wherein, is the current moment acting on the rigid body, i = 1 indicates the moment acting on the rigid body A, i = 2 indicates the moment acting on the rigid body B, indicates the current study of the spherical hinge constraint internal force; In step three, the method further comprises: Substitute formula (17), formula (18), formula (21)-(23) into formula (12), and the following formula is obtained: (24) wherein, is a coefficient matrix, denotes an increment of the constraint force between time steps and ; is a three-dimensional vector calculated by acting on the rigid body with the resultant force and the moment and the velocity and the displacement at the constraint node: is the inverse matrix of the inertia tensor of the rigid body A, is the inverse matrix of the inertia tensor of the rigid body B.

Images