Soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element

Abstract

The invention provides a soft tissue deformation simulation method based on a dual-mode hybrid dynamics finite element, and the method comprises the steps: dividing 3D grid data into a linear region and a non-linear region, carrying out the data pre-calculation of units in the linear region and the non-linear region, and obtaining the pre-calculation information of the units; the pre-calculation information of the unit is utilized, according to the area to which the unit belongs, a nonlinear finite element model is adopted to solve the stress and strain on the unit in the nonlinear area and the acting force of the unit on the node of the unit, and a linear finite element model is adopted to solve the stress and strain on the unit in the linear area and the acting force of the unit on the node of the unit; calculating the resultant force of the unit acting force borne by each node according to the acting force of all the units on the node; and the displacement of each node of the current moment unit is solved by using explicit time integration, the position of each node at the next moment is solved, and deformation information at the next moment is obtained. According to the method, the simulation accuracy is remarkably improved while the calculation efficiency is kept.

Description

Technical Field

[0001] The present invention relates to the technical field of surgical simulation, and in particular to a soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element. Background Art

[0002] With the in-depth application and penetration of virtual reality technology in the medical field, surgical simulation has become an important tool for surgeon training and an important means to improve surgeon skills and reduce surgical risks. The surgical simulation system can simulate the real surgical environment, allowing doctors to practice without actual patients, thereby improving the accuracy and safety of surgical operations. However, real-time and accurate simulation of soft tissue deformation has always been a major challenge in the field of surgical simulation.

[0003] Finite element method has become the mainstream method for soft tissue deformation simulation in surgical simulation. Due to its solid mathematical foundation, although the finite element method can obtain relatively reliable results and achieve certain effects in surgical simulation, the high computational cost of this method and the real-time requirements of the surgical simulation system have brought considerable challenges to the real-time deformation simulation of soft tissue. At present, the TLED (Total Lagrangian Explicit Dynamics) algorithm proposed by Miller et al. has done a relatively outstanding job in the field of deformation simulation. It uses the Lagrangian dynamics framework and explicit time integration to simulate the nonlinear deformation simulation of soft tissue in real time on a larger grid scale. Nonlinear finite elements are usually more accurate than linear finite elements when dealing with large deformations, and correspondingly, more computational effort will be required. When simulating a surgical scene, the requirements for simulation accuracy in each area of ​​the scene are different. Usually, accurate simulation can be performed on the object undergoing the surgical operation to ensure the accuracy of the surgical simulation. Reducing the simulation accuracy of the surrounding objects can reduce the computational burden. During the surgical process, surgical operations are often concentrated in a certain area of ​​the tissue or organ, where surgical operations such as cutting and detection will occur. This area requires a more complex model to enhance the authenticity of the surgical simulation. Tissues and organs far away from the surgical area do not deform much during surgery, so a simpler model can be used to improve calculation efficiency.

[0004] The invention patent with application number 201810566190.1 discloses a soft tissue deformation method based on the finite element model of the octree grid, including: drawing a three-dimensional model of the soft tissue, constructing multiple uniform hexahedral grids based on the AABB method, generating an octree grid based on the hexahedral grid based on the octree grid generation algorithm, modeling the hexahedral grid by the finite element method, solving the deformation process of the soft tissue, assembling the unit stiffness matrices on adjacent unit nodes into the total stiffness matrix of the discrete domain, numerically solving the dynamic equilibrium differential equation of each matrix by time integration under the condition of dynamic equilibrium, obtaining the node displacement that changes with time, and displaying the soft tissue force deformation scene by rendering the node displacement. The above invention realistically simulates the process of tensile deformation of the epidermis of soft tissue of arbitrary shape in virtual surgery, has high real-time performance, reduces the amount of calculation, and solves the problem that the number of finite element grids in the past is complex and cannot simulate the soft tissue deformation process in real time. In the past, the finite element method had the problems of complex number of grids, large amount of calculation, and difficulty in real-time simulation of the soft tissue deformation process when simulating soft tissue deformation. Especially when dealing with the situation where the epidermis of soft tissues of arbitrary shapes is deformed by tension, traditional methods often have imprecise grid division or low computational efficiency, resulting in simulation results that are not realistic enough and lack real-time performance. The soft tissue deformation method of the finite element model based on the octree grid proposed in the above invention draws a three-dimensional model of the soft tissue, and uses the AABB method to construct a uniform hexahedral grid as a basis, and then combines the octree grid generation algorithm to effectively optimize the grid structure, reduce the number of unnecessary grids, and maintain the accuracy of the simulation. In addition, during the finite element modeling and solution process, the above method assembles the unit stiffness matrices on adjacent unit nodes to form a total stiffness matrix, and uses time integration to solve the dynamic equilibrium differential equation, thereby realizing real-time simulation of the soft tissue deformation process. Therefore, it can solve the problems of complex grids, large amount of calculations, and difficulty in real-time simulation in the previous finite element method when simulating soft tissue deformation, and provides a more realistic and efficient soft tissue deformation simulation method. Summary of the invention

[0005] In order to solve the technical problems of large computational complexity and poor real-time performance of soft tissue deformation simulation methods, the present invention proposes a soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element, which combines linear and nonlinear finite element methods and directly calculates the force of the unit on each node connected to it through unit force analysis, thereby achieving accurate simulation of soft tissue deformation and improving the efficiency and accuracy of surgical simulation.

[0006] In order to achieve the above object, the technical solution of the present invention is implemented as follows: a soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element, the steps of which are as follows:

[0007] Step 1: dividing the tissue organ region into a linear region and a nonlinear region according to the 3D grid data, and performing data pre-calculation on the units in the linear region and the nonlinear region to obtain pre-calculation information of the units;

[0008] Step 2: Using the pre-calculated information of the unit, according to the area to which the unit belongs, a nonlinear finite element model is used to solve the stress, strain and the force of the unit on its nodes in the nonlinear area, and a linear finite element model is used to solve the stress, strain and the force of the unit on its nodes in the linear area;

[0009] Step 3: Calculate the resultant force of the unit forces on each node based on the forces of all units on their nodes;

[0010] Step 4: Use explicit time integration to solve the displacement of each node of the unit at the current moment, obtain the position of each node at the next moment, and obtain the deformation information at the next moment.

[0011] Preferably, the nonlinear area is a concentrated area for surgical operations, and the linear area is a part far away from the surgical operation area; the rectangular space is divided into a set of cubes, each cube is further divided into five tetrahedral units, and one tetrahedral unit is associated with four nodes; the units in the linear area are linear units, and the units in the nonlinear area are nonlinear units; the linear area and the nonlinear area in the surgical scene are connected by the intermediate nodes connecting the linear units and the nonlinear units to form a continuous surgical simulation scene; the force of the unit on each node connected to it is directly calculated through unit force analysis, and the resultant force of all the connected units on each node is solved, then the intermediate node is subjected to the force from the linear unit and the nonlinear unit, and the two areas can be connected through the intermediate node.

[0012] Preferably, the data obtained by precalculation includes the coefficient matrix of the initial coordinates 0 D X and the initial deformation submatrix of any node for calculating the strain-displacement matrix;

[0013] The deformation calculation of tetrahedral elements includes mesh unit information and the connection relationship between the unit and the node. The elements in the mesh are represented in the form of a set. The mesh contains three lists: tetrahedral list, surface polygon list and node list. The tetrahedral list stores the node number at the corresponding node associated with the tetrahedral element, and all information related to the node can be found through the number. The node number stored in the tetrahedral element is represented by an arrow, and the arrow points to the unit corresponding to the number; the surface polygon list maintains the connection information with the corresponding node;

[0014] Each element also includes information related to deformation calculation, including pre-calculated information in each unit and the initial volume of the unit; the node list stores the node mass, the initial position of the node, the resultant force of the unit and external forces acting on the point, the displacement of the node at the previous moment, and the node displacement information.

[0015] In the data pre-calculation, for each unit, the spatial derivative of its shape function and the linear strain-displacement matrix are calculated, the initial volume of the unit is calculated, the mass and damping of each node are calculated, and the topological information between the nodes, edges, units and surface polygons of the mesh is generated.

[0016] Preferably, the Cauchy strain tensor and the Cauchy stress tensor are used to describe the strain and stress of the linear unit respectively; the positions of the four nodes of the tetrahedral unit at the initial moment are The positions after deformation are The method for calculating the Cauchy strain tensor is: solve the displacement of each node in the tetrahedral unit, and then obtain the displacement gradient; after calculating the Cauchy strain tensor, obtain the Cauchy stress tensor σ;

[0017] In nonlinear units, the strain tensor uses Green's strain tensor, and the stress tensor often uses Second Piola-Kirchhoff stress tensor.

[0018] Preferably, the force exerted by the nonlinear unit on the node is

[0019]

[0020] in, is the vector form of the Second Piola-Kirchhoff stress tensor, represents the three normal stress components, represents the three shear stress components; is the strain-displacement matrix, 0 V is the initial volume of the unit;

[0021] In the linear unit, the resultant force f = (σ·n)A of the stress on each surface of the tetrahedral unit is evenly distributed to each node of the surface triangle, and the force of the tetrahedral unit on all nodes is obtained. Where σ is the Cauchy stress tensor, n is the normal vector of the tetrahedron surface triangle, and A is the area of ​​the tetrahedron surface triangle.

[0022] Preferably, the stress tensor The components of are obtained by differentiating the Green strain tensor through the Neo-hookean strain energy function:

[0023]

[0024] Where μ is the shear modulus, δ ij is a second-order tensor, λ is the first parameter of the Lame constant, t J is the Jacobian determinant, is the Right Cauchy-Green deformation tensor, and

[0025] The strain-displacement matrix Written as a collection of submatrices:

[0026] in, is the strain-displacement matrix The i-th submatrix of , and: in, is the deformation gradient of the nonlinear unit, the initial deformation submatrix of any node The calculation formula is:

[0027]

[0028] Solve for the initial deformation submatrix of any node Need to solve the derivative of the shape function with respect to the space coordinates ( 0 x, 0 y, 0 z) represents the initial spatial coordinates of any node; the spatial derivative of the shape function is solved by the isoparametric finite element formula,

[0029] According to the partial differential law, the shape function h of the node i The partial derivative relationship of and isoparameters is:

[0030]

[0031] in, 0 x1, 0 x2, 0 x3 all represent equal parameters;

[0032] The relationship between the partial derivatives of the shape function with respect to the isoparameters and the partial derivatives of the shape function with respect to the initial coordinates and the partial derivatives of the initial coordinates with respect to the isoparameters is solved by the chain rule. After a series of identity transformations, we get:

[0033]

[0034] Find the shape function h i Derivatives with respect to space:

[0035]

[0036] in, is the derivative of the shape function with respect to the isoparameter in the isoparametric element;

[0037] The method for calculating the Cauchy strain tensor is as follows: Assume that the initial coordinates of any node inside the tetrahedron unit are 0 x i , the position of the node after the tetrahedral element is deformed is t x i , the initial coordinates of any node are:

[0038]

[0039] Where b = (b1 b2 b3) T is the weight vector, the coefficient matrix 0 D X for

[0040]

[0041] Similarly, the position after deformation t x i for:

[0042]

[0043] Where b = (b1 b2 b3) T is the weight vector, the coefficient matrix t D X for:

[0044]

[0045] The final position after deformation is:

[0046]

[0047] in, is the mapping matrix of the position of any node inside the tetrahedron unit before and after deformation, and the displacement of any node 0 x k represents the coordinate component of the node in the k direction at the initial moment, It represents the deformation of the node in the k direction at time t;

[0048] Get the displacement gradient The calculation formula is:

[0049]

[0050] Where I is the identity matrix and the Cauchy strain tensor is expressed as

[0051]

[0052] The Cauchy stress tensor σ is: Among them, df is the force acting on any tiny plane in space, dA is the area of ​​the tiny plane, and n is the normal vector of the tiny plane.

[0053] Preferably, the method for calculating the resultant force of the unit forces acting on each node is as follows: for each node in the grid, there are a number of units connected to it, and the force acting on the node by the unit e connected to the node is t F e , then the resultant force F acting on the node v is v for If there is an external force acting on the node v, a force R v , then the resultant force on node v is R v -F v ;

[0054] Each node uses Newton's law to analyze the force. For a static system, the node is at rest and the equilibrium relationship is:

[0055]

[0056] For a dynamic system without damping effect, the dynamic equilibrium equation of the node is: Among them, m v is the mass of the node, u v is the displacement of the node, is the node acceleration, F v is the force exerted by the unit on the node, R v is the external force on the node;

[0057] For nodes with damping effects, the dynamic equilibrium equation is: Among them, d v is the damping factor of the node, is the node speed.

[0058] Preferably, in a dynamic system, the finite element mass matrix is ​​converted into a diagonal matrix by a direct mass concentration method. Assuming that the damping matrix D is proportional to the mass matrix D=αM, α is a constant, the damping matrix D is expressed as a diagonal matrix, and the stiffness term K(U)U of the dynamic equilibrium equation is:

[0059] Among them, e represents the unit in the set, Represents the force of the unit on the node, which is the force F e A vector; U represents a displacement vector;

[0060] Due to the deformation of the object, the resultant force of the internal stress on node i is F i , the mass of node i is M ii , damping is D ii , the current displacement is U i , the external force is R i , then the dynamic equilibrium equation of node i is: represents the velocity of node i, represents the acceleration of node i.

[0061] Preferably, the hybrid finite element model composed of the nonlinear finite element model and the linear finite element model is solved by the difference method to solve the dynamic equilibrium equation;

[0062] The central difference method uses difference instead of differentiation and solves the displacement U through iteration. i Function value at time t+Δt t+Δt U i , solve for the displacement U at time t i The first derivative of The method is:

[0063] For the unknown displacement U i The first and second derivatives at time t are

[0064]

[0065] Substitute into the dynamic equilibrium equation The equilibrium equations at the nodes are obtained:

[0066]

[0067] in, t F i and t R i They represent the unit force and external force on node i respectively;

[0068] Through the identity transformation we get:

[0069]

[0070] The displacement of each node at the next moment can be solved, and then the deformation information of the object in the dynamic system can be obtained.

[0071] Preferably, the step length Δt is smaller than the critical step length Δt cr ; Critical step length Δt in linear elastic analysis cr for: Among them, L e is the minimum characteristic unit length in the unit set, and c is the longitudinal wave velocity of the material.

[0072] Compared with the prior art, the present invention has the following beneficial effects: the present invention proposes an innovative dual-mode hybrid dynamic finite element algorithm, which combines the high efficiency of the linear finite element method and the accuracy of the nonlinear finite element method, and realizes fast and accurate simulation of soft tissue deformation by introducing a dual-mode hybrid strategy.

[0073] The present invention dynamically selects a linear or nonlinear finite element model for calculation according to the deformation degree of soft tissue, adopts a more accurate and higher nonlinear finite element model for the surgical operation area to improve the accuracy of surgical simulation, and adopts a linear finite element model with higher computational efficiency in the non-surgical area to obtain faster computational efficiency, uses explicit integration to simulate the dynamic system, integrates the hybrid model in the unit force analysis stage, and completes the real-time simulation of the hybrid finite element model, thereby significantly improving the accuracy of the simulation while maintaining the computational efficiency. Compared with the method of purely using a nonlinear finite element model, the present invention improves the computational efficiency of some units and improves the real-time performance of deformation simulation. Through experimental verification, the present invention can provide a more efficient and accurate soft tissue deformation simulation method for the surgical simulation field, provide a more realistic training environment for surgeons, thereby improving surgical skills and reducing surgical risks, and can become an important research direction in the field of surgical simulation in the future, providing new ideas and methods for scholars and researchers in related fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0075] Figure 1 It is a flow chart of the present invention.

[0076] Figure 2 Schematic diagram of two ways of dividing a rectangular unit into tetrahedral units.

[0077] Figure 3 An example of a generated tetrahedral mesh is shown below.

[0078] Figure 4 Schematic diagram of the data structure of the present invention.

[0079] Figure 5 Schematic diagram of a tetrahedral unit.

[0080] Figure 6 Schematic diagram of direct matrix concentration of rod elements.

[0081] Figure 7 Schematic diagram of direct mass concentration of tetrahedral elements.

[0082] Figure 8 Schematic diagram of solving the first-order derivative of displacement using the central difference method of the present invention.

[0083] Fig. 9 The diagram is an operating effect diagram of the dual-mode hybrid dynamic finite element deformation of the present invention, wherein all units in (a) are linear units, all units in (b) are nonlinear units, and the upper half in (c) is a linear unit and the lower half is a nonlinear unit.

[0084] Fig.10 This is the error analysis simulation scene diagram of the present invention. DETAILED DESCRIPTION

[0085] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0086] like Figure 1 As shown in Figure 1, a soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element is proposed.

[0087] Step 1: Divide the tissue and organ regions into linear regions and nonlinear regions according to the 3D grid data, and pre-calculate the data of the units in the linear region and the nonlinear region to obtain the required unit information.

[0088] The division of tissue and organ regions is mainly to determine which part is the linear region and which part is the nonlinear region. Specifically, by analyzing the concentrated area of ​​surgical operations, the area can be divided into nonlinear regions, and a more accurate nonlinear finite element model can be used for simulation to ensure the accuracy of surgical simulation. The part far away from the surgical operation area is divided into a linear region, and a linear finite element model with higher computational efficiency is used to reduce the burden of calculation. In the division process, factors such as the shape of the mesh model, the scale and position of the nonlinear finite element must also be considered to ensure the rationality and effectiveness of the division. At the same time, corresponding pre-calculated data is generated for each area to obtain the coefficient matrix of the initial coordinates. 0 D X and the initial deformation submatrix of any node for calculating the strain-displacement matrix , so as to carry out real-time simulation calculations later.

[0089] In order to facilitate the program to control the scale of the grid, the existing deformation simulation program is used to divide the regular space to generate a set of finite element units and nodes. First, the rectangular space is divided into a set of cubes, and each cube is further divided into five tetrahedrons, such as Figure 2 As shown in A and B, they represent two types of cube segmentation. Both types of units are necessary. If only one of them is used, such as type A, the edges on the interface of adjacent cubes will cross instead of overlap. However, if the type A unit is only coplanar with the type B unit, and the type B unit only has a common plane with the type A unit, then the generated mesh is correct. You can control the number of generated meshes, the shape and size of each unit, and the location of the mesh. These flexibilities bring convenience to deformation simulation. The program generates a simple three-dimensional mesh containing 120 nodes and 300 tetrahedral units, such as Figure 3 shown.

[0090] In order to use the hybrid finite element algorithm to solve the soft tissue deformation, the units in the model need to be divided into two categories: linear units and nonlinear units. For complex tissues and organs or specific surgical operations, the division of the surgical area may become a difficult task, and it is necessary to consider factors such as the shape of the mesh model, the scale and position of the nonlinear finite elements. Therefore, when the algorithm is actually used for surgical simulation, it is necessary to prepare a specific mesh model for a specific surgical scenario. In the simplest case, the units can be divided into levels, such as Figure 3 There are 5 layers of units in it, some of which can be set as linear units, and the rest as nonlinear finite elements. Assuming that the surgical interaction occurs on the upper surface of the cuboid-shaped soft tissue mesh, the upper part of the cuboid is set as a nonlinear unit and the lower unit is set as a linear unit.

[0091] The finite element calculation process mainly focuses on the mesh unit information and the connection relationship between the unit and the node. The graphics drawing process mainly focuses on the object surface mesh, so the data structure of the deformation simulation program is as follows: Figure 4As shown. The elements in the grid are represented in the form of a set. The grid mainly contains three lists: tetrahedron list, surface polygon (triangle) list and node list. There is direct association information in each list. Each tetrahedron is associated with 4 nodes. The system stores the node number in the corresponding node associated with the tetrahedron. All information related to the node can be found by numbering the tetrahedron. The node number stored in the tetrahedron unit is represented by an arrow, and the arrow points to the unit corresponding to the number. The triangle list also uses a similar method to maintain the connection information with the corresponding node. Each element also includes other information related to deformation calculation, such as pre-calculated information and unit volume at the initial time in each unit. The node mass, the initial position of the node, the resultant force on the node, the node displacement information, etc. must be saved in the node.

[0092] It is not enough to simulate deformation by only including node coordinates, tetrahedron and triangle topology information, so Figure 4 It is just a schematic diagram, and there is still a lot of information that is not shown. For real-time simulation, the calculation should be completed in the preprocessing stage as much as possible, and the data should be preprocessed directly at the simulation node. For each unit, use equations (18) and (22) to calculate the spatial derivative of its shape function and the linear strain-displacement matrix, calculate the initial volume of the unit, and calculate the mass and damping of each node. Generate topological information between mesh nodes, edges, units, and surface polygons. The unit information should also include the pre-calculated physical quantities in equation (18) and the pre-calculated data in equation (7) mentioned later. t D X And the initial volume of the unit; the node information should also include the node mass (obtained through a pre-made physical model or measurement), the initial position of the node (determined during modeling), the resultant force of the unit and external forces on the node (obtained by calculating the unit force on the node and the external force and summing them up), the displacement of the node at the previous moment (obtained by calculating or recording the previous time step), the displacement of the node at the current moment (obtained by calculating the deformation of the current time step) and the displacement of the node at the next moment (predicted by the explicit time integration method). The present invention uses tetrahedrons as units to divide the entire object area, and each unit contains four nodes. Consider that there is currently a tetrahedron unit with nodes P1, P2, P3, and P4. Figure 5 shown.

[0093] The nodes of the tetrahedron are P1, P2, P3, and P4 respectively, and the unit product of the tetrahedron is calculated as:

[0094] V=((P2-P1)×(P3-P1))·(P4-P1) / 6 (1)

[0095] Step 2: Unit force analysis: According to the area to which the finite element unit belongs, a nonlinear finite element model is used to solve the stress, strain and force of the unit on its nodes in the nonlinear area, and a linear finite element model is used to solve the stress, strain and force of the unit on its nodes in the linear area.

[0096] The hybrid finite element algorithm integrates the linear finite element model and the nonlinear finite element model in the unit analysis stage. The linear area and the nonlinear area in the surgical scene are connected by the intermediate nodes connecting the linear unit and the nonlinear unit to form a continuous surgical simulation scene. The unit force analysis directly calculates the force of the unit on each node connected to it, and then solves the resultant force of all the connected units on each node. The intermediate node will be subject to the force from the linear unit and the nonlinear unit, so the two areas can be connected through the intermediate node without assembling the overall stiffness matrix.

[0097] The simulation process of the finite element method is as follows: the continuum is idealized as a collection of several discrete units. The unit types are rectangles, triangles or cuboids, etc. The units are connected at the nodes. Each node is connected to several units. Each connected unit will exert a force on the node. The external forces are concentrated at these nodes. In the finite element force analysis process, the unit force analysis is first performed to calculate the force of each unit on its node, and then the node force analysis is performed to analyze the resultant force on each node. In the unit force analysis stage, the principle of virtual work is used to obtain the force of the unit on its node. In the node force analysis stage, the resultant force on each node is solved, and finally the displacement of the node in the unit is updated through the equilibrium equation.

[0098] During unit analysis, the dual-mode hybrid dynamics finite element deformation algorithm uses two finite element units and decides whether to use a linear model or a nonlinear model to solve the model based on the type of unit.

[0099] In the hybrid finite element model, the finite element units are divided into linear units and nonlinear units according to their locations, so the finite element formula can be written as:

[0100]

[0101] Among them, O represents the surgical area, N represents the non-operative area, and I represents the common area of ​​the two areas. The formulas of the two parts are similar, but the finite element formulas behind them are different. The following introduces the unit force analysis methods of linear finite element units and nonlinear finite element units.

[0102] For linear finite element units, the deformation of the object is assumed to be less than 10%, so the initial configuration can be used to approximate the configuration after deformation. The Cauchy strain tensor and Cauchy stress tensor are used to describe the strain and stress of the object respectively. The relationship between strain and stress in the linear model is fixed and does not change with the degree of deformation. Therefore, the computational efficiency of solving the linear model is higher than that of the nonlinear model, but the simulation accuracy is worse than that of the nonlinear model.

[0103] The tetrahedral unit with four nodes P1, P2, P3, P4 is as follows Figure 6 As shown, the initial positions of the four nodes are The positions after deformation are In order to calculate the Cauchy strain tensor, it is necessary to solve the displacement of each point in the tetrahedron and then find the displacement gradient. Assume that the initial coordinates of any point inside the tetrahedron are 0 x i , the position of this point after the unit is deformed is t x i , its weight relative to other nodes of the tetrahedron remains unchanged during the deformation process. The initial coordinates of any point are shown in the following formula:

[0104]

[0105] Where b = (b1 b2 b3) T is the weight vector, the coefficient matrix 0 D X for

[0106]

[0107] Similarly, the position after deformation t x i It can be expressed as:

[0108]

[0109] Where b = (b1 b2 b3) T is the weight vector, the coefficient matrix t D X for:

[0110]

[0111] The final position after deformation is:

[0112]

[0113] in, is a 3×3 matrix, which is the mapping matrix of the position of any point inside the tetrahedron before and after deformation. Then the displacement of any node is 0 x k represents the coordinate component of the node in the k direction at the initial moment, Represents the deformation of the node in the k direction at time t. Coefficient matrix It is constant during the simulation and can therefore be pre-calculated. This gives the displacement gradient The calculation formula is:

[0114]

[0115] Where I is the unit matrix, and the Cauchy strain tensor can be expressed as

[0116]

[0117] After calculating the strain tensor, the Cauchy stress σ is obtained. Since the Cauchy stress σ is defined as:

[0118]

[0119] Among them, df is the force on any small plane in space, dA is the area of ​​the plane, and n is the normal vector of the plane. Therefore, the resultant force of the stress on each surface of the tetrahedral unit can be obtained by the following formula:

[0120] f=(σ·n)A (12)

[0121] Among them, σ is the Cauchy stress, n is the normal vector of the tetrahedron surface triangle, and A is the triangular area of ​​the tetrahedron surface. Because the tetrahedron used is a linear finite element unit, the resultant force f of the force can be evenly distributed to each node of the surface triangle, and then the force of the tetrahedron on all nodes can be obtained: t F~ is the collection of stresses or forces at all nodes, which are calculated by distributing the resultant force equally to each node.

[0122] The nonlinear model of the dual-mode hybrid dynamic finite element deformation algorithm in the surgical area is described using the TLED algorithm. The main idea is: when the object undergoes large deformation or large distortion, linear strain tensors such as Cauchy strain cannot be used to solve it, because linear strain tensors are very inaccurate when the deformation is large. The commonly used nonlinear strain tensor is the Green strain tensor. When the strain tensor uses the Green strain tensor, the stress tensor often uses the Second Piola-Kirchhoff stress tensor. At the same time, the relationship between stress and strain of the material is not fixed. A relatively simple nonlinear model is the hyperelastic model, and many soft tissues can be approximately described as hyperelastic materials. The force exerted by the nonlinear unit on its nodes is

[0123]

[0124] in, is the vector form of the Second Piola-Kirchhoff stress tensor, represents the three normal stress components, represents the three shear stress components. is the strain-displacement matrix, 0 V is the initial volume of the unit. This formula performs force analysis at the unit level. All quantities in the formula can be obtained through local information on the unit (unit information calculated during preprocessing). This feature is very helpful for solving dynamic systems.

[0125] In order to solve the stress tensor in formula (13) Using the Neo-Hookean model, The components can be obtained by differentiating the Green strain tensor through the Neo-hookean strain energy function W as follows:

[0126]

[0127] Where μ is the shear modulus, δ ij is a second-order tensor representing the Kronecker delta function, λ is the first parameter of the Lamé constant, t J is the Jacobian determinant, is the Right Cauchy-Green deformation tensor, and

[0128]

[0129] Among them, the strain-displacement matrix It is directly related to the geometry of the element. In nonlinear analysis, the geometry of the element is constantly changing, so the strain-displacement matrix will also change accordingly. Since the tetrahedron contains 4 nodes, the strain-displacement matrix It can be written as a collection of its submatrices:

[0130]

[0131] in, is the strain-displacement matrix The i-th submatrix of corresponds to the i-th node of the unit, and the solution process is:

[0132]

[0133] in, is the deformation gradient of the nonlinear unit, the initial deformation submatrix of any node The calculation formula is:

[0134]

[0135] Solve for the initial deformation submatrix of any node Need to solve the derivative of the shape function with respect to the space coordinates ( 0 x, 0 y, 0 z) represents the initial spatial coordinates of any node. The spatial derivative of the shape function is usually solved by the isoparametric finite element formula. The main idea of ​​the isoparametric finite element formula is to connect the displacement of any point inside the unit with the node of the unit through the shape function. In the tetrahedron isoparametric unit, the shape function h on the four nodes P1, P2, P3, and P4 of the tetrahedron is i It can be expressed as h i (r,s,t), where the values ​​of the parameters r, s, and t are in the range [-1,1], and

[0136]

[0137] Among them, i represents the local label of the corresponding vertex in the corresponding unit.

[0138] According to the partial differential law, the partial derivative relationship between shape function and isoparameter is as follows:

[0139]

[0140] Here 0 x1, 0 x2, 0 x3 indicates equal parameters.

[0141] Formula (20) solves the relationship between the partial derivatives of the shape function with respect to the isoparameters, the partial derivatives of the shape function with respect to the initial coordinates, and the partial derivatives of the initial coordinates with respect to the isoparameters from the chain rule. After a series of identity transformations, we get:

[0142]

[0143] Therefore, the shape function h can be obtained i Derivatives with respect to space:

[0144]

[0145] Equation (22) shows how to compute the spatial derivative of the shape function on an isoparametric element. is the derivative of the shape function with respect to the isoparametric element. Since the isoparametric coordinates of each point in each element are constant, the function value of the derivative of the shape function with respect to the isoparametric element at the integration point remains constant throughout the simulation. The initial positions of the four nodes of each tetrahedral element are These positions are also unchanged during the deformation simulation process, so the spatial derivative of the shape function obtained by formula (22) is also unchanged during the entire deformation simulation process. Therefore, the initial deformation submatrix of any node can be obtained using formula (18): Then, the strain-displacement matrix of the element is obtained using formula (16): The force of the nonlinear element on its nodes can be obtained. The initial deformation submatrix of any node in each element This information remains constant throughout the deformation simulation and the initial deformation submatrix is ​​calculated It requires a large amount of calculation, so the initial deformation submatrix of any node in each unit can be The information is calculated in the pre-calculation stage and this pre-calculated data is directly used during the real-time simulation to accelerate the deformation simulation.

[0146] Step 3: Node force analysis: Calculate the resultant force of each node based on the forces of all units on their nodes.

[0147] After the linear and nonlinear region unit force analysis, all forces need to be applied to the nodes. The nodes between the two regions are called intermediate nodes. The units connected by the intermediate nodes are partly nonlinear units and partly linear units. Therefore, formulas (2) and (3) need to calculate the intermediate nodes independently, and solve the two sets of equations simultaneously. However, the present invention does not need to assemble the overall stiffness matrix, so a more direct method can be used to integrate the intermediate nodes into the algorithm. During the unit analysis, the force of each unit on the node to which it is connected will be calculated. This process does not restrict the nodes, so whether it is a node in the nonlinear region, a node in the linear region or an intermediate node, it will be affected by its adjacent units. After the unit analysis, the node force analysis is performed. This process calculates the resultant force received by each node, so all nodes are naturally included without the need for special processing.

[0148] After unit force analysis, the hybrid finite element algorithm has obtained the force of the tetrahedron unit on its nodes In the next step, we need to solve the resultant force of the unit forces on each node v. For each node in the grid, there are several units connected to it. Let the force of the unit e connected to it be t F e , then the resultant force F acting on the node is v for

[0149]

[0150] At this time, there may be an external force exerting a force R on the node v v , then the resultant force on the node is R v -F v .

[0151] The displacement of each point in the finite element unit is obtained by interpolating the displacement of the node. The resultant force of the external force on each point inside the unit is ultimately equivalent to the force of the node, and the force between adjacent units is also equivalent to the force of the unit on the node of the adjacent unit. Therefore, each node can be subjected to force analysis using Newton's law. For a static system, the node is at rest and the equilibrium relationship is:

[0152]

[0153] For a dynamic system without damping effect, the dynamic equilibrium equation of this node is

[0154]

[0155] Among them, m v is the mass of the node, u v is the displacement of the node, the node acceleration is the node displacement u v The second derivative with respect to time, F v is the force exerted by the unit on the node, R v is the external force acting on the node.

[0156] Similarly, for nodes with damping effects, the dynamic equilibrium equation is:

[0157]

[0158] The difference between formula (26) and formula (25) is that this formula has an additional damping term: where d v is the damping factor of the node, the node velocity is the node displacement u vThe first derivative of with respect to time.

[0159] The equation KU=R is a static equilibrium equation, which means that when an object is in a static equilibrium state, the product of the stiffness matrix and the displacement vector is equal to the external force vector, where K represents the stiffness matrix, U represents the displacement vector, and R represents the external force vector. All physical quantities in a static system are independent of time. Static equilibrium equations can also be applied to quasi-static scenarios. The state in a quasi-static scenario changes over time, but the change is so slow that the inertia and damping effects can be ignored. This type of system can be used in problems such as foundation settlement analysis and dam water injection, but tissue deformation and cutting operations in surgical scenarios will cause the system state to change rapidly over time. For example, tissues and organs will rebound after being pressed by surgical instruments, causing vibrations. Therefore, it is not appropriate to describe it using a static system equilibrium equation. When the displacement in the scene changes rapidly over time and the inertia effect cannot be ignored, it is appropriate to use a dynamic system to describe the scene. In this system, acceleration information will appear in the control equations, and the damping effect related to velocity should also be included. In the equilibrium equation, acceleration is the second-order derivative of displacement with respect to time, and velocity is the first-order derivative of displacement with respect to time. Compared with static systems, displacement variables in dynamic systems become more important because both velocity and acceleration are derivatives of displacement with respect to time. The basic control equations of dynamics are expressed by force balance, and the mechanical principle used is Newtonian mechanics.

[0160] In static analysis, an exact solution can usually be obtained, but in dynamic analysis, it is much more difficult to obtain an exact solution because the distribution of materials and damping effects must be considered. In dynamic finite element analysis, at least a mass matrix must be generated to be used in conjunction with the stiffness matrix. Usually, the assembly of the mass matrix M and the assembly of the stiffness matrix K can be performed in parallel, which can significantly improve the computational efficiency. The unit mass matrix is ​​generated in the local coordinate system and transformed to the global coordinate system. The method of assembling the unit mass matrix into a global mass matrix is ​​the same as the method of assembling the overall stiffness matrix. The biggest difference between the mass matrix and the stiffness matrix is ​​that the mass matrix can be converted into a diagonal matrix with the help of the direct lumping method. Because the solution of the inverse matrix of a diagonal matrix has a great advantage in calculation over the solution of a general inverse matrix, and the diagonal matrix brings convenience to the time integration of the unit level, the method of lumping the mass has a great advantage in the calculation process.

[0161] There are many methods to generate concentrated mass matrices. Here we only introduce the direct mass lumping method. The goal of direct mass lumping is to generate a diagonal concentrated mass matrix (DLMM) for a dynamic system. The method used is to directly distribute all the mass of the unit e to each node according to the degree of freedom of the node, without considering the influence between the nodes. The simplest example is a prismatic rod element with two nodes: length l, cross-sectional area A, mass density ρ, such as Figure 6 shown. Figure 6 The total mass of the one-dimensional unit in is ρAl, which is evenly divided into two parts and distributed to the two nodes of the unit. Similarly, this method can be used to concentrate the mass of a tetrahedral unit with a volume of V and a density of ρ on four nodes, such as Figure 7 shown.

[0162] The finite element mass matrix is ​​converted into a diagonal matrix by the direct mass concentration method. On this basis, it is assumed that the damping matrix D is proportional to the mass matrix, that is, D = αM, α is a constant, then D is also expressed as a diagonal matrix. The stiffness term K(U)U can be written as follows:

[0163]

[0164] Among them, e represents the unit in the set, Represents the force of the unit on the node, which is the force F e The vector of the internal forces. U represents the displacement vector. Equation (27) indicates that the internal forces acting on each node can be calculated by assembling the overall stiffness matrix K and multiplying the displacement vectors of all nodes, that is, the formula K(U)U; it can also be solved by performing unit stress analysis at the unit level, solving the unit forces on all its nodes, and then calculating the resultant force of the unit forces on each node. The two solution methods have the same calculation results, but the calculation methods of the two methods are very different. The first method has pre-calculated the information of the final overall influence between the nodes in the process of assembling the overall stiffness matrix. The second method first calculates the unit force on each node, and then solves the resultant force on each node based on the topological relationship. The first method is more efficient because the topological information is processed in the process of solving the overall stiffness matrix, while the second method is more flexible in dealing with topological changes, because the topological information is used in the process of calculating the node force, so it is more suitable for applications with mesh topology changes.

[0165] Because the mass matrix M and the damping matrix D are both diagonal matrices, that is, the matrix has non-zero elements only on the diagonal, and the elements at other positions are all 0, that is, the mass and damping information of each node will only affect the node itself and have nothing to do with other nodes. Assume that due to the deformation of the object, the resultant force of the internal stress on node i is Fi , whose mass is M ii , damping is D ii , the current displacement is U i , the external force is R i , then the dynamic equilibrium equation of node i is:

[0166]

[0167] Therefore, equations (27), (28) and the properties of the mass matrix and damping matrix can be used to solve the deformation simulation at the element and node level. represents the velocity of node i, represents the acceleration of node i.

[0168] Step 4: Solving the hybrid model: Use explicit time integration to solve the displacement of each node of the tissue organ (unit) at the current moment, and obtain the position of each vertex at the next moment, so as to obtain the deformation information at the next moment.

[0169] In addition to the aforementioned use of different models for different regions to reduce the amount of calculation in the real-time deformation simulation process, the real-time deformation of tissue and organ deformation in surgical simulation also requires the use of appropriate numerical integration strategies to solve the dynamic equilibrium equations. The integral in the time domain can be solved by explicit integration and implicit integration. The commonly used implicit integration method is unconditionally stable, but it requires an iterative method to solve a set of nonlinear algebraic equations, so the amount of calculation is very large. In contrast, the explicit integration method does not require iteration and has a small amount of calculation.

[0170] The central difference method is an explicit integration method for approximate solutions to second-order differential equations. The central difference method is used to solve the dynamic equilibrium equation of the hybrid finite element model. This method must be able to solve the information of the next state under the current state. Assuming that the initial simulation time is time 0, the current time is time t, and the time step is Δt, it can be known that the next simulation time is t+Δt. The unknown function in the dynamic equilibrium equation of formula (28) is U i , and other quantities have been calculated when solving the equation. The main idea of ​​the central difference method is to use differences instead of differentials and solve the unknown function U through iteration. i Function value at time t+Δt t+Δt U i , where the upper left subscript uniformly represents the moment when the physical quantity is measured, and solves U at time t i The first derivative of The principle is as follows Figure 8 shown.

[0171] Second-order derivative The principle is similar to this. For the unknown function U iThe first and second derivatives at time t are replaced using equations (29) and (30).

[0172]

[0173] Substituting formulas (29) and (30) into the dynamic equilibrium equation, we can obtain the equilibrium equation at the node:

[0174]

[0175] in, t F i and t R i They represent the unit force and external force acting on node i respectively.

[0176] Finally, we get the following through identity transformation:

[0177]

[0178] Formula (32) can efficiently solve the displacement of each node at the next moment, and then obtain the deformation information of the object in the dynamic system.

[0179] The implicit integration method is unconditionally stable, while the explicit integration method is only conditionally stable. The time step Δt must be less than the critical step Δt cr In linear elastic analysis, the critical step length Δt cr for:

[0180]

[0181] Among them, L e is the minimum characteristic unit length in the unit set, c is the longitudinal wave velocity of the material, and

[0182]

[0183] Where E is Young's modulus, ν is Poisson's ratio, and ρ is material density. When performing either material nonlinear or geometric nonlinear analysis, a smaller critical time is required than in linear analysis because the longitudinal wave velocity c usually becomes higher as the shape of the object changes, but in general, this can be reduced by using a critical step size Δt cr Measure the time step.

[0184] Although the explicit integration method is conditionally stable, it can be seen from formulas (34) and (35) that when the Young's modulus of the material is very small, that is, the material is very soft, the critical time step of explicit integration is relatively large, so it is possible to stably simulate the deformation of soft tissues at a longer time step. Human tissues and organs, especially internal organs, are relatively soft, so explicit integration can be used for simulation.

[0185] In order to verify the real-time performance of the hybrid finite element model, the present invention implements the simulation model of the hybrid finite element model and the simulation model of the TLED deformation algorithm on a standard PC, and compares the simulation results of the hybrid finite element model with the TLED deformation algorithm in terms of real-time performance. The nonlinear soft tissue constitutive relationship in all simulations adopts the new Hooke model, the Young's modulus is set to 5000, the Poisson's ratio is 0.35, and the material density is 1000kg / m 3 , the running effect is as follows Fig. 9 shown. Fig. 9 In (a), all units are linear units, in (b), all units are nonlinear units, and in (c), the upper half is linear units and the lower half is nonlinear units.

[0186] Since the refresh rate of the graphics card is fixed, for example, the frame rate of the simulation PC used is 75Hz, in order to visually observe the performance difference between the two algorithms, a larger grid scale is used in the simulation of the present invention to reduce the frame rate to a range that is convenient for observation and comparison. Otherwise, when the grid is small, the algorithm frame rate will be limited to the refresh frame rate of the graphics card. In the simulation, the soft tissue grid containing 3564 nodes and 14450 units is deformed using the dual-mode hybrid dynamic finite element deformation algorithm of the present invention to simulate the deformation, and the direct relationship between the proportion of linear units and the simulation frame rate is investigated. The frame rate changes obtained by simulation are shown in Table 1.

[0187] Table 1 Performance analysis of dual-mode hybrid dynamics finite element deformation algorithm

[0188]

[0189] It can be seen from Table 1 that as the proportion of linear units in the hybrid model continues to increase, the frame rate also gradually increases. This is mainly because a part of linear finite element units are mixed in the model, and the computational efficiency of linear finite element units is higher than that of nonlinear finite elements. When the proportion of linear units is 0, all units in the model are nonlinear units. It can be seen that the dual-mode hybrid dynamic finite element deformation algorithm is superior to the TLED algorithm in terms of computational efficiency, but the efficiency advantage of the algorithm comes at a price. The error analysis later shows that the error of this method is gradually increasing as the efficiency increases.

[0190] In the simulation, different displacements are applied to some nodes of the grid and the proportion of linear finite element units in the grid is gradually increased to verify the difference in simulation results between the dual-mode hybrid dynamic finite element deformation algorithm and the TLED algorithm under the corresponding state. The simulation scenario is as follows Fig.10 shown. Fig.10The middle rectangular block area is divided into a grid containing approximately 4000 tetrahedral units. The bottom surface of the cuboid is fixed, and the top surface is the surgical interaction area. The upper units of the cuboid use nonlinear units, and the edges are represented by light colors. The lower units use linear units, and the edges are represented by dark colors. The ratio of linear units to nonlinear units in the hybrid finite element model is recorded as r. The TLED algorithm also uses the same grid and applies the same displacement. The difference is that all units in the grid processed by the TLED algorithm are nonlinear. Therefore, the TLED algorithm can be regarded as a special case of the dual-mode hybrid dynamic finite element deformation algorithm r = 0. When a downward displacement is applied to the soft tissue at the center of the upper surface of the cuboid, a depression will appear on the upper surface. Fig.10 Shown is a simulation scenario for r=30% error analysis.

[0191] In order to measure the accuracy difference between the dual-mode hybrid dynamic finite element deformation algorithm and the TLED algorithm of the present invention, it is necessary to measure the "distance" between the deformed meshes processed by the two algorithms. The shape of the deformed mesh is directly related to the position of the mesh nodes. The difference in deformation of objects with the same position at the initial moment can be measured by the difference in displacement of the corresponding nodes, which is defined as the modulus of the difference in displacement of the corresponding nodes, and the average difference is used to measure the calculation error of the dual-mode hybrid dynamic finite element deformation algorithm and the TLED algorithm. Since the present invention considers a dynamic system including inertia and damping, the deformation of the model is time-series. Use two algorithms to apply displacement to the surface of the cuboid at the initial moment, and calculate the average error after 1000 iterations.

[0192] When a displacement vector with a length of 0.6 is applied to the upper surface of the cuboid, the average errors when the linear unit ratio is 10%, 20%, 30% and 40% are tested respectively, as shown in Table 2.

[0193] Table 2 Error analysis of hybrid finite element model

[0194]

[0195] Table 2 shows that as the proportion of linear units in the grid increases, the error of the algorithm gradually increases. This is mainly because the units in the grid will deform as the displacement is applied. If there are linear units that undergo large deformation at this time, the error will increase. Therefore, as the number of linear units in the model increases, the average error of the algorithm will become larger and larger.

[0196] When the displacement length applied to the upper surface of the cuboid is 1, the average error when the linear unit ratio is 10%, 20%, 30% and 40% is shown in Table 3.

[0197] Table 3 Error analysis of hybrid finite element model II

[0198]

[0199]

[0200] From the simulation results and analysis in Table 3, it can be seen that as the applied displacement increases, the error of the algorithm will also increase. This is mainly because the increase in displacement will increase the deformation of the rectangular soft tissue, causing more units to produce large deformations, so the error of the model will increase.

[0201] As the number of linear units in the hybrid finite element model increases, the computational efficiency will become higher and higher, while the computational accuracy will gradually decrease; as the displacement applied to the upper surface of the soft tissue increases, the error will increase accordingly, which is determined by the characteristics of the linear finite element, because the linear finite element processing is only suitable for small deformations, and increasing the applied displacement will increase the degree of deformation of the object, so the error will also increase. The simulation accuracy and efficiency of the hybrid finite element algorithm of the present invention are related to the degree of displacement of the nodes in the grid and the proportion of the linear units in the grid. Therefore, the computational efficiency and computational accuracy can be adjusted according to the specific simulation requirements to meet the needs of different simulations, which is more flexible than a completely linear finite element model or a completely nonlinear finite element model.

[0202] This paper deeply explores the application of dual-mode hybrid dynamic finite element algorithm in surgical simulation, aiming to improve the accuracy and efficiency of soft tissue deformation simulation. Through the organic combination of linear and nonlinear finite element methods, an innovative dual-mode hybrid strategy is successfully implemented, which can dynamically select the appropriate calculation model according to the deformation characteristics of soft tissue.

[0203] In terms of algorithm principle, the basic framework and implementation steps of the dual-mode hybrid dynamic finite element algorithm are elaborated in detail. Through key steps such as discretization of the continuum, node force analysis, and explicit time integral solution, the deformation behavior of soft tissue can be accurately simulated. At the same time, the specific application of the algorithm in the surgical simulation system is discussed, demonstrating its potential and value in actual surgical simulation.

[0204] In terms of experimental verification, multiple test cases were designed to comprehensively evaluate the performance of the dual-mode hybrid dynamic finite element algorithm. The experimental results show that the present invention can not only significantly improve the realism and accuracy of soft tissue deformation simulation while maintaining high efficiency, but also show good stability and adaptability. The algorithm has shown excellent performance when dealing with large-scale soft tissue deformation, complex boundary conditions and multiple physical properties.

[0205] Compared with other algorithms, the dual-mode hybrid dynamic finite element algorithm of the present invention shows obvious advantages in calculation time, accuracy and adaptability. It not only provides a more efficient and accurate soft tissue deformation simulation method for the field of surgical simulation, but also provides new ideas and methods for scholars and researchers in related fields.

[0206] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element, characterized in that: The steps are as follows: Step 1: dividing the tissue organ region into a linear region and a nonlinear region according to the 3D grid data, and performing data pre-calculation on the units in the linear region and the nonlinear region to obtain pre-calculation information of the units; Step 2: Using the pre-calculated information of the unit, according to the area to which the unit belongs, a nonlinear finite element model is used to solve the stress, strain and the force of the unit on its nodes in the nonlinear area, and a linear finite element model is used to solve the stress, strain and the force of the unit on its nodes in the linear area; Step 3: Calculate the resultant force of the unit forces on each node based on the forces of all units on their nodes; Step 4: Use explicit time integration to solve the displacement of each node of the unit at the current moment, obtain the position of each node at the next moment, and obtain the deformation information at the next moment.

2. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 1 is characterized in that: The nonlinear area is the concentrated area of ​​surgical operations, and the linear area is the part far away from the surgical operation area; the rectangular space is divided into a set of cubes, each cube is further divided into five tetrahedral units, and one tetrahedral unit is associated with four nodes; the units in the linear area are linear units, and the units in the nonlinear area are nonlinear units; the linear area and the nonlinear area in the surgical scene are connected by the intermediate nodes connecting the linear units and the nonlinear units to form a continuous surgical simulation scene; the force of the unit on each node connected to it is directly calculated through unit force analysis, and the resultant force of all connected units on each node is solved. The intermediate node is subjected to the force from the linear unit and the nonlinear unit, and the two areas can be connected through the intermediate node.

3. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 2 is characterized in that: The data obtained by precalculation includes the coefficient matrix of the initial coordinates 0 D X and the initial deformation submatrix of any node for calculating the strain-displacement matrix; The deformation calculation of tetrahedral elements includes mesh unit information and the connection relationship between the unit and the node. The elements in the mesh are represented in the form of a set. The mesh contains three lists: tetrahedral list, surface polygon list and node list. The tetrahedral list stores the node number at the corresponding node associated with the tetrahedral element, and all information related to the node can be found through the number. The node number stored in the tetrahedral element is represented by an arrow, and the arrow points to the unit corresponding to the number; the surface polygon list maintains the connection information with the corresponding node; Each element also contains information related to deformation calculation, including pre-calculated information in each unit and the initial volume of the unit; the node list stores the node mass, the initial position of the node, the resultant force of the unit and external forces acting on the point, the displacement of the node at the previous moment, and the node displacement information; In the data pre-calculation, for each unit, the spatial derivative of its shape function and the linear strain-displacement matrix are calculated, the initial volume of the unit is calculated, the mass and damping of each node are calculated, and the topological information between the nodes, edges, units and surface polygons of the mesh is generated.

4. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 2 or 3, characterized in that: The Cauchy strain tensor and Cauchy stress tensor are used to describe the strain and stress of the linear unit respectively; the initial positions of the four nodes of the tetrahedral unit are The positions after deformation are The method for calculating the Cauchy strain tensor is: solve the displacement of each node in the tetrahedral unit, and then obtain the displacement gradient; after calculating the Cauchy strain tensor, obtain the Cauchy stress tensor σ; In nonlinear units, the strain tensor uses Green's strain tensor, and the stress tensor often uses Second Piola-Kirchhoff stress tensor.

5. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 4 is characterized in that: The force exerted by the nonlinear unit on the node is: in, is the vector form of the Second Piola-Kirchhoff stress tensor, represents the three normal stress components, represents the three shear stress components; is the strain-displacement matrix, 0 V is the initial volume of the unit; In the linear unit, the resultant force f = (σ·n)A of the stress on each surface of the tetrahedral unit is evenly distributed to each node of the surface triangle, and the force of the tetrahedral unit on all nodes is obtained. Where σ is the Cauchy stress tensor, n is the normal vector of the tetrahedron surface triangle, and A is the area of ​​the tetrahedron surface triangle.

6. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 5 is characterized in that: The stress tensor The components of are obtained by differentiating the Green strain tensor through the Neo-hookean strain energy function: Where μ is the shear modulus, δ ij is a second-order tensor, λ is the first parameter of the Lame constant, t J is the Jacobian determinant, is the Right Cauchy-Green deformation tensor, and The strain-displacement matrix Written as a collection of submatrices: in, is the strain-displacement matrix The i-th submatrix of , and: in, is the deformation gradient of the nonlinear unit, the initial deformation submatrix of any node The calculation formula is: Solve the initial deformation submatrix of any node Need to solve the derivative of the shape function with respect to the space coordinates ( 0 x, 0 y, 0 z) represents the initial spatial coordinates of any node; the spatial derivative of the shape function is solved by the isoparametric finite element formula, According to the partial differential law, the shape function h of the node i The partial derivative relationship of and isoparameters is: in, 0 x1, 0 x2, 0 x3 all represent equal parameters; The relationship between the partial derivatives of the shape function with respect to the isoparameters and the partial derivatives of the shape function with respect to the initial coordinates and the partial derivatives of the initial coordinates with respect to the isoparameters is solved by the chain rule. After a series of identity transformations, we get: Obtain the shape function h i Derivatives with respect to space: in, is the derivative of the shape function with respect to the isoparameter in the isoparametric element; The method for calculating the Cauchy strain tensor is as follows: Assume that the initial coordinates of any node inside the tetrahedron unit are 0 x i , the position of the node after the tetrahedral element is deformed is t x i , the initial coordinates of any node are: Where b = (b1 b2 b3) T is the weight vector, the coefficient matrix 0 D X for Similarly, the position after deformation t x i for: Where b = (b1 b2 b3) T is the weight vector, the coefficient matrix t D X for: The final position after deformation is: in, is the mapping matrix of the position of any node inside the tetrahedron unit before and after deformation, and the displacement of any node 0 x k represents the coordinate component of the node in the k direction at the initial moment, It represents the deformation of the node in the k direction at time t; Get the displacement gradient The calculation formula is: Where I is the identity matrix and the Cauchy strain tensor is expressed as The Cauchy stress tensor σ is: Among them, df is the force acting on any tiny plane in space, dA is the area of ​​the tiny plane, and n is the normal vector of the tiny plane.

7. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 5 or 6, characterized in that: The method for calculating the resultant force of the unit forces acting on each node is as follows: for each node in the grid, there are a number of units connected to it. Suppose the force acting on the node by the unit e connected to the node is t F e , then the resultant force F acting on the node v is v for If there is an external force acting on the node v, a force R v , then the net force on node v is R v -F v ; Each node uses Newton's law to analyze the force. For a static system, the node is at rest and the equilibrium relationship is: For a dynamic system without damping effect, the dynamic equilibrium equation of the node is: Among them, m v is the mass of the node, u v is the displacement of the node, is the node acceleration, F v is the force exerted by the unit on the node, R v is the external force on the node; For nodes with damping effects, the dynamic equilibrium equation is: Among them, d v is the damping factor of the node, is the node speed.

8. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 7, characterized in that: In the dynamic system, the finite element mass matrix is ​​converted into a diagonal matrix by the direct mass concentration method. Assuming that the damping matrix D is proportional to the mass matrix D = αM, α is a constant, then the damping matrix D is expressed as a diagonal matrix, and the stiffness term K(U)U of the dynamic equilibrium equation is: Among them, e represents the unit in the set, Represents the force of the unit on the node, which is the force F e A vector; U represents a displacement vector; Due to the deformation of the object, the resultant force of the internal stress on node i is F i , the mass of node i is M ii , damping is D ii , the current displacement is U i , the external force is R i , then the dynamic equilibrium equation of node i is: represents the velocity of node i, represents the acceleration of node i.

9. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 8, characterized in that: The dynamic equilibrium equations are solved by using the central difference method to solve the hybrid finite element model consisting of the nonlinear finite element model and the linear finite element model; The central difference method uses difference instead of differentiation and solves the displacement U through iteration. i Function value at time t+Δt t+Δt U i , solve for the displacement U at time t i The first derivative of The method is: For the unknown displacement U i The first and second derivatives at time t are Substitute into the dynamic equilibrium equation The equilibrium equations at the nodes are obtained: in, t F i and t R i They represent the unit force and external force on node i respectively; Through the identity transformation we get: The displacement of each node at the next moment can be solved, and then the deformation information of the object in the dynamic system can be obtained.

10. The soft tissue deformation simulation method based on dual-mode hybrid dynamics finite element according to claim 9, characterized in that: The time step Δt is less than the critical step Δt cr ; Critical step length Δt in linear elastic analysis cr for: Among them, L e is the minimum characteristic unit length in the unit set, and c is the longitudinal wave velocity of the material.

Images