Soft body physics simulation method and system based on viscoelastic constraint position dynamics

Abstract

This application relates to the field of software simulation technology and discloses a software physical simulation method and system based on viscoelastic constrained position dynamics. The method includes: viscoelastic constitutive calibration by fusing a three-parameter model and the Ogden model to construct an organ-specific parameter set; and further establishing a viscoelastic constraint system through volume-preserving constraints and tensile constraints. Based on this system, strain rate detection and adaptive iterative solution are performed, and the memory effect is calculated by combining temperature field adjustment and historical strain buffer, ultimately obtaining an accurate deformation response sequence and evolution state. This application solves the problem that existing software physical simulation methods cannot accurately simulate temperature changes and fatigue accumulation effects.

Description

Technical Field

[0001] This application relates to the field of software simulation technology, and in particular to a software physical simulation method and system based on viscoelastic constrained position dynamics. Background Technology

[0002] Existing soft physics simulation techniques are mainly based on linear elasticity or mass-spring models, achieving physical simulation by simulating the deformation of objects. These methods exhibit good real-time performance in many applications, especially in simple soft physics simulations, where they can handle some basic tension and compression phenomena. However, these traditional methods have limitations in accurately simulating complex viscoelastic properties and the effects of high temperatures on material behavior. This results in poor performance in some high-precision applications, failing to truly reflect the physical properties of soft materials, particularly in terms of fatigue accumulation and deformation effects under long-term stress and temperature influences.

[0003] The shortcoming of existing technologies lies in the fact that most simulation methods fail to adequately consider the long-term effects of high-temperature effects and fatigue accumulation on materials. Traditional methods typically use fixed material models, neglecting the changes in the physical properties of soft materials under repeated deformation and high-temperature conditions over time. For example, when simulating the cutting of soft tissue by surgical instruments, existing methods often fail to reflect the viscoelastic changes in tissue caused by temperature increases, nor do they consider fatigue damage to tissue under prolonged contact or repeated loading. This deficiency limits their applicability and accuracy in complex application environments, especially in scenarios requiring precise simulation of tissue deformation, damage, and fatigue effects.

[0004] To address the shortcomings of existing technologies, the technical problems that current methods cannot accurately model have been further derived, mainly concerning the calculation of heat conduction, temperature-viscoelastic coupling, and cumulative fatigue effects. Traditional simulation methods lack dynamic responses to temperature changes and cannot simulate the profound impact of localized high-temperature regions caused by surgical instruments or environmental factors on tissue mechanical behavior. Furthermore, the lack of long-term dynamic tracking of fatigue accumulation effects makes it impossible to accurately simulate the damage process and deformation recovery behavior of materials under long-term repeated deformation. Therefore, developing a simulation method that comprehensively considers temperature field changes, thermodynamic corrections, and fatigue accumulation effects has become a key technical challenge in solving this problem. Summary of the Invention

[0005] This application provides a software physical simulation method and system based on viscoelastic constrained position dynamics, which solves the problem that existing software physical simulation methods cannot accurately simulate temperature changes and fatigue accumulation effects, improves the accuracy and reliability of the simulation process, and can more realistically reflect the physical behavior of materials, especially under high temperature and multiple loading conditions.

[0006] In a first aspect, this application provides a software physics simulation method based on viscoelastic constrained position dynamics, the software physics simulation method based on viscoelastic constrained position dynamics comprising:

[0007] Step S1: Viscoelastic constitutive calibration of soft tissue mechanical parameters is performed by fusing the three-parameter model and the Ogden model to obtain an organ-specific parameter set;

[0008] Step S2: Based on the organ-specific parameter set, construct a multi-level system of volume-preserving constraints and stretching constraints for the vertices of the tetrahedral mesh to obtain a viscoelastic constraint system;

[0009] Step S3: Based on the viscoelastic constraint system, the constraint projection is adaptively iteratively solved by strain rate detection to obtain the deformation response sequence;

[0010] Step S4: Based on the deformation response sequence, the viscosity coefficient is adjusted by temperature field coupling using the heat conduction equation to obtain thermodynamic correction parameters;

[0011] Step S5: Calculate the memory effect of the cumulative deformation by combining the thermodynamic correction parameters with the historical strain buffer to obtain the evolution state.

[0012] Secondly, this application provides a soft physics simulation system based on viscoelastic constrained position dynamics, the soft physics simulation system based on viscoelastic constrained position dynamics comprising:

[0013] The calibration module is used to perform viscoelastic constitutive calibration of soft tissue mechanical parameters by fusing the three-parameter model and the Ogden model, and obtain an organ-specific parameter set.

[0014] The construction module is used to construct a multi-level system of volume-preserving constraints and stretching constraints on the vertices of the tetrahedral mesh based on the organ-specific parameter set, thereby obtaining a viscoelastic constraint system.

[0015] The solution module is used to adaptively iteratively solve the constraint projection based on the viscoelastic constraint system through strain rate detection to obtain the deformation response sequence;

[0016] The adjustment module is used to adjust the viscosity coefficient by temperature field coupling according to the deformation response sequence using the heat conduction equation, so as to obtain thermodynamic correction parameters;

[0017] The calculation module is used to calculate the memory effect of cumulative deformation by combining the thermodynamic correction parameters with the historical strain buffer to obtain the evolution state.

[0018] Thirdly, a software physics simulation device based on viscoelastic constrained position dynamics is provided, comprising: a memory and at least one processor, wherein the memory stores instructions; the at least one processor invokes the instructions in the memory to cause the software physics simulation device based on viscoelastic constrained position dynamics to execute the aforementioned software physics simulation method based on viscoelastic constrained position dynamics.

[0019] Fourthly, a computer-readable storage medium is provided, wherein instructions are stored therein, which, when executed on a computer, cause the computer to perform the aforementioned software physical simulation method based on viscoelastic constrained position dynamics.

[0020] The technical solution provided in this application solves the problem of existing soft-body physical simulation technologies being unable to accurately simulate complex material properties by introducing a soft-body physical simulation method based on viscoelastic constrained position dynamics, particularly in applications involving high temperatures, fatigue accumulation, and deformation memory. By integrating a three-parameter model and the Ogden model, this application achieves a more refined calibration of soft tissue mechanics, enabling simulations to accurately reflect the elastic and viscous properties of soft materials under different deformation conditions. Especially in cases involving high temperatures or periodic loading, the simulation method in this application can not only handle instantaneous mechanical behavior but also accurately track the memory effect of materials under multiple deformations and temperature changes, thus overcoming the limitations of traditional simulation methods in adapting to these complex conditions. Through adaptive strain rate adjustment and temperature field coupling, the system can adjust the physical properties of the material in real time, ensuring the accuracy and stability of simulation results under different environments. This innovative approach provides strong technical support for applications involving soft body deformation and thermal effects in surgical simulation, medical image modeling, and virtual reality.

[0021] By accurately simulating the deformation and temperature effects of soft tissue under the action of surgical instruments, this application enables doctors to perform precise surgical operations in a virtual environment, thereby effectively reducing the risks of actual surgery. Furthermore, innovations in handling repeated loading and fatigue damage, and in modeling fatigue accumulation and material degradation, allow this technical solution to not only simulate conventional deformation but also provide dynamic adjustment functions based on material fatigue characteristics. This enables the system to more realistically reflect the performance degradation process of materials after long-term use. In these applications, the efficiency and accuracy of the algorithm directly affect the actual effect of the simulation system, driving the development of simulation technology towards greater complexity and precision, and broadening the application scope of simulation systems in practical operations. Therefore, the technical features of this application make the application of software simulation in the fields of medicine, industry, and even robotics more accurate and reliable, solving the high-precision requirements that traditional methods struggle to meet. Attached Figure Description

[0022] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This is a schematic diagram of an embodiment of the soft physics simulation method based on viscoelastic constrained position dynamics in this application.

[0024] Figure 2 This is a schematic diagram of one embodiment of the soft physical simulation system based on viscoelastic constrained position dynamics in this application.

[0025] Figure 3 This is a schematic block diagram of the structure of the soft physical simulation device based on viscoelastic constrained position dynamics in an embodiment of the present invention. Detailed Implementation

[0026] This application provides a soft physics simulation method and system based on viscoelastic constrained position dynamics. The terms "first," "second," "third," "fourth," etc. (if present) in the specification, claims, and accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms "comprising" or "having" and any variations thereof are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0027] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the soft physics simulation method based on viscoelastic constrained position dynamics in this application includes:

[0028] Step S1: Viscoelastic constitutive calibration of soft tissue mechanical parameters is performed by fusing the three-parameter model and the Ogden model to obtain an organ-specific parameter set;

[0029] Step S2: Based on the organ-specific parameter set, construct a multi-level system of volume-preserving constraints and stretching constraints for the vertices of the tetrahedral mesh to obtain a viscoelastic constraint system;

[0030] Step S3: Based on the viscoelastic constraint system, the constraint projection is adaptively iteratively solved by strain rate detection to obtain the deformation response sequence;

[0031] Step S4: Based on the deformation response sequence, the viscosity coefficient is adjusted by temperature field coupling using the heat conduction equation to obtain the thermodynamic correction parameters;

[0032] Step S5: Calculate the memory effect of cumulative deformation by combining thermodynamic correction parameters with historical strain buffer to obtain the evolution state.

[0033] It is understood that the executing entity of this application can be a software physical simulation system based on viscoelastic constrained position dynamics, or it can be a terminal or a server; the specific implementation is not limited here. This application's embodiments use a server as an example for illustration.

[0034] Specifically, viscoelastic constitutive calibration of soft tissue mechanical parameters was performed by fusing a three-parameter model and the Ogden model to obtain an organ-specific parameter set. First, a HY-0580 mechanical dynamic compression testing device was used to apply periodic compressive loads to soft tissue samples from the lungs, liver, kidneys, and heart. The measurement system acquired force-displacement response curves for each tissue sample at different strain rates. For the three-parameter model, a mechanical relationship was established involving a series-parallel structure of two springs and a damper. The first elastic constant characterizes the instantaneous elastic response, the second elastic constant characterizes the delayed elastic response, and the viscosity coefficient characterizes the viscous damping properties. By performing an inverse Laplace transform on the experimental data, the creep coefficient was calculated by numerically fitting the strain-time growth law under constant stress, while the relaxation coefficient was extracted by the stress-time decay law under constant strain. Subsequently, the Ogden model was introduced to describe the nonlinear characteristics of large deformation. This model defines the principal elongation of the material by the stretch ratio, and the nonlinear parameters control the stiffness amplitude and nonlinearity of the stress-strain curve. To achieve a smooth transition between the two models, a sigmoid weighting function is used for strain interval mapping. This function controls the transition rate from the linear to the nonlinear region through the transition slope, and the critical strain value defines the center position of the transition. During the least squares parameter fitting process, the sum of squared residuals between the theoretical predictions of the hybrid constitutive model and the measured force-deformation curves is minimized. Parameters such as the first elastic constant, the second elastic constant, the viscosity coefficient, the nonlinear parameter, the transition slope, and the critical strain value are iteratively adjusted until the residuals converge. For example, in lung tissue samples, the experimentally measured maximum force was 6.4 Newtons, average force was 2.74 Newtons, and initial stiffness was 1.8 Newtons per millimeter. The first elastic constant obtained by fitting with a three-parameter model is about 0.3 times that of liver tissue, reflecting the low stiffness of lung tissue. The viscosity coefficient is about 1.8 times higher than that of myocardial tissue, reflecting the high damping characteristics of alveolar structure. The small nonlinear parameter value of the Ogden model indicates that the stress increase of lung tissue is relatively slow under large deformation. The critical strain value of the sigmoid function is set at around 0.15, so that the weight is mainly dominated by the three-parameter model at small strain and gradually transitions to the Ogden model at large strain, ultimately forming an organ-specific parameter set containing all mechanical parameters of the organ.

[0035] Based on the organ-specific parameter set, a multi-level construction of volume-preserving and stretching constraints is performed on the vertices of the tetrahedral mesh to obtain a viscoelastic constraint system. First, the organ-specific parameter set obtained in step S1 is input into the Tetgen tetrahedral mesh generation tool. This tool reads the surface triangular mesh file of soft tissue and adaptively adjusts the internal node density according to the material stiffness information in the parameter set. For low-stiffness organs such as lung tissue, a sparser mesh is used, controlling the number of vertices to be between 745 and 2000. For medium-stiffness organs such as the liver, the number of vertices is set to 8000 to 15000. The number of vertices in heart muscle tissue is expanded to 29040. The volume-preserving constraint is defined as the first-level constraint relationship. It is calculated by cross-product of the vector obtained by subtracting the third vertex from the first vertex and the vector obtained by subtracting the fourth vertex from the first vertex, then by dot-product of the cross-product result with the vector obtained by subtracting the second vertex from the first vertex, and finally multiplied by a coefficient of one-sixth to obtain the current tetrahedral volume. This volume is compared with the initial volume to constitute the volume constraint deviation. The tension constraint is defined as the second layer of constraint relationship. For each edge of the tetrahedron, the difference between the current length and the initial edge length is calculated. Simultaneously, a mass weight is assigned to each vertex, with the mass weight inversely proportional to the number of tetrahedra connected to that vertex. Subsequently, the first elastic constant, the second elastic constant, and the viscosity coefficient from the organ-specific parameter set are embedded into the constraint stiffness coefficient. The dynamic calculation method for the constraint stiffness coefficient is the sum of the first and second elastic constants multiplied by an exponential decay function. The time constant of the exponential term is obtained by dividing the viscosity coefficient by the sum of the two elastic constants to obtain the relaxation time constant. For example, a tetrahedral mesh with 12,500 vertices is constructed for kidney tissue. The initial volume of a tetrahedron is 0.0082 cubic millimeters. The initial length of the longest side in the stretch constraint is 3.2 millimeters, and the length of the shortest side is 1.1 millimeters. Substituting the first elastic constant of 8,500 Pascals, the second elastic constant of 6,200 Pascals, and the viscosity coefficient of 420 Pascals into the kidney tissue parameter set into the calculation, the initial constraint stiffness is 14,700 Pascals. After 1 second, the stiffness decreases to 9,300 Pascals, which reflects the viscoelastic softening characteristics.

[0036] Based on a viscoelastic constraint system, the constraint projection is adaptively iteratively solved using strain rate detection to obtain the deformation response sequence. At the beginning of each time step, the positions of all vertices of the tetrahedral mesh are predicted. The prediction calculation is to add the current vertex position, velocity multiplied by the time step, external force divided by mass, and then multiply by the square of the time step. The time step is set to 0.01 seconds. After the prediction is completed, the initial deformation state is obtained, including the predicted position, current velocity, and external force vector of each vertex. The differential calculation process extracts the vertex position of the current time step from the initial deformation state, traverses each tetrahedral element to calculate the current length of the six edges and compares it with the initial length to obtain the strain. The strain is defined as the length change divided by the initial length. The strain values ​​of the same vertices and edges from the previous time step are read from the historical deformation record buffer. The difference between the two is calculated according to the vertex index and edge index. The difference result is divided by the time step of 0.01 seconds to obtain the strain rate parameter. The calculated strain rate parameters are compared with two preset thresholds. The first threshold is set at 0.5 per second as the criterion for high-speed deformation, and the second threshold is set at 0.1 per second as the criterion for low-speed deformation. When the strain rate in a certain region exceeds 0.5 per second, that region is marked as a stiffness reduction region, and the constraint stiffness needs to be reduced to 0.7 times the original value. Regions with strain rates below 0.1 per second are marked as stiffness preservation regions, maintaining the original stiffness. Regions between the two thresholds have their stiffness adjustment coefficients determined by linear interpolation, thus obtaining the constraint stiffness adjustment strategy. Based on the adjustment strategy, the Lagrange multiplier is calculated for each constraint in the viscoelastic constraint system. The Lagrange multiplier is calculated by dividing the negative value of the constraint deviation by the sum of the mass weights of all relevant vertices and the squares of the constraint gradient magnitudes, plus the regularization parameter, and then dividing by the square of the current adaptive stiffness. The regularization parameter is set to 0.01. After calculating the Lagrange multiplier, a position correction is applied to each vertex. The correction amount equals the vertex's mass weight multiplied by the Lagrange multiplier, then multiplied by the constraint gradient. The iteration process is executed a maximum of 10 times. After each iteration, the absolute value of the deviation of all constraints is checked to see if it is less than the convergence threshold of 0.0001 meters. For viscoelastic constraints, the projection correction amount is added to the conventional correction by multiplying the viscosity coefficient by the time step by the current position minus the position of the previous time step. This term introduces a velocity-dependent viscous damping effect through historical position information. After the iteration converges, the velocities of all vertices are updated. The velocity update is the projected position minus the predicted position divided by the time step. At the same time, the current strain value of each vertex is recorded and stored in the history buffer for use in the next time step.For example, in simulating liver tissue being pressed by surgical instruments, the initial vertical coordinate of a vertex is 20 mm, the external force is a negative 5 Newton force applied by the instrument, and the vertex mass is 0.002 kg. After position prediction, the vertical coordinate becomes 19.75 mm. The strain rates of the four edges connected to this vertex are calculated to be 1.0 s / s, 0.7 s / s, 0.8 s / s, and 0.7 s / s, respectively, all exceeding the 0.5 s / s threshold. Therefore, the stiffness is reduced to 0.7 times the original stiffness of 12000 Pascals, i.e., 8400 Pascals. After constrained projection iteration, the final position of the vertex is 19.82 mm, and the velocity is updated to -1.8 mm / s, forming a deformation response sequence that includes the position correction and velocity update value.

[0037] Based on the deformation response sequence, the viscosity coefficient is adjusted by temperature field coupling using the heat conduction equation to obtain thermodynamic correction parameters. The temperature field distribution within the simulation space is constructed based on the instrument contact position information recorded in the deformation response sequence. The initial temperature is set to 37 degrees Celsius to simulate normal human body temperature. When the surgical instrument contacts the tissue, the local heat source temperature is set according to the instrument type. The local temperature generated by the electrosurgical unit rises to 60-80 degrees Celsius, while that of the ultrasonic scalpel rises to 45-55 degrees Celsius. This yields spatial temperature data including the initial temperature and the local heat source temperature. The heat conduction equation is numerically solved based on the spatial temperature data. The heat conduction equation describes the temperature variation with time and space, where the partial derivative of temperature with respect to time is equal to the thermal conductivity multiplied by the Laplace operator of temperature. The thermal conductivity is set to 0.5 watts per meter per Kelvin. The partial differential equation is discretized and solved using the finite difference method. The temperature gradient of each grid node is calculated as the temperature difference between adjacent nodes divided by the node spacing, and the temperature diffusion rate is calculated as the thermal conductivity multiplied by the temperature gradient. This yields a transient temperature field including the temperature gradient and the temperature diffusion rate. The local temperature value in the transient temperature field is compared with the reference temperature of 37 degrees Celsius using an exponential function. The viscosity coefficient correction factor is calculated as the natural constant raised to the power of -0.03 times the temperature difference. This exponential function reflects the physical law that the viscosity of the tissue decreases with increasing temperature. When the temperature rises from 37 degrees Celsius to 60 degrees Celsius, the correction factor decreases to 0.54. The viscosity coefficient in the viscoelastic constraint system is adjusted multiplicatively based on the viscosity coefficient correction factor. The adjusted viscosity coefficient is equal to the original viscosity coefficient multiplied by the correction factor. Simultaneously, the elastic modulus of regions exceeding the damage threshold temperature is permanently reduced. The damage threshold is set to 70 degrees Celsius for more than 2 seconds. Regions meeting this condition are marked as solidification zones, and their first and second elastic constants are permanently reduced to 0.3 times their original values. This yields thermodynamic correction parameters that include the corrected viscosity coefficient and damage zone identification. For example, in simulating an electrosurgical unit cutting liver tissue, the temperature at the contact point is set to 75 degrees Celsius. Within a 5-millimeter radius around this point, the temperature diffuses through the heat conduction equation: 68 degrees Celsius at a distance of 1 mm, 52 degrees Celsius at a distance of 3 mm, and 42 degrees Celsius at a distance of 5 mm. The original viscosity coefficient at the contact point is 380 Pascals per second. The correction factor corresponding to the temperature of 75 degrees Celsius is 0.46, resulting in an adjusted viscosity coefficient of 175 Pascals per second. Since the temperature exceeds 70 degrees Celsius and lasts for 3 seconds, this area is marked as a damaged area. The original first elastic constant of 9500 Pascals is reduced to 2850 Pascals, and the second elastic constant of 7200 Pascals is reduced to 2160 Pascals, forming thermodynamic correction parameters for subsequent memory effect calculations.

[0038] The evolutionary state is obtained by calculating the memory effect of cumulative deformation using thermodynamic correction parameters combined with a historical strain buffer. A historical strain buffer is established for each tetrahedral mesh vertex based on the thermodynamic correction parameters. This buffer stores the strain values ​​of the most recent 100 time steps, with each time step interval being 0.01 seconds, corresponding to a total historical record of 1 second. The buffer uses a circular queue data structure; when a new strain value is added, the oldest strain value is overwritten, thus obtaining a deformation history sequence containing the strain values ​​of the most recent preset number of time steps. An exponential decay weighting is applied to the strain values ​​of each time step in the deformation history sequence. The weight is calculated as the negative time step index of the natural constant multiplied by the time step length divided by the power of the memory decay time constant. The memory decay time constant is set to 2 seconds, with strain values ​​further away from the current time having smaller weights. Then, a weighted summation operation is performed, multiplying the strain values ​​of all time steps by their corresponding weights and summing them to obtain the memory function value characterizing the degree of cumulative deformation. The memory function value is compared with a preset memory threshold of 0.8. For regions exceeding the threshold, the elastic modulus is reduced by 5% of the original elastic modulus. This process simulates the fatigue accumulation effect of the material, thus obtaining material degradation parameters that include fatigue accumulation indicators and the modulus reduction ratio. Fourier spectrum analysis is performed on the deformation history sequence based on the material degradation parameters. The Fourier transform converts the time-domain strain history signal into frequency components in the frequency domain. The presence of periodic loading characteristics is determined by identifying the dominant frequency peaks in the spectrum. When a dominant frequency is detected in the range of 0.5 to 2 Hz, it is determined to be cyclic loading. A fatigue acceleration factor is applied to the detected cyclic loading region, with the acceleration factor set to 1.5 times, increasing the fatigue accumulation rate in that region by 1.5 times. Simultaneously, the residual strain is calculated as 0.1 times the maximum strain value in the historical strain sequence. The strain recovery time constant is set to 10 seconds to describe the exponential decay of residual strain over time after the removal of external force, thus obtaining the evolution state including residual strain and strain recovery time constant. For example, simulating the periodic compression of heart tissue during surgery, the deformation history sequence of a certain vertex records the strain values ​​of the most recent 100 time steps. The strain values ​​of the first 50 time steps fluctuate between 0.05 and 0.15, and the strain values ​​of the last 50 time steps fluctuate between 0.02 and 0.08. An exponentially decaying weighted summation is performed on this sequence, with the weight of the most recent time step decreasing from 1.0 to the weight of the furthest time step at 0.37. The memory function value of the weighted summation is 0.92, exceeding the threshold of 0.8. Therefore, the first elastic constant of this vertex decreases by 5% from the original 11000 Pascals to 10450 Pascals. Fourier spectrum analysis detects a dominant frequency of 1.2 Hz, which is within the cyclic loading range. After applying a fatigue acceleration factor of 1.5 times, the fatigue accumulation rate of this vertex accelerates. The maximum value of the historical strain sequence is 0.18, corresponding to a residual strain of 0.018. When the external force is removed, this residual strain decays exponentially with a 10-second time constant, forming a complete evolutionary state describing the long-term mechanical behavior changes of the tissue.

[0039] In one specific embodiment, step S1 includes:

[0040] A three-parameter constitutive relation was constructed from the tissue mechanics data obtained by mechanical dynamic compression test to obtain the creep coefficient and relaxation coefficient, which include the first elastic constant, the second elastic constant and the viscosity coefficient.

[0041] The stress function of the Ogden model is weighted and fused based on the creep coefficient and relaxation coefficient to obtain a mixed stress relationship that includes the stretch ratio and nonlinear parameters.

[0042] By mapping the strain in the mixed stress relationship to the transition interval through the sigmoid weighting function, a mixed constitutive model containing the transition slope and critical strain value is obtained.

[0043] Based on the hybrid constitutive model, the least squares method was used to fit the measured force-deformation curves of lung, liver, kidney and heart tissues to obtain organ-specific parameter sets.

[0044] Specifically, when constructing the constitutive relationship of a three-parameter model for the tissue mechanics data obtained from the mechanical dynamic compression test, a cyclic load was first applied to the soft tissue sample using the HY-0580 mechanical dynamic compression test equipment. The equipment's sensors acquired force and displacement signals in real time, with a sampling frequency set to 1000 Hz. During the test, a compressive force ranging from 0 to 6.4 Newtons was applied to the lung tissue sample, and the corresponding displacement changes from 0 to 3.6 mm were recorded, resulting in the original force-displacement data sequence. The three-parameter model employs a combination structure of two springs and a damper. The first elastic constant represents the instantaneous elastic response in parallel with the first spring, the second elastic constant represents the delayed elastic response in series with the damper, and the viscosity coefficient represents the viscous drag characteristics of the damper. The calculation of the creep coefficient requires extracting the stage in the experimental data where the stress remains constant while the strain increases over time. A data segment is selected where the force value is stable at 4.2 Newtons for 2 seconds. The displacement is observed to gradually increase from an initial 2.3 mm to 2.8 mm. The displacement increment is divided by the initial length to obtain the strain increment. The relationship curve between the strain increment and time is the creep curve. Fitting this curve yields the creep coefficient, which includes an instantaneous term related to the reciprocal of the second elastic constant and a time-dependent term related to the first elastic constant and the viscosity coefficient. The calculation of the relaxation coefficient involves the process where the strain remains constant while the stress decays over time. In the experimental data, a data segment is selected where the displacement remains constant at 3.0 mm while the force value decays from an initial 6.0 Newtons to 4.5 Newtons. The stress decay curve is fitted using an exponential function. The fitting parameters include a combination of the two elastic constants and a time constant determined by the viscosity coefficient. This yields the creep coefficient and relaxation coefficient, which include the first elastic constant, the second elastic constant, and the viscosity coefficient. When weighting the stress function of the Ogden model based on the creep coefficient and relaxation coefficient, the Ogden model describes the nonlinear stress-strain relationship of the material under large deformation. Its stress function is based on the stretch ratio definition, which is the ratio of the deformed length to the initial length. The nonlinear parameter controls the degree of nonlinearity of stress variation with the stretch ratio. The core of weighting fusion is to determine the contribution ratio of the three-parameter model and the Ogden model in different strain ranges. The fusion process first converts the creep coefficient and relaxation coefficient into the stress prediction value of the three-parameter model at the current strain. This stress value is calculated by substituting the strain into the constitutive relation of the three-parameter model. At the same time, the current strain is converted into the stretch ratio, which is equal to 1 plus the strain value. The stretch ratio is then substituted into the stress function of the Ogden model to calculate the stress prediction value of the Ogden model. The weighting fusion adopts a weighted average method. The stress prediction value of the three-parameter model is multiplied by the weight factor 1 and the weight is subtracted. The stress prediction value of the Ogden model is multiplied by the weight factor weight. The two are added together to obtain the mixed stress value. This mixed stress value contains information on both the stretch ratio and the nonlinear parameter, forming a mixed stress relationship that includes the stretch ratio and the nonlinear parameter.

[0045] When mapping the strain in a mixed stress relationship to the transition interval using a sigmoid weighting function, the mathematical form of the sigmoid weighting function is 1 divided by 1 plus the negative transition slope of the natural constant multiplied by the strain minus the power of the critical strain value. The function's output value varies between 0 and 1. When the strain is less than the critical strain value, the function output is close to 0, indicating that the three-parameter model dominates; when the strain is greater than the critical strain value, the function output is close to 1, indicating that the Ogden model dominates. The transition slope controls the steepness of the transition from 0 to 1; a larger transition slope indicates a more abrupt transition, and a smaller transition slope indicates a smoother transition. The critical strain value defines the central location where the transition occurs. For lung tissue, the critical strain value is set to 0.15, and the transition slope is set to 5. When the strain increases from 0.1 to 0.15, the output of the sigmoid function increases from 0.07 to 0.5; when the strain increases from 0.15 to 0.2, the output of the sigmoid function increases from 0.5 to 0.93. The output value of the sigmoid function is used as a weight and substituted into the aforementioned weight fusion formula to replace the original fixed weights. This allows the mixed stress relationship to automatically adjust the contribution ratio of the three-parameter model and the Ogden model according to the current strain value. Thus, a mixed constitutive model containing the transition slope and the critical strain value is obtained. This model mainly reflects the viscoelastic characteristics of the three-parameter model at small strains and mainly reflects the nonlinear characteristics of the Ogden model at large strains, achieving a smooth transition near the critical strain value. When fitting the measured force-deformation curves of lung, liver, kidney, and heart tissues using the least squares method based on a hybrid constitutive model, the least squares method determines the optimal parameters by minimizing the sum of squared residuals between the theoretical predictions and the measured values. The fitting process first initializes the parameter set for each organ tissue, including the first elastic constant, the second elastic constant, the viscosity coefficient, the nonlinear parameter, the transition slope, and the critical strain value. Then, for each measured data point, its strain value is substituted into the hybrid constitutive model to calculate the theoretical stress value. The difference between the theoretical stress value and the measured stress value is squared, and the sum of the squared differences of all data points is used to obtain the sum of squared residuals. The gradient descent algorithm or the Levenberg-Marquardt algorithm is used to iteratively adjust the parameters. In each iteration, the partial derivatives of the sum of squared residuals with respect to each parameter are calculated, and the parameter values ​​are adjusted according to the direction of the partial derivatives to reduce the sum of squared residuals. The iteration process continues until the change in the sum of squared residuals is less than a set threshold of 0.0001 or the number of iterations exceeds 1000, thereby obtaining the organ-specific parameter set for each organ tissue.

[0046] In one specific embodiment, step S2 includes:

[0047] The organ-specific parameter set is input into the tetrahedral mesh generation tool for surface triangular mesh conversion, resulting in a set of tetrahedral elements containing vertex position information.

[0048] Define volume-preserving constraints for each tetrahedral element in the tetrahedral element set to obtain the first layer of constraint relationships based on the initial volume and vertex positions;

[0049] Based on the first layer of constraints, the edges of the tetrahedral element are stretched and constrained to obtain the second layer of constraints based on the initial edge length and mass weight.

[0050] Based on the first and second layer constraint relationships, the first elastic constant, the second elastic constant, and the viscosity coefficient in the organ-specific parameter set are embedded into the constraint stiffness coefficient to obtain a viscoelastic constraint system containing dynamic stiffness coefficient and relaxation time constant.

[0051] Specifically, when the organ-specific parameter set is input into the tetrahedral mesh generation tool for surface triangular mesh conversion, the Tetgen tetrahedral mesh generation tool first reads the surface triangular mesh file of the soft tissue. This file contains the three-dimensional coordinates and connection relationships of the surface vertices. Each triangular face is defined by three vertex indices. The first elastic constant in the organ-specific parameter set is used to determine the tissue stiffness level. The mesh density control parameter is adaptively set according to the stiffness level. For lung tissue, the first elastic constant is 3450 Pascals, which is low stiffness. The mesh density parameter is set to a larger value to control the generated tetrahedral size to be larger. For heart tissue, the first elastic constant is 11200 Pascals, which is high stiffness. The mesh density parameter is set to a smaller value to control the generated tetrahedral size to capture finer deformation details. The Tetgen tool executes the Delaunay tetrahedral meshing algorithm. This algorithm inserts internal nodes inside the surface triangular mesh and ensures that the circumsphere of any tetrahedron does not contain other nodes by using the circumsphere criterion. The meshing process extends layer by layer from the surface triangles inward. Each time a new node is inserted, it checks whether the circumsphere of the existing tetrahedron has been destroyed. If it has been destroyed, the relevant tetrahedron is deleted and reconnected to form a new tetrahedral element. After meshing, each tetrahedral element consists of four vertices, and each vertex stores three-dimensional spatial coordinates. The vertex position information includes three components: x-coordinate, y-coordinate, and z-coordinate. All tetrahedral elements constitute a tetrahedral element set, which is stored in an array structure. Each element of the array contains the indices of the four vertices pointing to the vertex coordinate array.

[0052] When defining volume-preserving constraints for each tetrahedral element in the tetrahedral element set, the volume-preserving constraints are based on the directed volume calculation of the tetrahedron. The calculation method is as follows: perform a cross product operation between the vector from the first vertex to the third vertex and the vector from the first vertex to the fourth vertex. The cross product result is a new vector. Perform a dot product operation between this new vector and the vector from the first vertex to the second vertex. The dot product result is a scalar. Multiply this scalar by one-sixth to obtain the directed volume of the tetrahedron. The sign of the directed volume depends on the vertex order; a positive value indicates a right-handed system, and a negative value indicates a left-handed system. The initial volume is calculated immediately after mesh generation. Substitute the initial coordinates of the four vertices of each tetrahedron into the above volume calculation formula to obtain the initial volume value. This initial volume value is stored in the constraint data structure as a reference for the constraint. The volume preservation constraint function is defined as the current volume minus the initial volume. The current volume is recalculated based on the real-time positions of the four vertices. When the tissue deforms, the vertex positions change, causing a deviation between the current volume and the initial volume. The output value of the constraint function is this deviation. A positive deviation indicates tetrahedral expansion, and a negative deviation indicates tetrahedral compression. The approximately incompressible nature of soft tissue requires this deviation to be as close to zero as possible. Therefore, the goal of the volume preservation constraint is to adjust the vertex positions to make the constraint function value approach zero. This yields the first layer of constraint relationship based on the initial volume and vertex positions. This constraint relationship uses the initial volume as the input parameter, the vertex positions as variables, and the constraint function value as the output.

[0053] When defining the stretch constraints for the edges of a tetrahedral element based on the first-level constraint relationships, each tetrahedral element contains six edges. These six edges connect the first vertex to the second vertex, the first vertex to the third vertex, the first vertex to the fourth vertex, the second vertex to the third vertex, the second vertex to the fourth vertex, and the third vertex to the fourth vertex, respectively. The definition of the stretch constraints relies on the vertex position information already determined in the first-level constraint relationships. An initial edge length is calculated for each edge, which is the Euclidean distance between the initial coordinates of the two endpoints. The calculation method is to take the square root of the sum of the squares of the differences in the x-coordinates, y-coordinates, and vertical coordinates of the two endpoints. The initial edge lengths of the six edges are calculated and stored separately. The allocation of mass weights is based on the number of tetrahedrons connected to each vertex. The entire tetrahedral cell set is traversed to count the number of times each vertex appears as a tetrahedron vertex. Vertices that appear more often are connected to more tetrahedrons and their mass weights are set to smaller values. Vertices that appear less often are assigned larger mass weights. Specifically, the mass weight is equal to the total mass divided by the number of tetrahedrons connected to that vertex. The total mass is calculated based on the tissue density and total volume. The soft tissue density is usually set to 1050 kg per cubic meter, which is close to the density of water. The tensile constraint function is defined as the current side length minus the initial side length. The current side length is recalculated based on the real-time positions of the two endpoints. When the tissue is deformed by force, the side length changes. The output value of the tensile constraint function is the change in side length. A positive value indicates tension, and a negative value indicates compression. The purpose of the tensile constraint is to limit excessive changes in the side length to maintain the continuity and shape of the material. When solving the constraint, the mass weights of the two endpoints need to be considered simultaneously. Vertices with larger mass weights move less during position correction, while vertices with smaller mass weights move more. This results in a second layer of constraint relationship based on the initial side length and mass weights. This constraint relationship uses the initial side length and mass weights as input parameters, the positions of the two endpoints of the side as variables, and the constraint function value as the output.

[0054] When the first elastic constant, second elastic constant, and viscosity coefficient of the organ-specific parameter set are embedded into the constraint stiffness coefficient based on the first and second layer constraint relationships, the constraint stiffness coefficient determines the strength of the constraint correction of the vertex position during the solution process. The larger the stiffness coefficient, the stronger the constraint and the larger the magnitude of the vertex position correction; the smaller the stiffness coefficient, the weaker the constraint and the smaller the magnitude of the vertex position correction. The formula for calculating the dynamic stiffness coefficient is the first elastic constant plus the second elastic constant multiplied by an exponential decay function. The independent variable of the exponential decay function is the negative current time divided by the relaxation time constant, which is the viscosity coefficient divided by the sum of the first and second elastic constants. This formula reflects the time-dependent characteristics of viscoelastic materials. At the initial moment when deformation begins, the current time is zero, the exponential term equals 1, and the dynamic stiffness coefficient equals the first elastic constant plus the second elastic constant. At this time, the constraint stiffness is at its maximum, reflecting the instantaneous elastic response of the material. As time progresses, the exponential term gradually decays, and the dynamic stiffness coefficient gradually decreases. When time approaches infinity, the exponential term approaches zero, and the dynamic stiffness coefficient approaches the first elastic constant. At this time, the constraint stiffness reaches its minimum steady-state value, reflecting the long-term elastic response of the material. The relaxation time constant controls the rate of stiffness decay; a large relaxation time constant results in slow decay, while a small relaxation time constant results in rapid decay. A larger viscosity coefficient and a larger relaxation time constant indicate stronger material viscosity, while a larger elastic constant and a smaller relaxation time constant indicate faster elastic recovery. Dynamic stiffness coefficients are assigned to the volume-preserving constraint of the first layer of constraints and the tensile constraint of the second layer of constraints. Each constraint has a different stiffness value at different times. The constraint solving algorithm calculates the Lagrange multiplier based on the current stiffness value and then calculates the vertex position correction, thus obtaining a viscoelastic constraint system containing dynamic stiffness coefficients and relaxation time constants. This constraint system transforms the mechanical parameters of the organ-specific parameter set into numerical parameters for constraint solving, establishing a bridge between the material constitutive model and the geometric constraint model.

[0055] In one specific embodiment, step S3 includes:

[0056] Based on a viscoelastic constraint system, the position of the vertices of the tetrahedral mesh is predicted to obtain the initial deformation state containing the predicted position, current velocity and external force information.

[0057] Based on the initial deformation state, the current strain and historical strain are differentially calculated to obtain the strain rate parameter characterizing the deformation rate;

[0058] The strain rate parameter is compared with the preset strain rate threshold to obtain a constraint stiffness adjustment strategy that includes a stiffness reduction region and a stiffness maintenance region.

[0059] Based on the constraint stiffness adjustment strategy, Lagrange multipliers are calculated and vertex positions are corrected for each constraint in the viscoelastic constraint system to obtain a deformation response sequence containing position correction and velocity update values.

[0060] Specifically, when predicting the position of vertices of a tetrahedral mesh based on a viscoelastic constraint system, the position prediction employs the explicit Euler integral method. This method calculates the state at the next moment based on the current state information. The calculation of the predicted position first reads the current position coordinates and current velocity vector of each vertex from the viscoelastic constraint system. The current position coordinates include three components: horizontal, vertical, and vertical coordinates. The current velocity vector also includes three directional components. Then, all external forces acting on the vertex are collected. The sources of external forces include gravity, the contact force of surgical instruments, and the interactive force applied by the user through a tactile device. Gravity is calculated as the vertex mass multiplied by the gravitational acceleration of 9.8 m / s², with the direction along the negative vertical axis. The contact force is determined by a collision detection algorithm to determine whether the instrument model intersects with the tissue surface mesh. When an intersection is detected, a force is applied at the intersection point in the direction of instrument movement. The repulsive force is proportional to the penetration depth, which is the shortest distance from the instrument surface to the tissue surface. The interaction force directly reads the force components in three directions from the force sensor of the tactile device. The resultant force of the vertex is obtained by vector addition of all external force vectors acting on the same vertex. The resultant force is divided by the vertex mass to obtain the acceleration vector. The formula for calculating the predicted position is the current position plus the current velocity multiplied by the time step plus the acceleration multiplied by the square of the time step. The time step is set to 0.01 seconds. This formula is applied to the three coordinate components to obtain the three coordinate values ​​of the predicted position. The predicted position, current velocity, and external force information together constitute the initial deformation state. The initial deformation state is stored in the vertex state array. Each element of the array contains nine values: the three-dimensional coordinates of the predicted position of a vertex, the three-dimensional vector of the current velocity, and the three-dimensional vector of the resultant external force.

[0061] When performing differential calculations on the current strain and historical strain based on the initial deformation state, the strain is defined as the relative degree of deformation of the material. For each edge of a tetrahedral element, the calculation of the current strain first extracts the predicted position coordinates of the two endpoints of the edge from the initial deformation state. The Euclidean distance between the predicted positions of the two endpoints is calculated as the current edge length. The change in edge length is obtained by subtracting the initial edge length from the current edge length. The change in edge length is divided by the initial edge length to obtain the current strain. A positive value indicates that the edge is stretched, a negative value indicates that the edge is compressed, and a zero value indicates that the edge length remains unchanged. The historical strain is read from the historical deformation record buffer. The historical deformation record buffer maintains a circular queue for each edge, and the queue stores the number of strains in the past several time steps. The queue length is set to 5, corresponding to the historical data of the past 0.05 seconds. The strain variable from the previous time step (0.01 seconds) is read from the queue as the historical strain variable. The difference calculation subtracts the historical strain variable from the current strain variable to obtain the strain increment. The strain increment reflects the strain change of the edge in the most recent time step. The strain increment is divided by the time step length of 0.01 seconds to obtain the strain rate parameter. The unit of the strain rate parameter is per second, which represents the deformation rate of the edge. A positive strain rate parameter indicates the tensile rate, and a negative strain rate indicates the compression rate. The larger the absolute value, the faster the deformation speed. A strain rate parameter value is calculated for each edge. The strain rate parameters of all edges constitute the strain rate parameter set. This set is used for subsequent threshold comparison and judgment.

[0062] When comparing the strain rate parameter with a preset strain rate threshold, the preset strain rate threshold includes two boundary values: the first threshold is a high-speed deformation threshold set to 0.5 seconds per second, and the second threshold is a low-speed deformation threshold set to 0.1 seconds per second. The comparison process iterates through each strain rate parameter value in the strain rate parameter set. First, the absolute value of the strain rate parameter is taken to eliminate the influence of the positive or negative sign. The absolute value is compared with the first threshold. If the absolute value is greater than 0.5 seconds per second, the region where the edge is located is determined to be a high-speed deformation region. In high-speed deformation regions, the constraint stiffness needs to be reduced to simulate the strain softening effect of the material. Strain softening refers to the phenomenon that the stiffness of the material decreases during rapid deformation. The edge is marked as a stiffness reduction region, and the stiffness adjustment coefficient is recorded as 0.7. The stiffness adjustment coefficient means that the adjusted stiffness is 0.7 times the original stiffness. If the absolute value is less than 0.1 seconds per second, it is determined to be... The region containing this edge is a low-speed deformation region, where the original constraint stiffness remains unchanged. This edge is marked as a stiffness-preserving region and the stiffness adjustment coefficient is recorded as 1.0. If the absolute value is between 0.1 seconds and 0.5 seconds, it is determined to be a medium-speed deformation region. The stiffness adjustment coefficient for the medium-speed deformation region is calculated by linear interpolation. The interpolation formula is: stiffness adjustment coefficient equals 1.0 minus the absolute value of strain rate minus 0.1, divided by 0.4, and then multiplied by 0.3. This interpolation formula ensures a smooth transition of the stiffness adjustment coefficient between 0.7 and 1.0. As the strain rate increases from 0.1 seconds to 0.5 seconds, the stiffness adjustment coefficient linearly decreases from 1.0 to 0.7. The marking information of all edges and the stiffness adjustment coefficient constitute the constraint stiffness adjustment strategy. This strategy is stored in the form of a data structure, which contains three fields: edge index, region label, and stiffness adjustment coefficient.

[0063] When calculating the Lagrange multiplier and correcting the vertex positions for each constraint in a viscoelastic constraint system based on a constraint stiffness adjustment strategy, the constraint solution employs a positional dynamics method. This method corrects the positions of violating vertices to satisfy the constraints through iterative projection. The maximum number of iterations is set to 10. Each iteration sequentially processes volume preservation constraints and tension constraints. For volume preservation constraints, the Lagrange multiplier calculation first calculates the constraint function value, i.e., the current volume minus the initial volume. Then, it calculates the gradient of the constraint function with respect to the four vertex positions. The gradient calculation involves partial derivative operations of vector cross products and dot products. The gradient of the second vertex is the vector of the first vertex minus the third vertex and the vector of the first vertex minus the fourth vertex. The gradient of the third vertex is the cross product of the vectors of the first and second vertices minus the second, and the vector of the first and fourth vertices minus the fourth, multiplied by one-sixth. The gradient of the fourth vertex is the cross product of the vectors of the first and second vertices minus the second, and the vector of the first and third vertices minus the third, multiplied by one-sixth. The gradient of the first vertex is the sum of the negative gradients of the other three vertices to ensure momentum conservation. The squared magnitude of the gradient vectors is calculated separately. The Lagrange multiplier equals the negative constraint function value divided by the mass weights of the four vertices multiplied by the sum of the squared magnitudes of the gradients, plus the regularization parameter divided by the square of the current constraint stiffness coefficient. The regularization parameter is set to 0.01 to prevent the denominator from being zero. The current constraint stiffness coefficient is determined according to the stiffness adjustment strategy. The adjustment coefficient is obtained by multiplying the original dynamic stiffness coefficient, which comes from the dynamic stiffness calculation of the viscoelastic constraint system. The adjusted stiffness coefficient reflects the influence of strain rate on material stiffness. After calculating the Lagrange multiplier, position corrections are applied to the four vertices. The position correction for each vertex is equal to the mass weight of that vertex multiplied by the Lagrange multiplier and then multiplied by the constraint gradient of that vertex. The position correction is a three-dimensional vector. The position correction is added to the current position of the vertex to obtain the corrected position. For tension constraints, the calculation of the Lagrange multiplier first calculates the constraint function value, which is the current side length minus the initial side length. The constraint gradient is a unit vector in the direction of the edge, and the gradient of the first endpoint of the edge is the distance from the first endpoint to the second endpoint. The unit vector is used, and the gradient at the second endpoint is the unit vector pointing from the second endpoint to the first endpoint. The Lagrange multiplier is equal to the negative constraint function value divided by the sum of the mass weights of the two endpoints and the squares of the gradient magnitudes, plus the regularization parameter divided by the square of the adjusted constraint stiffness coefficient. The position corrections at the two endpoints are calculated and added to the current position. Specifically for viscoelastic constraints, the position correction formula needs to add a viscous damping term. The viscous damping term is the viscosity coefficient divided by the time step, multiplied by the current position minus the position of the previous time step. This term introduces a velocity-dependent damping effect through historical position information, resulting in greater damping force during rapid deformation. During the iteration process, after each correction, it is checked whether the absolute value of the function values ​​of all constraints is less than the convergence threshold of 0.If all constraints meet the convergence condition, the iteration terminates early. If convergence is not achieved after the maximum number of iterations (10), the iteration is forcibly terminated, and the current correction result is used. After iteration, the positions of all vertices are corrected to approximately satisfy the constraints. The position correction is the corrected position minus the predicted position. The velocity update value is calculated by subtracting the true position from the corrected position and dividing by the time step size. The velocity update value reflects the actual velocity of the vertices under the constraints. The position correction and velocity update values ​​are stored together in the deformation response sequence. The deformation response sequence is time-series data, recording the position correction and velocity update values ​​of all vertices at each time step. This sequence is used to drive graphics rendering and force feedback output.

[0064] In one specific embodiment, the strain rate parameter characterizing the deformation rate is obtained by performing a difference calculation on the current strain and historical strain based on the initial deformation state, including:

[0065] Extract the vertex position information of the current time step from the initial deformation state, and measure the side length change of the tetrahedral element to obtain the strain at the current time step, including the stretching and compression directions.

[0066] Extract the position information of the corresponding vertex of the previous time step from the historical deformation record buffer, and measure and calculate the side length of the same tetrahedral element to obtain the historical moment strain including the initial side length and the deformed side length.

[0067] The strain change is obtained by calculating the element-by-element difference between the strain variable at the current moment and the strain variable at the historical moment according to the vertex index, which includes both positive and negative strain increments.

[0068] The strain rate parameter, which includes the instantaneous deformation rate and the average deformation rate, is obtained by performing a division operation based on the strain change and a preset time step, and by performing absolute value processing on the calculation result.

[0069] Specifically, the vertex position information for the current time step is extracted from the initial deformation state. These vertex positions include the x, y, and y coordinates in three-dimensional space, accurately reflecting the position of each vertex at the current moment. The position data of each vertex is stored in a vertex state array. Next, based on the three-dimensional coordinates of each vertex, the system calculates the edge length of each tetrahedral element. Specifically, the system obtains the current edge length by calculating the Euclidean distance connecting the two endpoints of the edge. For each edge of the tetrahedral element, the calculated current edge length is compared with the edge length in the initial deformation state to derive the strain of each edge. The strain reflects the degree of stretching or compression of the edge during the deformation process. A positive strain indicates stretching, a negative strain indicates compression, and a zero strain indicates no deformation.

[0070] The system extracts vertex position information from the previous time step from the historical deformation record buffer. This historical data reflects the deformation of the tetrahedral elements in the previous time step. In the historical deformation record buffer, the position and deformation record of each vertex are stored chronologically, and the system extracts the vertex position from the buffer for the previous time step. Based on these historical vertex position information, the system calculates the edge length at the historical moment. This edge length represents the length of each edge of the tetrahedral element in the previous time step. Then, the system calculates the strain at the historical moment by comparing it with the initial edge length. This strain represents the deformation of each edge in the previous time step, the degree of stretching or compression.

[0071] The strain at the current time step is compared with the strain at the historical time step using element-wise differences based on vertex indices. The difference between the current strain at the current time step and the strain at the historical time step is calculated for each edge, yielding the strain increment. The strain increment reflects the deformation change between the current and historical time steps. If the current strain is greater than the historical strain, it indicates that the edge has undergone further stretching; if the current strain is less than the historical strain, it indicates that the edge has experienced more compression. In this way, the system can obtain the specific amount of deformation change for each tetrahedral element between these two time steps, including the strain increment in both the stretching and compression directions.

[0072] The instantaneous deformation rate is calculated based on the strain change and a preset time step. The instantaneous deformation rate refers to the deformation rate per unit time, reflecting the deformation speed of the edge within the current time step. By dividing the strain increment by the time step, the system can calculate the deformation rate of each edge per second. Simultaneously, the system also calculates the average deformation rate by taking the strain increments over multiple time steps, which helps describe the deformation trend of the edge over a longer period. All calculated strain rate parameters are ultimately processed using absolute values ​​to ensure that the calculated deformation rate is positive regardless of the deformation direction. This processing allows the system to handle the deformation rate of all edges consistently, whether they are under tension or compression.

[0073] In one specific embodiment, step S4 includes:

[0074] Based on the instrument contact position information in the deformation response sequence, the temperature field distribution in the simulation space is constructed to obtain spatial temperature data including the initial temperature and the local heat source temperature.

[0075] Numerical solutions to the heat conduction equations based on space temperature data yield transient temperature fields containing temperature gradients and temperature diffusion rates.

[0076] The temperature-dependent viscosity coefficient correction factor is obtained by performing an exponential function operation between the local temperature value in the transient temperature field and the reference temperature value.

[0077] The viscosity coefficient in the viscoelastic constraint system is multiplicatively adjusted based on the viscosity coefficient correction factor, and the elastic modulus is permanently reduced in the region exceeding the damage threshold temperature, resulting in thermodynamic correction parameters that include the corrected viscosity coefficient and damage region identifiers.

[0078] Specifically, based on the instrument contact location information in the deformation response sequence, the system constructs the temperature field distribution within the simulation space. This contact location information indicates the specific location where the instrument contacts the tissue during the simulation. Whenever the instrument contacts the tissue, heat is generated at that contact point and in the surrounding area. The temperature distribution is not limited to the contact point but is also affected by heat conduction from the surrounding area. Therefore, the system uses this information to calculate the changes in the temperature field throughout the simulation space, obtaining spatial temperature data including the initial temperature and local heat source temperatures. The initial temperature is typically set to the tissue's ambient temperature, such as 37 degrees Celsius, while the local heat source temperature is adjusted according to the type and properties of different instruments; for example, devices such as electrosurgical units or ultrasonic scalpels generate locally high temperatures.

[0079] Based on spatial temperature data, the heat conduction equation is numerically solved to obtain the transient temperature field. The heat conduction equation describes the spatial and temporal variation of temperature. In this process, the system considers the propagation of heat within the tissue and the heat flow caused by temperature differences between different regions. Through numerical solutions, the system can obtain the temperature gradient (i.e., the rate of temperature change) and the temperature diffusion rate (i.e., the speed at which heat diffuses in space) at each moment. These two parameters are important components of the transient temperature field, helping the system accurately predict temperature changes. An exponential function operation is performed between the local temperature values ​​in the transient temperature field and the reference temperature value to obtain a temperature-dependent viscosity coefficient correction factor. The reference temperature is typically the normal body temperature of the tissue (e.g., 37 degrees Celsius). Increased temperature leads to decreased tissue viscosity; therefore, the correction factor calculation considers the effect of temperature changes on the viscosity coefficient. The higher the temperature, the smaller the correction factor, indicating a greater decrease in viscosity. Specifically, the correction factor is obtained by performing an exponential function operation on the difference between the current temperature and the reference temperature; this processing reflects the nonlinear characteristics of the effect of temperature on viscosity.

[0080] The viscosity coefficient in the viscoelastic constraint system is multiplicatively adjusted based on a viscosity coefficient correction factor. The adjusted viscosity coefficient more accurately reflects the viscoelastic properties of the tissue under high-temperature conditions. This adjustment ensures that the viscoelastic behavior of the tissue matches actual physical conditions during simulation. For regions exceeding the damage threshold temperature, the system permanently reduces the elastic modulus. The damage threshold temperature is typically set at 70 degrees Celsius; when the temperature exceeds this value and remains above it for a certain period, the elastic modulus of the tissue permanently decreases. This treatment reflects irreversible damage to the tissue under high-temperature conditions, thus yielding thermodynamic correction parameters that include the corrected viscosity coefficient and damage region identifiers.

[0081] In one specific embodiment, step S5 includes:

[0082] Based on thermodynamic correction parameters, a historical strain buffer is established for each tetrahedral mesh vertex to obtain a deformation history sequence containing the strain values ​​of the most recent preset number of time steps.

[0083] The strain values ​​at each time step in the deformation history sequence are assigned exponential decay weights and weighted summation is performed to obtain the memory function value that characterizes the degree of cumulative deformation.

[0084] The memory function value is compared with a preset memory threshold, and the elastic modulus of the vertex region exceeding the threshold is reduced to obtain material degradation parameters that include fatigue accumulation indicators and modulus reduction ratio.

[0085] Fourier spectrum analysis was performed on the deformation history sequence based on the material degradation parameters to identify periodic loading characteristics and apply fatigue acceleration factors to the regions where cyclic loading was detected, thus obtaining the evolution state including residual strain and strain recovery time constant.

[0086] Specifically, based on thermodynamic correction parameters, a historical strain buffer is established for each tetrahedral mesh vertex. This buffer stores the strain values ​​from the most recent preset number of time steps. Each vertex maintains a historical strain record, documenting its deformation over the past few time steps. This historical record allows the system to track the cumulative deformation of each vertex. This data is used for subsequent memory effect calculations, helping the system accurately reflect the long-term cumulative effect of deformation during simulation. An exponential decay weighting is applied to the strain values ​​at each time step in the deformation history sequence. Specifically, the system assigns a decay factor to each time step in the historical strain data. This factor is calculated based on the distance between time steps and a preset decay constant. Historical strain values ​​from further away have smaller weights, reflecting the smaller impact of these more distant strain values ​​on the current deformation state. Then, the system performs a weighted summation operation, multiplying the strain value of each time step by the corresponding decay factor and summing the results to obtain a memory function value characterizing the degree of cumulative deformation. This memory function value quantifies the overall effect of deformation over multiple time steps, representing the overall deformation of the material under long-term stress.

[0087] The system compares the memory function value with a preset memory threshold. If the memory function value exceeds the threshold, the system considers that significant fatigue accumulation has occurred in that region and performs an elastic modulus reduction process on that region. This means that the system reduces the elastic modulus of the region according to the degree of fatigue accumulation, simulating the phenomenon of materials gradually losing elasticity due to repeated stress. The reduction ratio of the elastic modulus is adjusted according to the magnitude of the memory function value, thus obtaining material degradation parameters that include fatigue accumulation indicators and the modulus reduction ratio. These parameters reflect the deterioration changes that occur in the microstructure under long-term repeated stress.

[0088] Based on material degradation parameters, the system performs Fourier spectrum analysis on the deformation history sequence. The Fourier transform converts the strain history signal in the time domain to the frequency domain, thereby identifying periodic loading characteristics. If the system detects periodic loading in certain regions, i.e., loading frequencies within a certain range (e.g., between 0.5 and 2 Hz), these regions are considered to be undergoing fatigue cycles. For these regions where cyclic loading is detected, the system applies a fatigue acceleration factor, typically set to 1.5 times. This acceleration factor speeds up the fatigue accumulation rate in that region, thus more realistically reflecting the long-term effects of high-frequency fatigue loading on the material. Through these processes, the system ultimately obtains the evolution state including residual strain and strain recovery time constant, indicating the residual deformation and recovery process of the microstructure after fatigue loading. The strain recovery time constant describes the rate of strain recovery over time after the removal of the external force, typically set to 10 seconds.

[0089] The above describes the soft physics simulation method based on viscoelastic constrained position dynamics in the embodiments of this application. The following describes the soft physics simulation system based on viscoelastic constrained position dynamics in the embodiments of this application. Please refer to [link / reference]. Figure 2 One embodiment of the software physics simulation system based on viscoelastic constrained position dynamics in this application includes:

[0090] The calibration module is used to perform viscoelastic constitutive calibration of soft tissue mechanical parameters by fusing the three-parameter model and the Ogden model, and obtain an organ-specific parameter set.

[0091] The construction module is used to construct a multi-level system of volume-preserving constraints and stretching constraints on the vertices of the tetrahedral mesh based on the organ-specific parameter set, thereby obtaining a viscoelastic constraint system.

[0092] The solution module is used to adaptively iteratively solve the constraint projection based on the viscoelastic constraint system through strain rate detection to obtain the deformation response sequence;

[0093] The adjustment module is used to adjust the viscosity coefficient by temperature field coupling according to the deformation response sequence using the heat conduction equation, so as to obtain thermodynamic correction parameters;

[0094] The calculation module is used to calculate the memory effect of cumulative deformation by combining the thermodynamic correction parameters with the historical strain buffer to obtain the evolution state.

[0095] above Figure 2 The software physics simulation system based on viscoelastic constraint position dynamics in this embodiment of the invention will be described in detail from the perspective of modular functional entities. The software physics simulation device based on viscoelastic constraint position dynamics in this embodiment of the invention will be described in detail from the perspective of hardware processing.

[0096] Reference Figure 3 This invention also provides a software physics simulation device based on viscoelastic constrained position dynamics. This device can be a server, and its internal structure can be as follows: Figure 3As shown, the software physics simulation device based on viscoelastic constrained position dynamics includes a processor, memory, display screen, input device, network interface, and database connected via a system bus. The processor in this computer design provides computational and control capabilities. The memory of the software physics simulation device based on viscoelastic constrained position dynamics includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system, computer programs, and database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The database of the software physics simulation device based on viscoelastic constrained position dynamics is used to store the data corresponding to this embodiment. The network interface of the software physics simulation device based on viscoelastic constrained position dynamics is used for communication with external terminals via a network connection. When the computer program is executed by the processor, it implements the above-described method.

[0097] Those skilled in the art will understand that Figure 3 The structure shown is merely a block diagram of a portion of the structure related to the present invention and does not constitute a limitation on the software physical simulation device based on viscoelastic constrained position dynamics on which the present invention is applied.

[0098] The present invention also provides a computer-readable storage medium, which can be a non-volatile computer-readable storage medium or a volatile computer-readable storage medium, wherein the computer-readable storage medium stores instructions that, when the instructions are executed on a computer, cause the computer to perform the steps of the software physical simulation method based on viscoelastic constrained position dynamics.

[0099] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0100] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a software physics simulation device based on viscoelastic constrained position dynamics (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0101] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A soft-body physical simulation method based on viscoelastic constrained position dynamics, characterized in that, The method includes: Step S1: Viscoelastic constitutive calibration of soft tissue mechanical parameters is performed by fusing the three-parameter model and the Ogden model to obtain an organ-specific parameter set; Step S2: Based on the organ-specific parameter set, construct a multi-level system of volume-preserving constraints and stretching constraints for the vertices of the tetrahedral mesh to obtain a viscoelastic constraint system; Step S3: Based on the viscoelastic constraint system, the constraint projection is adaptively iteratively solved by strain rate detection to obtain the deformation response sequence; Step S4: Adjust the viscosity coefficient using the thermal conduction equation based on the deformation response sequence to obtain thermodynamic correction parameters; including: constructing a temperature field distribution in the simulation space based on the instrument contact position information in the deformation response sequence to obtain spatial temperature data including the initial temperature and local heat source temperature; numerically solving the thermal conduction equation based on the spatial temperature data to obtain a transient temperature field including the temperature gradient and temperature diffusion rate; performing an exponential function operation on the local temperature value in the transient temperature field and the reference temperature value to obtain a temperature-related viscosity coefficient correction factor; multiplying the viscosity coefficient in the viscoelastic constraint system based on the viscosity coefficient correction factor, and permanently reducing the elastic modulus of regions exceeding the damage threshold temperature to obtain the thermodynamic correction parameters including the corrected viscosity coefficient and damage region identification. Step S5: Calculate the memory effect of cumulative deformation using the thermodynamic correction parameters combined with the historical strain buffer to obtain the evolution state. This includes: establishing a historical strain buffer for each tetrahedral mesh vertex based on the thermodynamic correction parameters to obtain a deformation history sequence containing the strain values ​​of the most recent preset number of time steps; assigning exponential decay weights to the strain values ​​of each time step in the deformation history sequence and performing a weighted summation operation to obtain a memory function value characterizing the degree of cumulative deformation; comparing the memory function value with a preset memory threshold and performing elastic modulus reduction processing on vertex regions exceeding the threshold to obtain material degradation parameters containing fatigue accumulation indicators and modulus reduction ratios; performing Fourier spectrum analysis on the deformation history sequence based on the material degradation parameters to identify periodic loading characteristics and applying fatigue acceleration factors to regions where cyclic loading is detected to obtain the evolution state containing residual strain and strain recovery time constant.

2. The soft physical simulation method based on viscoelastic constrained position dynamics according to claim 1, characterized in that, Step S1 includes: A three-parameter constitutive relation was constructed from the tissue mechanics data obtained by mechanical dynamic compression test to obtain the creep coefficient and relaxation coefficient, which include the first elastic constant, the second elastic constant and the viscosity coefficient. The stress function of the Ogden model is weighted and fused based on the creep coefficient and relaxation coefficient to obtain a mixed stress relationship that includes the stretch ratio and nonlinear parameters. The strain in the hybrid stress relationship is mapped to the transition interval through the sigmoid weighting function to obtain a hybrid constitutive model that includes the transition slope and the critical strain value. Based on the hybrid constitutive model, the least squares method was used to fit the measured force-deformation curves of lung, liver, kidney and heart tissues to obtain the organ-specific parameter set.

3. The soft physical simulation method based on viscoelastic constrained position dynamics according to claim 1, characterized in that, Step S2 includes: The organ-specific parameter set is input into a tetrahedral mesh generation tool for surface triangular mesh conversion, resulting in a tetrahedral cell set containing vertex position information. For each tetrahedral element in the tetrahedral element set, a volume-preserving constraint is defined to obtain the first layer of constraint relationships based on the initial volume and vertex position; Based on the first layer of constraint relationships, the edges of the tetrahedral elements are stretched and constrained to obtain the second layer of constraint relationships based on the initial edge length and mass weight. Based on the first layer of constraint relationships and the second layer of constraint relationships, the first elastic constant, the second elastic constant, and the viscosity coefficient in the organ-specific parameter set are embedded into the constraint stiffness coefficient to obtain the viscoelastic constraint system containing the dynamic stiffness coefficient and the relaxation time constant.

4. The soft physical simulation method based on viscoelastic constrained position dynamics according to claim 1, characterized in that, Step S3 includes: Based on the viscoelastic constraint system, the position of the tetrahedral mesh vertex is predicted to obtain the initial deformation state containing the predicted position, current velocity and external force information; Based on the initial deformation state, the current strain and historical strain are differentially calculated to obtain the strain rate parameter characterizing the deformation rate; The strain rate parameter is compared with a preset strain rate threshold to determine the constraint stiffness adjustment strategy, which includes a stiffness reduction region and a stiffness maintenance region. Based on the constraint stiffness adjustment strategy, Lagrange multipliers are calculated and vertex positions are corrected for each constraint in the viscoelastic constraint system to obtain the deformation response sequence containing position correction amounts and velocity update values.

5. The soft physical simulation method based on viscoelastic constrained position dynamics according to claim 4, characterized in that, The step of performing differential calculations on the current strain and historical strain based on the initial deformation state to obtain strain rate parameters characterizing the deformation rate includes: Extract the vertex position information of the current time step from the initial deformation state, and measure the change in the side length of the tetrahedral element to obtain the strain at the current time, which includes the stretching direction and the compression direction. Extract the position information of the corresponding vertex of the previous time step from the historical deformation record buffer, and measure and calculate the side length of the same tetrahedral element to obtain the historical moment strain including the initial side length and the deformed side length. The strain change is obtained by calculating the element-by-element difference between the strain at the current moment and the strain at the historical moment according to the vertex index, which includes both positive and negative strain increments. The strain rate parameter, which includes the instantaneous deformation rate and the average deformation rate, is obtained by performing a division operation based on the strain change and a preset time step, and by performing absolute value processing on the calculation result.

6. A soft physical simulation system based on viscoelastic constrained position dynamics, characterized in that, For implementing the soft physics simulation method based on viscoelastic constrained position dynamics as described in any one of claims 1-5, the soft physics simulation system based on viscoelastic constrained position dynamics comprises: The calibration module is used to perform viscoelastic constitutive calibration of soft tissue mechanical parameters by fusing the three-parameter model and the Ogden model, and obtain an organ-specific parameter set. The construction module is used to construct a multi-level system of volume-preserving constraints and stretching constraints on the vertices of the tetrahedral mesh based on the organ-specific parameter set, thereby obtaining a viscoelastic constraint system. The solution module is used to adaptively iteratively solve the constraint projection based on the viscoelastic constraint system through strain rate detection to obtain the deformation response sequence; The adjustment module is used to adjust the viscosity coefficient by temperature field coupling according to the deformation response sequence using the heat conduction equation, so as to obtain thermodynamic correction parameters; The calculation module is used to calculate the memory effect of cumulative deformation by combining the thermodynamic correction parameters with the historical strain buffer to obtain the evolution state.

7. A soft physical simulation device based on viscoelastic constrained position dynamics, characterized in that, Including memory and A processor, wherein the memory stores a computer program that can run on the processor, and the processor executes the computer program to implement the soft physical simulation method based on viscoelastic constrained position dynamics as described in any one of claims 1 to 5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is run by the processor, it causes the processor to execute the software physics simulation method based on viscoelastic constrained position dynamics as described in any one of claims 1 to 5.

Images