FEA Modeling of Soft Tissue Interaction for Active Needles with a Rotational Tip Joint

A finite element analysis (FEA) model was developed for investigating the design and the tissue interaction of actively steerable needles with a rotational tip joint in soft tissue, based on the coupled Eulerian-Lagrangian (CEL) mesh. The new model is algorithmically simple in that the proposed CEL strategy allows needle insertion simulations along non-predetermined paths in soft tissue without the re-mesh at every direction change of the needle tip. For the features, various steering motions can be simulated and studied for different needle designs in a simpler approach by changing the needle geometry and boundary conditions. The developed FEA model, using the thoroughly measured material properties, was validated for predicting insertion path and estimating the tissue interaction forces inside two different gelatin tissue phantoms. Further, using the validated model, the effect of tip geometry on tissue was briefly investigated. Given bevel angles, it was found that the ratio of tip length to the diameter dominates the tissue damage gradient. The results demonstrate that the proposed model effectively examines the steering and tissue insertion of actively steerable needles and investigates the tip design to minimize the tissue damage by the needle steering.


I. INTRODUCTION
In percutaneous needle insertion, precise targeting is not easy because the flexible needle is unexpectedly deflected due to its interaction with the surrounding tissue. Sometimes, clinicians need such a deflection to detour obstacles and organs prone to damage under excessive force. Still, it is not easy to precisely control the needle path. Steerable needles have been extensively studied to address these issues employing robotics technologies for almost two decades, categorized into the passive and the active depending on the steering principle [1]. The passive ones use the base manipulation only, whereas the conceptually newer, active ones have additional steering input in their tips for improved capability [2][3][4][5]. However, regardless of the principle, the steerability, i.e., changeability of insertion direction or paths, commonly relies on the mechanical interaction between the tip and tissue [6]. The needle insertion in soft tissue also results in unavoidable tissue damage because of the crack propagation ahead of the tip and the friction between the needle shaft and tissue [7]. Ironically, more-capable active needles may cause a higher damage gradient due to its localized high curvature motion near the tip unless it is appropriately designed and controlled. Thus, understanding the tip-tissue interaction is the key to achieving precise and safe needle insertion in soft tissue with better steerability.
Finite element analysis (FEA) techniques have been employed for studying the interaction because it is challenging to account for the large deformation in the viscoelastic tissue environment analytically. For example, several FEA models were developed for the passive needle, as a simple rod element [8] and for a notched body needle [9]. Advanced FEA models, based on the cohesive element (CE) theory, also have been developed and simulated to investigate crack propagation, the resulting compressive force, cutting forces required for penetrating a soft tissue, and tissue damage [7,[10][11][12]. The CE theory grants a reasonable estimation of the contact conditions between needles and tissue. However, since the element region must be predetermined, the insertion path should be pre-defined as a fixed boundary condition, limiting its application to the active needles that change their tip orientation and path during insertion responding to the steering input.
A recently proposed solution was applied to a unique active needle with prismatic tip joints [13], using the CE coupled with a dynamic-explicit FEM model implemented in ABAQUS (Dassault System, France). Although it is the first and the only working model for active needles, the analysis should be interrupted at every integration step to update the changed needle path and the mesh boundary conditions, unnecessarily increasing the computation time. The Arbitrary Lagrangian-Eulerian (ALE) formulation, studied for passive needles [14][15][16], can avoid the issue relevant to the CE region. Still, it also requires a similar mesh update from the Lagrangian mesh to a newly generated Eulerian mesh at every step.
The coupled Eulerian-Lagrangian (CEL) formulation represents an alternative solution, especially for more efficient FEA modeling of actively steerable needles in soft tissue, where the Eulerian and Lagrangian domains are defined separately. The Eulerian domain in the CEL strategy is fixed and defined as a boundary condition [17][18][19]. Then, the Lagrangian mesh can deform in the Eulerian domain according to the applied forces, displacements, or interaction with other rigid or deformable elements. Such a CEL method has been used for investigating the stress distribution in soft tissue as a sequence of the insertion of a passive biopsy needle without steering [19] and the effect of bevel geometry on the insertion path [20]. In these studies, Eulerian brick elements, applied to the tissue phantom domain, have a super fine (0.1mm edge length) mesh size to assure accuracy. Otherwise, when the mesh size is not properly optimized, the computation time would become significantly longer.
In addition to the mesh-related algorithmic issue of the previously developed models, many of them employed fully elastic materials for soft tissue, limiting the modeling accuracy. For example, a Eulerian hydrocode coupled with a Lagrangian mesh was developed to estimate the insertion force of passive needles [21]. The study highlighted how such FEA models underestimate the insertion force attributed to the hyperelastic reaction force, especially when deeply inserted.
This paper aims to address the issues above by developing an improved FEA model suitable for actively steerable needles, based on the CEL mesh strategy that doesn't require the re-mesh or mesh update at every step. Although an FEA model exists for the unique needle with prismatic tip joints [13], even if the re-mesh issue of the model is ignored, such joints are inferior to the rotational or bendable joints suitable for achieving tighter turns, i.e., better steerability. Therefore, this study presents, validates and, applies an algorithmically simpler FEA model, thanks to the CEL mesh strategy, to more-capable active needles with a rotational tip joint for the first time. Further, the developed model considers and is validated in the hyperelastic behavior of tissue environment for more accurate estimation, then investigates the effect of tip geometries of active needles other than the bevel angle on tissue damage for the first time.
The rest of this paper is organized as follows. First, the FEA model is presented in Section II, followed by its experimental validation with two different needle prototypes in two different gelatin tissue phantoms in Section III and results and discussion in Section IV. The validated model is then applied to a case study on the effect of tip geometry in Section V, followed by discussion and conclusions in Section VI.

II. FEA MODELING
This section describes the details of the FEA model and its implementation, carried out in the commercial software ABAQUS (Dassault Systèmes Simulia Corp, Providence, RI, USA) [22], including the environment design, meshes, contact and friction conditions, and boundary conditions, followed by the parameter values used for the simulation.

A. ENVIRONMENT DESIGN
The active needle studied in this research is composed of a needle tip with an asymmetric bevel profile connected to the needle body through a rotational joint (Fig. 1a). The dimensions included are those of the prototypes used for the experimental validation in Section 3, determined by the capability of the SLA 3D printer and a standard resin in the group (Form3, Formlabs, USA). In addition, a passive needle with a symmetric tip is included in this research for validating the FEA model and demonstrating the generality of the proposed solution (Fig. 1b). Although this section focuses on the implementation procedure for active needles, the FEA model for passive needles shares the same features and strategy, except for the steering boundary conditions. On the other hand, the tissue phantom is defined as a 70mm x 70mm x 35mm volume.

B. MESHES
The FEA model considered a 3D mesh for both the needle (fully Lagrangian mesh) and the tissue (CEL mesh) required by the CEL formulation. Still, only one-half symmetric geometry is selected for the reasonable computational time. The Eulerian domain defines not the geometry of the tissue but the boundary in the z-and x-directions for the tissue material to be able to move during the simulation (Fig. 2a) Then, the portion of tissue material is defined by the Eulerian volume fraction (EVF). An EVF=1 represents an element of the Eulerian domain (mesh) where the tissue geometry is fully present, whereas EVF=0 is an empty element that is still included in the boundary but neither crossed nor filled at the considered integration step. The calculation of the EVF parameter is carried out based on the Lagrangian mesh that defines the time-dependent position of the tissue phantom, accommodating the variation of the needle position.
At the top of the model where the needle penetrates the tissue, the Eulerian domain has been set 5 mm higher than the Lagrangian (tissue) domain for tissue deformation along t he negative y-direction during the insertion (Fig. 2b). The progressive penetration of the needle creates the crack propagation ahead of the needle tip and the elastic deformation on the tissue surrounding the needle body. At this stage, a fine mesh region having a 20 mm x 10 mm cross-section and a global element length of 0.5 mm has been defined in the global mesh size of 3 mm and represents the region of the model where the needle is expected to move during the insertion. Both regions of the tissue geometry have been meshed with the 8-node linear cubical EC3D8R element of ABAQUS [22]. Besides, the needle assembly has been meshed with the hexahedral C3D8 element of ABAQUS [22] with a global size of 0.5 mm and a minimum size of 0.1 mm at the tip.

C. CONTACTS AND FRICTIONS
The contacts between bodies or meshes have been defined considering a general contact algorithm coupled with the penalty contact method, where the Lagrangian domain of the needle assembly and tissue are defined as the master and the slave, respectively. Also, five mesh strategies have been employed to achieve the high accuracy and more realistic tissue deformation around the tip, considering three coarser and one finer one. For this reason, the selected mesh sizes and strategy would compromise the computational time and accuracy because no significant variations have been identified in the reaction forces during the validation summarized in Section III.
Regarding the rotational tip joint, a multi-point constraint (MPC) is added, where the needle body and the tip are constrained to translate as rigid entities while they freely rotate around their centers of rotation. This constraint helps avoid calculating the contact stress in the joint socket, reducing the computation burden. Still, it allows for the rotation of the needle tip based on the applied steering boundary conditions.
The friction between the needle body and tissue is assumed to be Coulomb friction. Although viscous friction can be involved, the tissue phantoms used in this study show low dependency on the velocity (refer to Fig. 4). Also, the insertion force for a passive needle with a symmetric tip, the sum of cutting and friction forces, linearly increases after the punction as the insertion depth increases (refer to Fig. 6(a)). Since the cutting force is consistent for the constant speed and straight

Figure. 3 The proposed CEL mesh and tip steering strategies: (a) Mesh definition, (b) Update procedure for the CEL mesh for the tissue phantom and the fully Lagrangian mesh of the needle assembly, (c) tip steering for the active needle for three different conditions, straight (I), upward (II), and downward (III)
insertion into the homogeneous phantoms, the friction can be modeled as a simple linear function of the insertion depth only, i.e., Coulomb friction. However, it is still hard to directly estimate the coefficient during insertion experimentally. Thus, five candidates (µ=0.01, 0.05, 0.1, 0.18, and 0.25) were applied to the simulation model, whose results were compared with the experiment regarding the end-point force at the needle base. As presented in Fig. 6 (a), the coefficient of 0.1 resulted in the minimum difference of the end-point forces for both phantoms and was thus, chosen and applied to all the simulations in this study.

D. BOUNDARY CONDITIONS
The tissue phantom is encapsulated in a fixed region, called the encastre boundary condition in ABAQUS, where only vertical displacement opposite to the needle insertion direction (positive y-direction) is allowed. The insertion speed of the needle assembly is controlled through a velocity boundary condition applied to the end-point of the needle body. The tip steering is carried out by a separate application of displacement boundary conditions on the needle tip, Fig. 2c. The needle insertion, defined as the rigid movement of the needle assembly, starts from the initial step of the simulation, and progresses until the final integration step. For the FEA model, penetration and tip steering are performed as nonsimultaneous and sequential operations to avoid contact issues between the needle tip and the tissue during the onset of their interaction. After the initial penetration is made, i.e., the needle tip has entered the tissue, steering is initiated, following the medical scenarios of steerable needles. The displacement applied to the two sides starts not from zero but a pre-tensioned configuration keeping the joint kinematics, defined as in Fig. 2c. In the experimental validation, the pre-tension is removed from the signals measured on the tendon force sensors.

E. OTHER CONSIDERATIONS AND PARAMETERS
All the FEA models for the active needle start from the straight shape (Condition I) and, by changing the displacement equilibrium on the two tendons, upwards (Condition II) and downward (Condition III) rotations can be achieved and utilized to steer the needle within the tissue phantom, Fig. 2c. The input values for the steering are the linear displacements applied to two tendons, one for each side, and the measured parameters are the reaction force on the tendons and the relevant rotation of the needle tip. Furthermore, considering the time step increment in the FEA models, the angular velocity of the needle tip also can be estimated.
The hyper-elastic material properties of tissue phantoms are experimentally measured and inserted into the FEA model using the Ogden strain energy potential model constants, detailed in Section III. In addition, the mechanical properties of the needles have been characterized by tensile tests. All the property parameters measured and used for the simulation are summarized in Table I, as well as the geometry and control parameters.

III. EXPERIMENT
This section describes the experiments for validating the developed FEA model, which uses a 3D printed active needle prototype inserted into two different tissue phantoms simulating the hyperelastic behavior of healthy and cancerous tissues, which is represented by softer and harder tissue, respectively.

A. NEEDLE PROTOTYPES
The needles, shown in Fig. 5, were manufactured using the 3D stereolithography additive manufacturing process and the Acrylate polymer (MD-R001CR, ApplyLabWork, Torrance, CA, USA). The deposition layer thickness was set to 0.025 mm and cured with the UV light of 385~405nm for 120 minutes at 60°C.
Since the needle diameter can influence the resulting mechanical properties, normal for 3D printed parts, tensile tests have been carried out for the tensile specimens having a cylindrical calibrated region with a diameter ranging from 1mm to 5mm. The thickness-dependent Young's moduli for the considered resin is shown in Fig. 3, and the 3mm diameter is selected for its consistent printing quality and reasonable deflections in the tissue phantoms fabricated. The tip is steered by pulling and releasing two tendons of Aramid Sewing Thread 1200D (Jiangsu Yuhuan New Material Technology Co., Ltd., Taizhou, Jiangsu, PRC) with Young's modulus of 112.4GPa. Compared to the remaining parts of the experimental assembly, the high rigidity of the tendon elements allows us to consider it as a rigid element that does not affect the overall stiffness of the system. For the same reason, the tendons have not been included in the FEA model but are provided by boundary conditions.
A Teflon heat shrink cover has been applied on the whole outer surface of the needle body, except for the joint and the needle tip part, separating the tendons and tissue during the insertion. This configuration may minimize any variation in the smooth friction behavior between the tendons and the needle body and the tendon kinematics.
Then, the two tendons are connected to the needle tip on one side and to the tension-force sensing structures on the other side (Fig. 5b), where four strain gauges (GFLAB-3-50-5LJC-F, Tokyo Measuring Instruments Lab, Tokyo, Japan) are installed in a Wheatstone bridge configuration. The sensor structures are fabricated with the same material, deposition strategy, and curing conditions employed for the needle.

B. TISSUE PHANTOMS
A bovine gelatin powder (Sammi Gelatin, SAMMI INDUSTRY CO., Republic of Korea) is used to fabricate the two tissue phantoms. The ratio of the powder to the water controls the material softness, and 11.11% and 16.67% are selected for the softer (TP1) and the harder (TP2), respectively. For each ratio, a cylindrical phantom was fabricated and used for compression tests, characterizing its properties, while a parallelepiped block was prepared and used for the needle insertion experiments.
As demonstrated in various recent studies [6,9,14], tissue phantoms are an efficient yet accurate proxy alternative to biological materials. Furthermore, a transparent phantom allows real-time monitoring of the tip pose and needle body shape during the insertion, an important aspect of this study. Thus, it is an essential compromise that does not affect the overall accuracy or reliability of the proposed investigation procedure and developed simulation model.
The compression tests, utilized for the determination of the mechanical properties of the two tissue phantoms employed in this research, were performed at 2 mm/min, 10 mm/min, and 20 mm/min speeds (Fig. 4a). For each speed, five tests were repeated on the Instron 5969 Dual Column Testing System (Instron, Norwood, MA, USA), where the compression load has been recorded with a 2580 load cell (Instron, Norwood, MA, USA) having a maximum load of 1kN and a 0.1% resolution. The stroke was measured by the vertical columns encoder of the testing system. The support plate was cleaned with a rubbing alcohol solution between tests to remove possible contamination that might influence the friction behavior during the test. According to the compression test results reported in Fig. 4b and 4c, the two selected tissue phantom material did not show a strong strain rate dependency thus, strain rate dependency has been neglected to simplify the setting for the simulation model, and a hyperelastic formulation has been employed.
The Ogden strain energy potential formulation, Eq. (1), was used for characterizing the material behavior, already employed for the constitutive modeling of hyperelastic materials, such as the case of tissue phantoms and biological materials, in finite element modeling [23]. In Eq. (1), are the deviatoric principal stretches, calculated according to Eq.
(2) based on the principal stretches ( ) and is the total volume ratio, standing for the change in volume of the considered geometry concerning the initial geometry and accounts for elastic deformations only.

(a) Compression tests of the fabricated tissue phantom and the velocity-dependent experimental engineering stress-strain curves with the Ogden strain energy potential model fitting curves for (b) TP1 and (c) TP2 tissue phantom materials.
The model constants have been estimated by iterative regression analysis to minimize the difference between the engineering stress-curves area integrals relevant for the experimental and constitutive modeling curves. Because the two tissue phantoms have a slight strain rate-dependent behavior, the material constants were calculated for the three compression speeds and then averaged to obtain a single set of model constants in the FEA simulations, as reported in Table  I. The engineering stress-strain curves are shown with the fitting in Fig. 4b (TP1) and Fig. 4c (TP2) and show good agreement for all three experimental results.

C. EXPERIMENTAL SETUP
The experimental set-up, Fig. 5a, has a step motor (PV6 NEMA23 2-phase hybrid stepper model, HANPOSE, Guangzhou, PRC) coupled with an linear guide (SFU1204 model, HANPOSE, Guangzhou, PRC), which moves the whole needle system forward and backward with respect to the tissue phantom. The needle system includes two linear motors pulling and releasing the tendons, tension-force sensing structures, and the needle assembly. The needle is secured into position through a connecting plate, fastened with four screws and, on its back, a small-scale load cell (FX29K0-100A-0010-L, TE Connectivity, Schaffhausen, Switzerland), with a resolution equal to 50N ± 0.5%, is installed to record the compression force at the needle base during the tests.
The steering is applied utilizing two linear motors connected to the sensor structures that transfer the motion to the tendons and the needle tip (Fig. 5b). Note that for the passive needle, these motors were not used. The correlation between the digital signal from the Arduino interface connected to the Wheatstone bridge and the applied load has been established through pre-tests. The calibrated loads, from 0.2N to 7.1N, have been connected to the tendons, resulting in the calibration chart of Fig. 5c. The time-dependent position of the needle has been recorded through a high-resolution camera with a sampling rate of 50Hz. The recordings were analyzed using a feature recognition algorithm, implemented in MATLAB (R2021b, MathWorks, MA, USA) [24], allowing the calculation of the angle of the needle tip based on the tension applied on the tendons by the linear motors.

D. TEST PROCEDURE
Three different sets of experiments were done, one with the passive needle (E1) and the other two with the active needle (E2, E3). The E1 is to find the friction coefficient between the needle and the tissue, required for the FEM simulations and validate the model based on the generated stress field ahead of the tip and the resultant compressive force measured at the needle base by the end-point load cell. The E2 measures the global friction between the tendons and the tendon grooves, and the tip joint because these frictions are not included in the FEA model but the measured values by the tension-force sensors. The experiments were done with the tip deflected to a similar level of the tissue experiments and steered in the air. Although the deflection varies during insertion, the E2 with a fixed deflection will provide a rough but good estimate of the value for the relatively small tip deflection during the experiments. The E3 is the needle insertion experiment using the two tissue phantoms in two sub-steps. First, a pretension displacement is applied to both tendons to ensure alignment with the needle body, followed by the straight insertion into the tissue phantoms (Condition I). Second, the tip is steered upward (Condition II) and downward (Condition III), from a static position, without further insertion. The estimated global

Fig. 6 Experimental setup for validation: (a) Entire set-up with details of the sensor devices and needle hinge joint region, (b) assembly between the tendon force sensor and the linear motor, (c) calibration between the digital output of the tendon force sensors, and (d) simple insertion in soft tissue phantom (TP1)
VOLUME XX, 2017 9 friction in E2 and the pre-tension forces are subtracted from the measured tension force in E3, resulting in the pure steering force reacting to the tissue interaction at the tip of the active needles. The friction-decoupled steering forces are compared with the FEA simulation results to evaluate the model accuracy. Regarding the steering forces in the FEA simulation, the contact pressure on the whole top and bottom surface of the needle tip has exported from the results and multiplied by its relevant half areas calculated from the CAD design (SolidWorks 2021, Dassault System, France) [25].
The experimental images were also compared with the FEA contour plot of the von Mises equivalent stress gradient resulting from the relevant FEA model. For both E1 and E3, the insertion speed in the experiments and simulations was set to 1mm/s.

IV. RESULTS AND DISCUSSION
The results of E1 and E3 after the removal of global friction estimated in the E2 and the pre-tension forces, are presented in Fig 6 and 7, respectively, which are discussed below for the friction measurement, the passive needle insertion, and the active needle steering.

A. FRICTION COEFFICIENT
According to the E1 result (Fig. 7a), the compressive reaction force at the needle base matches well the FEA simulation with the friction coefficient of 0.1 among the five values tested (0.01, 0.05, 0.1, 0.18, 0.25), regardless of the tissue phantoms. With the value, the maximum deviation of the compressive force between experiments and the simulation are 8.7% and 8.9% for TP1 and TP2, respectively, and thus, selected for all the simulations hereafter. So naturally, an increase in the friction coefficient results in a higher compressive reaction force on the needle body, but the overall shapes are pretty similar. The first load drop is caused by the initial elastic deformation of the outer tissue surface during the initial contact with the needle tip. After the needle tip creates the first crack on the tissue surface, the following penetration is characterized by a relaxation of the tissue, represented by the faint drop in the load. This result indicates that the friction interaction between the needle body and the cavity created within the tissue is a linear correlation, where the friction coefficient acts as a proportional factor for the compressive base force. Although the selected 0.1 is smaller than the 0.18 used in [19], further decrease is desired to minimize tissue damage, for which the developed FEM model would be useful.

B. PASSIVE NEEDLE DEFLECTION AND POTENTIAL FRACTURE PROPAGATION
The passive needle experiences a slight downward deflection, primarily because the initial orientation of the tip was not perfectly normal to the outside surface of the tissue phantom. On the other, no boundary conditions limiting the z-direction (vertical in the figures) displacement have been introduced in the implementation of the FEA model. Thus, the slightly asymmetric stress field around the needle would be caused by a small deflection of the needle body, partly attributed to the contact conditions and meshes used. It is also strongly related to the absence of a predetermined crack path in the developed model, ensured by utilizing the CEL strategy. Nevertheless, the deflection difference is pretty small, resulting in the almost same y-position of the tip between the experiments and FEM simulation for three different insertion depths (l1, l2, and l3 in Fig. 7b-7d), whose maximum deviation is only 0.81% with respect to the experiment results.
The deformation field around the needle during the insertion shows good agreement with the von Mises equivalent stress

Fig. 6 Passive needle experiments and FEA simulations at the insertion speed of 1mm/s: (a) friction coefficient estimation based on the base force measurement, (b)-(d) insertion experiments in TP1 and von Mises equivalent stress fields obtained by FEA simulation for different insertion depth, e) base force comparison between the FEA and the experiments for two different phantoms.
gradient calculated by the implemented FEA model (Fig. 7b-7d), where a slightly higher stress field at the deflected side (bottom in the figure).
Two different scales for the stress are presented: the maximum values equal to the failure stress of TP1 (Top) and a lower maximum scale value showing a comprehensive plot of the stress distribution around the needle body (Bottom). It is also interesting that the asymmetric gradient of the von Mises stress field seems to be a function of the insertion length, a fact that strengthens the slight deflection of the needle body during the insertion.
This passive needle result validates the generality of the proposed FEA model. The model is developed for active needles with a rotational tip joint, but it can be applied to passive and active needles by simply switching the steering boundary conditions and the needle body geometry. In addition, although a detailed investigation is not the scope of this study, an intermediate type (e.g., an active needle with a symmetric tip) would yield an intermediate stress field between those of a passive needle with a symmetric tip (refer to Fig. 6) and those of an active needle with an asymmetric bevel tip (refer to Fig. 7), another use scenario of the model.

C. ACTIVE NEEDLE DEFLECTION AND STEERING FORCE COMPARISON
The results of the active needle experiments are accompanied with the post-processed images using the feature recognition, tip-angle estimation included, and the von Mises stress fields around the tip. Although the needle body and tip contours are identifiable, a few white residual pixels could not be filtered out during the post-process. However, they do not affect the penetration length and steering angle calculations. Fig. 7a shows the needle inserted with only the pre-tension applied to both tendons, and thus, the deflection is caused only by the asymmetric shape of the needle tip, detail of Fig. 7b. Fig. 7c-7d present the subsequential steps in the order, the upward steering and the downward steering (Refer to Fig. 2c).
First of all, the overall tip deviation shows a relatively large error: the experimental result of 6.05mm while the FEA results of 2.58mm. One of the error sources would be the Teflon cover that was not accounted for in the FEM model but creates a concave surface near the joint at the bottom, as clearly shown in the post-processed images. Such a surface increases the drag, especially at the lower part of the needle tip, further enhancing the effect of the asymmetric tip and the upward deflection. On the other hand, the stress is relieved near the joint compared to the tip in the FEA simulation, implying no such increased drag. Another source is the unavoidable mechanical flexibility existing at every component from the steering motor to the tip joint, including the tendon, which allows undesired tip rotation counterclockwise due to the tissue reaction force during the insertion. As a result, the joint angle is not a fixed zero during insertion. However, it is likely to be a positive value that is not possible in the FEA model where a theoretical geometry is considered. The reduced tip angle in tissue (Fig. 7c-7d) compared to those in the air, upward (from 20° to 17.1°) and downward (from 23° to 20.9°), also proves the mechanical flexibility. Further investigations are necessary to reduce the gap, especially in the needle design of both the FEA model and the prototype.
On the other hand, as clearly shown in Fig. 7e, the FEApredicted steering forces accurately match the experimental results during the entire motion even with the rough estimation of global friction, regardless of the direction and tissue

Fig. 7 Active needle experiments and FEA simulations at the insertion speed of 1mm/s: (a) full insertion TP1 and tip deviation in TP1, (b)-(d) full insertion and tip orientation with the von Mises stress field around the tip, (e) steering force comparisons between the experiments and the FEA simulation.
phantoms. The maximum differences for upward (condition II, Fig. 7) and downward (condition III, Fig. 7) steering are 11.7% for TP1 and 12.9% for TP2, respectively. The direction-dependent steering force is due to the asymmetric tip geometry that dominates the interaction, well predicted in the FEA simulation. It also should be noted that the aforementioned concave surface and flexibility issues in the overall tip deviation are not the case for this steering force.
This active needle result validates the effectiveness of the FEA model for accurately investigating the steering force, i.e., tip-tissue interaction of active needles. The validated model can advance current steerable needle research by allowing quantitative and controlled investigation of the other aspects relevant with the interaction in an efficient way. As an example, a case study is conducted and presented in the following section.

V. A CASE STUDY: TIP GEOMETRY AND TISSUE DAMAGE
In steerable needle research, a model-based approach is efficient and essential for investigating the effects of design parameters on tissue, especially when there is no appropriate sensor directly measuring the effect. Tissue damage assessment during insertion is one of such situations because, without a model, it can only be inferred from steering forces experimentally measured. Thus, the effect of the tip geometry on tissue damage is investigated as a case study to demonstrate the usefulness of the proposed and validated model toward tip design optimization in the future, using the geometry and steering parameters summarized in Table II of Appendix A. Although the effect of the bevel angle has already been investigated for passive needles [8], the other geometry parameters have not been studied and for active needles. Since the hyperelastic behavior is normally characterized by sudden failure stress remarkably close to the ultimate tensile strength of the material for biological materials, the ratio of the von Mises equivalent stress in the tissue to the ultimate strength of the material is used as a metric of tissue damage. One should note that the stress can not be measured experimentally.
The results relevant for both minimum and maximum damages are shown in Fig. 8a-8b for both TP1 and TP2, as well as the 3D plots comprising the summary of the results and the interpolation planes in Fig. 8c-8d. Since all analyses considered the penetration depth of 30 mm at the same insertion speed of 1 mm/s for the same upward steering, the results can be directly compared. The only variables involved are the needle tip geometry and the hyperelastic tissue properties.
The overall results for both TP1 and TP2 share the same trend for the ratio of needle tip length to the diameters and the bevel angle. Also, the influence of the ratio is far more pronounced than that of the bevel angle. The higher the ratio, the more severe tissue damage is, probably due to the two factors below.

Fig. 8 Case study results: (a) maximum and minimum damage cases in TP1, (b) Maximum and minimum damage cases in TP2, (c) tissue damage in TP1 as a function of the ratio and the bevel angle, (d) tissue damage in TP2 as a function of the ratio and the bevel angle.
A long and blunt needle increases the relevant contact area and the force transferred to the tissue, making the needle tip have to displace more tissue to achieve the desired steering angle. Further, the low sharpness of the needle tip edge makes it difficult for the crack to reorient itself and follow the new propagation trajectory imposed by the steering. Combining these two factors creates an uneven force distribution in the tissue, resulting in a stress gradient and, ultimately, a higher damage distribution in the soft tissue cavity. Since TP1 is softer than TP2, the effect of this uneven stress gradient results in more significant damage in TP2. In addition, such a long and blunt tip leads to a significant and continuous variation of the contact conditions between the needle and the tissue during the steering operation, thus resulting in a higher damage gradient.
It is also interesting to highlight that the 30° bevel angle allows a clear drop in tissue damage. The crack ahead of the needle tip primarily develops according to a mode-I propagation; the 30° bevel angle allows for a smoother change in the crack propagation direction, leading to lower tissue damage. For the same reason, contrary to the results of the other bevel angles of 45° and 60°, the combination of a 30° bevel angle with a short tip allows minimizing the maximum damage.

VI. CONCLUSIONS
The research presented in this paper detailed a finite element analysis (FEA) model to investigate the physical interaction between steerable needles and soft tissues during medical examinations. Since the model is based on the CEL method that does not require predetermined insertion paths and remesh steps, the proposed model and approach are suitable for studying active needles that change the shape according to the steering inputs. However, they are also applicable to various needles by simply updating the boundary conditions, even passive needles without tip steering. An experimental set-up was designed, implemented, and utilized to validate the developed FEA model for both passive and active needles, showing good accuracy and reliability for both designs in two different tissue phantoms, especially for the steering forces. The validated model was then applied to the research on the tip geometry, concluding that the ratio between the needle tip length and its diameter is the dominant factor in controlling the damage gradient on the tissue. Although the potential has been demonstrated for estimating key parameters for passive and active needles, such as design optimization and defining the control parameters, there are still room for improvements in the future. Especially, the model will be validated in a more realistic tissue environment that may have non-negligible viscoelastic deformation and viscous friction.