Dynamics Simulation Method, System and Storage Medium for Multi-Body System Based on Lie Group to Describe Rotation

Abstract

The present invention provides a multi-body system dynamics simulation method, system and storage medium based on Lie group description of rotation, including Lie group generalized-α integration method and Lie group backward difference method. The beneficial effects of the present invention are as follows: The rotation described by Lie group avoids the singularity of the three rotation parameters, making the integration process continuous and without frequent restart; it avoids the additional constraint equations of Euler quaternions and overcomes the rotation speed saturation problem of Euler quaternions at the same time; it is more suitable for large-scale calculations compared with the three rotation parameters and Euler quaternions.

Description

Technical Field

[0001] The present invention relates to the technical field of multi-body system dynamics simulation, and particularly to a multi-body system dynamics simulation method, system and storage medium based on Lie group to describe rotation. Background Art

[0002] In order to solve differential-algebraic equations, a suitable numerical integration algorithm needs to be adopted. General numerical integrators need to meet basic conditions such as accuracy, efficiency and robustness. According to whether iteration is required in the calculation process, the integrators can be divided into two types: implicit and explicit. In multi-body system dynamics simulation, the mainstream general integrators are generally in implicit format, and the most widely used in practical applications are the integrator families of the generalized α method and the integrator families of backward difference formulas.

[0003] The above integrators can solve the dynamic equations described by quaternions. The quaternion method can solve most multi-body dynamics simulation problems. The advantage is that it avoids rotation singularities, and the disadvantage is that additional constraint equations are introduced in the calculation, increasing the number of equations to be solved. Research shows that global parameters with three components, such as Euler angles and rotation vectors, will lead to rotation singularities. To solve such singularity problems, special strategies need to be adopted. For example, for Euler angles, re-parameterization can be carried out by switching variables to avoid singular points; this is also the strategy adopted by the commercial software ADAMS. In addition, the Hilber-Hughes-Taylor method and the Newmark method based on local rotation vectors can also be adopted, where the local rotation vector represents the increment of each step, that is, the rotation parameters of each time step are used to describe the rotation of the current rotation frame relative to the rotation frame of the previous time step.

[0004] All methods of describing rotation with three parameters have singular configurations and are very likely to encounter problems in large-scale multi-body system dynamics simulation, especially for finite element flexible bodies. Although the singularity can be avoided by replacing two sets of parameters with each other, this processing method is not perfect and is not suitable for high-performance computing.

[0005] On the other hand, quaternions themselves will add a set of constraint equations to the dynamic equations; on the other hand, if the description method of quaternions is adopted, saturation problems will occur in high-speed rotation problems. For example, as Figure 1 shown, a constant spatial torque acts on a rigid body sphere. Theoretically, the angular velocity of the rigid body should increase linearly. However, a strange angular velocity saturation phenomenon will occur in the actual calculation results, as Figure 2As shown. The simulation found that when calculating using the Lagrange equation derived from quaternions, no matter how the integrator is selected and how the error control is selected, the angular velocity will reach saturation at a certain moment. This saturation phenomenon can be solved by using the Euler equation, but due to limitations, its integration time step is generally very small, and the simulation calculation efficiency is not high. Summary of the Invention

[0006] The present invention provides a multi-body system dynamics simulation method based on Lie group to describe rotation, including the following steps:

[0007] Step 1: Input parameters into the integrator. The input parameters include the number of output times n, the initial time step h, the initial time t(1), and the initial values. The initial values include the initial value q of displacement n , the initial value of velocity , the initial value of acceleration , the initial value λ of the Lagrange multiplier n , and the initial value of the acceleration auxiliary variable a n ;

[0008] Step 2: Estimate the predicted value of displacement at the next time step , the predicted value of velocity , the predicted value of acceleration , and the predicted value of the Lagrange multiplier , and the predicted value of the acceleration auxiliary variable

[0009] Step 3: Solve the Jacobian matrix jacobi(n + 1) and solve the residual residual(n + 1);

[0010] Step 4: Newton iteration, correct the estimated value, and solve the tolerance value to1;

[0011] Step 5: Determine whether the result converges. If so, then execute Step 6. Otherwise, return to execute Step 3;

[0012] Step 6: Solve the error value err;

[0013] Step 7: Determine whether the step size meets the acceptance criteria. If so, then execute Step 8. Otherwise, update the step size and then return to execute Step 2;

[0014] Step 8: Output the displacement q at the next time step n+1 , velocity , acceleration a n+1 , and the Lagrange multiplier λ n+1 ;

[0015] Step 9: Determine whether the time t is the last step. If so, then end. Otherwise, execute Step 10;

[0016] Step 10: n = n + 1, update the step size, and then return to execute Step 2.

[0017] As a further improvement of the present invention, in the said Step 3, the Jacobian matrix is where η, β, γ are the integral parameters of the generalized-α integral, M is the mass matrix of the kinetic equation, and F q and are the partial derivatives of the generalized force with respect to displacement and velocity, respectively.

[0018] The present invention also provides a multi-body system dynamics simulation system based on Lie group description of rotation, including: a memory, a processor, and a computer program stored on the memory, and the computer program is configured to implement the steps of the multi-body system dynamics simulation method of the present invention when called by the processor.

[0019] The present invention also provides a computer-readable storage medium, and the computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of the multi-body system dynamics simulation method of the present invention when called by a processor.

[0020] The beneficial effects of the present invention are: the rotation described by the Lie group avoids the singularity of the three rotation parameters, making the integration process continuous and without the need for frequent restart; it avoids the additional constraint equations of the Euler quaternion and overcomes the problem of rotational speed saturation of the Euler quaternion; it is more suitable for large-scale calculations compared with the three rotation parameters and the Euler quaternion. Description of the Drawings

[0021] Figure 1 is a schematic diagram of a rotating rigid body;

[0022] Figure 2 is a schematic diagram of the rotation saturation phenomenon;

[0023] Figure 3 is a flowchart of the Lie group generalized-α integration method;

[0024] Figure 4 is a flowchart of the Lie group backward difference algorithm. Detailed Embodiments

[0025] The present invention utilizes the Lie group characteristics of rotational motion to convert the relevant calculations of the rotation matrix into operations near the identity matrix, thereby avoiding the phenomenon of numerical singularities.

[0026] In the dynamics of multi-body systems, the choice of rotation parameters directly affects the design of integrator algorithms. In engineering practice, the most widely used rotation parameters are Euler angles, rotation vectors, and quaternions. To achieve the avoidance of redundant constraints and singularities in dynamic simulations, a feasible solution is to use Lie group methods to describe rotations. Different from integrators based on global rotation parameters such as quaternions or Euler angles, Lie group integrators use local rotation frames.

[0027] Now consider a matrix R that satisfies

[0028] R(t) = A T (t n )A(t)

[0029] It is easy to know that R ∈ SO(3), that is, R lies on the Lie group manifold SO(3) composed of the set of rotation matrices. The matrix A can be parameterized using the local rotation vector θ as

[0030]

[0031] where θ is the local rotation vector, and is the skew-symmetric matrix of the local rotation vector, such as

[0032]

[0033] The rotation matrix derived from the above formula satisfies the orthogonality requirement. This method can not only ensure the orthogonality of the rotation matrix at each step but also avoid the singularities of rotations without introducing redundant constraints.

[0034] Utilizing the Lie group properties of rotational motion, the related calculations of the rotation matrix are transformed into operations near the identity matrix, thus avoiding the phenomenon of numerical singularities. This idea can be directly applied to the implementation of the generalized α integrator and the BDF integrator.

[0035] Lie group generalized α integration method

[0036] For the Lie group generalized α integrator, assume that at time t n the rotation matrix is A n , the angular velocity is ω n , the angular acceleration is and the auxiliary angular acceleration is a n ; and assume that the next predicted time step h is known. Next, solve for the state variables at t n+1 = t n + h. Assume that the predicted value of a n+1 is and the correction is Then, from the Lie group generalized α scheme:

[0037]

[0038] where θ n+1 represents the rotation vector in the local coordinates at time n+1, is θ n+1 The estimated value of ω n+1 is the body angular velocity at time n+1, Yes n+1 estimated value of; is the body angular velocity at time n+1, yes A n+1 is the direction cosine matrix at time n+1, and β, γ, and η are the generalized α integrator parameters.

[0039] The above prediction values are:

[0040]

[0041] where α m and α f is the generalized α integrator parameter.

[0042] Substituting the formula for the next time step into the dynamic equations to solve them is the Lie group generalized alpha integral format. In Lie group generalized alpha integrals, there's no requirement that the norm of θ must be less than 2π, meaning that the rotation angle in a single time step can exceed one full rotation. This significantly improves simulation efficiency for high-speed systems.

[0043] Lie Group Integrator Backward Difference Method

[0044] In the Lie group integrator backward difference method, the time t n is the reference node, except for the rotation vector θ at time n n =0, other time nodes t need to be calculated n+1 (j=1,2,3,…,k+1) relative to t n The local rotation variable θ at the nj moment n-j . Given a rotation matrix A n-j The size of A n-j satisfy

[0045]

[0046] The definition is shown in formula (1.1), from which θ can be calculated n-j .

[0047] Record the time step vector from 1 to k+1 steps

[0048]

[0049] in, Assume that the rotation vector at time n is the rotation vector variable θ(t) on the uniform grid t n ,t n -h,t n -2h,…,t n -kh value is obtained by iterating the following formula

[0050]

[0051] Assume that x n+1 , x n is any variable at time n+1 and time n, the post-difference operator The definition is as follows

[0052]

[0053] Therefore, the difference matrix Can be expressed as

[0054]

[0055] The expression of the difference matrix can be derived

[0056]

[0057] in is the coefficient matrix of the difference matrix.

[0058] The local rotation vector and its first-order derivative with respect to time can be written as

[0059]

[0060] Among them, their corresponding estimated values are

[0061]

[0062] Calculate θ n+1 ,Can Calculate the rotation matrix.

[0063] The present invention discloses a multi-body system dynamics simulation method based on Lie group description of rotation, while the translational motion variables are still described by Cartesian coordinates, which is no different from a general integrator.

[0064] like Figure 3 As shown in Figure 2, the generalized-α integral method for multi-body dynamics based on Lie group description of rotation includes the following steps:

[0065] Step 1: Input parameters into the integrator, including the number of outputs n, the initial step length h (the initial step length h is the initial time step length h, which will change as the calculation progresses), and the initial time t n , and initial values, including the initial value q of the displacement n , the initial value of the speed Initial value of acceleration Initial value λ of the Lagrange multiplier n , and the acceleration auxiliary variable a n The initial value of .

[0066] Step 2: Estimate the displacement at the next time step Speed estimate Acceleration estimate and Lagrange multiplier estimates and acceleration auxiliary variable estimates

[0067] Step 3: Solve the Jacobian matrix jacobi(n+1), solve the residual(n+1), the Jacobian matrix is where η, β, and γ are the integral parameters of the generalized-α integral, M is the mass matrix of the kinetic equation, and F q and are the partial derivatives of the generalized force with respect to displacement and velocity, respectively.

[0068] Step 4: Newton iteration, correct the estimated value, solve the tolerance value to1, and use Newton iteration to obtain the variable x.

[0069] Step 5: Determine whether the result converges. If so, proceed to step 6; otherwise, return to step 3.

[0070] Step 6: Solve the error value err.

[0071] Step 7: Determine whether the step size meets the acceptance criteria. If so, proceed to step 8. Otherwise, update the step size and then return to step 2.

[0072] Step 8: Output the displacement q of the next time step n+1 ,speed acceleration a n+1 and the Lagrange multiplier λ n+1 and other parameters.

[0073] Step 9: Determine whether time t is the last step (i.e., determine whether the integration is completed). If so, then end; otherwise, execute step 10.

[0074] Step 10: n=n+1, update the step size, and then return to step 2.

[0075] In step 2, at each time step t n introduce α m and α f as well as variable a n and a n+1 , such that they are related to and satisfy the following relationship

[0076]

[0077] Similarly, when calculating the generalized coordinates q and velocity , the acceleration terms used are interpolated by introducing γ and β

[0078]

[0079] Introduce a correction quantity x such that and let the predicted value be the predicted value of a n+1 . Substitute the expression of a n+1 into the above formula and simplify to obtain the expressions of acceleration, velocity and position vector

[0080]

[0081] In the Lie group integrator, the physical quantities describing translation and deformation are described by the above formulas; but the description of rotation is by the Lie group method, which means that the rotation variables used in each step of calculation are local rotation variables relative to the reference time node. Assume that at time t n the rotation matrix is A n , the angular velocity is ω n , and the angular acceleration is a n+1 the predicted value is the correction quantity is Taking time t n as the reference time node, for the rotation variables, q n = θ n = 0, then at t n+1 = t n + h, the state variables are

[0082]

[0083] where the angular velocity ω n+1 and its predicted value the angular acceleration term and its predicted value the rotation vector θ n+1 and its predicted value

[0084] The rotation matrix for each time step is calculated by the following formula

[0085]

[0086] Substituting the above into the system of equations to be solved, a non-linear system of equations about the unknown small quantity x is obtained. At each step of the calculation, θ n+1 is obtained, and the rotation matrix is given by Equation (1.5). In the program implementation, for each time step, the part of the rotation variables corresponding to Equation (1.4) needs to be used, and the rotation matrix is updated by Equation (1.5).

[0087] Backward difference method for describing rotation based on Lie group method:

[0088] In the Lie group backward difference method, except for the parameters used to describe rotation, the discretization of all other generalized coordinates is obtained using the conventional backward difference format; for the rotation parameters, at each time step, differences need to be taken using the rotation parameters at multiple known time steps, and they are all local quantities relative to the reference configuration at the current moment. As Figure 4 shown, the following introduces a backward difference integrator for describing rotation based on Lie group method:

[0089] Near the time t n , using local coordinates as the generalized coordinates for describing rotation, satisfying

[0090]

[0091] where A n is the direction cosine matrix at time n. A(t) and θ(t) are the direction cosine matrix and the local rotation vector at time t, respectively.

[0092] Define again

[0093]

[0094] In the above formula, A represents the rotation matrix, θ represents the local rotation vector, the subscript represents the time step, n is the time step, and n - j is the time step.

[0095] The next time step can be calculated using θ n-1 , θ n-2 , θ n-3 , …, θ n-k-1 as follows

[0096]

[0097] For the formula description, see Equations (1.2, 1.3). Among them can be obtained from the time data and calculated using two algorithms:

[0098] The first is the algorithm of the quasi-fixed step size formula:

[0099] Step (1): Normalize and take where h = t n -t n-1 , is the current time step, is historical time data;

[0100] Step (2): Calculate the Vandermonde matrix P k =V(x k );

[0101] Step (3): Generate matrix O k , whose (i,j) component is

[0102] Step (4): Calculation Replace k with k+1 to get the desired value.

[0103] The second is the algorithm of variable step size formula:

[0104] Step a: Calculation

[0105] Step b: Take and

[0106] Step c: Calculate Q = (1 k o T ) / O

[0107] Step d: Calculation Replacing k with k+1 gives the desired value, which is the upper triangular matrix of the matrix calculated in step c.

[0108] The above steps only change the parameters such as displacement, velocity, acceleration and Lagrange multiplier of the next time step estimated in step 2 of the generalized α-Lie group integrator. The other steps are the same as above.

[0109] Figure 4 The Lie group backward difference algorithm is a Figure 3 Improvements to the Lie group generalized-α integral method. The specific improvements include:

[0110] In step 1, the parameters also include the initial order K, and the initial value of μ(1) is used to replace a n Initial value of; after step 1, perform the judgment step;

[0111] Determination step: Determine whether n is greater than or equal to the specified number (the specified number is 3). If so, solve the difference matrix Dnk+1 of the k+1 step and then execute step 2. Otherwise, solve the difference matrix Dnk of the k step and then execute step 2.

[0112] In the step 2, calculate t(n + 1) = t(n) + h to estimate the displacement prediction value at the next time step. Velocity prediction value Acceleration prediction value And Lagrange multiplier prediction value And

[0113] In the step 4, calculate the Newton iteration residual newnrm without solving for to1.

[0114] In the step 6, solve for two types of error parameters err and err1.

[0115] In the step 7, determine whether the step size meets the acceptance criteria. If so, then execute step 8; otherwise, update the step size, order K, update Dnk+1, and then return to execute step 2.

[0116] In the step 8, first calculate the order-changing elements e1, e2, e3, and then output the displacement q at the next time step n+1 , velocity acceleration Lagrange multiplier λ n+1 , μ(n + 1);

[0117] In the step 10, n = n + 1, and solve for three types of error indicators error1, error2, and error3.

[0118] The multi-body system dynamics simulation method further includes:

[0119] Step 11: Determine whether n is greater than the specified value (the specified value is 3). If so, then execute step 12; otherwise, Dnk = e3, and then return to execute the determination step.

[0120] Step 12: Solve for the relevant error indicators, step size, and order: err1, err2, err3, t1, h1, k1.

[0121] Step 13: Update the order k and the step size h, and then return to execute the determination step.

[0122] The present invention has the following beneficial effects:

[0123] 1. The rotation described by the Lie group avoids the singularity of the three rotation parameters, making the integration process continuous and eliminating the need for frequent restarts.

[0124] 2. It avoids the additional constraint equations of the Euler quaternion and overcomes the rotation speed saturation problem of the Euler quaternion.

[0125] 3. Compared with the three-parameter rotation and Euler quaternion, it is more suitable for large-scale calculations.

[0126] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For ordinary technicians in the technical field to which the present invention belongs, several simple deductions or substitutions can be made without departing from the concept of the present invention, which should be regarded as falling within the protection scope of the present invention.

Claims

1. A dynamic simulation method for multi-body systems based on Lie group to describe rotation, characterized in that, It includes the following steps: Step 1: Input parameters into the integrator. The input parameters include the number of output times n, the initial time step h, the initial time t(1), and the initial values, where the initial values include the initial value of displacement q n , the initial value of velocity , the initial value of acceleration , the initial value of the Lagrange multiplier λ n , and the initial value of the acceleration auxiliary variable a n ; Step 2: Estimate the displacement prediction value for the next time step Velocity prediction value Acceleration prediction value and Lagrange multiplier prediction value and acceleration auxiliary variable prediction value Step 3: Solve the Jacobian matrix jacobi(n+1) and solve the residual residual(n+1); Step 4: Perform Newton iteration, correct the estimated value, and solve the tolerance value to1; Step 5: Determine whether the result converges. If so, then execute Step 6; otherwise, return to execute Step 3; Step 6: Solve the error value err; Step 7: Determine whether the step size meets the acceptance criterion. If so, then execute Step 8; otherwise, update the step size, and then return to execute Step 2; Step 8: Output the displacement q at the next time step n+1 , velocity , acceleration a n+1 and the Lagrange multiplier λ n+1 ; Step 9: Determine whether the time t is the last step. If so, then end; otherwise, execute Step 10; Step 10: n=n+1, update the step size, and then return to execute Step 2; In step 2, the rotation is described using the Lie group method, and θ is calculated at each step. n+1 After that, the rotation matrix is updated by formula (1.5); in the program implementation, at each time step, the corresponding rotation variable part of formula (1.4) is required, and the rotation matrix is updated by formula (1.5). Among them, A n is a rotation matrix, β and γ are generalized-α integral parameters, is the skew-symmetric matrix of the local rotation vector, and θ n+1 represents the rotation vector in the local coordinate system at the (n + 1)-th moment; in the said step 1, the parameters further include the initial order K, and the initial value of μ(1) is used to replace the initial value of a n ; after step 1, a judgment step is executed. Judgment step: Determine whether n is greater than or equal to the specified number. If so, then solve the difference matrix Dnk+1 for K+1 steps, and then execute Step 2; otherwise, solve the difference matrix Dnk for K steps, and then execute Step 2; In step 2, calculate t(n + 1) = t(n) + h to estimate the displacement prediction value at the next time step Velocity prediction value Acceleration prediction value And Lagrange multiplier prediction value And In step 4, calculate the Newton iteration residual newnrm without solving for to1; In the said Step 6, solve two types of error parameters err and err1; In the said Step 7, determine whether the step size meets the acceptance criterion. If so, then execute Step 8; otherwise, update the step size, the order K, update Dnk+1, and then return to execute Step 2; In step 8, first calculate the order-changing elements e1, e2, and e3, and then output the displacement q, n+1 velocity, acceleration, Lagrange multiplier λ, n+1 and μ(n + 1); In the said Step 10, n=n+1, solve three types of error indicators error1, error2, and error3; This multi-body system dynamics simulation method further includes: Step 11: Determine whether n is greater than the specified value. If so, then execute Step 12; otherwise, Dnk=e3, and then return to execute the judgment step; Step 12: Solve the relevant error indicators, step size, and order: err1, err2, err3, t1, h1, k1; Step 13: Update the order k and the step size h, and then return to execute the judgment step.

2. The multi-body system dynamics simulation method according to claim 1, wherein: In the said step 3, the Jacobian matrix is where η, β, γ are the integral parameters of the generalized-α integration, M is the mass matrix of the dynamic equation, F q and are the partial derivatives of the generalized force with respect to displacement and velocity respectively.

3. The multi-body system dynamics simulation method according to claim 1, characterized in that: In the said step 2, the physical quantities describing the translation and deformation are expressed by the following formula: where x is the correction amount, and η, β, and γ are the generalized-α integration parameters.

4. The multi-body system dynamics simulation method according to claim 1, characterized in that: In the said judgment step, the specified number is 3; in the said Step 11, the specified value is 3.

5. The multi-body system dynamics simulation method according to claim 1, wherein: In the said step 2, the next time step is calculated using θ n-1 , θ n-2 , θ n-3 , …, θ n-k-1 and so on Among them, is calculated through the following steps: Step (1): Normalize and take where h = t n -t n-1 , which is the current time step, and is the historical time data; Step (2): Calculate the Vandermonde matrix P k = V(x k ); Step (3): Generate matrix O k , whose (i, j) component is Step (4): Calculate Replacing k with k + 1 gives the desired value.

6. The multi-body system dynamics simulation method according to claim 1, characterized in that: In the step 2, the next time step is calculated using θ n-1 , θ n-2 , θ n-3 , …, θ n-k-1 and so on Among them, is calculated through the following steps: Step a: Calculate Step b: Take and Step c: Calculate Q = (1 k o T ) / O Step d: Calculate Replacing k with k + 1 gives the desired value.

7. A multi-body system dynamics simulation system based on Lie group description of rotation, characterized in that, It includes: A memory, a processor, and a computer program stored on the memory, the computer program being configured to implement the steps of the multi-body system dynamics simulation method according to any one of claims 1-6 when called by the processor.

8. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, the computer program being configured to implement the steps of the multi-body system dynamics simulation method according to any one of claims 1-6 when called by a processor.

Images