A semi-implicit iterative simulation method for hyperelastic bodies based on peridynamics

Abstract

The present invention discloses a semi-implicit iterative simulation method for hyperelastic bodies based on peridynamics, which is a simulation method with high precision, high robustness, and high performance. The peridynamics theory is used to model the energy of hyperelastic materials, and the dimensionality of the particle neighborhood of codimension objects is extended. The semi-implicit iterative method is used for elastic solution, which can ensure convergence and stably handle the overshoot problem in the iterative process through an adaptive step size. The present invention can simulate the dynamic processes and interaction behaviors of various codimension materials such as hyperelastic materials, hyperelastic fabrics, and thin shell-like objects with high precision, high robustness, and in real time.

Description

Technical Field

[0001] The present invention belongs to the field of computer graphics, and particularly relates to a semi-implicit iterative simulation method for hyperelastic bodies based on peridynamics. Background Art

[0002] Hyperelastic objects have very unique material properties, and their stress-strain relationships and strain-geometric relationships often have non-linear characteristics. Therefore, efficiently and stably simulating three-dimensional hyperelastic objects, fabrics, and thin shells, etc., which are co-dimensional hyperelastic bodies, is an important issue in the field of physical simulation in graphics. Although the finite element method based on the elastic stress method can provide accurate simulation solutions for hyperelastic bodies, due to its dependence on specific linearization methods, the assembly of the global stiffness matrix, and the solution of sparse linear systems, it cannot meet the requirements of efficient and real-time physical simulation in graphics. The position-based simulation method re-describes the dynamic general constraints using particle positions and uses an iterative method to solve the global constraint linear system, which has high efficiency and real-time performance. However, since it focuses on the constraint description of the dynamic equation, its control equation does not have real physical meaning, so it cannot meet the accurate simulation of specific material parameters.

[0003] The reasons for the inability to combine the accuracy and efficiency of traditional hyperelastic body simulations mainly include two aspects: (1) The implicit time integration method based on traditional finite elements involves the assembly of the global stiffness matrix and the solution of large-scale sparse linear systems. Due to the existence of global quantities, it is difficult to directly parallelize, and the linearization method based on Newton's method involves the calculation of the Hessian matrix, which is extremely time-consuming, so this method does not have high efficiency and real-time performance; (2) The position-based method uses particle position constraints to describe the control equation and uses Jacobi iteration or Gauss-Seidel iteration to solve it. It does not require stress calculation and linearization and can perform simulations efficiently. However, since the general constraints are directly assembled into a linear system for solution, the inconsistency of the stiffness between the constraints will cause the sparse linear system to be ill-conditioned, so for some material parameters, it cannot guarantee convergence or has a large error. Summary of the Invention

[0004] The technical problem solved by the present invention: Aiming at the fact that traditional simulation methods cannot efficiently and accurately simulate three-dimensional and two-dimensional co-dimensional hyperelastic objects, a high-precision, high-robustness, and high-performance semi-implicit iterative simulation method for hyperelastic bodies based on peridynamics is provided. Since it belongs to the particle method, it can avoid the problem of inefficient operation of the finite element method based on grid discretization; use a semi-implicit solver that guarantees convergence, and at the same time use an adaptive step size to solve the overshoot problem. The invention provides a unified solution for three-dimensional objects and two-dimensional co-dimensional objects based on hyperelastic energy modeling and iterative solution, and finally can achieve efficient, accurate, and high-robustness simulation of hyperelastic bodies.

[0005] The present invention is a simulation method with high precision, high robustness, and high performance. It uses the peridynamics theory to perform energy modeling of hyperelastic materials. For codimension objects, the vertex normal is used to introduce virtual particle bonds to expand the dimension of the particle neighborhood, solving the degeneracy problems of the particle shape matrix and deformation gradient matrix in the simulation of codimension objects, thereby providing a unified solution for two-dimensional and three-dimensional objects. The present invention uses a semi-implicit iterative method for elastic solution, decomposes the Piola-Kirchhoff stress tensor into a positive definite term and a negative definite term to construct an iterative formula that ensures convergence, and introduces an adaptive step size to solve the overshoot problem in the global iterative convergence process. The present invention can efficiently and real-time simulate the dynamic processes and interaction behaviors of various hyperelastic materials, hyperelastic fabrics, and thin-shell-like objects, etc., which are codimension materials. Since it does not involve the assembly of a global linear system and the operation of inner product, it has good parallelization characteristics.

[0006] The specific technical solution of the present invention is as follows:

[0007] A semi-implicit iterative simulation method for hyperelastic bodies based on peridynamics, as Figure 1 shown, includes the following steps:

[0008] Step 1: Read in the mesh parameters of the three-dimensional geometric model of the simulation object, initialize the vertex positions of the mesh model as the particle positions of the peridynamics hyperelastic body particles, initialize the particle velocity to zero, establish the particle neighborhood relationship according to the mesh topology, and save the particle positions at the initial moment as the rest shape;

[0009] Step 2: If the current simulation time step is not the initial simulation time step, use the results of the hyperelastic solver and the time integrator from the previous simulation time step to update the particle positions and velocities. The obtained results are the results of the previous simulation time step and can be used for output or visualization; use the time integrator to perform time integration on the external force field (such as gravity, load, etc.) applied to the simulation object specified by the user, update the current particle velocity and position, and obtain the time integrator result of the current simulation time step as the input to the hyperelastic solver;

[0010] Step 3: Calculate the particle shape matrix and deformation gradient matrix according to the current particle positions, particle velocities, and the particle neighborhood relationship and rest shape stored in Step 1;

[0011] Step 4: According to the specified constitutive model of the hyperelastic body material, decompose the Piola-Kirchhoff stress tensor into a positive definite term and a non-positive definite term, and calculate the positive definite term and non-positive definite term of the above tensor according to the particle shape matrix and deformation gradient matrix obtained in Step 3;

[0012] Step 5: Calculate the coefficient matrix and source term of the sparse linear system in the form of Jacobi iteration for implicit time integration according to the results obtained in Step 3 and Step 4;

[0013] Step 6: Perform Jacobi iteration for a set number of steps based on the coefficient matrix and source term obtained in Step 5 to solve the relevant sparse linear system, and calculate the particle positions after this iteration.

[0014] Step 7: Calculate the particle position gradient based on the difference between the particle positions after one iteration in Step 6 and the particle positions before iteration, and calculate an adaptive iteration step size that can avoid overshoot problems using the one-dimensional elastic energy of the hyperelastic material. Update the particle positions using this adaptive step size and the position gradient.

[0015] Step 8: Repeat Steps 3 to 7 and increment the iteration step number until the position gradient residual is less than the set error or the algorithm reaches the specified iteration step number. The result obtained is the result of the hyperelastic solver.

[0016] Step 9: Increment the simulation time step and repeat the operations between Steps 2 - 8 until the termination time or the set simulation termination condition is reached, and complete the entire simulation process.

[0017] In Step 3, if the object to be simulated is a two-dimensional codimension object, perform dimension expansion on the hyperelastic body particles as follows:

[0018] ① During initialization, pre-compute and store the angle-weighted average normal of the vertex i where the particle is located where J represents the serial number of the triangle adjacent to vertex i, ξ represents the thickness from the outer boundary layer to the neutral layer of the two-dimensional codimension object, n J represents the outer normal of triangle J, and θ J represents the interior angle of triangle J at vertex i;

[0019] ② At each simulation time step, calculate the angle-weighted average normal of the vertex i where the particle is located after deformation where the superscript * represents the corresponding quantity in the deformed configuration;

[0020] ③ The shape matrix expression of particle i is The deformation gradient matrix expression is Let x represent the position of the particle in the undeformed initial configuration, and x i represent the position of particle i in the undeformed initial configuration; let represent the position of particle i in the deformed configuration at iteration step k; j represents the adjacent particle of the particle; ω ij represents the scalar weight between particles i and j; represents the Kronecker product.

[0021] In Step 4, the decomposition format of the Piola - Kirchhoff stress tensor is: where To deform the gradient matrix through standard SVD decomposition The obtained diagonal matrix, U i and V i Are the left and right rotation matrices corresponding to the above SVD decomposition; Is the diagonal part of the Piola - Kirchhoff stress tensor calculated using the deformed gradient diagonal matrix ; Is the positive definite part of the above stress tensor diagonal matrix; Is the negative definite part of the above stress tensor diagonal matrix.

[0022] In step 5, the operations are as follows:

[0023] ① The Jacobi iteration form of implicit time integration is expressed as:

[0024] ② The coefficient matrix of the sparse linear system is expressed as:

[0025] ③ The source term of the sparse linear system is expressed as:

[0026] In the above equations, y i , y j Are the positions of particles i and j to be solved in the iteration form; Represents the deformed configuration position of particle i at iteration step k; Represents the deformed configuration position of particle i at simulation time step t; Represents the velocity of particle i at simulation time step t; m i Is the mass of particle i; h is the time step; the scalar V i , V j Are the volumes of particles i and j respectively; Where r i Is the spherical neighborhood radius of particle i; I represents the identity matrix.

[0027] In step 6, the Jacobi iteration process operates as follows:

[0028] ① The sparse linear system solved in the Jacobi iteration process is: for any particle i,

[0029] ② Solve using the iterative method, where the initialization The iteration process is Until the residual converges or m reaches the specified maximum number of iteration steps.

[0030] In step 7, the calculation of the adaptive step size and the update operation of the particle position are as follows:

[0031] ① Calculate the adaptive step size to avoid overshoot problems: where E is the energy density function;

[0032] ② Update the particle position using the adaptive step size:

[0033] The present invention also provides a server, which includes a memory and a processor. The memory stores a computer program, and the computer program is configured to be executed by the processor. The computer program includes instructions for executing the steps in the above method.

[0034] The present invention also provides a computer-readable storage medium, on which a computer program is stored. The computer program, when executed by a processor, implements the steps of the above method.

[0035] Compared with the existing methods, the present invention has the following advantages:

[0036] (1) The present invention uses the peridynamics theory to model hyperelastic materials, and introduces virtual bonds for two-dimensional co-dimension objects to expand the dimension of the particle neighborhood, solving the problem of matrix degradation in the simulation of co-dimension objects, thereby providing a unified solution for three-dimensional and two-dimensional co-dimension hyperelastic objects;

[0037] (2) The present invention proposes a semi-implicit iterative solver. By decomposing the stress matrix and assembling the corresponding iterative formula, global convergence can be guaranteed. Since it does not rely on other linearization methods and the calculation of global quantities, it has the characteristics of high efficiency and real-time;

[0038] (3) The overall solution proposed by the present invention has very high stability. By using the adaptive step size, the overshoot problem in the iterative convergence process is solved. It can simulate hyperelastic materials with a wide range of stiffnesses, and can be widely applied to various hyperelastic material constitutive models and user-defined materials, so as to simulate the required material properties;

[0039] (4) Each step in the present invention has high parallel characteristics, so it is suitable for modern GPU architectures. Brief Description of the Drawings

[0040] Figure 1 is the implementation flowchart of the present invention.

[0041] Figure 2 is the schematic diagram of the unified solver. Detailed Embodiments

[0042] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be described in detail below through specific implementations and drawings, but it does not constitute a limitation to the present invention.

[0043] 1. The hardware platform of the method of the present invention uses an Intel i9-10900K octa-core CPU with a main frequency of 3.5 GHz, an NVIDIA GeForce RTX A4000 graphics card with a video memory of 12 GB. The system program is written in C++. For the parallel computing part, the CUDA language is used for acceleration, and the program is compiled and executed with the help of Microsoft Visual Studio 2019. Open source libraries such as OpenGL and Eigen are used during the development process.

[0044] 2. The present invention is based on peridynamics theory, and its specific form is as follows:

[0045] First, consider the implicit time integration formula:

[0046]

[0047] where represents the position of particle i at time step t after deformation, represents the velocity of particle i at time step t, represents the indirect velocity, M i is the mass matrix of the particle, and h is the time step. Given the reference position x i of particle i and the neighborhood {j}, in the peridynamics framework, the internal force f between particles is determined by the particle bond:

[0048]

[0049] where T is the stress state basis function, V is the volume, and the peridynamics constitutive model is expressed as:

[0050]

[0051] where P is the first Piola-Kirchhoff stress tensor, which is a function of the deformation gradient matrix F, and F is expressed as:

[0052]

[0053] where

[0054]

[0055] Therefore, the implicit integration formula can be transformed from the above formula to:

[0056]

[0057] where the external force is irrelevant to the hyperelastic solution and can be pre-computed at the time integrator. Therefore, after transformation, the hyperelastic solver finally uses the following equation as the implicit time integration formula to be solved:

[0058]

[0059] 3. Figure 2 is the flowchart of the semi-implicit iterative simulation solver for hyperelastic bodies based on peridynamics. In the initialization stage before simulation, first read the three-dimensional geometric model, establish peridynamic hyperelastic body particles one by one corresponding to the model grid vertices, and store the position x i of the initial reference configuration of particle i. According to the topological constraints, establish and store the adjacency relationship {j} for each particle i, preprocess the topological relationship, including the relationship of the voxels (triangular mesh or tetrahedral mesh) where the corresponding vertices of the particles are located and the relationship of the vertices contained in the voxels, and set the initial time step t = 0.

[0060] 4. After the simulation starts, if it is not the initial time step t = 0, first update the particle positions solved by the hyperelastic solver in the previous step to the current positions, use the position difference between the hyperelastic solver and the time integrator in the previous step divided by the time step as the velocity increment in the previous step to update the current particle velocities. The above variables satisfy the implicit time integration equation and can be further used for output or visualization; perform time integration: Before the start of hyperelastic solution in each frame, use gravity as the external force to update the velocities of the particles at the next moment, and use these velocities to update the positions of the particles at the next moment. This pre-computed velocity and position are used as the inputs to the hyperelastic iterative solver. The above process can be expressed as:

[0061] ① If t≠0, update the current positions to the results solved by the hyperelastic solver:

[0062] ② Update the intermediate velocities and intermediate positions pre-computed by the time integrator: Increment t.

[0063] 5. Perform hyperelastic solution. In the initial stage, set the iteration step k = 0. For any particle i, set the initial position to the result at time step t, that is, initialize In each subsequent iteration step k, calculate the corresponding strain measures:

[0064] ① For three-dimensional hyperelastic bodies, for each particle i, calculate the shape matrix and the deformation gradient matrix:

[0065] ②For two-dimensional co-dimension objects, for each particle i, first calculate the reference vertex normal of the vertex corresponding to the particle and the current vertex normal: where the triangle where the vertex corresponding to the particle is located is preprocessed into a topological relationship, so it can be queried separately for each particle i, and then calculate the non-degenerate form of the shape matrix and the deformation gradient matrix after the particle neighborhood dimension is extended: Calculate the deformation gradient F i After that, perform SVD decomposition on it:

[0066] 6. Use the diagonal form of the deformation gradient Calculate the decomposed form of the Piola - Kirchhoff stress tensor:

[0067] ①For three-dimensional hyperelastic bodies, use the neo-Hookean hyperelastic model The diagonal form of the Piola stress tensor is expressed as: where In the above formula, I1 = tr(F T F), I2 = F T F:F T F, 13 = det(F T F);

[0068] ②For two-dimensional co-dimension objects, use the tensile strain energy and the bending strain energy to model the constitutive model respectively. The tensile strain energy is expressed as: where And select n = 3. The diagonal form of the Piola stress tensor is expressed as: where The basis function of the bending stress state is expressed as which will directly contribute to the coefficient matrix of the implicit time integration formula.

[0069] 7. Calculate the coefficient matrix and the source term of the sparse linear system for Jacobi iteration:

[0070] ①For three-dimensional hyperelastic bodies, for any particle i, the coefficient matrix is calculated as follows:

[0071]

[0072] The source term is expressed as:

[0073] ②For two-dimensional co-dimension objects, for any particle i, the coefficient matrix is calculated as follows:

[0074]

[0075] The source term is expressed as:

[0076] 8. For each particle i, perform Jacobi iteration for the specified maximum number of iterations, with the initial conditions being: The iteration formula is:

[0077]

[0078] 9. For each particle i, calculate the step size that can avoid overshooting problems: Update the particle position using the above step size as: Increment k.

[0079] 10. Wait until the hyperelastic iterative solver converges or the iteration step k reaches the maximum number of iterations, then return the result to the time integrator to re-enter the loop, reset the temporary variables, start stepping, and enter the calculation of the next frame until the simulation ends.

[0080] In the present invention, the calculation process and iteration process of the relevant quantities of particle i can both be performed in parallel on the GPU. The hyperelastic constitutive model can be replaced with the corresponding constitutive according to actual needs and calculated according to the decomposition method of the present invention.

[0081] The above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit them. Those skilled in the art can modify or equivalently replace the technical solutions of the present invention without departing from the spirit and scope of the present invention. The protection scope of the present invention shall be defined by the claims.

Claims

1. A semi-implicit iterative simulation method for hyperelastic bodies based on peridynamics, the steps of which include: 1) Read in the three-dimensional geometric model of the object to be simulated. The three-dimensional geometric model is a mesh model. Take the vertex positions of the three-dimensional geometric model as the particle positions of the peridynamics hyperelastic body particles, initialize the particle velocities to zero, establish the particle neighborhood relationship according to the mesh topology of the three-dimensional geometric model, and save the particle positions at the initial moment as the static shape; 2) If the current simulation time step is not the initial simulation time step, update the particle positions and velocities using the simulation results of the previous simulation time step; then use the time integrator to perform time integration on the external force field applied to the simulation object specified by the user, and update the particle velocities and positions; 3) Calculate the shape matrix and deformation gradient matrix of the particles according to the current particle positions, particle velocities, and the particle neighborhood relationship and static shape stored in step 1); 4) According to the specified constitutive model of the hyperelastic material, decompose the Piola-Kirchhoff stress tensor into a positive definite term and a non-positive definite term, and calculate the positive definite term and the non-positive definite term according to the shape matrix and deformation gradient matrix obtained in step 3); 5) Calculate the coefficient matrix and source term of the sparse linear system in the form of Jacobi iteration for implicit time integration according to the results obtained in steps 3) and 4); 6) Perform Jacobi iteration to solve the sparse linear system according to the coefficient matrix and source term obtained in step 5), and calculate the particle positions after this iteration; 7) Calculate the position gradient of the particles according to the difference between the particle positions obtained by iteration in step 6) and the particle positions before iteration, and calculate the adaptive iteration step size using the one-dimensional elastic energy of the hyperelastic material. Update the particle positions using the adaptive iteration step size and the position gradient; 8) Iteratively execute steps 3) to 7) until the gradient residual of the position gradient is less than the set error or the number of iterations reaches the specified number of iteration steps; 9) Increment the simulation time step, and repeat steps 2) to 8) until the termination time is reached or the set simulation termination condition is met, and the simulation process is completed.

2. The method according to claim 1, wherein If the object to be simulated is a two-dimensional codimension object, the dimensions of the hyperelastic body particles are extended. The specific method is as follows: during initialization, pre-compute and store the angle-weighted average normal of the vertex i where the particle is located. At each simulation time step after initialization, calculate the angle-weighted average normal of the vertex i where the deformed particle is located. Express the shape matrix of particle i as Express the deformation gradient matrix as x represents the position of the particle in the undeformed initial configuration, x i represents the position of particle i in the undeformed initial configuration; represents the position of particle i in the deformed configuration at simulation time step k; j represents the adjacent particle of particle i; ω ij represents the scalar weight between particles i and j; represents the Kronecker product.

3. The method according to claim 2, wherein The decomposition format of the Piola-Kirchhoff stress tensor is: Among them, is the diagonal matrix obtained by performing SVD decomposition on the deformation gradient matrix , and U i and V i are the left and right rotation matrices corresponding to the SVD decomposition; is the diagonal matrix of the Piola - Kirchhoff stress tensor calculated using the diagonal matrix of the deformation gradient ; is the positive definite part of the diagonal matrix of the above stress tensor; is the negative definite part of the diagonal matrix of the above stress tensor.

4. The method according to claim 3, characterized in that, The specific implementation method of step 5 is as follows: Express the Jacobi iteration form of implicit time integration as: Then the coefficient matrix of the sparse linear system The source term of the sparse linear system y i 、y j are the positions of particles i and j to be solved in the iteration form; represents the deformed configuration position of particle i at simulation time step t; represents the velocity of particle i at simulation time step t; m i is the mass of particle i; h is the time step; the scalar V i ,V j are the volumes of particles i and j, respectively; where r i is the spherical neighborhood radius of particle i; I represents the identity matrix.

5. The method according to claim 1 or 2 or 3, characterized in that, The adaptive step size of particle i at iteration step k+1 is where E is the energy density function, represents the deformed configuration position of particle i at iteration step k, and y i is the position of particle i to be solved; the particle position of particle i after updating with the adaptive step size at iteration step k+1 is 6. The method according to claim 1 or 2 or 3, characterized in that, The simulation results include the results of the hyperelastic solver and the results of the time integrator; the results of the time integrator are to perform time integration on the external force field specified by the user using the time integrator to update the current particle velocities and positions; taking the results of the time integrator as the input of the hyperelastic solver and executing steps 3) to 7) to obtain the results of the hyperelastic solver.

7. A server, characterized in that, It includes a memory and a processor. The memory stores a computer program, and the computer program is configured to be executed by the processor. The computer program includes instructions for executing each step in any one of claims 1 to 6.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, the steps of any one of claims 1 to 6 are implemented.

Images