3D Cohesive Finite Element Minimum Invasive Surgery Simulation Based on Kelvin-Voigt Model

Minimally invasive surgery is an important technique used for cytopathological examination. Recently, multiple studies have been conducted on a three-dimensional (3D) puncture simulation model as it can reveal the internal deformation state of the tissue at the micro level. In this study, a viscoelastic constitutive equation suitable for muscle tissue was derived. Additionally, a method was developed to define the fracture characteristics of muscle tissue material during the simulation process. The fracture of the muscle tissue in contact with the puncture needle was simulated using the cohesive zone model and a 3D puncture finite element model was established to analyze the deformation of the muscle tissue. The stress nephogram and reaction force under different parameters were compared and analyzed to study the deformation of the biological soft tissue and guide the actual operation process and reduce pain.


Introduction
Minimally invasive surgery refers to the procedure of surgery guided by medical images, which is performed using a puncture probe or catheter inserted through the skin. This technique can be used to deliver a specific drug to the lesions in a tissue for targeted treatment or to remove the affected part of a tissue for subsequent pathological examination to aid in the diagnosis of the disease. Puncture surgery has multiple advantages such as a small wound, mild pain, rapid postoperative recovery, and lesser complications compared with conventional treatment methods. Therefore, minimally invasive surgery is widely used for medical procedures such as pleural puncture, lumbar puncture, and tumor ablation [1][2][3][4]. The continuous advancement of modern medicine and computer technology has led to the development of virtual surgery, which can simulate surgery procedures by combining simulation modeling and virtual reality technology. Virtual surgery has enabled surgeons to conduct multiple simulation surgeries without complications. Additionally, it has provided technical support to improve the success rate of surgery and to train interns. It can be observed that the effective means to perform virtual surgery is to study the constitutive modeling and simulation of biological soft tissue.
Over the past few decades, several scholars have performed the modeling and simulation of the puncture process, and have adopted complex models or finite element analysis to predict and improve the offset of the target during the puncture process. The compression tests of soft tissue have been used in multiple fields such as medical science, surgery, biology, and material science. Several constitutive models have been developed for soft biological tissues owing to their complex mechanical characteristics of nonlinear stress-strain behavior [5], large deformation [6], tissue anisotropy [7], and viscoelasticity [8]. Various nonlinear strain energy functions have been developed for the soft tissue, that are mostly imitate from hyperelastic models of Neo-Hooke model, Mooney-Rivlin model, Ogden model and Yeoh model [9]. Among all these constitutive models, the Yeoh model is more accurate for modelling the uniaxial extension of hyperelastic material [10]. In addition, some other constitutive models have been available, which are based on the above constitutive models and own special range of application. For example, Nolan et al. [11] proposed a robust anisotropic hyperelastic formulation for the modelling of soft tissue; Fan et al. [12] studied the hyperelastic constitutive model for EPDM soft PSD of dual-pulse motors, and Chanda et al. [13] developed a novel soft composite based material framework to model tissue anisotropy. Eto et al. [14] designed a puncture experiment to observe the mechanical properties of soft tissue, compared the changes of mechanical injury in vivo and in vitro, and analyzed its mechanical properties. Although a majority of the abovementioned strain energy functions have been used for the calculation and simulation of soft tissue as constitutive functions, soft tissues are different from hyperelastic materials owing to their anisotropy and viscoelasticity characteristics, whereby the above constitutive functions have limited application. Additionally, the constitutive models neglect the increase in the viscoelastic characteristic and anisotropy during calculation. Puncture mechanics modeling is also a hot issue in recent years. The classical puncture theoretical model was proposed by Okamura in 2004, and then gradually became mature after a series of development and improvement [15]. Tan et al. [16] studied the effect of vibration frequency on biopsy needle insertion force. Yu et al. [17] established a needle-shaped ultrathin piezoelectric microsystem to inject onto ventional biopsy needles and real-time measurements of variation in tissue modulus. Liu et al. [18] studied mechanics of tissue rupture during needle insertion in soft tissue and Liao et al. [19] studied the mechanics of cutting bones. These studies are based on experiments.
The simulation study of soft tissues has gradually emerged with the development of computer technology. Presently, several studies have been conducted using the finite element theory to develop two-dimensional (2D) or three-dimensional (3D) simulation models for soft tissue simulation. Jiang et al. [20] studied the deformation of the flexible needle during puncture and discussed the insertion mechanism under different puncture parameters. Noble et al. [21] studied the dissection of arterial layers using a stiff, planar, penetrating external body, and formulated a novel model of the process by means of cohesive zone formalism. Tagawa et al. [22] built 2D cervical cord model to study the material mechanical properties for simulating soft tissue under large deformation. Gao et al. [23] studied on the influence of a biocompatible hydrophilic needle surface coating on a puncture biopsy process for biomedical applications by building 2D simulation model. Some algorithm and methods of soft tissue FEM were proposed to simulate soft tissue more precise and efficient. Xie et al. [24] studied the simulation of soft tissue deformation in real-time using Kalman Filter finite element method, and Mendizabal et al. [25] proposed a simulation method of hyperelastic materials in real-time using deep learning. These methods are different explorations of the existing simulation directions.
The research of 3D simulation started late, but it has the advantages of high accuracy and can reflect the real puncture situation. Now more and more simulation focused on 3D simulation. For example, 3D injured blood-perfused skin was simulated by FEM [26], and under low-frequency ultrasonic tomography a simulation of soft tissue was studied to obtain its mechanical properties [27]. Lin et al. [28] proposed a computational approach to investigate optimal cutting speed configurations in rotational needle biopsy cutting soft tissue. Wang et al. [29] conducted a simulation of coupling process of flexible needle insertion into soft tissue based on ABAQUS and pointed out the importance of 3D simulation.
It was observed from the abovementioned studies that it is difficult to perform the simulation of puncture surgery. Therefore, in this study, a viscoelastic constitutive equation suitable for muscle tissue was derived. Additionally, a method was developed to define the fracture characteristics of muscle tissue material during the simulation process. The fracture of the muscle tissue in contact with the puncture needle was simulated using the dual cohesive force model. A 3D puncture finite element model was established to analyze the deformation of the muscle tissue. Stress nephograms and reaction forces under different parameters were compared and analyzed to investigate the deformation of biological soft tissue. This study can provide guidance during actual surgical procedures and reduce the pain experienced by patients.

Modeling
The pork tenderloin muscle tissue is similar to human muscle tissue in foundation compositions. Hence, it was selected as the experimental specimen. In microstructure of muscle tissue, myofibrils are rod-like contraction proteins, which traverse the entire length of the muscle fiber. Interstitial cells are endomysium, perimysium and connective tissue, which maintain certain hardness and flexibility function. Generally, myofibrils are seen as the combination of hyperelastic characteristic and viscoelasticity characteristic. Interstitial cells are defined by a single hyperelastic body owing to its hardness and flexibility. Taking the uncertainty in the distribution of two main types of tissue into consideration, in following discussion, the muscle tissue is assumed as hyperelastic of unidirectional homogeneity, which involves some viscoelastic characteristics. Figure 1 presents a set of steps for parameter identification to construct an appropriate constitutive model for the muscle tissue. Initially, the material density was calculated using the size and mass of the specimen. Subsequently, the experimental data of the time-force curve were averaged and separated into two parts, which provided the hyperelastic and viscoelastic parameters. The compression stage was used to calculate the hyperelastic parameters through the uniaxial compression test as shown by the red dotted line. The entire process of compression and relaxation was used to determine the viscoelastic parameters as shown by the blue dotted line. After that, the test data is matched with material parameters, which means the hyperelastic parameters are hyperelastic material constants (C ij ) and viscoelastic parameters are relaxation time ( τ i ) and Prony series ( g i ).
Meanwhile, the parameters of hyperelastic and viscoelastic are coefficient into the form of Kelvin-Voigt model. Finally, the procedure automatically iterates to find the optimal solution.

Experimental Setup and Specimen Preparation
A compression apparatus was designed to measure the force, displacement and time during the soft tissue deformation. An aluminum alloy section frame set was constructed to stabilize the equipment, as shown in Figure 2. A pressure cell with an operational range of 0.01-50 N and sample frequency of 1000 Hz was used to measure the force. The pressure cell was supported by a sliding rail, which was driven using a 47 mm stepping motor in the vertical direction. Additionally, a displacement transducer was installed at the frame to record the change in the real-time displacement.
In this study, the tenderloin muscle tissue specimen was obtained from porkets (1-1.5 years) with an average body weight of 30 kg. Figure 2 shows the tenderloin muscle tissue growing in the sirloin butt to the clod above the vertebra. These samples were acquired and tested within three hours post-mortem from the abattoir. The samples were cut into a cube of 30 mm × 30 mm × 30 mm using a scalpel. The direction of compression was parallel to that of the muscle fiber during the experiment, and the temperature and humidity were 18 ºC and 27%, respectively.
In the available literatures, many classical models, such as Neo-Hooke model (N-H model), Mooney-Rivlin model (M-R model), and Yeoh model to present the hyperelastic properties. Additionally, different from these models deduced by the strain energy function. In this paper, comparison of these models is firstly performed to determine the optimal model at current cases presented here, which will be used in further analysis.
To obtain the accurate hyperelastic and viscoelastic material parameters, parameter matching process is divided into pre-matching process and formal matching process. The pre-matching process is used to determine an appropriate hyperelastic model, which is suitable for compression process simulation, and the formal matching process defines final material parameters, which combines both hyperelastic and viscoelastic by Kelvin-Voigt model of soft tissue. To accomplish the pre-matching process, the test related to the strain-stress relationship is carried out. In order to obtain the stress-strain relationship during deformation, uniaxial tensile tests were carried out on fiber oriented and vertical fiber-oriented muscle tissue samples. Due to the soft texture of muscle tissue, too fast stretching speed would produce rapid tension of muscle tissue and destroy its internal connective tissue, making it lose its buffering effect on muscle cells, which is quite different from the normal deformation state. Therefore, the tensile speed is selected as 5 mm/min to ensure that the internal connective tissue will not be damaged rapidly due to rapid tension. The sample was fixed on the instrument and stretched until the tissue was completely broken. The data were saved, and the stretching experiment was repeated. Eight groups of experiments were carried out along the fiber direction and vertical fiber direction respectively, and the standard deviation analysis and average ultimate strength stress were solved. The force-displacement curve obtained from the tests is transformed into strainstress curve, which is used to estimate the hyperelastic parameters. Figure 3 presents the results of tension and compression stress-strain curves of the porcine tenderloin muscle tissues, where the load is along the fiber direction and perpendicular to the fiber direction respectively. The strain-stress curve shows an obvious rubber-like mechanical characteristic, whose linear region at the small engineering strain of the curve is followed by a non-linear response.
The average strain-stress curve is used to match the four representative hyperelastic models, and the involved material constants are calculated by the FEM software program package. Using different hyperelastic constitutive equations to fit the stress-strain curve, the results of Figure 4 are obtained.
Based on the accurate deformation model of soft tissue, the third linear stage is abandoned, and the hyperelastic model of linear nonlinear stage is fitted to obtain the material constant of hyperelastic stage. At the same time, to reduce the calculation time, the fitting results are shown in Table 1 under the condition of ensuring the accuracy of the first two stages.
Generally, the puncture needle is made of stainless steel. The flexible needle made of nickel titanium alloy has the characteristics of shape memory, biocompatibility, corrosion resistance and excellent elastic deformation ability, which meets the basic requirements of biomedical puncture needle. Therefore, in this paper, the material property of the puncture needle is defined as nickel titanium alloy. Since the puncture needle is defined as a flexible body in the analysis, the deformation caused by contact can be considered, and the deformation of the puncture needle in the soft tissue can be simulated. The specific material parameters are shown in Table 2.

Muscle Constitutive Model
For the moderate deformation, the hyperelastic models presented above cannot completely describe the deformation characteristics of the soft tissue because the viscoelastic characteristics are obvious and significant. Thus, it is necessary to include the viscoelastic properties so as to efficiently present the moderate deformation process of soft tissue. In this paper, as shown in Figure 5, the modified Kelvin-Voigt model is employed to deal with these terms. The Kelvin-Voigt material model consists of some dashpots with viscosity η and some springs of modulus E , which is changed by hyperelastic body having Cauchy stress tensor σ . It is assumed that the strain ( ε h ) experienced by hyperelastic body is the same as that experienced by the dashpot ( ε d ) and the applied external stress equals the sum of the stress within the hyperelastic body ( σ h ) and dashpot ele- Then the viscoelastic part of stress can be expressed by Prony series as following where g P i is the Prony series parameters, and τ i is the relaxation time. Therefore, from Eq. (3), the solution of the strain of soft tissue can be obtained as Due to the complex composition of muscle tissue, the above equation derivations are based on the assumptions of isotropic and uniformity, which means both myofibrils and interstitial cells are seen as the combination of hyperelastic and viscoelastic body.

Finite Element Modeling of Insertion
Human tissue is complex and arbitrary. Hence, it is difficult to achieve high-quality modeling, which affects the simulation accuracy. A complex model requires excessive calculation and time. Therefore, the explicit solver was used and the constitutive equation model was improved to establish a regular 3D model in ABAQUS 6.14 to simulate the puncture process of the biopsy needle. This model can conserve the calculation time and reduce the amount of calculation, and the material property parameters defined by the experimental data ensure the calculation accuracy. Figure 6 shows the definition of geometric data for a 3D model. Figure 6(a) represents the geometric modeling of the biopsy needle, which is a hollow cylinder with angle of inclination. In detail, the length of the needle tube is 50 mm, the initial definition of the outer diameter of the needle tube is 2.77 mm, the inner diameter is 2.16 mm, and the bevel angle of the needle tip is 12° as shown in Figure 6(c). Figure 6(b) presented the geometric modeling of biological soft tissue, in which the size of the tissue model is 50 mm (height) × 50 mm (width) × 35 mm (thick). According to existing research, muscle soft tissue is a relatively soft material, and its equivalent Young's modulus is about 20 GPa [30]. Therefore, the hyperelastic model can better describe its deformation behavior than the linear elastic model, especially when the deformation is large. According to the fitting results of the parameters using the constitutive equation, the improved constitutive model is selected because it is closest to the deformation behavior of muscle tissue and can provide stability for all strain ranges. For the incompressibility of soft tissue, the Poisson's ratio of tissue model is set to 0.495. The puncture needle was initially placed 2 mm away from the cut plane of the tissue phantom, and at the same time, the movement of translation (insertion) was carried out to cut into the tissue phantom. The insertion speed is fixed, set as 2 mm / s, SPR is 0-5.

Definition of Contact and Boundary Conditions
The general contact algorithm was set to allow objects to contact each other. In addition, the contact of each critical area was further defined by the surface contact. The contact between the needle wall and tissue phantom had a friction coefficient of 0.3, while the contact between the needle tip surface and new surface due to fracture was assumed to be frictionless to prevent excessive deformation of the element. The tissue model was fixed on the distal surface to limit its axial movement in the puncture needle tube. The joint surface was inserted into both sides of the predetermined crack path, which is the direction of the thin cylindrical sleeve along the axial direction, as shown in Figure 7. A custom constraint was constructed to dynamically determine the X and Y displacements of the cohesive node at each increment to ensure that the needle tip cut along the cohesive surface. A slight adjustment was performed when a cohesive node was observed outside the predetermined crack path. If fixed constraints are used, it might result in excessive response to the human force of the system. Therefore, the mesh on the tissue model was specially arranged, that is, the material inside the puncture needle was meshed to ensure numerical accuracy and efficiency.
The central area was represented by a 6 mm × 6 mm red dashed square with a mesh size set to 0.25 mm by the thickness of the tissue model. The mesh size was gradually increased to 3 mm from the edge of the central area to the boundary of the tissue model. This arrangement resulted in 147120 elements. A global grid of 0.68 mm was used to generate 39515 elements. The unit type of the two parts was C3D8R. The bilinear cohesive force model was used to simulate the 3D fracture of the tissue in contact with the puncture needle.

Definition of Fracture Criterion
It is crucial to select an appropriate fracture model to accurately simulate the deformation states of the muscle and soft tissue. The cohesive force model was used to effectively describe the mechanical behavior of an interface. It was originally proposed to solve the elasticplastic fracture problem of ductile metal materials. In this model, one part of the crack surface is not subject to stress, while the other part is subject to stress, which is called "cohesive force". Therefore, the mechanical behavior of the complex deformation region near the crack tip can be described by the equivalent relationship between the opening force acting on the crack surface and the opening displacement. The cohesive force model was introduced into the deformation simulation of the puncture needle to effectively simulate the crack propagation behavior at the contact between the puncture needle and tissue. The bilinear tension displacement rule is a simple and effective cohesive force rule. The bilinear finite element model was used to provide theoretical guidance for the puncture process. There are two types of two-line cohesion models. The first is based on the description of traction separation, whereas the second is based on the description of continuum. Among them, the method based on the description of trace separation is widely used. Considering the methods based on the description of traction separation, the most commonly used constitutive model is the bilinear constitutive model shown in Figure 8.
It provides the linear elastic section before the material reaches the strength limit and the linear decreasing softening stage after the material reaches the strength limit. Note that the ordinate is stress and the abscissa is displacement, so the slope of the linear elastic segment represents the stiffness of the cohesive element. The area under the curve is the energy release rate when the material ruptures. Therefore, when defining the mechanical properties of cohesive zone, it is necessary to determine the specific shape of the above constitutive model, including stiffness, ultimate strength, critical fracture energy release rate, or the displacement of the element at the final failure. The commonly used definition method is to determine the cohesive zone constitutive model by giving the first three of the above parameters. Cohesive zone element can be understood as a kind of quasi two-dimensional element. It can be regarded as two faces separated by a thickness, which are connected with other solid elements respectively. The interface might be significantly thin (some interfaces have a thickness of 0.001 mm or lower) in practical problems. In the case of some problems, the unit contact of the cohesive zone might have a thickness of 0 mm (such as the crack problem), while the overall model might be large. If these two thicknesses are not introduced, a small interface will need to be created in a large model, which is significantly difficult. The significantly small interface thickness can be replaced by a finite thickness while modeling including these two thicknesses.
The cohesive zone element only considers out of plane forces, including normal stress and shear stress in XZ and YZ directions. The governing equation is as follows.
where, σ is the normal stress value,σ max is the maximum normal stress, and the corresponding crack interface opening displacement is δ 0 n . When the stress begins to decrease to zero after reaching its maximum value, the crack cracking is completed, and the corresponding displacement value is the final cracking displacement value δ f n . The critical value of fracture energy for each item is φ c n . The calculation equation is as follows: The two-line cohesive force model is simple and effective, which can be effectively calculated using the finite element method with high accuracy.
Interstitial cells are composed of myointima, fascicular membrane and connective tissue, which maintain a certain hardness and elastic function of muscle tissue. The uniaxial tensile tests are completed by QX-W200 electronic universal testing machine to determine the stiffness and ultimate strength of muscle and soft tissue, and the specific tensile test is in Section 2.1.
To study the critical fracture energy release rate in tissue and the cutting performance of muscle tissue under fiber cutting force, we built a cutting experimental platform and carried out five groups of cutting experiments from fiber to muscle tissue. Fracture toughness is one of the important parameters to measure energy release rate, and it is also an important mechanical parameter to evaluate muscle tissue. When the scalpel moves for a period of displacement after contacting the muscle tissue, the fracture toughness can be defined as: where, F is the cutting force, dɅ is the internal energy stored in the deformed material, dΓ represents the dissipated energy in the plastic flow process, and JdA is the newly generated cutting surface of the material in the cutting process. The fracture energy release (Je) rate is obtained by recording different experimental conditions, as shown in the Table 3.

Finite Element Modeling of Insertion
The results of the model were shown in detail through ABAQUS simulation. Figure 9 shows the change in stress at the outer surface during the puncture process. It can be observed from Figure 7 that the change in stress is mainly concentrated near the contact between the biopsy needle and tissue during the puncture process. According to the scale, the change was approximately in the range of 50-700 MPa. The change in stress of other parts of the soft tissue was not distinct and the stress at the edge remained constant, which is consistent with the actual situation of puncture operation. It can be observed from Figure 7(c) that although the puncture needle was subjected to stress changes, it was insignificant. Considering the top view, it can be observed that the change in stress had a circular distribution; the stress concentration (red color) was at the contact between the tissue and interior, and the stress in the interior was greater than that at the edge of the tissue due to the deformation extrusion.

Effects of the Needle Diameter and Experiment Verification
In order to study the effect of different inserting parameters on inserting force and stress change during puncture, simulation is carried out under different diameters, velocities, needle shapes and rotation angles. The specific parameter settings are shown in Table 4. To verify the correctness and accuracy of the simulation results, design the verification experiment as shown in Figure 10. The vertical inserting experimental platform was set up, the needles used in the experiment were 7#, 9# and 12# peritoneal puncture needle, and the shape of the tip was single-bevel. The driving force of the puncture test is provided by the precision motion platform manufactured by Beijing Paidiwei Instrument co., LTD., and the repeatability positioning accuracy is 5 μm. A 6-DOF force/torque sensor was used to measure the forces and torques when the needle was inserted into the tissue. The force sensor used was an K6D27 50 N/1 N•m with a range of 50 N and a theoretical accuracy of 0.25 N•m. The analogue signal from the force/torque sensor was converted into a digital signal by DAQ Card with a sample rate of 60 Hz. Inserting experiments were performed at 0.5 mm/s, and compared with the simulation results. All experiments with the same parameters were repeated five times to gather force data. Figure 11 shows the change in the stress nephogram inside the puncture needle during the puncture process compared with that of the experiment. The three processes represent the state when the needle contacts the soft tissue, when it enters 3 mm inside the tissue, and when it enters 8 mm inside the tissue.
It can be observed from Figure 11 that in the whole process of puncture, both the experimental results and the simulation results have the same deformation trend outside the puncture needle. The deformation of soft tissue inside the puncture needle is more severe than that outside the puncture needle, and because of the symmetry of the puncture needle geometry, the change of stress nephogram also presents a symmetrical distribution according to the central axis. The bulge phenomenon will appear outside and inside the puncture needle, which is determined by the mechanical properties of soft tissue hyperelasticity, and because of the  small space inside the needle tube, the internal tissue will be squeezed by the inner wall of the whole needle tube to the center, and the deformation will be more severe. When the puncture reaches a certain depth, the friction between the tissue and the inner wall of the needle tube will separate the internal biological soft tissue from the original whole, so as to achieve the purpose of biopsy sampling. The force exerted on the muscle tissue is difficult to obtain, and the knife reaction force is consistent with the force on muscle tissue. Therefore, the force exerted on the muscle tissue was determined by analyzing the reaction force of the biopsy needle. The reaction force at the reference point of the puncture needle end was used to compare the stress state of muscle tissue under different parameter variables. Figure 12 shows a comparison between the puncture and reaction forces at the reference point of the puncture needle during the puncture process. The red line represents the change in stress in the simulation result, and the blue line represents the change at the needle end measured by the force sensor. According to the analysis in Figure 12, the maximum range is about 4 N within the fluctuation range of puncture force, and the force collected in the experiment is slightly higher than the simulation result. The two curves analysis demonstrated that the stress of the soft tissue gradually increased with the progression of the puncture process. The rebound phenomenon of soft tissue occurred at point A because the tissue did not break after the puncture needle entered the soft tissue due to the stiffness force of the biological valve. Therefore, the force at this point gradually increased. The stiffness force was absent when the puncture needle Only the puncture force and friction between the puncture needle and inner and outer surfaces of the tissue were critical. Therefore, a sharp decline was observed. It can be concluded by an analysis of the puncture process that this simulation model is relatively consistent with the real operation situation and is reasonable. Figure 13 shows the stress nephogram under different needle diameters. Figure 13 represents needles with an outer diameter of 2, 1.77 and 1.3 mm, respectively. The three stress nephograms demonstrated a similar puncture depth. It can be observed from the stress nephogram that the stress variation and deformation characteristics of the soft tissue with different needle diameters were essentially the same, and the stress concentration was mainly distributed at the contact position between the outer wall of the needle tube and soft tissue. However, a difference was observed in the order of magnitude of stress. The stress concentration decreased with an increase in the needle diameter. A decrease in the needle diameter resulted in an increase in the range of stress change, and the stress concentration phenomenon increased.

Effect of Needle Diameter on Tissue Deformation and Stress
Considering the deformation state by observing Figure 13(a) and (b) in the puncture state, the cutting surface of the muscle tissue after puncture was smooth, and the upper chip part was relatively complete. In the case of Figure 13(c), it was observed that the cutting surface exhibited an undulation state, and the stress at the undulating position was significantly higher than that at the same position in other models.
The above analysis demonstrated that using a puncture needle within a certain range with a small diameter can cut the muscle tissue during the puncture process with a smooth cutting surface, complete chip, and without tearing. However, considerable damage was caused to the muscle tissue when the diameter of the puncture needle was significantly large. Therefore, it is crucial to select a smaller diameter puncture needle under a certain puncture speed to reduce tissue damage. Figure 14 shows the change in the reaction force during the puncture process when different diameter puncture needles were used. The green, blue and dotted lines represent biopsy needles with a diameter of 2, 1.77 and 1.3 mm, respectively.
It can be observed from Figure 14 that an increase in the diameter of the puncture needle increased the reaction force. This was because the contact surface between the tissue and the puncture needle increased when the diameter of the puncture needle was increased, and the sliding friction increased as well. The change trend of the three curves was similar. The reaction force increased slowly, continuously and smoothly with an increase in the puncture displacement, followed by a decrease in the reaction force. The reaction force fluctuated until the end of the puncture process. The points A, B and C in Figure 14 represent the critical state of tissue fracture when the stiffness force was the largest. The tissue valve easily  fractured in the case of needles with a small diameter. This phenomenon can be explained by the energy change of the system. The smaller diameter of the puncture needle increases the stress concentration, which enables the rapid accumulation of the muscle tissue and internal energy. This results in the advancement of the fracture point, that is, the displacement point of the first obvious reduction of reaction force decreases. Therefore, selecting the appropriate diameter of the needle can effectively control the tissue deformation and the change in the force exerted on the tissue. Figure 15 shows the stress change nephogram of the muscle and soft tissue for the same puncture needle depth and different puncture speeds. Figure 13 represents the outer diameter of the puncture needle when the speed was 15, 10, 5 and 2 mm/s, respectively. It can be observed from the stress nephogram that the stress change in the muscle and soft tissue differed at different puncture speeds. The deformation of the soft tissue on the inner side of the needle wall decreased with a decrease in the speed, and the degree of fit between the tissue and inner wall of the needle increased. In the case of the outer side of the puncture needle wall, the stress change and the degree of fit between the tissue and outer wall increased with a decrease in the velocity.

Effect of Puncture Speed on Tissue Deformation and Stress
It can be concluded from the above analysis that maintaining a higher puncture speed within a certain range can reduce the effect of tissue damage during the process of puncture, which is conducive for wound healing. Figure 16 shows the change of reaction force during the puncture process with different diameters of puncture needle. The pink dotted line represents the puncture speed of 15 mm/s, the blue dotted line represents the puncture speed of 10 mm/s, the green solid line represents the puncture speed of 5 mm/s, and the black dotted line represents the puncture speed of 2 mm/s. It can be observed from Figure 16 that the effect of puncture speed on the tissue force was minor compared with that of the diameter of the puncture needle, and the variation range of force was approximately the same for different puncture speeds. The fluctuation of the velocity curve was relatively higher compared with that of the needle diameter. This was because with the increase in the cutting displacement, the muscle tissue that accumulated deformation finally ruptured and released its internal energy, resulting in an instantaneous decrease in the reaction force. The reaction force fluctuated with a further increase in the cutting displacement, which is consistent with the results of the cutting process in the experiment. The muscle tissue was continuously cut with an increase in the displacement. The continuous fracture and tearing of the muscle tissue gradually destroyed the internal fiber network. With the continuous release of the internal energy of muscle tissue, the muscle tissue at the front end of cutting increasingly accumulated its internal energy until rupture. The continuous accumulation at the front end of the muscle tissue and the rupture of the muscle tissue at the contact point caused the repeated increase and decrease in the internal energy, leading to the continuous fluctuation of the reaction force of the puncture needle. Therefore, selecting a higher puncture speed can effectively control the soft tissue deformation during puncture biopsy. However, the change in the force exerted on the tissue was relatively mild when the speed was low, which can reduce the pain.

Optimization of Puncture Process for Pain Relief
It is crucial that surgeons minimize the pain experienced by patients during the process of puncture operation. In the field of mechanical engineering, pressure and stress are used to define pain. However, the extent of pain is unclear and this ambiguity is due to the diversity of different substances. In this section, pain indicators were proposed to compare the severity of pain. The following conditions should be fulfilled before deriving a pain index. First, the index should be applicable to different muscle tissues. Second, it should involve the puncture process. Third, it should be intuitive and clear. Based on the above analysis, we defined the pain index (PDI). It is generally believed that when the muscle stress increases to 0.3 gradually, the pressure begins to increase rapidly and shows obvious nonlinear characteristics. Here, it means that when the muscle is compressed to a certain extent by the puncture needle, the stress of 0.3 strain will be determined as the sensory basic pressure. The stress at a strain of 0.3 is defined as the basic stress, which is defined by the σ * . Then, the real-time stress during puncture is expressed as where R is the outer diameter of the puncture needle and r is the inner diameter of the puncture needle. Therefore, the definition of PDI is given, where PDI is the multiple of pain basis, which is used to determine the degree of pain. Moreover, according to the simulation results in this study, the degree of pain under different parameters can be obtained.
The simulation of muscle tissue puncture under different puncture diameters, puncture depth and puncture speeds, and the comparative analysis of stress nephogram and reaction force under different parameters are carried out to study the deformation of biological soft tissue in order to reduce pain. The deformation and stress of the contact between puncture needle and soft tissue directly affect the patient's pain and postoperative healing.
Through the reasonable 3D finite element model of puncture, from the perspective of muscle tissue stress and deformation, the puncture parameters have an important impact on empathy. For example, if the puncture speed is fast, the contact time between the needle and tissue can be shortened. However, an excessively high puncture speed can easily rupture the tissue. Therefore, the reasonable puncture parameters optimized according to the simulation model can effectively control the deformation degree of the muscle tissue and the force, which is helpful for wound healing and to reduce the pain.

Conclusions
Research on the inserting performance of muscle soft tissue is beneficial for surgical design, surgical simulation, tissue additive manufacturing, and bionic tissue design. The inserting performance of tissue refers to the mechanical properties of the inserting process after the surgical needle comes in contact with the soft tissue and subsequently undergoes extrusion deformation until the inserting leads to the destruction of the soft tissue. The cutting performance in the experimental study has certain limitations due to the irregular shape of the muscle tissue and differences between individuals. Therefore, it is crucial to establish accurate cutting simulation models to investigate the inserting performance of soft tissue. In this study, the following conclusions were obtained from the soft tissue inserting simulation.
(1) A viscoelastic constitutive equation suitable for muscle tissue was derived. The viscoelastic characteristic equations were compared and selected by considering the structural characteristics of the muscle tissue, and the Kelvin-Voigt model was improved. The viscoelastic characteristic was combined with the hyperelastic constitutive equation to obtain the constitutive equation suitable for muscle tissue. (2) Considering the tissue deformation of puncture biopsy, a method was developed to define the fracture characteristics of muscle tissue during the simulation process. The fracture of muscle tissue and needle was simulated using the two-line cohesion model. A 3D puncture finite element model was established to analyze the deformation of the muscle tissue. Additionally, the model was experimentally verified. (3) According to the actual shape of the puncture needle and muscle tissue, the 3D puncture finite element model was defined. Additionally, the contact characteristics, fracture criteria, and mesh division between the tool and muscle tissue were determined. The muscle tissue puncture simulation was performed under different puncture diameters, depths, and speeds. The stress nephograms were compared and analyzed to determine the reaction forces under different parameters. According to the simulation results, selecting a needle with a smaller diameter and a higher puncture speed can effectively control the soft tissue deformation during the puncture biopsy.