Affine multi-body system real-time simulation method based on projection near field dynamics
By adopting projected near-field dynamics framework and semi-implicit continuous iteration method in multi-body dynamics simulation, the problem of difficulty in paralleling traditional methods is solved, and efficient parallel multi-body dynamics simulation is achieved, which improves the running speed and supports real-time simulation.
Patent Information
- Application Number
- CN202411506813.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2024-10-28
- Publication Date
- 2025-03-07
AI Technical Summary
Traditional constraint multibody dynamics solution methods based on implicit Euler are difficult to parallelize, which makes it difficult to achieve efficient simulation on GPUs.
The projected near-field dynamics framework is adopted, and global solutions are performed through dual form constraints and semi-implicit continuous iteration method, and local solutions are performed in combination with Newton's method to achieve efficient parallel multi-body dynamics simulation.
It realizes efficient parallel multi-body dynamics simulation on the GPU, improves the running speed, realizes real-time simulation, and simplifies the calculation of potential energy such as collision and friction.
Smart Images
Figure CN119578027A_ABST
Abstract
Description
Technical Field
[0001] The invention belongs to the field of computer graphics, and in particular relates to a real-time simulation method of an affine multi-body system based on projected peridynamics. Background Art
[0002] In multi-body dynamics simulation, the implicit Euler-based method can provide relatively stable and accurate simulation results. In the solution process, the Hessian matrix must be established for various nonlinear constraints, the problem is converted into a nonlinear optimization problem, and the Newton method is often used to solve the matrix. However, since the contact often changes during the simulation, it is necessary to frequently and time-consumingly update the global Hessian matrix. Moreover, due to the impact of contact changes on the values of non-zero non-diagonal elements, the storage method of the Hessian matrix cannot be unified. This means that the implicit Euler solution method using Newton's method is difficult to implement efficiently on the GPU. Summary of the invention
[0003] The purpose of the present invention is to address the problem that the traditional implicit Euler-based constrained multibody dynamics solution method is difficult to parallelize, and to use the projected peridynamics framework to remodel the problem to achieve efficient parallel multibody dynamics simulation.
[0004] The technical solution adopted by the present invention to achieve its purpose is as follows:
[0005] A real-time simulation method for an affine multi-body system based on projected peridynamics comprises the following steps:
[0006] (1) Construct dual constraints for the colliding objects in the multi-body system, including the constraint relationships between the dual action points, such as the constraint potential energy, collision, friction, and joints;
[0007] (2) A global energy equation is modeled for all the constructed constraint relationships, and the nonlinear constraints in the equation are globally solved using a semi-implicit continuous iteration method. The position of each action point is calculated through iterative updates;
[0008] (3) Use Newton's method to perform local solution, calculate the object's affine matrix and position coordinates, and update them;
[0009] (4) Repeat the global solution and local solution process until the relative error converges to a preset value or reaches the maximum number of iterations;
[0010] (5) Update the object's position, rotation state, linear velocity, and angular velocity.
[0011] Furthermore, in step (1), the energy constraint relationship is constructed for the constraint potential between the dual action points:
[0012]
[0013] Where Ψ represents the constraint potential energy, y represents the vector position of the action point, l represents the distance between the dual action points, ∈ represents the constraint strength, and C represents a linear function about y and l. The implicit continuous iteration of the energy constraint relationship is:
[0014]
[0015] Where i, j represent two colliding objects, y represents the position of the object, ξ represents the constraint, and k represents the kth iteration. Represents the distance between the dual action points at the kth iteration.
[0016] Furthermore, in step (1), the log-barrier function is used to construct the constraint relationship between the collision and friction between the dual action points:
[0017]
[0018] Among them B N Refers to the N-order polynomial barrier function, d represents the distance between the dual action points, Represents the maximum distance between the dual application points of the energy constraint.
[0019] Furthermore, three joint constraints are constructed for the joints in step (1):
[0020] Construct spherical joint constraints: The rotation of a spherical joint in three-dimensional space is fixed by applying three distance constraints on the three axes of the joint coordinate system. The distance constraint of each axis calculates the coordinates of two action points on the object, and the rotation matrix of the joint anchor point is used to ensure that the rotation of the object complies with the constraints of the spherical joint.
[0021] Construct hinge joint constraints: The hinge joint uses the distance constraint between two points to fix the direction of the rotation axis. An additional distance constraint is applied to the position between the two objects to control the opening and closing angle of the hinge.
[0022] Constructing sliding joint constraints: The sliding joint uses two distance constraints to limit the relative translation of the objects along the sliding axis; for torsion, an additional distance constraint is applied to the set position between the two objects to fix the rotation axis.
[0023] Furthermore, the global energy equation modeled in step (2) is:
[0024]
[0025] Where y represents the coordinates of all action points, is the time derivative of y, and E represents the y represents the momentum potential energy, Ψ represents the constraint potential energy, A irepresents the affine matrix of object i, r iξ represents the local coordinates of the point of action of constraint ξ on object i, p i Represents the center coordinates of object i, y iξ Represents the global coordinates of the point of action of constraint ξ on object i.
[0026] Furthermore, when the semi-implicit continuous iteration method is used in step (2) to globally solve the nonlinear constraints in the equation, the iterative formula for the kth iteration is:
[0027]
[0028] Among them Ψ ξ represents the constraint potential energy of constraint ξ, represents the derivative of Ψ, express The sign of the positive term is The negative sign of d ξ Represents the distance between the dual action points.
[0029] Furthermore, in step (2), the iterative formula of the kth iteration is substituted into the global energy equation for global solution, and the result is:
[0030]
[0031] in represents the coordinates of object i in constraint ξ at the kth iteration, is a constant, and the expression is as follows:
[0032]
[0033] Where h is the simulation time step, m i is the mass of object i.
[0034] Furthermore, the specific steps of step (3) include:
[0035] Perform an inverse transformation on the coordinates of all the action points of the object according to the rotation matrix of the current object, and rotate the object back to its initial state;
[0036] The local rigid constraints of the object are optimized and solved using Newton's method. The optimization goal is to maximize the rigidity of the object while minimizing the distance error between all action points on the object, and solve the updated affine matrix.
[0037] A two-step strategy is adopted for updating, including: the first step is to update the affine matrix and adjust the rotation state of the object; the second step is to keep the affine matrix unchanged and use Newton's method to iteratively update the center position of the object according to the dual action point in the global solution.
[0038] Furthermore, the optimization objective function of the local solution is:
[0039]
[0040] in‖·‖ F represents the Frobenius norm, q i represents the position coordinates of object i before update, p i represents the center position of object i, ω iξ represents the weight of the potential energy term at the dual action point, r iξ represents the local coordinates of the point of action of constraint ξ on object i, δ represents the penalty coefficient, and A i represents an affine matrix and I is the third-order identity matrix.
[0041] Furthermore, the relative error in step (4) is:
[0042]
[0043] where εk 表 is the relative error of the kth iteration, Ψ(x 0 ) represents the initial position x 0 The potential energy, Ψ(x k ) represents the kth iteration position x k The potential energy, Ψ(x N ) represents the Nth iteration position x N The potential energy is the potential energy after final convergence.
[0044] Compared with the prior art, the present invention has the following advantages:
[0045] (1) The present invention uses a projected peridynamics framework and a semi-implicit continuous iteration method in the global solution. Each iteration step is highly parallel and there is no need to assemble a global Hessian matrix. Therefore, it is suitable for the existing GPU architecture and can use hardware acceleration to greatly improve the running speed of large-scale complex multi-body physics simulation situations, and can achieve real-time simulation.
[0046] (2) The present invention uses affine bodies to model rigid bodies and ensures the rigidity of objects in local solutions. Compared with the traditional method of using rigid body modeling, it is simpler and faster in the calculation of potential energy such as collision and friction.
[0047] (3) This method uses a two-step strategy in local solution, which greatly improves the convergence and accuracy of local object affine matrix and position calculation.
[0048] (4) The present invention redesigns the joint constraints in dual form, which is more intuitive and concise than the traditional constraint implementation form of limiting the degree of freedom, and can realize the combination of degrees of freedom of different joints. BRIEF DESCRIPTION OF THE DRAWINGS
[0049] Figure 1 It is a schematic diagram of constraint modeling using traditional methods.
[0050] Figure 2 It is a schematic diagram of the dual form constraint modeling of the present invention.
[0051] Figure 3A These are graphs of polynomial barrier functions and log-barrier functions of different highest degrees.
[0052] Figure 3B is different Polynomial barrier function curve graph when taking values.
[0053] Figure 4 It is a schematic diagram of the two-step strategy in the local solution process.
[0054] Figure 5 It is a collision effect diagram simulated by the Bullet engine and the method of the present invention.
[0055] Figure 6 This is a comparison chart of the time step cost under different numbers of objects. DETAILED DESCRIPTION
[0056] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is described in detail below through specific implementations and drawings, but it does not constitute a limitation of the present invention.
[0057] The present invention specifically discloses an affine multibody dynamics simulation method based on projective peridynamics, which adopts the projective peridynamics framework to globally solve the linear system and locally solve all nonlinear terms. Projective peridynamics is a method that can be easily implemented in parallel to solve nonlinear optimization in constrained multibody dynamics. This method can solve constrained multibody dynamics problems with nonlinear constraints based on projective peridynamics. This method models the constraints as peridynamic bonds, which helps to integrate semi-implicit continuous iteration methods to deal with nonlinear terms in the global step size. Traditional methods regard the object as an intangible rigid body, and solve each degree of freedom of the rigid body according to the constraints, such as Figure 1 Different from the traditional method, this method reformulates the original constrained multibody dynamics into a dual form, considering the system as a continuum subject to rigid constraints, such as Figure 2As shown. This method uses an affine matrix instead of a rotation matrix to represent the rotation state of an object, which helps to achieve coupling by imposing linear or nonlinear constraints between rigid and flexible bodies. In the global solution, this method uses a semi-implicit continuous iteration method (SISSM) to solve nonlinear constraints, which expands the scope of projected peridynamics in dealing with nonlinear problems and is very suitable for parallel implementation on GPUs; in the local solution, Newton's method is used to solve the rigid constraints of the object, and a two-step strategy is adopted to ensure the stable translation and rotation of the object.
[0058] The specific steps processed by this method within one time step include:
[0059] (1) Construct dual constraints for the colliding object pairs in the multi-body system (i.e., constraints between the dual action points in the coordinate systems of the two objects).
[0060] The soft kinematic constraints are constructed for the constraint potential between the dual action points, which is similar to the quadratic energy function of a spring oscillator. The energy constraint relationship is:
[0061]
[0062] Where Ψ represents the constraint potential, y represents the spatial coordinates of the action point, l represents the distance between the action points, ∈ represents the constraint strength, and C represents a linear function about y and l. The implicit continuous iteration of the constraint relationship is:
[0063]
[0064] Where i, j represent two colliding objects, y represents the coordinates of the objects, ξ represents the constraint, and k represents the kth iteration. Represents the distance between the dual action points at the kth iteration.
[0065] For collision and friction, hard kinematic constraints are constructed between the dual action points. The approximate polynomial used to construct the hard kinematic constraints using the log-barrier function is:
[0066]
[0067] Among them B N Refers to the N-order polynomial barrier function. The highest order of the polynomial barrier function is 2N. d represents the distance between the dual action points. The maximum distance between the dual action points of the energy constraint, the default value is 0.001. Unlike the original method, the approximate polynomial method accepts smaller penetration, Figure 3A The energy-distance curves of the polynomial barrier function and the log-barrier function are shown for different N values. Figure 3B Shows a different Force-distance curve of the polynomial barrier function.
[0068] For joints, the following three joint constraints are constructed:
[0069] i. Spherical joint
[0070] For an unconstrained spherical joint, a distance constraint is directly used to connect two points on two objects. For a spherical joint rotation angle constraint in three-dimensional space, three distance constraints are applied on the three axes of the joint frame to fix the spherical joint rotation angle. The expression of the distance constraint point on one axis in the rotation angle constraint of the spherical joint (i.e. the point where the constraint is established) is as follows:
[0071] y i1 =y i0 +R i r1
[0072] y j1 =y j0 +R(θ)R j r1
[0073] where y i0 represents the coordinates of the action point 0 on object i in the world coordinate system, r1 represents the position of the action point in the joint anchor point coordinate system, and R i represents the rotation matrix of object i, and R(θ) is the rotation matrix corresponding to the axial rotation angle.
[0074] ii. Hinge joint
[0075] For unconstrained hinge joints, use the distance constraint between two object points to fix the hinge rotation axis. For the hinge opening and closing angle constraint, an additional distance constraint is applied at a specific position between the two objects to fix the hinge joint opening and closing angle. The distance constraint point expression for the hinge joint opening and closing angle constraint is as follows:
[0076]
[0077] Where R(θ|a) is the rotation matrix corresponding to the rotation angle θ along the rotation axis a, and r0 represents the position of the joint action point on the object in the local coordinate system.
[0078] iii. Sliding joints
[0079] The sliding joint needs to dynamically generate two distance constraints to limit the relative translation of the object along the sliding axis. Project two points on object j along the sliding axis onto the sliding axis of object i, and establish two distance constraints between the projected points on object i and the two points on object j to realize a sliding joint without torsion restriction. For the case of torsion restriction, similar to a fixed hinge joint, an additional distance constraint is dynamically applied to a specific position between the two objects to achieve it. The expression of the sliding joint projection point is as follows:
[0080] y′ i =y i
[0081] y′ j =y j +λn j
[0082] where y i is the position of the joint action point on object i, y′ i Yes i The projection point of the sliding axis line from point j to the center of object j, where n j is the unit direction vector of the sliding axis, and λ measures the relative translation of the sliding joint.
[0083] (2) Perform a global solution for all constructed constraints and use a semi-implicit continuous iteration method to solve the position of the constraint action point after being affected by the constraints.
[0084] First, the global energy equation is modeled for the constraint energy between the dual action points to solve the optimization problem, which is expressed as:
[0085]
[0086] Where y represents the coordinates of all action points, is the time derivative of y, and E represents the y is the momentum potential energy, ψ is the sum of various constraint potential energy terms, A i represents the affine matrix of object i, r iξ represents the position of the action point of constraint ξ on object i in the coordinate system of object i. Constraint ξ includes elasticity, contact and other constraints. p i Indicates the center coordinates of object i, y iξ Represents the global coordinates of the point of action of constraint ξ on object i.
[0087] Then, the semi-implicit continuous iteration method is used to solve the nonlinear global energy constraint. At the kth iteration, the iteration formula is expressed as:
[0088]
[0089] Among them Ψ ξ represents the constraint potential energy of constraint ξ, represents the derivative of Ψ, express The sign of the positive term is The negative sign of d ξ =‖y jξ -y iξ ‖ represents the distance of constraint ξ. Substituting into the global energy equation, we can get the following formula:
[0090]
[0091] in represents the coordinates of object i in constraint ξ at the kth iteration, Can be regarded as a constant, the specific expression is as follows:
[0092]
[0093] Where h is the simulation time step, m i is the mass of object i.
[0094] (3) Rotate all constraint variables to their initial state (i.e., multiply the coordinates of all action points in the constraint coordinate system by the inverse of the current object rotation matrix), keeping p i The affine matrix is unchanged to avoid the influence of the initial value of the affine matrix on the convergence speed. Then, the Newton method is used for local solution, that is, the local rigid constraints of each object are solved to ensure the maximum stiffness of each object and minimize the distance of all constraint points on each object. The affine matrix of each object is calculated and updated, which requires the calculation of the inverse of the 9×9 matrix. In the local solution, the potential energy to ensure the position of the constraint point and the orthogonal potential energy Q to ensure the rigidity of the object are solved. The local optimization objective function is:
[0095]
[0096] in‖·‖ F represents the Frobenius norm. The Frobenius norm of a matrix is defined as the root of the sum of the squares of the absolute values of the elements of the matrix. i represents the position coordinates of object i before updating, ω iξ represents the weight of the potential energy term at the constraint point, δ represents the penalty coefficient, and I is the third-order unit matrix.
[0097] Finally keep A i The Newton method is used to iteratively update the center coordinates of the object and the position of the object according to the calculated affine matrix and the constraint points in the global solution.
[0098] The above update adopts a two-step strategy, such as Figure 4 shown.
[0099] (4) Repeat step (3) until the relative error converges to a preset value or a maximum number of iterations.
[0100] Iterative convergence can be measured by relative error, and its specific expression is:
[0101]
[0102] where ε k represents the relative error of the kth iteration, x 0 represents the initial position, x k represents the kth iteration position, x N Indicates the Nth iteration position.
[0103] (5) Update the solved position and rotation state of each object, and update the linear velocity and angular velocity of each object.
[0104] Experimental test:
[0105] The hardware platform used in this experiment uses an Intel i9-12900H sixteen-core CPU with a main frequency of 3.8GHz and an NVIDIA GeForce GTX 3080 graphics card with 16GB of memory. The system program is written in C++, and the CUDA language is used to accelerate the parallel computing part. The program is compiled and executed with the help of Microsoft Visual Studio 2019. Open source libraries such as Peridyno, OpenGL, and Qt are used in the development process. Linear BVH is used to accelerate collision detection.
[0106] This experiment tests the rigid body stacking collision effect at different orders of magnitude and compares it with the existing physics engine Bullet. The method of the present invention is more accurate than 10 4 In the case of large-scale objects, the simulation can be maintained at 40fps stably and robustly. Figure 5 Figure (a) shows the collision effect of 1000 objects under the Bullet engine, and Figure (b) shows the collision effect of 1000 objects achieved by the method of the present invention. Figure 6 The average time spent per time step for different numbers of objects is compared in the figure. It can be seen that when the number of objects reaches 10 5 After that, it can reach 3 orders of magnitude acceleration. In addition, this experiment also carried out unit tests on three basic joints, and constructed free joints, driven joints, and fixed joints respectively. The simulation effect is stable and normal. It can be concluded that the affine multi-body simulation using the method of the present invention is effective and stable, and the global loading calculation of the Hessian matrix is avoided, which is suitable for parallel implementation of modern GPUs and can improve the calculation efficiency.
[0107] The above implementation examples are only used to illustrate the technical solution of the present invention rather than to limit it. Ordinary technicians in the field can modify or replace the technical solution of the present invention with equivalents without departing from the spirit and scope of the present invention. The scope of protection of the present invention shall be based on the claims.
Claims
1. A real-time simulation method for an affine multi-body system based on projected peridynamics, comprising the following steps: (1) Construct dual constraints for the colliding objects in the multi-body system, including the constraint relationships between the dual action points, such as the constraint potential energy, collision, friction, and joints; (2) A global energy equation is modeled for all the constructed constraint relationships, and the nonlinear constraints in the equation are globally solved using a semi-implicit continuous iteration method. The position of each action point is calculated through iterative updates; (3) Use Newton's method to perform local solution, calculate the object's affine matrix and position coordinates, and update them; (4) Repeat the global solution and local solution process until the relative error converges to a preset value or reaches the maximum number of iterations; (5) Update the object's position, rotation state, linear velocity, and angular velocity.
2. The method according to claim 1, characterized in that In step (1), the energy constraint relationship is constructed for the constraint potential energy of the dual action point: Where Ψ represents the constraint potential energy, y represents the vector position of the action point, l represents the distance between the dual action points, ∈ represents the constraint strength, and C represents a linear function about y and l. The implicit continuous iteration of the energy constraint relationship is: Where i and j represent two colliding objects, y represents the position of the object, C represents the constraint, and k represents the kth iteration. Represents the distance between the dual action points at the kth iteration.
3. The method according to claim 1, characterized in that In step (1), the log-barrier function is used to construct the constraint relationship between the collision and friction between the dual action points: Among them B N Refers to the N-order polynomial barrier function, d represents the distance between the dual action points, Represents the maximum distance between the dual action points of the energy constraint.
4. The method according to claim 1, characterized in that In step (1), three joint constraints are constructed for the joints: Construct spherical joint constraints: The rotation of a spherical joint in three-dimensional space is fixed by applying three distance constraints on the three axes of the joint coordinate system. The distance constraint of each axis calculates the coordinates of two action points on the object, and the rotation matrix of the joint anchor point is used to ensure that the rotation of the object complies with the constraints of the spherical joint. Construct hinge joint constraints: The hinge joint uses a distance constraint between two points to fix the direction of the rotation axis, and an additional distance constraint is applied to the position set between the two objects to control the opening and closing angle of the hinge; Constructing sliding joint constraints: The sliding joint uses two distance constraints to limit the relative translation of the objects along the sliding axis; for torsion, an additional distance constraint is applied to the set position between the two objects to fix the rotation axis.
5. The method according to claim 2, characterized in that The global energy equation modeled in step (2) is: Where y represents the coordinates of all action points, is the time derivative of y, and E represents the y represents the momentum potential energy, Ψ represents the constraint potential energy, A i represents the affine matrix of object i, r iξ represents the local coordinates of the point of action of constraint ξ on object i, p i Indicates the center coordinates of object i, y iξ Represents the global coordinates of the point of action of constraint ξ on object i.
6. The method according to claim 5, characterized in that When the semi-implicit continuous iteration method is used in step (2) to globally solve the nonlinear constraints in the equation, the iterative formula for the kth iteration is: where ψ ξ represents the constraint potential energy of constraint ξ, represents the derivative of Ψ, express The sign of the positive term is The negative sign of d ξ Represents the distance between the dual action points.
7. The method according to claim 6, characterized in that In step (2), the iterative formula of the kth iteration is substituted into the global energy equation for global solution, and the result is: in represents the coordinates of object i in constraint ξ at the kth iteration, is a constant, and the expression is as follows: Where h is the simulation time step, m i is the mass of object i.
8. The method according to claim 1, characterized in that The specific steps of step (3) include: Perform an inverse transformation on the coordinates of all the action points of the object according to the rotation matrix of the current object, and rotate the object back to its initial state; The local rigid constraints of the object are optimized and solved using Newton's method. The optimization goal is to maximize the rigidity of the object while minimizing the distance error between all action points on the object, and solve the updated affine matrix. A two-step strategy is adopted for updating, including: the first step is to update the affine matrix and adjust the rotation state of the object; the second step is to keep the affine matrix unchanged and use Newton's method to iteratively update the center position of the object according to the dual action point in the global solution.
9. The method according to claim 1, characterized in that The optimization objective function of the local solution is: where ||·|| F represents the Frobenius norm, q i represents the position coordinates of object i before updating, p i represents the center position of object i, ω iξ represents the weight of the potential energy term at the dual action point, r iξ represents the local coordinates of the point of action of constraint ξ on object i, δ represents the penalty coefficient, A i represents an affine matrix and I is the third-order identity matrix.
10. The method according to claim 2, characterized in that The relative error in step (4) is: where ε k represents the relative error of the kth iteration, Ψ(x 0 ) represents the initial position x 0 The potential energy, Ψ(x k ) represents the kth iteration position x k The potential energy, Ψ(x N ) represents the Nth iteration position x N The potential energy is the potential energy after final convergence.