Three-dimensional Modeling and Simulation of Muscle Tissue Puncture Process

Needle biopsy is an essential part of modern clinical medicine. The puncture accuracy and sampling success rate of puncture surgery can be effectively improved through virtual surgery. There are few three-dimensional puncture (3D) models, which have little significance for surgical guidance under complicated conditions and restrict the development of virtual surgery. In this paper, a 3D simulation of the muscle tissue puncture process is studied. Firstly, the mechanical properties of muscle tissue are measured. The Mooney-Rivlin (M-R) model is selected by considering the fitting accuracy and calculation speed. Subsequently, an accurate 3D dynamic puncture model is established. The failure criterion is used to define the breaking characteristics of the muscle, and the bilinear cohesion model defines the breaking process. Experiments with different puncture speeds are carried out through the built in vitro puncture platform. The experimental results are compared with the simulation results. The experimental and simulated reaction force curves are highly consistent, which verifies the accuracy of the model. Finally, the model under different parameters is studied. The simulation results of varying puncture depths and puncture speeds are analyzed. The 3D puncture model can provide more accurate model support for virtual surgery and help improve the success rate of puncture surgery.


Introduction
Minimally invasive puncture is playing an increasingly important role in clinical treatment. Compared with traditional medicine, puncture biopsy has the advantages of minor trauma and faster recovery [1][2][3]. Virtual surgery can realize the simulation of the surgical process by combining simulation modeling and virtual reality surgery. With the continuous development of modern medicine and computer technology, virtual surgery technology has increasingly applied to puncture technology. Virtual surgery is beneficial for surgeons to increase the experience of surgical operations, which provides technical support for trainee training, and improve the success rate and efficiency of surgical procedures [4]. It is the key to establishing a precise surgery simulation model in virtual surgery. An accurate model can improve the accuracy of the contact deformation between the instrument and the organ during the operation, adapt to the operating environment in a variety of circumstances, then enhance the quality of surgery, reduce the biological soft tissue trauma and achieve accurate localization.
In recent years, many scholars have devoted themselves to the research on the modeling and simulation of the puncture process. The constitutive model of soft materials is crucial for understanding the complex behavior of biological soft tissues and organs [5]. Gao et al. [6] studied the effect of the coating on the puncture process; Eto et al. [7] used the puncture device to conduct puncture experiments to observe the mechanical properties changes of in-vivo and in-vitro tissues in the mechanical properties measurement experiment. Jerg et al. [8] proposed a new meshless needle insertion simulation method. But this method cannot capture the deformation component perpendicular to the needle bar. Wang et al. studied the tangential cutting characteristics of the porcine loin tissue surface [9], revealing the relationship between stress and pain. Jiang et al. studied the changes of the needle during puncture and discussed the insertion under different puncture parameters mechanism [10]. Xu et al. [11] established a deflection model based on mechanical principles when the needle is inserted into soft tissue to predict the bending of the needle tip inserted into the soft tissue. Dong et al. [12] established a fiber bending stage around the needle. The mechanical model of the puncture needle, the functional relationship between the curve parameters of the needle tip, the structural size of the puncture needle, and the elastic modulus of the material were studied. Gao et al. [13] proposed a needle-tissue coupling model based on the improved local constraint method and analyzed the tissue the change of force and displacement of each element node in the unit. Liu et al. [14] established a cutting force model based on the solid mechanics theory of unidirectional fiber-reinforced composites to predict the stiffness force and cutting force when a needle was inserted into a transversely isotropic tissue. Gzaiel et al. [15,16] established a cutting force model. A finite element model of bladed-cut neoprene was used to analyze shear stress, soft tissue stress state at different speeds and material failure characteristics, which determines the stress state during blade insertion. Jushiddi et al. [17] developed a three-dimensional FE model that simulated the interaction of biopsy needles when it's inserted into soft tissue. Liu et al. [18] established a mechanics-based simulation model of needle deflection in laterally isotropic tissue to predict and compensate for the deflection error of acupuncture. Zhang et al. [19] introduced a VS-based puncture model composed of three layers of soft tissue, which simulated the surface recovery phenomenon when the tissue surface was damaged. Fu et al. [20] established the cantilever beam model of flexible needle and analyzed the influence of the puncture process on tissue deformation. Cui et al. [21] proposed a liver model based on the improved mass-spring model to simulate the deformation of the liver during puncture. Zhang et al. [22] established the relationship between the force model and the force during the puncture process. It proves that the power of the puncture needle with different materials will be different, which results in calculation errors. Lin et al. [23] established a model to study the puncture force during puncture. It's indicated that the cutting speed (rotation and insertion) affects the tissue's reaction force and cutting force. Adam et al. [24] proposed a new kinematicbased modeling approach that uses a meshless complete Lagrangian explicit dynamics algorithm to handle extensive deformations and strains. Xiao [25] established a puncture needle and motion function model based on the MLESAC algorithm to locate the needle position and predict the target movement. Adagolodjo et al. [26] proposed a numerical method for applying a finite element model to realize the closed-loop control of robots. The high-frequency case of needle insertion and complex nonlinear phenomena are also studied. Oldfield et al. [27] proposed a three-dimensional finite element quasi-static simulation model. It is found that compared with the traditional linear needle feeding, the puncture method of multi-body needle sub-feeding or the puncturing method reciprocating needle feeding has a smaller displacement error of the target point. Wang et al. [28] and Jiang et al. [29] studied the puncture force and target point error and established a mathematical and finite element model of the puncture process.
With the development of virtual surgery and the advances in modern medicine, performing sound preoperative planning can reduce the risks and increase the success rate. A 3D dynamic puncture model with a high degree of accuracy can provide more opportunities for trial and error for the surgeon and more realistic practice for the trainee. Current research on 3D puncture models is sparse, and muscle tissue models are less accurate in material parameters, making it difficult to guide surgery in complex working conditions. A more accurate 3D puncture model is established in this paper. The mechanical properties of muscle tissue through tensile and compression tests are analyzed. A more suitable constitutive model is selected by fitting. In establishing the 3D model, the fracture criterion and the bilinear cohesion model are used. An in vitro puncture platform is set up to perform in vitro puncture experiments. The model's accuracy is verified by comparing experimental results and simulation results-research on 3D models under different puncture speeds and puncture depths. The results are analyzed and the law of changes is summarized for getting better puncture parameters.

Stress Analysis of Soft Tissue Fibers
Currently, soft muscle tissue is generally considered as hyperelastic material due to its high similarity to rubberlike materials. Compared with muscle tissue and rubber, they are both in a linear structure, and the fiber and macromolecule chains are similar to the curling stage, which can be found in Figure 1(a), without an external force acting. Under external tensile force, the fiber and macromolecule chains become in order and show straightened linear stages like Figure 1(b). Along with the strain increases, the stress of muscle tissue firstly shows the linear increase and then becomes nonlinear increase until the muscle tissue completely straightens. The difference between muscle tissue and rubber is that in the next stage of tensile, shown in Figure 1(c), the muscle fiber will stop deformation due to the protection of the myolemma. To protect muscle fiber, accompanied by strain deformation, the myolemma will prior to starting plastic change and fracture than the muscle fiber as the red connected tissue is shown in Figure 1(c). Thus, the stress-strain curve of muscle tissue shows a linear stage when rubber material still shows nonlinear, and its slope of the curve continues to increase. This is why when using the hyperelastic constitutive model, matching muscle tissue always results in a less precise stress-strain curve.

Muscle Tissue Modeling and Tool Modeling
In this section, a novel method for the constitutive model of muscle tissue is proposed by dividing the stress-strain curve into three parts and defining their mechanical characteristics. The assumption of three parts is according to the characteristic of the constitution of the muscle tissue and the characteristic of the stress-strain curve [30]. The first part of the curve is treated as hyperelastic material from the beginning of strain changing to the end of the nonlinear stage. The second part of the curve is regarded as plastic material due to the plastic deformation of the myolemma. The third part uses the fracture model to identify its mechanical characteristic.
The first part of the constitutive model of muscle tissue is a hyperelastic model. Most of the constitutive models of biological tissues are based on the continuum mechanics of isotropic and hyperelastic materials. The strain energy density function depends on the strain invariants I 1 , I 2 and J in the three right Cauchy-Green deformation tensors. If the elastomeric material is approximately incompressible, then I 3 is considered constant and equal to 1, which means that I 3 does not contribute to the strain energy. The most classical form of strain energy can be expressed as a polynomial of invariants.
where W is the strain energy function. I 1 and I 2 are the first two right Cauchy-Green deformation tensors, respectively. J is the deformation gradient tensor. C ij and D i are the material constants.
Second Piola-Kirchhoff stress tensor S can be derived by expressing the partial derivative of W concerning Cauchy-Green strain tensor E or right Cauchy-Green tensor based on the chain rule, The Cauchy stress tensor σ can be obtained by pushing forward S.
(1) In the case of simple tests, such as uniaxial compression, biaxial compression, and pure shear, the Cauchy stress tensor can be expressed by substituting deformation gradient F=diag(λ 1 , λ 2 , λ 3 ) and with simplification.
where the invariants are given by Eq.
(2) can be written as the following direct stress components by Eq. (3).
where, λ i is the principal stretch in the deformation phase, and the soft tissue is assumed as isotropy during the deformation process, which means I 3 =1.
Here, Mooney-Rivlin (M-R) is chosen to regard as a hyperelastic constitutive model to describe the first part of deformation.
The M-R model retains two polynomials because the shear stress is proportional to the shear angle when the hyperelastic material is subjected to a simple shear load, which makes the model more accurate in tensile and shear experiments.
The second part of the constitutive model is the plastic model. The plastic stage refers to the deformation of the material under the action of an external force. Because of the particularity of superelastic material structure, there are elastic deformation and plastic deformation in the deformation process. To solve the plastic stage, some people calculate the closed area formed by the soft tissue stress deformation curve and the soft tissue recovery curve without stress. Limited by the experimental environment, this paper conducted two compression experiments on the same soft tissue to explore the area formed by the stress-strain curves under the two experiments. As shown in Figure 6, the area of the closed graph formed by the two curves is the energy lost by plastic deformation. The formula is derived as follows.
As usual, the elastic-plastic model deformation can be divided into elastic deformation and plastic deformation. The general form is as follows, where F is the total deformation gradient, F el is the total elastic gradient, F pl is the total plastic gradient. This formula can be rewritten when applied to hyperelasticity.
where F hl is the total hyperelastic gradient. Decomposition of the deformation gradient in the above formula yields the decomposition rate in strain rate as follows.
where ε is the total strain rate, ε hl is the total superelastic strain rate, and ε pl is the total plastic strain rate.
The strain energy density equation can solve the hyperelastic strain rate. In addition, in non-elastic flows, it is assumed that the yield limit is set to zero. In the nonelastic flow of material, the inelastic deformation part, according to the definition of the flow rule, can be expressed as Eq. (11). where g i (σ, θ, H i, a ) is the potential energy, dλ i is the time change rate. According to the above method, the plastic stage stress and strain data are fitted to obtain the plastic stage model.
The third part of the constitutive model is the fracture model. After the plastic stage, muscle tissue begins to break. Due to the muscle tissue exhibiting overall shear damage while breaking, so shear failure is selected in the form of fracture damage. The stress-strain relationship in the fracture model is as Eq. (12).
where D 0 hl indicates the hyperelastic stiffness at the beginning. D hl is the hyperelastic degenerate stiffness, d is the stiffness degradation variable, which ranges in size from 0 (undamaged material) to 1 (destroyed material). Failure modes combined with failure criteria lead to a sharp decrease in the elastic stiffness, and the stiffness degradation is isomorphic and determined only by a single degenerate variable d. The available stress from continuous fracture mechanics is as Eq. (13). Finally, Cauchy stress and finite stress are available through a degenerative relationship,

Uniaxial Tensile Test
The universal testing machine is used to carry out uniaxial tensile tests on the fibrous and vertical fibrous muscle tissue samples to obtain the stress-strain relationship in the deformation process. To avoid the great difference between the normal deformation and the excessive stretching speed, a stretching speed of 5 mm/ min is adopted, which ensures that the internal connective tissue does not break down rapidly due to rapid tension during the experiment. Tensile tests are carried out on the experimental tissue, and 8 groups of experiments are carried out along the fiber direction and the vertical fiber direction, respectively. Standard deviation analysis and average ultimate strength stress calculation are carried out on the experimental data results.
The following equation can express the solution of the stress-strain relationship.
where σ is nominal stress; F is tensile force; A is the crosssectional area of the sample; ε is nominal strain; L is the length after deformation. L 0 is the initial length of the sample.
The average stress-strain curve obtained after processing the experimental data from multiple experiments is shown in Figure 2. Figure 3 shows the muscle tissue tensile states, fiber direction, and vertical fiber direction at stress maximum, respectively.
It can be seen from Figure 2 that the maximum stress of muscle tissue stretching along the fiber direction is much greater than that of muscle tissue stretching along the vertical fiber direction. The strain at fracture is slightly larger in the vertical fiber direction than in the fiber direction. Taking the average value of the ultimate stress measured from the 8 groups of experiments, it can be obtained that the average value of the ultimate stress is 0.085 MPa along the fiber direction and 0.026 MPa along the vertical fiber direction. The ultimate stress along the  fiber direction is larger than that along the vertical fiber direction, which is about 3.2 times of it.

Stress Relaxation Test
The stress relaxation experiment of muscle tissue is carried out using a universal testing machine. The experimental material is pork from the back of the pig and stored in a biologically active solution. The material is stretched at a speed of 200 mm/min. When the load reached 8 N, the experiment is stopped, and the keep distance between the ends of the clamp is the same. The change of the tensile force on the force sensor within 300 s can then be obtained as shown in Figure 4.

Uniaxial Compression Test
The universal testing machine is used to carry out the compression test. The test samples similar to the tensile test are taken, except the compression fixture is replaced, and 6 groups of the same compression tests are carried out. The compression speed is set to 5 mm/min. The force and displacement data obtained from the experiments are processed by Eq. (15), and the stress-strain curve as shown in Figure 5 is obtained. As shown in Figure 5, due to the individual differences of samples during the experiment, the results of the 6 groups of experiments are averaged, and the error line is made using mean ± variance for the stress at the same strain.
According to the stress-strain curve, the compression process can be divided into three stages. The first stage is the stage of small deformation, in which the stressstrain curve is linear, and the stress increases slowly with the increase of strain. The stress-strain curves of the second stage show the nonlinear law. In the third stage, the stress-strain curve shows linear rule again, and the stress increases rapidly with the increase of strain.
The universal testing machine was used to compress the same sample twice at a 5 mm/min loading rate. Six experiments were repeated to measure the area formed by two stress-strain curves. As shown in Figure 6, the shaded areas are the energy required for plastic deformation.

Parameter Fitting in the Hyperelastic Stage
In the finite element analysis of soft biological combinations, they are generally regarded as hyperelastic bodies. To better select suitable superelastic constitutive equations, different superelastic constitutive equations are used to fit the stress-strain curves in Figure 6. The fitting results of different hyperelastic constitutive equations are shown in Figure 7.   It can be seen from Figure 7 that the linearly, non-linearly and linear stress-strain curves obtained by experiments are greatly different from those obtained by fitting. Each hyperelastic model is not accurate enough in the third stage of the stress-strain curve. Although Yeoh model is closer to the actual curve, there are still some deviations. Considering that the soft tissue strain in the puncture process is minor, the third nonlinear stage [10] is abandoned. It's noted that the composition of the two-order material constant of the M-R model performs well in simulating the slight deformation of hyperelastic materials, especially in the shear experiment, and the operation time of the M-R model is less than that of the Yeoh model and Odgen model. Under careful consideration, the M-R model is selected for curve fitting under the condition of ensuring the accuracy of the first two stages.
In the simulation of muscle tissue puncture, the stiffness degradation variable is set as 0.95, and the failure displacement is set as 4×10 −006 . The displacement criterion is used when setting fault evolution. Linear softening is used for the softening model, and maximum recession is set for recession mode. The initial linear slope is 1140, and the final linear slope is 20000.
In this work, the material constants are matched by the above method, and the final result can be shown in Table 1. The last three parameters are the material constants of the scalpel.

Preparation for In Vitro Puncture Experiments
To verify the accuracy of the 3D model, in-vitro puncture experiments are required to compare with the simulation. The biological tissue with muscle fibers and connective tissue is used in the in vitro puncture test, namely live pig muscle. Epithelial cells and large blood vessels should be avoided in the selected test samples. Large pieces of muscle tissue are cut into 100 mm×100 mm×100 mm cubes, and the surface of the test specimen should be kept as flat as possible.
As soft tissues are susceptible to force degeneration and are not easily fixed, a specific in-vitro puncture test platform is built to ensure that the puncture tests are carried out smoothly. As shown in Figure 8, the puncture motion platform consists of a precision motion platform and a fixed base. A 6-axis mechanics sensor is mounted in the center of the limited base, and a puncture needle base connects the puncture needle to the six-axis mechanic's sensor.
This allows the movement of the sliding table to be controlled by the controller, which drives the needle to perform the puncture. The soft tissue fixation box holds the soft tissue in place and prevents the soft tissue from moving during the puncture process. Due to the small force on the puncture needle, a signal amplifier is required to amplify the signal from the sensor. During the puncturing process, the controller is used to control the running speed of the sliding table, and then the puncturing speed is controlled.

Geometric Models and Material Properties
The geometric model of the puncture needle and soft tissue is established, in which the length of the needle tube is 100 mm, the outer diameter is 2.7 mm, the inner diameter is 2.1 mm, and the angle of the needle tip is 17°, the size of the soft tissue is 100 mm×100 mm× 100 mm. The performance of tissues is complex and diverse, and there are many changes in muscle tissues during the puncture process. The complex changes are one of the important factors affecting the simulation accuracy, which also leads to the difficulty in establishing a high-quality model. The 3D simulation model of the needle biopsy process is established in ABAQUS by using the display solver. This model can save computation time and reduce computation amount. In the process of setting material parameters, the soft tissue characteristic parameters measured in the experiment are input into the model, which ensures accuracy to a certain extent.
To describe the deformation behavior better, the M-R hyperelastic model is selected for muscle tissue. For some of its material property parameters, the equivalent Young's  modulus is set as 20 GPa, and the Poisson's ratio is set as 0.495. The needle is placed above the model of muscle tissue. The puncture needle is given a speed of 2 mm/s to penetrate the muscle tissue model vertically at a constant speed.

Contact and Boundary Conditions
Although the stress curves of muscle tissue are different in the fiber direction and vertical fiber direction, the mechanical properties of muscle tissue in the fiber direction and vertical fiber direction are the same. In this study, only the mechanical properties of fibrous soft tissue were considered. Therefore, soft tissue can be defined as isotropic material. In the definition of the contact algorithm, the whole set is to allow the objects to contact each other. Surface contact is used to define some critical areas of contact more accurately. The muscle tissue model is set as fixed coupling to avoid the influence of muscle tissue movement in the direction of the puncture needle tube on the puncture process. To avoid excessive component deformation, the friction coefficient between the needle and the muscle tissue is defined as 0.3, and the friction at the tip of the needle is set as no friction. The mesh division of the organization model is arranged to ensure the accuracy and efficiency of operation. As shown in Figure 9, in the area where the muscle tissue is in contact with the puncture needle, the approximate global size of the mesh in this area is set to 0.8. From this area to the edge, the mesh size gradually increases. The approximate global size of the puncture needle model is set as 1. The unit types of the model are all C3D48. A bilinear cohesion model is also introduced to simulate the contact fracture between the tissue and the puncture needle more accurately.

Analysis of 3D Simulation Results
The puncturing process of the 3D puncture simulation is shown in Figure 10. As shown in Figure 10, the stress change of puncture needles and tissue is relatively insignificant. The stress variation area is mainly in the contact area between the puncture needle and the tissue, and the stress variation is in a circular distribution. Stress concentration occurred in the contact between tissue and puncture needle, and the stress in tissue is greater than that in tissue edge due to the extrusion of deformation. According to the stress scale, the change of stress is smaller beyond the tissue edge, which is consistent with the actual situation of puncture surgery.

3D Simulation Results of Different Puncture Depths
To analyze the stress of different puncture depths and puncture needles, three stages of complete puncture simulation are intercepted. The three stages are the first contact between the puncture needle and soft tissue, the tip of the needle entering the tissue, and the puncture needle penetrating the soft tissue 5 mm. As shown in Figure 11, the mechanical properties of soft tissue hyperelasticity make the puncture needle protrude both externally and internally. The small space inside the needle leads to the extrusion pressure from the inner wall of the needle to the center of the needle in the process of puncture. The deformation of soft tissue inside the needle is more severe than that outside the needle. When the puncture reaches a certain depth, there is enough friction between the inner wall of the needle and the internal tissue. When the puncture needle exits the tissue, it takes the internal tissue to separate from the original whole, to achieve the purpose of puncture biopsy. The reaction force experienced by the puncture needle is the same as the piercing force. We analyze the reaction force experienced by the puncture needle. Figure 12 shows the reaction force of the puncture needle during puncture. The tip of the needle touching the outer membrane of the soft tissue is counted as the starting point of the puncture process. Before point X, the needle tip has not been able to penetrate the soft tissue due to the rigidity of the biological valve. The puncture needle increases as the puncture depth increase until it punctures the biological valve. After point X, the puncture needle pierces the outer membrane of the soft tissue, and the reaction force drops suddenly. At this time, the needle tip will damage the soft tissues, making the needle tip more stressful. The part behind the needle tip rubs against the soft tissue, making the stress of this part smaller than the needle tip. The reaction force is regarded as the combined force of the cutting force of the needle tip on the soft tissue and the friction force received by the outer wall of the needle. As the puncture depth increases, the frictional force increases and the reaction force also increases.
The stress curve in Figure 12 is not smooth. As the puncture depth increases, multiple sudden drops appear on the stress curve. The ideal curve should be linear. It is preliminarily speculated that this phenomenon is due to the unevenness of the muscle and soft tissue. If a uniform tissue prosthesis is punctured, there are only two sudden drops of force, one at the moment of piercing the soft tissue and the other at the moment of piercing the tissue.
The following discussion of different puncture speeds will be further analyzed through experiments.

3D Simulation Results of Different Puncture Speeds
To analyze the stress changes at different puncture speeds, 4 groups of simulations with puncture depths of 3 mm, 6 mm, 9 mm, and 12 mm are designed in Figure 13.
The puncture depth of each group is the same, and the puncture velocities from left to right are 2 mm/s, 5 mm/s, and 10 mm/s, respectively. As shown in Figure 13, with the increase of puncture speed, the fitting degree between the outer wall tissue and the puncture needle gradually decreases, and the stress change increases. The fitting degree of the internal tissue of the puncture needle and the puncture needle decreases, and the strain increases.
The higher the fit of the puncture needle to the tissue, the greater the puncture resistance and the more likely it is that tissue adhesions will occur, which can cause secondary wound injury. Therefore, maintaining a high puncture speed within a certain range during puncture is beneficial in reducing tissue damage and increasing the rate of wound healing.
To verify the accuracy of the model, the built in-vitro puncture platform is used for in-vitro puncture experiments, the feed speed of the sliding table is controlled, the puncture stress at three speeds of 2 mm/s, 5 mm/s, and 10 mm/s is counted, and the puncture stress is calculated. The results are compared with the simulation results. Figure  It can be seen from Figure 14 that during the puncture process, the stress curve has multiple sudden drops. Through experimental analysis, the reasons for this phenomenon can be divided into two types. One is the vibration of the puncture needle during the puncture process, and the other is caused by the interaction between the tissue and the puncture needle. The degree of sudden drop caused by the interaction between the tissue and the puncture needle is greater than the degree of sudden drop caused by needle vibration. Actually, the interaction between the tissue and the puncture needle is analogous to cutting a workpiece with a tool. During the cutting process of the tool, some high-temperature chips will form a built-up edge on the tool edge. During the puncture process, due to the stickiness of the soft tissue, part of the soft tissue adheres to the needle tip, forming a soft tissue "built-up edge".
Due to a "built-up edge" in soft tissues, the reaction force curve shows a nonlinear growth. When the "builtup edge" disappears, there will be a sudden drop point. During the puncture process, the "built-up edge" phenomenon also continuously produced, disappeared, and reproduced. The frequency of occurrence determines the number of sudden drop points during the puncture process. With the increase of the puncture speed, the frequency of the occurrence and disappearance of "built-up edge" also increased. The sudden drop points also increased, and the puncture stability of the needle decreased.
As shown in area Z in Figure 14, there is a certain deviation between the experimental data and the simulation results. During the experiment, the puncture needle can be seen as a cantilever structure fixed at one end, and it's affected by gravity. The puncture platform inevitably produces small amplitude vibration during the puncture process. Due to the influence of equipment limitations and experimental environmental factors, certain errors appear between simulation and experimental data. The fixation method of soft tissue in in-vitro puncture experiments will also bring some errors. However, it is not as influential as equipment constraints and experimental environmental factors. It can be seen from Figure 14 that the experiment and simulation curves have a high degree of coincidence, and the above-mentioned error influence is small and can be ignored.

Conclusions
This paper studies finite element simulation of the interaction between flexible puncture needle and muscle tissue and establishes a high-precision 3D dynamic puncture model for characteristics of muscle tissue. The main contributions of this study are listed as follows.
(1) An efficient method of modeling muscle tissue composition is proposed for the properties of muscle tissue. The stretching, compression, and stress relaxation experiments of muscle tissue are carried out and the complex mechanical properties of muscle tissue are analyzed. By fitting multiple constitutive models, the M-R model is adopted to simulate and determine the characteristic parameters of each stage of muscle tissue. The results show that the M-R model considers both fitting accuracy and calculation speed. (2) An accurate 3D dynamic puncture model is established. The experimental data is used to set soft tissue parameters during building the 3D model. Bilinear cohesive force model is used to define the fracture during puncture contact. This can improve the calculation speed of the model on the premise of ensuring accuracy. The failure criterion is adopted for the fracture characteristics of muscle tissue materials. The model's accuracy is verified by comparing the results of the in-vitro puncture experiment with the simulation results. The accuracy and practical application value of the 3D model is much higher than that of the 2D model.