Real-time interactive high-precision cable simulation method and system

Abstract

A real-time interactive high-precision rope simulation method and system, including an initialization stage and a simulation stage, the initialization stage is based on the initial state of the rope to be simulated, the continuous rope is discretely processed and stored, and the simulation parameters are set; the simulation stage uses an improved projection mechanics algorithm, based on the motion constraints derived based on the Cosserat rope model, uses matrix decomposition and cache acceleration, and outputs the motion state of each logical frame of the rope on the basis of real-time, realizes the physical precision enhancement of the rope simulation. The improved projection mechanics algorithm is used in the application, the heading and angular velocity are added to the original projection mechanics framework, and the matrix decomposition and cache acceleration are used, so that the running state of the rope after stress and collision of various materials can be simulated in real time, and the accuracy close to that of the offline simulation of the commercial finite element software can be ensured, and the reality of the simulation of the rope stretching, shearing, bending and twisting phenomena in the real-time interactive application is improved.

Description

Technical Field

[0001] This invention relates to a technology in the field of materials science, specifically a real-time interactive method and system for simulating cables with an average error of less than three percent. Background Technology

[0002] Ropes and cables are among the most common devices in daily life. Simulating the motion of ropes and cables using computers has always been a subject of considerable interest. In related technologies, simulation algorithms and techniques with high physical accuracy for ropes and cables have made significant progress; for example, the commercial software Abaqus employs the finite element method (FEM). A common characteristic of these classic algorithms is that they derive the equations of motion based on forces and solve these nonlinear partial differential equations using the finite element method (FEM) or the finite difference method (FDM), typically yielding physically accurate results. However, the main drawback of these force-based methods is the low computational efficiency in solving rigid differential equations, making real-time calculations impossible.

[0003] Some research focusing on real-time interactive applications in the graphics domain has proposed fast simulation methods that sacrifice physical accuracy for higher performance in order to reduce computational overhead, such as position-based dynamics (PBD). These methods are commonly used to simulate cables in real-time graphics applications such as games, offering high computational efficiency. However, their drawbacks include unclear physical meanings of parameters and room for improvement in physical accuracy. Therefore, performing high-precision real-time cable physics simulations on PCs with limited computing resources remains a challenge. Summary of the Invention

[0004] This invention addresses the problem of insufficient physical accuracy in rope simulation in existing real-time graphics applications by proposing a real-time interactive high-precision rope simulation method and system. Based on an improved projection mechanics algorithm, orientation and angular velocity are incorporated into the original projection mechanics framework, and matrix decomposition and caching are used for acceleration. This method can simulate the stress and post-collision running state of ropes made of various materials in real time, while ensuring accuracy similar to offline simulations using commercial finite element software. This improves the realism of simulating rope tension, shearing, bending, and torsion phenomena in real-time interactive applications.

[0005] This invention is achieved through the following technical solution:

[0006] This invention relates to a real-time interactive high-precision rope simulation method, comprising an initialization stage and a simulation stage, wherein: In the initialization stage, the continuous rope is discretized and stored according to the initial state of the rope to be simulated, and then the simulation parameters are set; In the simulation stage, an improved projection mechanics algorithm is used, based on the motion constraints derived from the Cosserat rope model, matrix decomposition and caching are used to accelerate the process, and the motion state of each logical frame of the rope is output in real time, thereby enhancing the physical accuracy of the rope simulation.

[0007] This invention relates to a real-time interactive rope simulation system for implementing the above-described method, comprising: an interaction module, a simulation module, a rendering module, an input device, and a display device, wherein: the interaction module acquires user input in real time and changes the motion state of the rope or other objects in the virtual scene based on the input; the simulation module simulates the motion state of the rope in real time according to the rope simulation method; the rendering module reconstructs the geometric model of the rope in the virtual scene in real time based on the rope motion state and parameters and outputs a rendered image; the input device is used to receive user input; and the display device outputs the rendered rope motion animation.

[0008] Technical effect

[0009] This invention utilizes an improved projection mechanics algorithm to decompose the constraints affecting cable motion into potential energy constraints derived from the Cosserat model and dynamically generated constraints such as collision and animation constraints. It also generalizes the performance bottleneck in solving potential energy constraints in projection mechanics—solving high-dimensional linear equations—to solving a problem of sparse matrices and vector products. Thanks to the unchanged sparse matrix after the decomposition, the solution process is greatly accelerated by decomposing the sparse matrix and caching the decomposition results, thus achieving real-time performance. A position-based dynamics algorithm is used to solve the remaining constraints. The technical results are: achieving accuracy close to that of non-real-time simulations in commercial finite element software, with an average error of less than 3%, while improving performance by an order of magnitude, thus meeting real-time requirements; compared with simulation methods commonly used in real-time graphics applications such as games, the accuracy is greatly improved, and the parameters and weights have realistic physical meaning. Attached Figure Description

[0010] Figure 1 This is a flowchart of the method of the present invention;

[0011] Figure 2 This is a flowchart of the initialization process of the present invention;

[0012] Figure 3 This is a schematic diagram of the reconstructed rope model of the present invention;

[0013] Figure 4 This is a flowchart of the simulation process of the present invention;

[0014] Figure 5 This is a system diagram of the present invention;

[0015] Figures 6-9 This is a schematic diagram illustrating the effect of an example. Detailed Implementation

[0016] like Figure 1As shown in the figure, this embodiment relates to a real-time interactive high-precision rope simulation method, which includes an initialization stage and a simulation stage. In the initialization stage, the continuous rope is discretized and stored according to the initial state of the rope to be simulated, and then the simulation parameters are set. In the simulation stage, an improved projection mechanics algorithm is adopted, the orientation of discrete rope segments is introduced, and the potential energy constraints derived from the Cosserat rope rod model are used. Matrix decomposition and caching are used to accelerate the simulation, and the motion state of each logical frame of the rope is output in real time, thereby enhancing the physical accuracy of the rope simulation.

[0017] like Figure 2 As shown, the initialization phase specifically includes:

[0018] The first step is to obtain the control points input by the user to define the center line of the rope: the user inputs the rope radius and control points according to the initial shape of the rope in their mind; after obtaining the control points in three-dimensional space, in this embodiment, a Catmull-Rom spline curve with C1 continuity that passes through all control points is selected as the center line of the rope.

[0019] The second step is to obtain the sampling frequency N of the user input: According to the definition of the Cosserat theoretical model, let the centerline of the rope of length L be... By uniformly sampling N points along the center line according to the arc coordinates s: r(0), r(L / N-1), r(2L / N-1)...r(L), the cable is divided into N-1 segments, each segment consisting of two vertices {x}. i x i+1} and a quaternion u representing rotation i Definition. Using u i After rotating the world coordinate system, we obtain the local coordinate system in the Cosserat rope model. Then, the N vertices are particles, and the rope is divided into N-1 segments, which are discrete rope segments, serving as data structures stored in the computer.

[0020] The particle has mass, position, and velocity attributes, with the initial position set at the sampling point and the initial velocity being 0. The discrete rope segment has orientation, angular velocity, geometric moment of inertia, and initial length attributes, with the initial orientation set as the tangent direction of the rope centerline at the starting particle of the discrete rope segment, the initial angular velocity being 0, and the initial length being the distance between the two corresponding vertices of the discrete rope segment.

[0021] The third step involves the user inputting the geometric and material parameters of the cable: such as... Figure 3 As shown, based on the cross-sectional radius input by the user, a circular surface with the corresponding radius is drawn at each particle. The normal direction of this circular surface is the tangential direction of the center line at that point. The points on these circular surfaces are connected in sequence to form a triangular mesh, thereby restoring the geometric model of the rope.

[0022] The geometric parameters include the initial length and the cross-sectional radius.

[0023] The material parameters mentioned include density, Young's modulus, and Poisson's ratio.

[0024] The geometric moment of inertia of the discrete rope segment is obtained by multiplying the rope density input by the user by the average length of the discrete rope segment with the particle as the starting and ending point, to obtain the mass attribute of each particle, and then the geometric moment of inertia of the discrete rope segment can be calculated based on the cross-sectional radius.

[0025] The fourth step is to accept the simulation parameters input by the user, including: the total duration of the simulated motion, the time step of each change in the motion state of the simulated cable, and the number of iterations for each state change.

[0026] The fifth step is to store the geometric parameters and material parameters of the rope, as well as information on particles and discrete rope segments.

[0027] like Figure 4 As shown, the specific steps of the simulation phase are as follows:

[0028] Step 1) Divide the entire process of simulating cable movement into several logic frames to be calculated according to the settings.

[0029] The time step set for the interval between adjacent logical frames.

[0030] Step 2) Within each logical frame, input the motion state of the cable in the previous frame, and use the improved projection mechanics algorithm to calculate the motion state of the cable in that frame, specifically including:

[0031] 2.1) Predicting the position and orientation of the rope: Based on the particle's velocity and position in the previous frame, and the net external force acting on it in this frame, predict the particle's position in this frame, specifically: in: For the predicted positions of all particles in this frame, x (t) The positions of all particles in the previous frame are given by h, where h is the set time step, and v is the position of all particles in the previous frame. (t) Let f be the velocity of all particles in the previous frame, M be the diagonal matrix consisting of the masses of all particles, and f be the velocity of all particles in the previous frame. ext The net external force acting on the particle.

[0032] 2.2) Based on the angular velocity and orientation of the discrete rope segment in the previous frame and the net external torque acting on it in this frame, predict the orientation of the discrete rope segment in this frame, specifically: in: For the orientation of this frame of the discrete rope segment, u (t) For the orientation of the discrete rope segment in this frame, ω (t)ω represents the angular velocity of the discrete rope segment, obtained from the motion state of the previous frame; τ represents the net external torque, obtained from user input, collision handling, or animation design. (Diagonal matrix) J represents the geometric moment of inertia of the rope segment along the center line. Since the quaternion is treated as a four-dimensional vector in the linear solution process, and based on the definition of a quaternion, we know that... n =l n ρdiag(0, J1, J2, J3), where l n Let J1, J2, and J3 be the initial length of the rope segment. Since the entire rope is uniform, the same density ρ is used throughout. J1, J2, and J3 are the geometric moments of inertia of the cross-sections; when the cross-section is circular, this can be determined by... Figure it out.

[0033] 2.3) Initialize Constraints: Constraints include potential energy constraints derived from the Cosserat model, as well as collision constraints and animation-driven constraints added to enhance the realism and interactivity of the physics. The potential energy constraints remain unchanged after the cable parameters are determined and only need to be initialized when the cable parameters are modified. Collision constraints and animation-driven constraints change with the movement of the cable and user interaction, and need to be re-initialized in each frame. Specifically, this includes: i) Initializing the constraints derived from the stretching and shear potential energy of the Cosserat model. ii) Initialize the constraints derived from the bending and torsional potential energies of the Cosserat model as follows: iii) Initialize collision constraints as and iv) Initialize animation constraints as Among them: W sEi (q, p) i ) represents the current particle p i The difference in tensile and shear potential energy between the current state and the target state, where min represents the constraint minimizing this potential energy difference; Γ n The strain, which measures the degree of tension and shear, is equal to the direction vector of the particle on the tangent of the center line, minus the orientation direction vector of the discrete rope segments starting from that particle.

[0034] The direction vector on the tangent of the centerline is approximately equal to the position difference between adjacent particles within the discrete unit divided by the initial length. therefore Where: Im(q) represents the imaginary part of the four elements q; W BTi (q, p) i ) represents the current discrete rope segment p i The bending and torsional potential energy difference between the current state and the target state; min represents the constraint being to minimize this potential energy difference; Ω nThe strain used to measure the degree of tension and shear is equal to the rate of change of curvature between the discrete rope segment and its adjacent discrete rope segments, i.e. c j The strength of the collision constraint; δ j It is an indicator function, which is 1 when the particle's collider is still in contact with other colliders, and 0 when they are no longer in contact; R is the radius of the cable's cross-section. Let j be the position of particle j in three-dimensional space. c represents the position of the particle projected onto the surface of other colliders. j Strength of the animation constraint The position set for the particle within this animation frame.

[0035] 2.4) Use the position-based dynamics (PBD) method to solve animation constraints: determine the animation keyframe time, keyframe state information and interpolation method, and calculate the state information of each logical frame at the time point.

[0036] The state information includes: the position of the particle controlled by the animation and the orientation of the discrete rope segment. The solution method is similar to Macklin's solution of fixed point constraints in the Extended Position Dynamics Method (XPBD), except that the target fixed point position is set as the position of the particle controlled by the animation in the state information.

[0037] 2.5) For all constraints derived from potential energy, solve for the projected position of the constraint, i.e., the position that satisfies constraint i and has the minimum generalized distance from the current motion state q. Specifically, solve for the position p under constraint i that minimizes the distance between the current motion state q and a certain position p. i generalized distance The minimum value pi is the projection position of the constraint, where A, B, and S are constant matrices constructed based on the constraint, and w i To constrain the corresponding non-negative weights, δ C The indicator function is defined if and only if the projection position p i It is 0 when the constraint is satisfied, and infinite when the constraint is not satisfied.

[0038] The solution method for the projection position under tensile and shear constraints is as follows: at the projection point pi that satisfies the constraints of tensile and shear potential energy, the strain, which measures the degree of tensile and shear stress of the discrete element, is 0, i.e. Where: particle position change These represent the positions of the two particles in the discrete element during this iteration, and the normal to the cross-section of the discrete element. Since the particle's position and the rope segment's orientation are independent, the problem can be further divided into two optimization problems: optimizing the particle's position. And optimization of the local orientation of the rope segment: The solution for particle position optimization is in Δx n When d3 is reached, the positional difference of the particles is parallel to the normal of the cross-section, and the distance between the particles is the same as in the initial state, meaning the tensile deformation of this discrete element is 0. The solution for optimizing the local orientation of the rope segment is... When the shear deformation of the discrete element is 0, it is obtained that... To represent from Δx n The quaternion of d3 rotation. Therefore, the solution for the tension and shear potential energy constraints is...

[0039] The solution method for the projection position under bending and torsional constraints is as follows: the projection p that satisfies the bending and torsional potential energy constraints. i At this point, the rate of change of curvature Ω measures the degree of bending and twisting of the discrete element. n It should be 0, when measuring the system state q projected onto p. i When the distance is u, i.e. n u n+1 and When calculating the distance between locations, the quaternion is simply treated as a four-dimensional vector, and the Frobenius norm, the same norm used for calculating positional distances, is employed. The optimal solution can be obtained as follows: and Therefore, the solution to the curvature and torsional potential energy constraints is A. i =I8, B i =I8,

[0040] 2.6) Based on the constraints derived from the potential energy, the motion state of the current frame is obtained, i.e., the vector q representing the motion state, so that it achieves the minimum projection position in step 2.5, specifically: Where: q represents the current frame motion state, including the positions of all particles and the orientations of all discrete rope segments, i.e., q = [x1, ..., x2]. N u1, ..., u N-1 ], a diagonal matrix M, where the elements on the diagonal are the masses of all particles and the moments of inertia of all discrete rope segments. The remaining symbols are the same as in step 2.5.

[0041] Preferably, solving for the motion state q of the current frame is equivalent to solving a linear equation. When the particle mass, time step, number of constraints, and parameters are fixed, the above left-hand matrix... The result remains unchanged. Before the first calculation, the matrix on the left side of the above equation is decomposed using LU, and the decomposition result is stored. The calculation is performed using the upper and lower triangular matrices obtained after LU decomposition, which significantly improves speed.

[0042] The current frame motion state is equivalent to solving the linear equation system A. * x = b, where: sparse matrix When the linear equation system A is solved directly by using the elimination method * The time complexity of x = b is O(n^2). 3 ), where n is a matrix A * The dimensions. Since position and the four elements representing orientation are represented using three-dimensional and four-dimensional vectors respectively, A * The dimension is the sum of 3 times the number of particles and 4 times the number of discrete rope segments, which is a high-dimensional and sparse matrix. Directly solving the linear equation system will incur a considerable computational cost.

[0043] This embodiment employs an improved solution method, specifically including:

[0044] i) For storing sparse matrices, while traditional array storage methods can randomly retrieve elements in the matrix in constant time, they consume a lot of unnecessary memory and are not conducive to accelerating computation by leveraging sparsity. A better approach is to use the Compressed Sparse Column (CSC) format, which efficiently retrieves non-zero elements from a row and supports sparse matrix arithmetic operations. This format stores a sparse matrix by storing non-zero values, column ranges, and row indices.

[0045] ii) Although vector b changes in each iteration due to the projection result in step 2.5, as long as the selection matrix S of any constraint i remains constant... i constant matrix A i Weight w i And all particle masses, cable moments of inertia, and time steps remain unchanged, with a sparse matrix A. * It also remains unchanged. Therefore, the constraints generated in real time based on the motion state and input are separated out, while the potential energy constraints derived from the Cosserat model are retained. The weights and constant matrices of this part of the constraints remain unchanged; therefore, as long as the geometric parameters, material parameters, and simulation parameters of the cable remain unchanged, the sparse matrix A remains unchanged. * It will remain unchanged.

[0046] iii) Use the LU decomposition method to decompose the invariant sparse matrix A * The construction process proves that the matrix is ​​a positive definite symmetric matrix, therefore matrix A can be... *The unique product of a lower triangular matrix L and an upper triangular matrix U, i.e., PA * Q = LU, where P and Q are permutation matrices obtained using methods such as minimum degree sorting, which are beneficial for subsequent decomposition of the sparse matrix A. * If the filling of non-zero elements is minimized during the process, then the linear equation system A can be solved. * x = b is transformed into solving Ly = b and LUx = b separately. Due to the trigonometric properties of L and U, the elimination process can be completed in O(n log n) time. 2 Completed within )

[0047] iv) Cache the decomposition results, the lower triangular matrix L and the upper triangular matrix U, and reuse them in subsequent simulations. Although the decomposition itself requires O(n) time... 3 However, thanks to the method described above, the sparse matrix A... * Since it remains unchanged, it only needs to be decomposed once during cable initialization. The technical advantage is that the time complexity of solving for the vector q representing the motion state in each logical frame is reduced by an order of magnitude compared to methods like elimination that directly solve linear equations.

[0048] v) For constraints other than the potential energy constraint derived from the Cosserat model after decomposition, position-based dynamics (PBD) is used for preprocessing or post-processing, as in steps 2.4 and 2.7. After decomposing the constraints, the above sparse matrix A * The condition for invariance is that the mass of all particles, the moment of inertia of the cable, and the time step remain unchanged. Therefore, when the user changes the geometric parameters, material parameters, or simulation parameters of the simulation algorithm, the sparse matrix A... * This will be reconstructed and decomposed. Most of these constraints are related to scene interactions (such as collision constraints) and user interactions (such as animation constraints), and are unrelated to phenomena caused by internal stresses such as tension, shearing, bending, and torsion of the rope itself, thus having a minimal impact on the physical accuracy of the rope simulation.

[0049] 2.7) Use the Position-Based Dynamics (PBD) method to solve collision constraints to resolve self-collisions of the rope; collisions between the rope and other objects in the scene are usually detected and handled by the graphics engine after the triangular mesh is reconstructed, using, but not limited to, the method in Macklin's Extended Position Dynamics (XPBD).

[0050] 2.8) Using the motion state q obtained above, update the particle's velocity and the angular velocity of the discrete rope segment, including: the particle's velocity is set to the difference between the new and old positions, divided by the time step. The angular velocity of the discrete rope segment is set to the difference between the new and old orientations, divided by the time step.

[0051] 2.9) Output simulation results including the position and velocity of each particle, as well as the orientation and angular velocity of each discrete rope segment.

[0052] like Figure 5 As shown in the figure, this embodiment relates to a real-time interactive rope simulation system that implements the above-described method. The system includes: an interaction module, a simulation module, a rendering module, an input device, and a display device. The interaction module acquires user input in real time and changes the motion state of the rope or other objects in the virtual scene based on the input. The simulation module simulates the motion state of the rope in real time according to the rope simulation method. The rendering module recreates the geometric model of the rope in the virtual scene in real time based on the rope's motion state and parameters and outputs a rendered image. The input device receives user input. The display device outputs the rendered rope motion animation.

[0053] The interactive module acquires user input and modifies the geometric and material parameters of the rope based on the user input, or changes the rope's motion state, or creates animations for the rope's motion. Since the rope physics simulation method in this embodiment is real-time, the interactive module includes: a geometric parameter modification unit, a material parameter modification unit, a simulation parameter modification unit, a scene control unit, and an animation creation unit. Specifically, the geometric parameter modification unit re-initializes the rope based on the user-adjusted geometric parameters in the visual editor, such as the rope's centerline control points, cross-sectional radius, and sampling frequency. This includes generating Catmull-Rom spline curves based on control points, uniformly sampling arc coordinates according to the sampling frequency to obtain sampling points and discrete rope segments, and reconstructing the cross-section and connected triangular meshes based on the cross-sectional radius. The result is an updated geometric model of the rope, updating the discretized information of the rope stored in the computer. The material parameter modification unit updates the parameters in the simulation algorithm based on the user-adjusted material parameters in the visual editor, such as rope density, Young's modulus, and Poisson's ratio. This includes: updating the particle mass matrix, updating the discrete rope segment moment of inertia, and updating the weights of tensile and shear, bending and torsional potential energy constraints; the simulation parameter modification unit updates the operation of the rope simulation algorithm based on the simulation algorithm parameters adjusted by the user in the visual editor, such as time step and number of iterations. This includes: modifying the time between each updated motion state of the rope and the previous motion state, and solving the minimum weighted sum of the predicted motion state to the projection distance of each constraint by the number of iterations; the scene control unit updates the simulated scene in the computer based on user input information, such as key operations and mouse operations, including: moving characters in the scene, directly modifying the motion state of the rope in the scene, etc., to obtain the updated virtual scene; the animation creation unit creates corresponding animations and animation constraints in the simulation algorithm based on the user's keyframe time and keyframe state. The state of particles and discrete rope segments controlled by the animation in the current logical frame can be obtained by linear interpolation between keyframes, thus obtaining the rope motion result controlled by the animation in part or all of the time.

[0054] The simulation module includes a motion prediction unit, a projection position calculation unit, a projection distance optimization unit, and a motion state update unit. Specifically: the motion prediction unit processes the motion information from the previous frame, including the positions, velocities, and net external forces acting on all particles, as well as the orientation, angular velocities, and net external torques acting on all discrete rope segments, according to the equations of motion in classical mechanics, to obtain the motion result of the rope neglecting internal stress; the projection position calculation unit calculates the projections of the predicted rope motion state onto each constraint in parallel based on the predicted motion state; the projection distance optimization unit solves a system of linear equations based on the predicted motion state and the projection information onto all constraints to optimize the weighted sum of distances to all projections, obtaining a new predicted motion state result; and the motion state update unit updates the rope motion state based on the predicted motion state that satisfies the iteration count or the error tolerance, obtaining the motion result for this frame, and transmits it to the rendering module.

[0055] The rendering module calculates a circular surface with a preset rope radius at each particle's location, representing the rope's cross-section. The normal direction of the cross-section is the same as the tangent direction of the discrete rope segment where the particle is located. A certain number of points are sampled on each particle's cross-section, and connecting these points forms a triangular mesh, representing the rope's surface. Using a preset rope material, this triangular mesh is rendered, thus recreating the moving rope in the virtual scene. This rendering module is integrated into the Unity graphics engine.

[0056] Through specific practical experiments, under the hardware environment settings in Table 1 and the software environment settings in Table 2, the above-mentioned rope simulation algorithm was launched with the parameters in Table 3. Before conducting specific experiments and measuring the physical accuracy of the rope simulation, an accurate solution as a reference is needed. Since most rope motion scenarios are difficult to calculate the accurate value of the integral and therefore do not have a mathematically meaningful analytical solution, the same scenario was recreated in the tried-and-tested commercial finite element simulation software ABAQUS, and the simulation was performed using extremely small time steps to obtain a highly accurate numerical solution as a reference.

[0057] hardware devices Configuration CPU Amd5800 Memory 16 GB harddisk Western Digital SN550 graphics card NVIDIA 3070

[0058] Software Name Version Unity3D 2020.2.3 Abaqus 2021

[0059] Image ID Figure 7 Figure 8 ， Figure 9 R(mm) 10 6 <![CDATA[ρ(g / m 3 )]]> 1.3 1.3 E(MPa) * 10 v 0.25 0.25 L(m) 2 2 h(ms) 5 5 <![CDATA[g y (m / s 2 )]]> 9.861 9.861

[0060] In Table 3, R represents the radius of the cross-sectional area of ​​the rope (mm), and ρ represents the density of the rope (g / m³). 3 E is Young's modulus (MPa), v is Poisson's ratio, L is the total length of the rope (m), h is the time step (ms), and g is the total length of the rope. y The component of gravitational acceleration in the Y-axis direction (m / s²) 2The asterisk (*) represents parameters that vary in the experiment. In Abaqus, which serves as the reference solution, only the highest possible accuracy is considered, not performance; therefore, some parameters are inconsistent, such as extremely small time steps: h. Abaqus =1ms.

[0061] like Figure 6 As shown, this simulates a scenario where the two ends of different ropes are fixed and allowed to hang naturally under gravity. The ropes are initially placed horizontally with their natural length (without any elastic deformation). The comparison method is the PBD method, widely used in real-time graphics applications, implemented using ObiRope, the most popular rope simulation plugin in Unity. The Y-axis of the image corresponds to the position of each discrete unit of the rope in the Y-direction. Since, except for highly elastic ropes, the deformation in the Y-direction is smaller than the original length of the rope under gravity alone, the X-axis of the image uses the particle index x∈[1,N], rather than the particle position component on the X-axis, where N is the number of discrete units in Abaqus, which is a fixed value of 600 in this experiment. Since the output of Abaqus is relative deformation, it is added to the initial height of the rope. Figure 6 The experiment simulates ropes made of different materials, with each row corresponding to one material. (a) is a rubber rope with an E = 1 MPa, (b) is a rubber rope with an E = 10 MPa, and (c) is a hemp rope with an E = 100 MPa. It can be seen that the weights in the PBD are derived from empirical parameter tuning and do not possess physical meaning. In this experiment, the weight parameters of the PBD were carefully adjusted to better simulate a rubber rope with a Young's modulus of 10 MPa. However, it can be seen that modifications to the physical parameters are not reflected in the PBD parameters, therefore it cannot simulate ropes made of other materials. The algorithm of [other algorithm], on the other hand, simulates the performance of ropes made of different materials better, with an average error of less than 3%.

[0062] like Figure 7 As shown, to verify that this method can correctly drive and handle collisions through animation, and achieve relatively accurate results with fewer iterations, the experiment was designed with a fixed starting point of the rope and the following animation for the ending point: two keyframes. In the first keyframe, the endpoint is at its natural length from the starting point of the rope; in the second keyframe, the endpoint coincides with the starting point of the rope. The interval between keyframes is 3 seconds. In Abaqus, which serves as the reference solution, this interval is set to 30 seconds to obtain higher accuracy. The horizontal and vertical axes in the figure represent the X and Y components of the rope's position in world space, respectively, in meters. Similar to the previous experiment, since the output of Abaqus is relative, the initial position of the rope is added to it through post-processing. Figure 7The effects of different iteration numbers are illustrated in (a) 4 iterations, (b) 8 iterations, and (c) 16 iterations. It can be seen that the algorithm converges with fewer iterations, and the constraint strength is determined by the physical parameters, independent of the iteration number. However, in the PBD method, the stiffness of each constraint increases with the number of iterations, but this does not necessarily improve physical accuracy. In this experiment, the rope exhibited excessive elasticity with fewer iterations, while exhibiting excessive rigidity with higher iterations. Therefore, it has the drawback of requiring empirical parameter tuning to obtain a relatively accurate visual representation. Furthermore, parameter tuning in PBD becomes more difficult for scenarios where multiple physical parameters simultaneously determine the physical performance. For example, in this experimental scenario, Young's modulus and Poisson's ratio jointly determine the rope's flexibility and lowest descent position, which the PBD method struggles to accurately simulate simultaneously.

[0063] like Figure 8 As shown, in Figure 7 The scene demonstrates how collision constraints and animation-driven constraints can correctly control the movement of the cable and handle collisions caused by that movement, especially when the cable end point moves and collides with itself. The collision constraints generated by detecting particles and the implicit surface determined by the centerline of the Cosserat theoretical model also handle self-collisions better than PBD. Figure 9 In this paper, some scenarios provided by Obi Rope are added, and the algorithm used simulates rope breakage, combined force, and torsion phenomena. It is proven that the algorithm can adequately reproduce the physical behavior of ropes.

[0064] Table 4 summarizes the performance data from the previous two experiments. The PD algorithm framework (Standard-PD) with only rope orientation and angular momentum included, the PD algorithm with matrix decomposition and caching (Cache-PD) used in this invention, and the PBD algorithm (implemented by the Obi Rope team) were tested under the same scenario. Only the performance of updating the rope was statistically analyzed, the degrees of freedom of the system state and the number of iterations were standardized, and other parameters were the same as in Table 2. A real-time requirement of 30 frames per second was set, meaning each computation must take less than 33ms. This demonstrates that the algorithm has achieved real-time performance, with a significant improvement over the previous version.

[0065] Scene N Standard-PD This invention PBD Figure 6 280 45.8ms 17.46ms 0.92ms Figure 7 280 61.2ms 18.70ms 1.30ms

[0066] Table 5 shows the performance data of cable simulation methods with similar physical accuracy to the previous method in recent years. Among them: SOLER's method is also based on projection mechanics. Although it is faster in calculation, it does not add animation and collision constraints, which reduces interactivity and realism; WEN's method automatically downsamples the cable in real time, sacrificing some accuracy for real-time performance. Therefore, the actual accuracy of WEN's method should be reduced; GAZZOLA uses classical mechanics to derive the equations of motion after discretization and calculates using the explicit Euler method. However, it cannot guarantee the stability of the algorithm when the cable height is deformed.

[0067] Table 5

[0068]

[0069]

[0070] Compared with existing technologies, this method achieves an excellent balance between physical accuracy and computational efficiency in cable simulation. The physical accuracy is similar to that of non-real-time finite element methods, and it represents a significant improvement over the PBD method used in current real-time graphics applications. The computational efficiency is also improved compared to other cable simulation algorithms in recent years. This invention adds angular momentum to the original projection mechanics algorithm and decomposes the original projection mechanics constraints into fixed cable constraints and dynamically generated collision and animation constraints. By using matrix decomposition and caching to accelerate the solution of linear equations, it can maintain the original accuracy while achieving real-time performance, thus overcoming the problem of insufficient physical accuracy in simulating cables in current ordinary real-time graphics applications and enhancing the performance of cable simulation. Using the tensile, shear, bending, and torsional potential energy of the cable as constraints, the projection force algorithm can be extended to simulate various cable phenomena, achieving physical accuracy similar to that of commercial finite element software.

[0071] The above-described specific implementations can be partially adjusted by those skilled in the art in different ways without departing from the principles and purpose of the present invention. The scope of protection of the present invention is defined by the claims and is not limited to the above-described specific implementations. All implementation schemes within the scope of the claims are bound by the present invention.

Claims

1. A real-time interactive high-precision rope and cable simulation method, characterized in that, It includes an initialization phase and a simulation phase. In the initialization phase, the continuous rope is discretized and stored according to the initial state of the rope to be simulated, and then the simulation parameters are set. In the simulation phase, an improved projection mechanics algorithm is used. Based on the motion constraints derived from the Cosserat rope model, matrix decomposition and caching are used to accelerate the simulation. The motion state of each logical frame of the rope is output in real time, thereby enhancing the physical accuracy of the rope simulation. The specific steps of the simulation phase are as follows: Step 1) Divide the entire process of simulating cable movement into several logical frames to be calculated according to the settings; Step 2) Within each logical frame, input the motion state of the cable in the previous frame, and use the improved projection mechanics algorithm to calculate the motion state of the cable in that frame, specifically including: 2.1) Predicting the position and orientation of the rope: Based on the particle's velocity and position in the previous frame and the net external force acting on it in this frame, predict the particle's position in this frame, specifically: ,in: For the predicted positions of all particles in this frame, This represents the positions of all particles in the previous frame. For the set time step, The velocity of all particles in the previous frame. It is a diagonal matrix composed of the masses of all particles. The net external force acting on the particle; 2.2) Based on the angular velocity and orientation of the discrete rope segment in the previous frame and the net external torque acting on it in this frame, predict the orientation of the discrete rope segment in this frame, specifically: ,in: For the orientation of this frame of discrete rope segments, For the orientation of this frame of discrete rope segments, The angular velocity of the discrete rope segment is obtained from the motion state of the previous frame. The resultant external torque is obtained from user input, collision handling, or animation design; diagonal matrix The geometric moment of inertia represents the rope segment along the center line. Since the quaternion is treated as a four-dimensional vector in the linear solution process, and based on the definition of a quaternion, we know... ,in This is the initial length of the rope segment; since the entire rope is uniform, the same density is used throughout. Let be the geometric moment of inertia of the cross-section. When the cross-section is circular, it can be obtained from... , work out; 2.3) Initialize Constraints: Constraints include potential energy constraints derived from the Cosserat model, as well as collision constraints and animation-driven constraints added to enhance the realism and interactivity of the physics. The potential energy constraints remain unchanged after the cable parameters are determined and only need to be initialized when the cable parameters are modified. Collision constraints and animation-driven constraints change with the movement of the cable and user interaction, and need to be re-initialized in each frame. Specifically, this includes: i) Initializing the constraints derived from the stretching and shear potential energy of the Cosserat model... ii) Initialize the constraints derived from the bending and torsional potential energies of the Cosserat model as follows: iii) Initialize the collision constraints as and iv) initialize animation constraints, where: For the current particle The difference in tensile and shear potential energy between the current state and the target state. This means the constraint is to minimize the potential energy difference; To measure the degree of tension and shear, strain is equal to the direction vector of the particle on the tangent to the center line, minus the orientation direction vector of the discrete rope segments originating from that particle. The strength of the collision constraint; It is an indicator function, which is 1 when the particle's collider is still in contact with other colliders, and 0 when they are no longer in contact; R is the radius of the cable's cross-section. Let j be the position of particle j in three-dimensional space. The position of the particle projected onto the surface of other colliders; Strength of the animation constraint The position set for the particle within this animation frame; 2.4) Solve animation constraints using the position-based dynamics (PBD) method: determine the animation keyframe time, keyframe state information and interpolation method, and calculate the state information of each logical frame at the time point; The state information includes: the position of the particle controlled by the animation and the orientation of the discrete rope segment. The solution method is similar to Macklin's solution of fixed point constraints in the Extended Position Dynamics Method (XPBD), except that the target fixed point position is set as the position of the particle controlled by the animation in the state information. 2.5) For all constraints derived from potential energy, solve for the projected position of the constraint, i.e., satisfying constraint i and relating it to the current motion state. The point where the generalized distance is minimized is specifically: finding the position where the current motion state q is minimized under constraint i. generalized distance Get the minimum value That is, the projection position of the constraint, where A, B, and S are constant matrices constructed based on the constraint. To constrain the corresponding non-negative weights, The indicator function is defined if and only if the projection position is... The value is 0 when the constraint is satisfied, and infinite when the constraint is not satisfied. 2.6) Based on the constraints derived from the potential energy, the motion state of the current frame is obtained, i.e., the vector representing the motion state. To achieve the minimum projection position in step 2.5, specifically: Where: q is the desired motion state of the current frame, including the positions of all particles and the orientations of all discrete rope segments, i.e. The diagonal matrix M contains the elements on the diagonal as the mass of all particles and the moment of inertia of all discrete rope segments; the remaining symbols are the same as in step 2.

5. 2.7) Use the Position-Based Dynamics (PBD) method to solve collision constraints to resolve self-collisions of the rope; collisions between the rope and other objects in the scene are usually detected and handled by the graphics engine after the triangular mesh is reconstructed, using, but not limited to, the method in Macklin's Extended Position Dynamics (XPBD) method. 2.8) Using the motion state obtained above The velocity of the particle and the angular velocity of the discrete rope segment are updated, including: the velocity of the particle is set as the difference between the new position and the old position, divided by the time step; the angular velocity of the discrete rope segment is set as the difference between the new orientation and the old orientation, divided by the time step. 2.9) Output simulation results including the position and velocity of each particle, as well as the orientation and angular velocity of each discrete rope segment.

2. The real-time interactive high-precision cable simulation method according to claim 1, characterized in that, The initialization phase specifically includes: The first step is to obtain the control points used by the user to define the center line of the rope: the user inputs the rope radius and control points; after obtaining the above control points in three-dimensional space, the Catmull-Rom spline curve with C1 continuity that passes through all control points is selected as the center line of the rope. The second step is to obtain the sampling frequency N of the user input: According to the definition of the Cosserat theoretical model, let the centerline of the rope of length L be... ;Sampling N points uniformly along the center line based on arc coordinates s: The rope is then divided into N-1 segments, each consisting of two vertices. And a four-element representing rotation Definition; Adoption After rotating the world coordinate system, we obtain the local coordinate system in the Cosserat rope model. Then, the N vertices are particles, and the rope is divided into N-1 segments, which are discrete rope segments, serving as data structures stored in the computer. The third step involves the user inputting the geometric and material parameters of the rope: based on the cross-sectional radius input by the user, a circular surface with the corresponding radius is drawn at each particle, and the normal direction of the circular surface is the tangential direction of the center line at that point; the points on these circular surfaces are connected in sequence to form a triangular mesh, thereby restoring the geometric model of the rope; The fourth step is to accept the simulation parameters input by the user, including: the total duration of the simulated motion, the time step of each change in the motion state of the simulated cable, and the number of iterations for each state change. The fifth step is to store the geometric parameters, material parameters, and information about particles and discrete rope segments of the cable.

3. The real-time interactive high-precision cable simulation method according to claim 2, characterized in that, The particle has mass, position, and velocity attributes, with the initial position set at the sampling point and the initial velocity being 0. The discrete rope segment has orientation, angular velocity, geometric moment of inertia, and initial length attributes, with the initial orientation set as the tangent direction of the rope centerline at the starting particle of the discrete rope segment, the initial angular velocity being 0, and the initial length being the distance between the two corresponding vertices of the discrete rope segment. The geometric parameters include the initial length and the cross-sectional radius; The material parameters mentioned include density, Young's modulus, and Poisson's ratio; The geometric moment of inertia of the discrete rope segment is obtained by multiplying the rope density input by the user by the average length of the discrete rope segment with the particle as the starting and ending point, to obtain the mass attribute of each particle, and then the geometric moment of inertia of the discrete rope segment can be calculated based on the cross-sectional radius.

4. The real-time interactive high-precision rope simulation method according to claim 1, characterized in that, The direction vector on the tangent of the center line is approximately equal to the position difference between adjacent particles within the discrete unit divided by the initial length. ,therefore ,in: This represents the imaginary part of the four elements q; For the current discrete rope segment The difference in bending and torsional potential energy between the current state and the target state. This means the constraint is to minimize the potential energy difference; The strain used to measure the degree of tension and shear is equal to the rate of change of curvature between the discrete rope segment and its adjacent discrete rope segments, i.e. .

5. The real-time interactive high-precision cable simulation method according to claim 1, characterized in that, The solution method for the projection position under tensile and shear constraints is specifically as follows: the projection position that satisfies the constraints of tensile and shear potential energy. At this point, the strain that measures the degree of tension and shear in a discrete element should be 0, i.e. , These represent the positions of the two particles in the discrete element during this iteration, and the normal to the cross-section of the discrete element. Since the particle's position and the rope segment's orientation are independent, the problem can be further divided into two optimization problems: optimizing the particle's position. And optimization of the local orientation of the rope segment: The solution for particle position optimization is in When the positional difference of the particles is parallel to the normal of the cross-section, and the distance between the particles is the same as in the initial state, the tensile deformation of the discrete unit is 0; the solution for optimizing the local orientation of the rope segment is obtained in When the shear deformation of the discrete element is 0, it is obtained that... As a representative from arrive The quaternion of rotation; from this, the solution for the tension and shear potential constraints is obtained as follows: ; The solution method for the projection position under bending and torsional constraints is specifically as follows: the projection position that satisfies the bending and torsional potential energy constraints. At this point, the rate of change of curvature measures the degree of bending and twisting of the discrete element. It should be 0, in measuring the system state. To projection When the distance is, that is and When calculating the distance between them, the quaternion is simply treated as a four-dimensional vector, and the Frobenius norm, which is the same as the one used to calculate the distance between positions, is employed. The optimal solution can be obtained as follows: as well as Therefore, the solutions for the curvature and torsional potential energy constraints are obtained as follows: .

6. The real-time interactive high-precision rope simulation method according to claim 1, characterized in that, Solving for the motion state q in the current frame is equivalent to solving a linear equation. When the particle mass, time step, number of constraints, and parameters are fixed, the matrix on the left side of the above equation... The result remains unchanged; before the first calculation, the matrix on the left side of the above equation is decomposed using LU, and the decomposition result is stored. Using the upper and lower triangular matrices after LU decomposition for calculation can significantly improve speed; The current frame motion state is equivalent to solving a system of linear equations. , where: sparse matrix When the linear equation system is solved directly by using the elimination method The time complexity is n is a matrix The dimensions; since position and the four elements representing orientation are represented using three-dimensional and four-dimensional vectors respectively, therefore The dimension is the sum of 3 times the number of particles and 4 times the number of discrete rope segments, which is a high-dimensional and sparse matrix. Directly solving the linear equation system will incur a considerable computational cost.

7. The real-time interactive high-precision rope simulation method according to claim 6, characterized in that, The method for solving the motion state of the current frame specifically includes: i) Use a compressed sparse column format to store sparse matrices, storing non-zero values, column ranges, and row indices to store a sparse matrix; ii) Decompose the constraints generated in real time based on the motion state and input, while retaining the potential energy constraints derived from the Cosserat model; iii) Use LU decomposition to decompose the invariant sparse matrix. : Matrix The unique product of a lower triangular matrix L and an upper triangular matrix U, i.e. Where P and Q are permutation matrices obtained using the minimum degree sorting method, the system of linear equations will be solved. Transform into solving separately and ; iv) Cache the decomposition results, the lower triangular matrix L and the upper triangular matrix U, and reuse them in subsequent simulations; v) For constraints other than the potential energy constraint derived from the Cosserat model after the splitting in step ii), a position-based dynamics method is used for preprocessing or postprocessing.

8. A real-time interactive cable simulation system for implementing the method of any one of claims 1 to 7, characterized in that, include: The system includes an interaction module, a simulation module, a rendering module, an input device, and a display device. The interaction module acquires user input in real time and changes the motion state of the rope or other objects in the virtual scene based on the input. The simulation module simulates the motion state of the rope in real time according to the rope simulation method. The rendering module recreates the geometric model of the cable in the virtual scene in real time based on the cable's motion state and parameters, and outputs the rendered image; the input device is used to receive user input; and the display device outputs the rendered cable motion animation.

9. The real-time interactive cable simulation system according to claim 8, characterized in that, The interactive module includes: a geometric parameter modification unit, a material parameter modification unit, a simulation parameter modification unit, a scene control unit, and an animation creation unit. Specifically: the geometric parameter modification unit reinitializes and updates the geometric model of the rope based on the user-adjusted rope geometric parameters in the visual editor, and updates the discretized rope information stored in the computer; the material parameter modification unit updates the parameters in the simulation algorithm based on the user-adjusted rope material parameters in the visual editor; the simulation parameter modification unit updates the operation of the rope simulation algorithm based on the user-adjusted simulation algorithm parameters in the visual editor; the scene control unit updates the simulated scene in the computer based on user input information; and the animation creation unit creates corresponding animations and animation constraints in the simulation algorithm based on the user's keyframe times and keyframe states. Linear interpolation between keyframes yields the state of particles and discrete rope segments controlled by the animation within the current logical frame, resulting in the rope motion outcome controlled partially or entirely by the animation.

Images