Simulation of Coupling Process of Flexible Needle Insertion into Soft Tissue Based on ABAQUS

: In order to get to the desired target inside the body, it is essential to investigate the needle-tissue coupling process and calculate the tissue deformation. A cantilever beam model is presented to predicting the deflection and bending angle of flexible needle by analyzing the distribution of the force on needle shaft during the procedure of needle insertion into soft tissue. Furthermore, a finite element (FE) coupling model is proposed to simulate the needle-tissue interactive process. The plane and spatial models are created to relate the needle and tissue nodes. Combined with the cantilever beam model and the finite element needle-tissue coupling model, the simulation of needle-tissue interaction was carried out by the ABAQUS software. The comparing experiments are designed to understand the needle-tissue interactions, by which the same points in the experiments and simulation are compared and analyzed. The results show that the displacements in x and z directions in the simulation can accord with the experiments, and the deformation inside the tissue mainly occurs in the axial direction. The study is beneficial to the robot-assisted and virtual needle insertion procedure, and to help the physicians to predict the inside tissue deformation during the treatments.


Introduction
At present, the minimally invasive and local treatments have become the major development trend in the medical field, especially the application of minimally invasive surgery is more and more extensive. As a typical minimally invasive surgery, needle insertion is mainly used in the biopsy of organs, local anesthesia and brachytherapy [Abolhassani, Patel and Moallem (2007)]. In the actual insertion, the causes of insertion errors include the resolution of organ image, inaccurate focus location, tissue deformations, doctors' experience among others. Carr et al. [Carr, Hemler, Halford et al. (2001)]. The success rate of the operation is directly affected by the accuracy of insertion, and imprecise insertion can lead to serious complications. The majority of insertion operations is done by the rigid needles, while the rigid needles have little deformation during the insertion process, and the obstacles cannot be avoided. In many cases, multiple insertions are required to carry out the diagnosis and treatment procedure, so the researchers are replacing the rigid needles with flexible needles [Webster, Memisevic and Okamura (2005) ;Webster, Kim, Cowan et al. (2006)]. However, the insertion and retraction of the flexible needles can also cause the deformations and displacements of the surrounding soft tissue, for which the needle tip fail to reach target positions. In recent years, many scholars have devoted themselves to the research of insertion modelling and insertion process. The research focuses on two aspects of the insertion process, one is needle deformation, and the other is tissue deformation. By using nonholonomic kinematics, Webster et al. [Webster, Kim, Cowan et al. (2006)] controlled the needle to reach the target position according to the planned path. The motion model can effectively improve the accuracy of the insertion targeting and promote the therapeutic effect. Abayazid et al. [Abayazid, Roesthuis, Reilink et al. (2012); Glozman and Shoham (2004)] proposed a cantilever beam model of the flexible needle with a virtual spring. Using this model, they calculated the forward and inverse kinematics of the needle tip, and proposed a method to simulate and predict the needle's path in real-time. However, when the deflection is small, this model is more effective for predicting the needle trajectory. When the bending degree of the needle is large, the error is large. ] analyzed the insertion trajectory and tissue deformation during the percutaneous needle insertion procedure. They described a method for estimating the force distribution of the needle during the insertion process. Using the discrete finite element method, they analyzed the interaction force between the needle and the tissue, and calculated the displacement and deformation of planar tissue. Cotin et al. [Cotin, Delingette and Ayache (1999); Yu and Dong (2019)] established a volume model of medical images and an elastic model of tissue deformation, and performed real-time calculation of liver deformation by finite element method. Alterovitz et al. [Alterovitz, Goldberg, Pouliot et al. (2003); Alterovitz, Goldberg and Okamura (2005)] developed a 2D dynamic FE model for the needle-tissue interactions. They simulated the effects of physician-controlled parameters such as needle friction, the displacement of needle tip and deformation of soft tissue. Zhang et al. [Zhang, Wang, Sun et al. (2018)] proposed a physical model of virtual soft tissue, which is a twist model based on the Kriging interpolation and membrane analogy. Gao et al. [Gao, Lei, Lian et al. (2016)] proposed a needle-tissue coupling model based on the improved local constraint method, and studied the analysis method of force and displacement of each element node in the soft tissue. Gao et al. [Gao, Zheng, Yao et al. (2010)] established a virtual spring model of soft tissue inserted by the flexible needle in two-dimensional space, on the basis of analyzing the force and deflection of the flexible needle, and obtained the mathematical model of the needle shaft and the trajectory of the tip. Chentanez et al. [Chentanez, Alterovitz, Ritchie et al. (2009)] proposed a new coupling algorithm and optimization technique to improve the efficiency of solution and calculation, and a coupling simulator was used to predict the deformation of the tissue. Kshrisagar et al. [Kshrisagar, Francis, Yee et al. (2019)] presented the node based smoothed finite element method for the linear and nonlinear problems. Yamaguchi et al. [Yamaguchi, Tsutsui, Satake et al. (2014)] used the hydrogel as a bionic soft tissue for the insertion experiments. The FE tools, like ADINA and ABAQUS, has been used for simulating the liquid-structure and thermal coupling procedures [Peng, Chen, Yu et al. (2018); Gao and Wang (2018)]. The finite element simulation was carried out in threedimensional, which can dynamically analyze the needle deflection and the tissue deformation during the needle insertion process. However, they didn't consider the interaction between the flexible needle and soft tissue in three-dimensional space, and those models didn't reflect the evolution of the flexible needle and soft tissue during the insertion procedure. In this paper, a method for analyzing the force of the needle shaft during inserting soft tissue in two-dimensional and three-dimensional has been presented. The position of the needle tip is predicted through this bending model of the needle in real time. Two coupling models are proposed to simulate the interaction force between the needle and the soft tissue. The bending model of the needle and the coupling model are applied to the two-dimensional (2D) and three-dimensional (3D) insertion simulations. At last, a series of experiments have been presented to prove the feasibility of the above method. The paper is organized as follows. Firstly, the forces of the flexible needle in the insertion process are analyzed, and the prediction model of the needle trajectory is established. Secondly, the quadrilateral and hexagon mesh types of plane and spatial tissue are constructed with the finite element method, and the coupling models of needle-tissue are created in the two-dimensional and three-dimensional spaces. In Section 4, the insertion simulation and experiment are carried out, and the simulation results are compared with the experimental results. Finally, the conclusions and future work are given.

Force analysis of the needle
In this section, a force analyzation and a deflection model of insertion needle are presented to predict needle trajectory. During the insertion process, the needle will be inclined to the direction of bevel because of the asymmetry pressure [Alterovitz, Goldberg and Okamura (2005) ;Misra, Reed, Douglas et al. (2008); Glozman and Shoham (2006)]. The following study assumes that the needle and soft tissue are linearly deformed in two-dimensions. The trajectory of the needle shaft is consistent with that of the tip. During the insertion process, the needle is subjected to the interaction forces between the needle and the tissue. The force F C at the bevel tip of the flexible needle is the main cause of the bending. The part of the needle that enters into the tissue is subjected to the friction Ff. The direction of the friction is opposite to the insertion direction, and it is parallel to the needle shaft. Similarly, the needle tip is subjected to friction parallel to the oblique plane Ft. Assuming the radial pressure of the needle shaft is inserted into the soft tissue as a triangular distribution load q0. The force analysis of the needle is shown in Fig. 1.

F = F -F sinθ -F cosθ -F = 0 q (L -L ) F = -F sinθ -+ F cosθ = 0 2 q (2L -3L L + L ) M = M + -F Lcosθ = 0 6(L -L )
(1) Using Eq. (1), the balance equation for the force in the x direction and the y direction that describes the force state at any time are presented. The bending and path of the needle can be determined by the moment equation of needle, which will be described in later. The curvature equation is used to calculate the needle deflection and the tip angle deviates from the direction of the needle in two-dimensional space.

Needle track prediction model
Setting the needle base as the coordinate origin of the x direction, the point on the needle is x (L0≤x≤L). Building the relationship between the deflection W and the coordinates in x direction as W(x), so the relationship between needle deflection and needle curvature is shown in Eq. (2).
The curvature equation for flexible needles in the online elastic range and pure bending is shown in the Eq. (2), ρ is the bending curvature of the needle, and M is the total moment force equation of the beam model, E is the elastic modulus of the flexible needle, I is the moment of inertia of the needle. So W(x) is shown in the Eq. (3).

W(x) = -( Mdx )dx EI
(3) Because Ft is so small compared to C F and f F , the effect of friction on the tip is negligible. Bending moment equation can be calculated as Eq. (4).
[ ] According to the force balance Eqs. (1) and (4), the deflection equation and the angular equation of the needle shaft can be obtained as shown in the Eqs. (5) and (6).
[ ] (1). Therefore the deflection and the deflection angle at the needle tip are determined by the insertion force in the x direction and the moment of the needle base in the x-y plane.
3 Needle-tissue coupling model By analysing the force relationship between the needle and soft tissue, the finite element method of the quadrilateral and hexahedral mesh type is proposed for the tissue in the plane and space. The coupling analysis model of the needle insertion into soft tissue is established under different conditions. The quasi-static analysis of the insertion process is regarded as the process of inserting the entire needle shaft into multiple discrete needle segments. The needle segments are connected by nodes, each needle segment is analysed. The needle and tissue are analysed separately. The interaction between the flexible needle and the tissue is translated into the interaction between the needle shaft nodes and the soft tissue nodes.

The 2D needle-tissue coupling model
The finite element method is often used to analyse the process of flexible needle insert into soft tissue. It is necessary to study the interaction relation between the needle and soft tissue when analysing the deformation of soft tissue. The coupling model of the force between the needle and soft tissue is established. The flexible needle is seemed as a whole that is consisted of a certain number of discrete needle segments, and one node is connected between each segment. The interaction force between the needle and the tissue is transformed into the effect of the needle nodes on the nodes of the soft tissue elements in the finite element method. The position of the needle nodes during the insertion process can be divided into two categories, and the following two situations are explained.
As the needle deepens, the needle tip will gradually deviate from the initial direction, and the position of the needle will constantly change. It is assumed that when the distance between the needle node and the closest tissue node is less than one tenth of the grid size, it is considered as the nodes superposition state. In other words, when the needle node and soft tissue node are approximately coincident in a certain precision range, as shown in Fig. 2. Defining the set of overlapping nodes P={P1, P2, P3, P4, ···, Pn}.   Figure 3: Coupling model under general cases If the needle nodes are not overlapped with the soft tissue element nodes, the force between the needle and the soft tissue can be analysed by local constraint method. As is shown in Fig. 3, defining a point set P={P12, P23, P34, P45, ···, Pn(n+1)}, N={N1, N2, N3, N4, ···, Nn} and M={M1, M2, M3, M4, ···, Mn}. The distance between P12 and N1 node in the x direction is 12 1 P N X , and the distance between P12 and N1 node in the y direction is When analysing the whole tissue, the forces of each node of the tissue element are calculated as shown in Eq. (7).

-λ )δ F + (-
In the above equation, When analysing the interaction between the needle and the tissue in the three-dimensional space, similar to the plane, the network division is performed using the element type of the regular hexahedron for the tissue in the space. When the flexible needle is not rotated, the needle trajectory will always remain in the same plane. When the plane of the needle insertion trajectory is not parallel to any plane of the grid element within the tissue, then the spatial process of the needle insert soft tissue needs to be analysed.

The 3D needle-tissue coupling model
Studying the process of insertion in three-dimensional space, the relative positions of the flexible needle and the soft tissue are also divided into two types, the special position and the general position. The insertion coupling model of the two relative positions is analysed below. As shown in Fig. 4, define the coincidence node set P={P1, P2, P3, P4, ···, Pn}. Since the needle nodes coincide with the tissue nodes, when the forces of the tissue element are studied in three-dimensional space, the mutual friction force is directly applied to the coincident node of the tissue element in the form of a segment. If the needle trajectory does not pass through the nodes of the tissue element, that is, when the needle nodes are located inside the tissue finite element. The interaction between the needle and the tissue needs to be applied to the node of the tissue element through the needle node. Using a local constraint method to analyse the forces of a regular hexahedron element. As shown in Fig. 5, define the node set P={P12, P23, P34, P45, ···, Pn(n+1)}, and the node set of the tissue element A={A1, A2, A3, A4, ···, An}, B={B1, B2, B3, B4, ···, Bn}, C={C1, C2, C3, C4, ···, Cn} and D={D1, D2, D3, D4, ···, Dn}. When the needle node P12 is inside the element, it is defined that the distance from the A1 node in the x direction is 12 1 P A X , the distance in the y direction is 12 1 P A Y , the distance in the z direction is 12 1 P A Z , and the distance ratio variable When analysing the tissue element on the insertion trajectory as a whole, when i is greater than 1, the force distribution relationship between the needle node and each node of the tissue element is as shown in the following Eq. (8).

Result analysis of simulations and experiments 4.1 Result analysis of 2D FE model
The markers are used to study the influence of the needle movement on tissue deformation. In the ABAQUS simulation, establishing the tissue model with the size of 0.1×0.1 m 2 in the part module in two-dimensional. The model of 21 G type of needle is established, with the tip angle of 20°, the length of 0.2 m and the diameter of 0.8 mm.
The material properties of needles and soft tissue are given in the property module. The displacement freedom of the three edges are constrained fully as 0, and the insertion edge is free. In the mesh the overall grid size is defined as 0.01 m at first. Select Quadrilateral element shape and the type of free to mesh the tissue model. Local seed is carried out in the nearby of the coupling to refine the mesh, defining the size of local seed as 0.001 m in the following simulations. So we can simulate the coupling state of arbitrary insertion time. The whole simulation process is divided into several steps according to the needle nodes' displacement, the static mechanical analysis is set in the definition analysis. The load of each node on the needle track in the tissue model is determined by the insertion step analysis. The coupling simulation process of insertion depths 0.01 m, 0.04 m and 0.07 m are shown in Fig. 6. Five markers are set in the route from the initial insertion direction of 5 mm and parallel to that direction, with the distance between each mark point is 0.02 mm. The displacement of tissue elements and the mark points are obtained through the simulation process as shown. The displacement of the tissue model continues to be bigger with the increase of the insertion depth. It is also obvious that the tissue closes to the needle shaft has a larger displacement than the distant. In the model, the displacement of the mark point near the pierce point is also greater than that the distant. The needle of 21 G and the same size of prosthesis as simulations are used for the insertion experiments. The component ratio of the prosthesis used in the experiment is polyvinyl alcohol (g): volume of distilled water (ml): volume of dimethyl sulfoxide (ml) =8:40:60 [Zhu, Gao and Zhao (2015)]. Because the mechanical properties of the Polyvinyl Alcohol (PVA) prosthesis are most similar to those of liver tissue. And it is highly transparent, conducive to the acquisition of images. Adjust the insertion angle of flexible needle to ensure that the bending of the needle is as far as possible in the plane of markers. A CCD camera is used to absorb all the images of tissue during insertion. Then it follows up processing those images to obtain the coordinates of all the markers in the prosthesis at different stages. The acquired images and processed images are shown in Fig. 7.  Fig. 8. The displacement of the markers in the x direction and the y direction is shown, respectively. It can be seen from the comparison diagram that the displacement of each mark point is increasing continuously in the x direction during the experiment. But the increase of displacement is unstable, because the dislocation vibrations between the needle and the tissue during insertion. In the experiment, the displacement of each mark point is suddenly decreased in the front part of the x direction. Because the tissue's wallet was inserted and the tissue rebounded. The displacement trend of each marker is consistent, and the maximum error of the No. 1 to No. 4 markers between simulation and experiment is less than 0.2 mm. In the simulation, the displacement of the No. 5 marker in the x direction is much smaller than in the experiment. Because the movement of the tissue is completely constrained in the simulation, but there is a gap between tissue and its container due to the irregularity of their surfaces in the experiment. It can also be seen from the diagram that the displacement of tissue is relatively small in y direction. As can be seen from the Fig. 8(b), the displacements are less than 0.1 mm, and the maximum error is less than 0.05 mm in the y direction, which occurs in Marker 1. Therefore, the errors are less than 0.5 mm, and the 2D FE model can reflect the procedure of flexible needle insertion into soft tissue.

Result Analysis of 3D FE model
Since the study of tissue deformation can only be obtained in the form of a planar image, the deflection and tip angle of the flexible needle during the insertion process are always on the same plane. When the plane of the needle deflection does not coincide with the plane of image capture. At this point, it is necessary to carry out insertion research in three-dimensional space. In order to simulate the process of flexible needle inserting soft tissue in three-dimensional space, the soft tissue is accurately meshed in real time. The predicted trajectory of the needle is applied to the finite element model of soft tissue to simulate the effect of insertion on tissue deformation. Fig. 9 shows the path of the needle and the displacement cloud of the piercing point section when the insertion depths are 10, 40, 70 mm. As shown in Fig. 9, as the insertion depth increases, the trajectory of the needle gradually deviates from the insertion direction. It can be seen from the cross-sectional displacement cloud diagrams of (a) and (b) in Fig. 9 that the deepening of the needle insertion gradually shifts the deep soft tissue and increases the displacement of the soft tissue near the trajectory. Comparing (a) and (b), it can be seen that the main deformation of the soft tissue caused by the flexible needle bending is no longer a vertical plane, but a plane in the direction of deflection. Make a bionic prosthetic soft tissue of the same scale as the model, and place it in the prosthetic container. To measure the tissue deformation, markers are placed on the upper surface and a side surface of the soft tissue. The two CCD cameras are respectively disposed on the two faces, vertically aligned with the plane. When the needle is inserted from the other side, the above two cameras are used to perform image acquisition of the tissue deformation process. The insertion and image acquisition system, the shootings of the piercing surfaces are shown in Fig. 10. The coordinate system of the insertion is shown in Fig. 10(b). In two perpendicular free faces, the markers are set corresponding to the equipotential position of the piercing point. Five markers are sequentially arranged in the piercing direction at intervals of 20 mm in the two faces. Five markers are respectively extracted as the research object of the soft tissue deformation. The setting of the markers in the simulation is shown in Fig. 10(c). The displacements of the ten markers are obtained through the simulations and experiments. The maximum error between the simulations and experiments is 0.3500 mm in the x direction, which occurs in Marker 10. The displacements are less than 0.05 mm in the y and z directions, which is much smaller than that in the x direction. Therefore, the simulation model agrees to the experiments very well, and the 3D FE model is valid for the needle-tissue interactive procedure. As shown in Figs. 11 and 12, with the deepening of the needle insertion during the insertion process, the mark displacement on the soft tissue in the x direction continues to increase. While the mark displacement in the y direction and the z direction changes, but they are smaller than the displacement in the x direction. As can be seen in the figure, the displacement of upper surface marker caused by the insertion process is relatively larger than the side surface. This is because the sides of the tissue are not completely free and are also partially constrained by the soft tissue container.

Conclusions and future work
In this paper, a novel needle-tissue coupling model is developed for simulating the needle-tissue interactions. The deflection and rotation angle of flexible needle are calculated with the cantilever beam model. The needle-tissue interaction is translated into the interaction between the needle shaft node and the tissue node, which implement the force transmission between the nodes of flexible needle and soft tissue. The 2D and 3D FE models are constructed in the ABAQUS software. The experiments are designed to verify the validation of the needle-tissue coupling model. The markers inside the PVA phantoms are set to be compared with the mark of the simulation model, and the displacements in the x and y directions in the simulation can agree to that in the experiments well. The results show that the needle-tissue coupling model can reflect the mechanism of needle insertion into soft, and the two and three-dimensional models can accurately simulate the deformations of flexible needle and soft tissue. Future work needs to investigate the parameters of real tissue for applying the needle-tissue coupling model into the inhomogeneous tissue environment. Furthermore, the novel coupling model can provide the accurate calculation for the virtual needle insertion, clinical diagnosis and treatment, which can help the unexperienced and experienced physicians to improve the skills and intraoperative information.