Modelling and real-time dynamic simulation of flexible needles for prostate biopsy and brachytherapy

ABSTRACT

Percutaneous needle insertion constitutes a widely adopted technique for performing minimally invasive operations. Robot-assisted needle placement and virtual surgical training platforms have the potential to significantly improve the accuracy of these operations. For this, the development of mathematical models that provide a complete characterization of the underlying dynamics of medical needles is considered of paramount importance. In this paper, we develop two three-dimensional nonlinear rigid/flexible dynamic models of brachytherapy and local anaesthetic transperineal biopsy (LATP) needles. The proposed models relax the assumptions of previous investigations, quantify the vibrational behaviour and the rigid-body dynamics of medical needles and allow for real-time haptic and visual feedback information. Their accuracy and computational efficiency are assessed and validated using commercial software. The results show that, among the examined methods, the Rigid Finite Element Method provides the most accurate and numerically efficient solution for capturing the dynamics of flexible medical needles.

KEYWORDS:

• Surgery simulation
• needle insertion
• flexible multibody dynamics
• rigid finite element method

1. Introduction

In the last few decades, minimally invasive surgery (MIS) and localized therapy have become integral parts of modern medical practices as they offer significant advantages over conventional open surgery, such as decreased recovery time, lower risk of infection and reduced patient discomfort [1]. One of the most widely used MIS procedures is percutaneous needle insertion, with a plethora of diagnostic and therapeutic applications, such as neurosurgery, deep brain stimulation, epidural anaesthesia, tissue biopsy and prostate brachytherapy [2].

The success of the aforementioned operations is highly dependent on the accuracy of needle placement, while imprecise targeting often leads to severe complications, such as false negatives in biopsy or ablation of healthy tissue [1]. However, accurate percutaneous needle placement is a highly challenging task. The limited visual feedback during the operation, combined with factors such as tissue anisotropy, heterogeneity and variability in anatomical structures among patients make navigation through tissue particularly complicated and decrease the procedure’s overall accuracy [3].

In recent years, there has been significant effort to provide robust solutions for increasing the accuracy of needle placement procedures. Research has mainly focused on the development of robotic systems that allow autonomous or semi-autonomous navigation and accurate needle placement to targeted locations inside soft tissue [4,5]. Furthermore, the development of high-fidelity visual/haptic surgical simulators has allowed the training of junior surgeons in a variety of surgical scenarios and the preoperative planning of complex procedures by experienced doctors, leading to an overall improvement of the needle’s placement accuracy [2,6,7].

A key requirement for the implementation of robotic guidance or simulated training solutions is the formulation of mathematical models that can accurately describe the dynamics of needle insertion. A common practice found in the literature is the division of the needle insertion modelling problem into three separate components [8]: [a)] the modelling of the soft tissue [9–11], the characterization of the interaction forces between the needle and the surrounding medium [1,12,13] and the modelling of the flexible needle [8,14,15]. This division allows the development of decoupled mathematical models for the needle and the tissue, while needle–tissue interaction is often modelled based on experimentally derived generalized force profiles at the needle’s base during needle insertion procedures [1].

This work focuses on the formulation of accurate and computationally efficient models of flexible medical needles. These are of paramount importance as they not only provide an accurate description of the needle’s behaviour but also allow for the minimization of the stochasticity that is present during the aforementioned experimentally derived models of needle–tissue interaction. As a result, accurate needle models can significantly improve the identification accuracy of needle–tissue interaction dynamics, and, thus, constitute a key prerequisite for the systematic modelling of needle insertion procedures.

1.1. Related work

A complete characterization of flexible medical needles aims to identify the relationship between the spatio-temporal behaviour of the needle’s geometry and the specified motion trajectory at its base (e.g. insertion velocity, axial and lateral rotation), under the effect of some arbitrary loading conditions acting on its shaft. Furthermore, complete needle models provide information about the reaction forces at the needle’s base, due to the combined effect of the input driving constraints and the input loading conditions. The following paragraphs provide a review of various computational methods that have been developed to tackle some or all of these modelling requirements.

In the study by Webster et al. [16], a novel non-holonomic model, based on a variation of the kinematic bicycle model, was proposed. The model allowed for a computationally efficient estimation of the needle’s deflection given the insertion and rotation speed of its base. However, the model did not provide any information about the vibrational behaviour of the needle and the profiles of the associated reaction forces. Furthermore, due to its various modelling assumptions, the method led to significant deviations when compared with experimental studies [14]. Extensions of this model have been used in the literature for the development of kineto-dynamic path planning algorithms [17,18].

The need for more accurate needle models and information about the dynamic properties of needle insertion has led to a transition from kinematic to mechanics-based modelling approaches. Researchers [15,19,20] developed quasi-static needle models, based on the linear Euler–Bernoulli beam theory, to describe needle deflection during percutaneous insertion procedures. While this modelling approach offered simplicity and computational efficiency, its linearity assumption led to invalid solutions when large deflections were considered [8,21]. To account for the geometric nonlinearity and minimize the targeting error of linear modelling approaches, DiMaio and Salcudean [22] and Goksel et al. [8] developed quasi-static methods based on nonlinear 2D and 3D Finite Element Method (FEM). Similarly, Adagolodjo et al. [23] presented a geometrically nonlinear quasi-static model of flexible needles based on FEM and the co-rotational formulation. These models provided both computational efficiency and high accuracy when compared with experimental results for describing the steady-state responses of both large and small needle’s deflections.

Both linear and nonlinear models described above estimate needle’s deflection assuming a quasi-static needle insertion. In other words, it is assumed that the needle is inserted at a rate that does not allow for any oscillatory transients to occur and deflection is approximated by a sequence of steady-state responses. This assumption, however, does not provide any information about the vibrational behaviour of the needle and does not account for the correlation between its deflection and the system’s input parameters, such as its insertion velocity and the rate of axial rotation. This poses significant limitations as needle insertion is a highly dynamic procedure [24]. As shown in the study by Podder et al. [25,26,], during brachytherapy operations, input parameters, such as the insertion velocity, span a wide range of values and considerably affect the needle’s deflection and vibration. More specifically, few studies [20,27–30] have proven the strong coupling of parameters such as the needle insertion speed, vibration and axial rotation to the needle–tissue interaction force and, subsequently, to its deflection. Furthermore, axial needle rotation and vibration can be used for minimizing needle deflection during needle insertion procedures [20,31–33]. Incorporation of such parameters naturally requires a dynamic model of flexible needles and not a quasi-static one. Therefore, accounting for the underlying needle dynamics is crucial for the complete characterization of needle insertion and, subsequently, for the development of accurate simulation solutions and precise control algorithms [24].

Due to these limitations, it is necessary to develop modelling methodologies that incorporate the underlying needle dynamics. Khadem et al. [14] proposed a novel dynamic model of a rigid/flexible 2D needle based on the Rayleigh-Ritz approximation technique and the linear Euler–Bernoulli beam theory. This work was one of the first to introduce the importance of dynamic modelling in needle insertion procedures, but it had significant limitations due to modelling assumptions, such as the beam’s linearity and unmodelled torsional dynamics. In the study by Martsopoulos et al. [34], both the Rayleigh-Ritz and the FEM approaches were employed for capturing the three-dimensional dynamics of a rigid/flexible model of brachytherapy needles. However, the proposed models were based on infinitesimal strain theory and, thus, they did not account for geometric nonlinearities and large needle deflections. Similarly, Haddadi and Hashtrudi-Zaad [24] extended the static angular spring FEM approach of Goksel et al. [8]: the model captured the dynamics of needle insertion, while accounting for geometric nonlinearities due to large deflection, but it was only limited to 2D applications.

1.2. Contributions

Extending the work of Khadem et al. [14,24], Haddadi and Hashtrudi-Zaad [] and Martsopoulos et al. [34], this paper develops and compares two different modelling approaches that aim to provide a complete characterization of the three-dimensional dynamics of flexible medical needles. For this, we employ the theory of flexible multibody dynamics (FMD), as it provides a plethora of computational tools that allow a systematic characterization of the needle’s spatio-temporal behaviour, under arbitrary driving constraints and loading conditions. The two formulations are examined based on three criteria: accuracy, computational efficiency and suitability to provide robust and stable solutions for the particular properties of the examined system (e.g. properties of numerical stiffness and shear locking due to the thinness of the needle are considered). The proposed methods are derived and solved using a systematic approach that allows direct comparisons to be performed. Special attention is paid to the computational efficiency of the methods, with the development of efficient solution strategies and numerical integration schemes, while concurrent execution and parallel programming concepts are also considered. Three simulation scenarios are developed and compared with commercial software solutions.

To the best of the authors’ knowledge, this is the first work to investigate and employ the theory of FMD for modelling the dynamics of flexible medical needles, and to provide a systematic comparison among different modelling techniques, while accounting for the challenges that are unique to needle insertion procedures. It is also the first to investigate both large rigid body motion of the needle’s base and large needle’s deflections and relax the modelling assumptions of previous research works (e.g. linearity assumptions, unmodelled dynamics, 2D models). Furthermore, to the best of the authors’ knowledge, this work is the first to apply the DDM and RFEM methods for modelling the dynamics of flexible needles. Even though approaches similar to the RFEM method have been used by Goksel et al. [8] and Haddadi and Hashtrudi-Zaa [24], there was no explicit reference to the RFEM. Instead, in both cases, the authors proposed the modelling of a flexible needle with a set of rigid bodies interconnected with angular springs. The meshing methodology of the RFEM was not taken into account, while the parameters of the angular springs were based on experimental data and not on the theory of continuum mechanics, as proposed by Wittbrodt et al. [35]. Furthermore, contrary to the RFEM, the springs had only rotational and no translational degrees of freedom.

We believe that this work provides a well-rounded modelling approach without the need for any simplifications that compromise the overall accuracy of the system. Our aim is that the flexibility that this model provides will be used for finding further correlations between the movement of the needle’s base and the forces acting on it during its interaction with the soft tissue. Furthermore, we hope that the computational efficiency of this model will allow its integration with needle–tissue interaction models and the derivation of model-based control algorithms for automating needle insertion procedures.

1.3. Medical needles and flexible multibody dynamics

This section provides an overview of different FMD formulations and presents some essential considerations for choosing the appropriate approach for modelling the dynamics of flexible medical needles. Shabana [36] a detailed review of available techniques of FMD systems theory is presented. The key difference between the methods lies in the selection of the generalized coordinates and the associated coordinate frames that are chosen to describe the system’s configuration. This allows the separation of the available formulations into two main categories: the local and the absolute frame approaches [37].

In the local frame approaches, the deformation field of a flexible body is described with the help of nodal coordinates that are measured with respect to some local frame of reference, while its rigid body motion is defined based on the inertial description of the frame’s position and orientation. Some of the most widely used local frame approaches are the floating frame of reference (FFR), the convected coordinate system (CCS) and the large rotation vector (LRV) [36]. As discussed by Shabana [36], the CCS formulation results in inexact modelling of the rigid body dynamics of flexible structures, while the LRV formulation leads to numerical singularities and locking problems when analysing slender flexible bodies, such as the ones considered in this application.

In the majority of the research in the literature, the FFR formulation is combined with classical linear finite-element methodologies or other linear implementations of the assumed modes method, such as the Rayleigh-Ritz or the Galerkin’s methods. Even though these approaches allow the description of arbitrarily large rigid body motions, they are only valid when small deformations of the flexible components are considered [38, pp. 309]. To address this, there have been various studies on the combination of the FFR formulation with the finite strain theory, to allow the analysis of geometric nonlinearities and large deflections [39–43]. The nonlinearity is added to the system with the addition of the higher-order terms on the flexible body’s elastic potential energy. Even though this methodology allows the straightforward incorporation of the system’s geometric nonlinearity and a computationally efficient solution, it suffers from the numerical problem of element locking, as described by Reddy [44]. While numerical remedies have been proposed to address this problem, such as the reduced-order Gaussian integration [45, pp. 279–280], due to their heuristic nature, they introduce uncertainty to the problem and, thus, are not considered in our work.

To allow the description of both large rigid body motion and large deformation, a new set of methods, known as absolute frame approaches, were developed [36]. Using these techniques, the pose and the deformation field of a flexible body is described using inertially defined nodal coordinates. The most widely used technique of this category is the absolute nodal coordinate formulation (ANCF), proposed by Shabana [46]. For beam elements, the choice of nodal coordinates in ANCF relaxes the assumptions of Euler–Bernoulli and Timoshenko beam theories and allows large deflections and cross-sectional deformations to be considered [47]. While the method requires a larger number of generalized coordinates and leads to highly nonlinear expressions for the system’s generalized elastic forces, it results in a constant mass matrix which significantly simplifies the associated equations of motion and accelerates computations [48].

Whilst ANCF seems to meet the requirements of flexible needle modelling, as defined in Section 1, the method has reportedly led to significant numerical and convergence problems when dealing with thin and stiff structures, such as the ones considered in our work [49–51]. This phenomenon, known as Poisson locking, stems from the coupling between the high-frequency cross-sectional deformation modes (found in thin and stiff elements) and the other deformations (e.g. axial and bending) [49]. To address this, different approaches have been proposed in the literature, such as the development of higher-order thin elements [52], the uncoupling between the shear deformation and antisymmetric bending [50] and the redefinition of the system’s degrees of freedom (DOFs) [53]. While some of the proposed techniques have the potential to alleviate the problem of Poisson locking, they introduce additional modelling complexity [50] or remove the initial benefits of the ANCF, such as the formulations of Gruber et al. [53] and von Dombrowski [54] that lead to a non-constant mass matrix.

As an alternative to the ANCF, our work proposes the modelling of flexible needles with the use of the discrete deformation mode (DDM) formulation, presented by Jonker and Meijaard [55]. A comparative analysis between DDM and a modified (locking-free) version of ANCF in the study by Schwab and Meijaard [50] shows that the DDM (which also belongs to the family of the absolute frame approaches) allows for an accurate and computationally efficient description of large reference motion and deformation of thin flexible elements without introducing further modelling complexity. Furthermore, we also propose the modelling of flexible medical needles with the help of the rigid finite element method (RFEM). This technique, which also belongs to the absolute frame approaches, is well suited for the current application, as it has successfully been applied for capturing the nonlinear dynamics of slender structures [56,57] without leading to the aforementioned numerical problems.

The importance of modelling both large rigid body motion and large needle deformation for prostate biopsy and brachytherapy procedures can be clarified with the help of Figure 1. In these procedures, and especially LATP biopsies that use a precision point transperineal access system [58], the surgeon’s hand performs large rigid body motions in order to guide the needle in the desired location. As illustrated in Figure 1, this can cause significant overall needle deflections, while needle deflections inside the tissue could still remain small.

Figure 1. Large deflections and large rigid body motion in prostate brachytherapy and biopsy.

As a result, it is crucial to provide a mathematical model that captures both large rigid body motion and large needle deflections. The reason for this is twofold:

• One of the possible applications of a high-fidelity needle model is the development of haptic simulators for training junior surgeons in LATP biopsies. It is, thus, important to provide a mathematical model that captures the dynamics of the needle insertion in its entirety without introducing any limiting assumptions. This would give trainees the ability to test even the ‘extreme’ regions of their input space as they would do in a real training environment.

• A highly accurate dynamic model of needle insertion can provide insights into the identification of needle–tissue interaction forces, as it captures the correlations between the needle’s state (deflection) and the imposed trajectory at its base. Introduction of modelling assumptions, such as the quasi-static needle insertion or small needle deflections, does not allow the investigation of these correlations, as they restrict the system’s input and state space. Thus, to understand the real dynamics of needle insertion procedures, it is critical to remove these assumptions and investigate the true correlation between input (imposed trajectory) and output (deflection/forces). By improving our understanding of these correlations, there is the potential to build improved algorithms for robot-assisted needle insertion procedures (e.g. by building better model-based control algorithms).

1.4. Needle–tissue interaction

This section analyses existing needle–tissue interaction models and formulations. As discussed by Rossa and Tavakoli [59], during needle insertion the deformed needle shaft compresses the surrounding tissue which, in turn, applies forces to the former influencing its deflection. As a result, needle deflection and tissue deformation are coupled effects that influence each other [59].

There are two main approaches for identifying these interactions, which, in this work, are referred to as constrained and unconstrained approaches. In the constrained approach, the needle and tissue are studied as a unified system. Needle–tissue interactions are imposed as kinematic constraints between the two mediums, and the system’s equations for both the needle and tissue are solved simultaneously [22,60,61]. The reaction forces between the two bodies are then obtained with the help of LaGrange multipliers. However, this approach presents some shortcomings. First, in cases where the tissue is modelled with the help of FEM, this technique requires that the tissue mesh is adapted at each time step to match the discretization of the needle and allow the enforcement of the kinematic constraints on the corresponding nodes of the tissue and the ones of the needle [22]. The computational cost of this remeshing makes the application of this technique unsuitable for interactive applications. One possible solution is using a meshless discretization, as in the study by Xu et al. [62], or using the technique of local constraints presented by Wang and Hirai [63]. Second, the model requires knowledge of the exact mechanical properties of the tissue on its interface with the needle. However, monitoring the medium’s properties on this interface is not usually plausible or applicable in real-time applications or in vivo tests [64]. Finally, coupling the two flexible mediums leads to a high-dimensional set of equations that is challenging to solve in real-time.

In the unconstrained approach, the needle–tissue interaction is modelled with the help of interface force models. These forces, which are functions of both the needle’s and tissue’s state, are imposed simultaneously but separately on the tissue and the needle, allowing the mathematical decoupling of the two. Furthermore, these models are dependent on a small number of parameters. The unconstrained modelling method is quite popular as it is both simple and computationally efficient, allowing its application to interactive surgery simulation and control. However, it suffers from limited accuracy as the identification of these force profiles is usually based on semi-empirical approaches and experimental observations. A plausible remedy to this is the use of online system identification algorithms for the tuning the model’s parameters [65].

One of the most widely used unconstrained models is the one proposed by Okamura et al. [1]. Experimental studies reported in the study Okamura et al. [1] showed that the axial forces acting on the needle’s base are the summation of stiffness, friction and cutting forces. These forces depend on the relative state (position and velocity) of the tissue and the needle and are distributed along the needle’s shaft. Since the work of Okamura et al. [1], significant research has been devoted to the definition of mathematical models that can accurately capture these force profiles. For example, whilst Okamura et al. [1] proposed a modified Karnopp friction model, different friction modelling approaches were then adopted by Asadian et al. [66] and Yang et al. [67]. Similar, simplified interaction force models have been developed for modelling the distributed lateral loads acting on the needle’s shaft as is inserted inside the tissue. The most popular approach is modelling the lateral tissue forces as virtual springs (both linear and nonlinear) of varying mechanical properties, similar to the Winkler foundation model [19,68–70].

This section provided an overview of existing modelling approaches in needle–tissue interaction. In Section 2, we develop mathematical needle models based on two formulations [a)] the DDM (Section 2.1) and the RFEM (Section 2.2). A simulation algorithm for obtaining the needle’s state is analysed in Section 2.3 while needle–tissue interaction models are formulated in Section 2.4. Real-time considerations are discussed in Section 2.5. Section 3 presents the results from three different simulation scenarios, the interpretation of which is provided in Section 4. Concluding remarks are given in Section 5.

2. Materials and methods

In this work, we focus on moderately flexible medical needles, such as those used in prostate brachytherapy and local anaesthetic transperineal prostate (LATP) biopsy, shown in Figure 2. We propose a model, shown in Figure 2, that considers a) the base of the needle as a rigid body that follows any arbitrary spatial trajectory, imposed by a surgeon or a robotic arm, and the shaft of the needle as a flexible body that is rigidly attached to the base and deforms elastically under a general state/time-dependent field of forces, which corresponds to the interactions between the needle and the surrounding tissue.

Figure 2. Brachytherapy/LATP needle model.

To describe the motion of the system shown in Figure 2, we introduce the following frames of reference: let $F$ be the inertial frame of reference and $f_c$ a frame rigidly attached to the rigid body’s centre of mass $C$. Furthermore, we introduce a floating frame of reference (FFR) $f$ which is rigidly attached to point $A$ of the rigid body. As shown in Figure 2, the orientation of the frames $f_c$ and $f$ is identical and fixed, i.e. there is no relative motion between them.

In the following formulations, the needle is considered as a slender and homogeneous beam with circular cross-section and constant mechanical and geometrical properties along its span. The material of the needle is isotropic and linear. Furthermore, we introduce the notation ${\underset{\_}{r}}_{AP/f_1}^{f_2}$ which describes the position of a vector $\overleftarrow{AP}$ with respect to frame $f_1$ expressed in frame $f_2$ [71, pp. 13–16]. The following sections present the kinematic and dynamic description of the aforementioned needle model based on the two proposed FMD formulations, namely, the DDM and the RFEM approaches. In both approaches, we first introduce the kinematic description of the model and then derive expressions for the virtual works of the forces acting on the structure. Next, using D’Alembert’s form of the principle of virtual work, we express the dynamic equilibrium of the system and derive the equations of motions. These steps are thoroughly presented in Sections 2.1 and 2.2.

2.1. Equations of motion: DDM approach

2.1.1. Kinematic description

The DDM method employs a spatial discretization of the needle’s geometry into $n_e$ beam elements, each of which is described using inertially measured nodal coordinates. More specifically, as shown in Figure 3, a complete characterization of the configuration of the beam element $j\in [1,n_e]$ can be acquired based on the inertially measured cross-sectional position and orientation at its nodes $A_j$ and $B_j$, while intermediate points are defined with the help of interpolation functions. The element’s vector of generalized coordinates can be defined as

Figure 3. Spatial DDM element at reference and deformed state.

(1) ${\underset{\_}{e}}_j={\left[\begin{array}{cccc}{\underset{\_}{r}}_{A_j}^T& {\underset{\_}{\theta }}_{A_j}^T& {\underset{\_}{r}}_{B_j}^T& {\underset{\_}{\theta }}_{B_j}^T\end{array}\right]}^T,$(1)

where ${\underset{\_}{r}}_{A_j}={\underset{\_}{r}}_{O{A}_j/F}^F$ and ${\underset{\_}{r}}_{B_j}={\underset{\_}{r}}_{O{B}_j/F}^F$ are the inertial positions of the nodal cross-sections’ midpoints $A_j$ and $B_j$, while ${\underset{\_}{\theta }}_{A_j}$ and ${\underset{\_}{\theta }}_{B_j}$ are the Euler angles, based on the sequence of intrinsic rotations $z_0$-$y{ }^{\mathrm{\prime }}$-$x{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}$, that define the cross-sectional orientation at nodes $A_j$ and $B_j$. These are defined as

${\underset{\_}{\theta }}_{D_j}={\left[\begin{array}{ccc}{\varphi }_{D_j}& {\theta }_{D_j}& {\psi }_{D_j}\end{array}\right]}^T\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}D=\{A,B\}.$

With the help of the Euler angles ${\underset{\_}{\theta }}_{A_j}$ and ${\underset{\_}{\theta }}_{B_j}$, the set of triad unit vectors, shown in Figure 3, are defined as

${\underset{\_}{e}}_{r_{D_j}}^F=R_{f_{D_j}}^F({\underset{\_}{\theta }}_{D_j}){\underset{\_}{e}}_{r_{D_j}}^{f_{D_j}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}r=\{x,y,z\}\mathrm{a}\mathrm{n}\mathrm{d}D=\{A,B\}.$

The rotation matrices $R_{f_{A_j}}^F$ and $R_{f_{B_j}}^F$ describe the orientation of the local frames $f_{A_j}$ and $f_{B_j}$ with respect to the inertial frame of reference $F$. The ${\underset{\_}{e}}_{x_{A_j}}^{f_{A_j}},{\underset{\_}{e}}_{y_{A_j}}^{f_{A_j}},{\underset{\_}{e}}_{z_{A_j}}^{f_{A_j}}$ and ${\underset{\_}{e}}_{x_{B_j}}^{f_{B_j}},{\underset{\_}{e}}_{y_{B_j}}^{f_{B_j}},{\underset{\_}{e}}_{z_{B_j}}^{f_{B_j}}$ are constant unit vectors pointing in the associated directions $x,y$ and $z$. The inertial position of point $P_j$ of element $j$, illustrated in Figure 3, can be obtained with the help of the generalized element coordinates and interpolation functions ${\varphi }_i(x_j)$ as

(2) ${\underset{\_}{r}}_{O{P}_j/F}^F={\varphi }_1(x_j){\underset{\_}{r}}_{A_j}+{\varphi }_2(x_j){\underset{\_}{e}}_{x_{A_j}}^F+{\varphi }_3(x_j){\underset{\_}{r}}_{B_j}+{\varphi }_4(x_j){\underset{\_}{e}}_{x_{B_j}}^F.$(2)

Expressions for the ${\varphi }_i(x_j),i=\{1,\dots 4\}$, can be chosen as in the study by Jonker and Meijaard [55], while the scalar $x_j$ is the point’s distance along the element’s elastic line measured using the co-rotational frame introduced by Jonker and Meijaard [55]. Based on equation (2)(2) ${\underset{\_}{r}}_{O{P}_j/F}^F={\varphi }_1(x_j){\underset{\_}{r}}_{A_j}+{\varphi }_2(x_j){\underset{\_}{e}}_{x_{A_j}}^F+{\varphi }_3(x_j){\underset{\_}{r}}_{B_j}+{\varphi }_4(x_j){\underset{\_}{e}}_{x_{B_j}}^F.$(2) , the inertial velocity of point $P_j$ can be expressed as

(3) ${\underset{\_}{\dot{r}}}_{O{P}_j/F}^F=Z_j(x_j,{\underset{\_}{e}}_j){\underset{\_}{\dot{e}}}_j,$(3)

where

${\mathit{Z}}_j(x_j,{\underset{\_}{e}}_j)=\left[\begin{array}{cccc}{\varphi }_1(x_j){\mathit{I}}_3& -{\varphi }_2(x_j){\mathit{R}}_{f_{A_j}}^F\mathit{S}({\underset{\_}{e}}_{x_{A_j}}^{f_{A_j}}){\mathit{G}}_{A_j}& {\varphi }_3(x_j){\mathit{I}}_3& -{\varphi }_4(x_j){\mathit{R}}_{f_{B_j}}^F\mathit{S}({\underset{\_}{e}}_{x_{B_j}}^{f_{B_j}}){\mathit{G}}_{B_j}\end{array}\right].$

The matrices $G_{A_j}=G({\underset{\_}{\theta }}_{A_j})$ and $G_{B_j}=G({\underset{\_}{\theta }}_{B_j})$ in the above expression map the derivatives of the Euler angles ${\underset{\_}{\dot{\theta }}}_{A_j}$ and ${\underset{\_}{\dot{\theta }}}_{B_j}$ to the respective frame rotational velocities, such that

${\underset{\_}{\omega }}_{f_{A_j}/F}^{f_{A_j}}={\mathit{G}}_{A_j}{\underset{\_}{\dot{\theta }}}_{A_j}\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{\omega }}_{f_{B_j}/F}^{f_{B_j}}={\mathit{G}}_{B_j}{\underset{\_}{\dot{\theta }}}_{B_j}.$

Based on equation (3)(3) ${\underset{\_}{\dot{r}}}_{O{P}_j/F}^F=Z_j(x_j,{\underset{\_}{e}}_j){\underset{\_}{\dot{e}}}_j,$(3) the inertial virtual displacement of point $P_j$ is given as

(4) $\delta {\underset{\_}{r}}_{O{P}_j/F}^F=\frac{\mathrm{\partial }{\underset{\_}{r}}_{O{P}_j/F}^F}{\mathrm{\partial }{\underset{\_}{e}}_j}\delta {\underset{\_}{e}}_j={\mathit{Z}}_j(x_j,{\underset{\_}{e}}_j)\delta {\underset{\_}{e}}_j,$(4)

while the acceleration of point $P_j$ takes the form

(5) ${\underset{\_}{\ddot{r}}}_{O{P}_j/F}^F={\mathit{Z}}_j(x_j,{\underset{\_}{e}}_j){\underset{\_}{\ddot{e}}}_j+{\underset{\_}{a}}_{v_{f_j}}(x_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j).$(5)

In the above expression, the Coriolis/centrifugal component of acceleration is given as

${\underset{\_}{a}}_{v_{f_j}}(x_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)={\varphi }_2(x_j){\underset{\_}{a}}_{v_f}^{A_j}+{\varphi }_4(x_j){\underset{\_}{a}}_{v_f}^{B_j},$

where

${\underset{\_}{a}}_{v_f}^{D_j}={\mathit{R}}_{f_{D_j}}^F\left({\mathit{S}}^2({\underset{\_}{\omega }}_{f_{D_j}/F}^{f_{D_j}}){\underset{\_}{e}}_{x_{D_j}}^{f_{D_j}}-\mathit{S}({\underset{\_}{e}}_{x_{D_j}}^{f_{D_j}}){\dot{\mathit{G}}}_{D_j}{\underset{\_}{\dot{\theta }}}_{D_j}\right)\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}D=\{A,B\}.$

Similarly, the orientation of the cross-sections along the element’s length can be written as

${\underset{\_}{\theta }}_j=\left[\begin{array}{c}{\varphi }_j\\ {\theta }_j\\ {\psi }_j\end{array}\right]=\left[\begin{array}{c}{\underset{\_}{s}}_{\bar{\varphi }}(x_j){\underset{\_}{\bar{\varphi }}}_j\\ {\underset{\_}{s}}_{\bar{\theta }}(x_j){\underset{\_}{\bar{\theta }}}_j\\ {\underset{\_}{s}}_{\bar{\psi }}(x_j){\underset{\_}{\bar{\psi }}}_j\end{array}\right]\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\left\{\begin{array}{c}{\underset{\_}{\bar{\varphi }}}_j={\left[{\varphi }_{A_j},{\varphi }_{B_j}\right]}^T\\ {\underset{\_}{\bar{\theta }}}_j={\left[{\theta }_{A_j},{\theta }_{B_j}\right]}^T\\ {\underset{\_}{\bar{\psi }}}_j={\left[{\psi }_{A_j},{\psi }_{B_j}\right]}^T\end{array}\right.,$

where expressions for the interpolation vector functions ${\underset{\_}{s}}_{\bar{\varphi }}(x_j)$, ${\underset{\_}{s}}_{\bar{\theta }}(x_j)$ and ${\underset{\_}{s}}_{\bar{\psi }}(x_j)$ can be found in the study by Schwab and Meijaard [50]. If the mapping from the nodal components of the Euler angles to the generalized element’s coordinates is defined as ${\bar{r}}_j=L_{\bar{r}}{\underset{\_}{e}}_j$, with $r=\{\varphi ,\theta ,\psi \}$, then the above expressions take the form

(6) ${\underset{\_}{\theta }}_j={\mathit{S}}_{\theta }(x_j){\underset{\_}{e}}_j\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{S}}_{\theta }(x_j)=\left[\begin{array}{c}{\underset{\_}{s}}_{\bar{\varphi }}(x_j){\mathit{L}}_{\bar{\varphi }}\\ {\underset{\_}{s}}_{\bar{\theta }}(x_j){\mathit{L}}_{\bar{\theta }}\\ {\underset{\_}{s}}_{\bar{\psi }}(x_j){\mathit{L}}_{\bar{\psi }}\end{array}\right].$(6)

2.1.2. Virtual work

Next, we derive the virtual work of the inertial forces of element $j$. To simplify the associated expressions, it is assumed that the inertial forces caused by the transversal shear deformation of the element’s cross-section are negligible due to the needle’s thickness. Based on this assumption, the virtual work of the inertial forces of element $j$ can be defined as

$\delta W_{f_j}^{(in)}={\int }_{m_{f_j}}({\underset{\_}{\ddot{r}}}_{O{P}_j/F}^F{)}^T\delta {\underset{\_}{r}}_{O{P}_j/F}^F\mathrm{d}{m}_{f_j}+{\int }_0^{l_{e_j}}\rho I_{p_j}{\ddot{\varphi }}_j\delta {\varphi }_j\mathrm{d}{x}_j,$

where $m_{f_j}$ is the mass, $l_{ej}$ the length, $I_{p_j}$ the polar moment of area and $\rho$ the density of element $j$. The second term of the above equations employs the first component of the sequence of Euler angles ${\varphi }_j$ to describe the inertial forces due to the element’s torsional dynamics. Using equations (4)(4) $\delta {\underset{\_}{r}}_{O{P}_j/F}^F=\frac{\mathrm{\partial }{\underset{\_}{r}}_{O{P}_j/F}^F}{\mathrm{\partial }{\underset{\_}{e}}_j}\delta {\underset{\_}{e}}_j={\mathit{Z}}_j(x_j,{\underset{\_}{e}}_j)\delta {\underset{\_}{e}}_j,$(4) , (5(5) ${\underset{\_}{\ddot{r}}}_{O{P}_j/F}^F={\mathit{Z}}_j(x_j,{\underset{\_}{e}}_j){\underset{\_}{\ddot{e}}}_j+{\underset{\_}{a}}_{v_{f_j}}(x_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j).$(5) ) and (6(6) ${\underset{\_}{\theta }}_j={\mathit{S}}_{\theta }(x_j){\underset{\_}{e}}_j\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{S}}_{\theta }(x_j)=\left[\begin{array}{c}{\underset{\_}{s}}_{\bar{\varphi }}(x_j){\mathit{L}}_{\bar{\varphi }}\\ {\underset{\_}{s}}_{\bar{\theta }}(x_j){\mathit{L}}_{\bar{\theta }}\\ {\underset{\_}{s}}_{\bar{\psi }}(x_j){\mathit{L}}_{\bar{\psi }}\end{array}\right].$(6) ), the above expression can be rewritten as

(7) $\delta W_{f_j}^{(in)}=({\mathit{M}}_{f_j}{\underset{\_}{\ddot{e}}}_j-{\underset{\_}{f}}_{v_{f_j}}{)}^T\delta {\underset{\_}{e}}_j,$(7)

where ${\mathit{M}}_{f_j}={\mathit{M}}_{f_j}^t+{\mathit{M}}_{f_j}^r$. The mass matrices ${\mathit{M}}_{f_j}^t$ and ${\mathit{M}}_{f_j}^r$ are defined as

${\mathit{M}}_{f_j}^t({\underset{\_}{e}}_j)={\int }_{m_{f_j}}{\mathit{Z}}_j^T(x_j,{\underset{\_}{e}}_j){\mathit{Z}}_j(x_j,{\underset{\_}{e}}_j)\mathrm{d}{m}_{f_j}\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{M}}_{f_j}^r={\mathit{L}}_{\bar{\varphi }}^T\rho I_{p_j}\left({\int }_0^{l_{e_j}}{\underset{\_}{s}}_{\bar{\varphi }}^T{\underset{\_}{s}}_{\bar{\varphi }}\mathrm{d}{x}_j\right){\mathit{L}}_{\bar{\varphi }}.$

The generalized Coriolis/centrifugal vector takes the form

${\underset{\_}{f}}_{v_{f_j}}({\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)=-{\int }_{m_{f_j}}{\mathit{Z}}_j^T(x_j,{\underset{\_}{e}}_j){\underset{\_}{a}}_{v_{f_j}}(x_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)\mathrm{d}{m}_{f_j}.$

The key component of the DDM method is the deformation modes, which provide a complete characterization of the beam element’s deformation field. These modes can be specified using six independent generalized strain expressions, which are functions of the element’s inertially defined generalized coordinates [72]. The generalized strains, which remain invariant under arbitrary rigid body motions and orthogonal coordinate transformations [72], take the form

(8) ${\underset{\_}{\overline{ϵ}}}_j=\left[{\overline{ϵ}}_{1_j}{\overline{ϵ}}_{2_j}{\overline{ϵ}}_{3_j}{\overline{ϵ}}_{4_j}{\overline{ϵ}}_{5_j}{\overline{ϵ}}_{6_j}\right]=\left[{\overline{l}}_j-l_{e_j}\left({({\underset{0.3em\$\$}{\underset{\_}{e}}}_{z_{A_j}}^F)}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{y_{B_j}}^F}-{({\underset{0.3em\$\$}{\underset{\_}{e}}}_{y_{A_j}}^F)}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{z_{B_j}}^F}\right)/2-l_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{z_{A_j}}^F}{l}_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{z_{B_j}}^F}{l}_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{y_{A_j}}^F}-l_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{y_{B_j}}^F}\right],$(8)

where ${\bar{l}}_j=\parallel {\underset{\_}{r}}_{B_j}-{\underset{\_}{r}}_{A_j}\parallel$ and ${\underset{\_}{e}}_l^j=({\underset{\_}{r}}_{B_j}-{\underset{\_}{r}}_{A_j})/{\bar{l}}_j$. In equation (8)(8) ${\underset{\_}{\overline{ϵ}}}_j=\left[{\overline{ϵ}}_{1_j}{\overline{ϵ}}_{2_j}{\overline{ϵ}}_{3_j}{\overline{ϵ}}_{4_j}{\overline{ϵ}}_{5_j}{\overline{ϵ}}_{6_j}\right]=\left[{\overline{l}}_j-l_{e_j}\left({({\underset{0.3em\$\$}{\underset{\_}{e}}}_{z_{A_j}}^F)}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{y_{B_j}}^F}-{({\underset{0.3em\$\$}{\underset{\_}{e}}}_{y_{A_j}}^F)}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{z_{B_j}}^F}\right)/2-l_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{z_{A_j}}^F}{l}_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{z_{B_j}}^F}{l}_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{y_{A_j}}^F}-l_{e_j}{(\underset{0.3em\$\$}{{\underset{\_}{e}}_i^j})}^T\underset{0.3em\$\$}{{\underset{\_}{e}}_{y_{B_j}}^F}\right],$(8) , the generalized strain ${\overline{ϵ}}_{1_j},$ is used to describe the beam’s axial elongation, ${\overline{ϵ}}_{2_j}$ is for torsion, while ${\overline{ϵ}}_{3_j},\dots ,{\overline{ϵ}}_{6_j}$ describe transversal bending. To account for geometric nonlinearities, second-order terms can be added in (8). The modified second-order generalized strains can be calculated as

(9) ${\underset{\_}{ϵ}}_j={\underset{\_}{\overline{ϵ}}}_j({\underset{\_}{e}}_j)+\tilde{\omega }({\overline{ϵ}}_j({\underset{\_}{e}}_j))={\underset{\_}{\eta }}_j({\underset{\_}{e}}_j,$(9)

where expressions for the function $\tilde{\omega }({\overline{ϵ}}_j({\underset{\_}{e}}_j))$, which incorporates the second-order terms, and the function ${\underset{\_}{\eta }}_j({\underset{\_}{e}}_j)$, which maps the generalized coordinates to generalized strain measures, are given by Jonker and Meijaard [55]. Furthermore, if it is assumed that the beam has been divided into a sufficient number of elements so that the strains remain small, then the generalized strains can be associated with the generalized stresses (Voigt notation) using the constitutive law written as [50]

(10) ${\underset{\_}{\sigma }}_j=S{ }^j{\underset{\_}{ϵ}}_j,$(10)

where $\mathit{S}{}^j$ is the element’s constant symmetric stiffness matrix, defined by Schwab and Meijaard [50]. Based on equation (9)(9) ${\underset{\_}{ϵ}}_j={\underset{\_}{\overline{ϵ}}}_j({\underset{\_}{e}}_j)+\tilde{\omega }({\overline{ϵ}}_j({\underset{\_}{e}}_j))={\underset{\_}{\eta }}_j({\underset{\_}{e}}_j,$(9) , the infinitesimal change of the generalized strains can be calculated as

(11) $\delta \underset{0.3em\$\$}{{\underset{\_}{ϵ}}_j}=J_{{\eta }_j}\left({\underset{\_}{e}}_j\right)\delta {\underset{\_}{e}}_j\text{where}{J}_{{\eta }_j}{\underset{\_}{e}}_j=\frac{\partial {\underset{0.3em\$\$}{\underset{\_}{\eta }}}_j}{\partial {\underset{0.3em\$\$}{\underset{\_}{e}}}_j}.$(11)

With the help of equations (9)(9) ${\underset{\_}{ϵ}}_j={\underset{\_}{\overline{ϵ}}}_j({\underset{\_}{e}}_j)+\tilde{\omega }({\overline{ϵ}}_j({\underset{\_}{e}}_j))={\underset{\_}{\eta }}_j({\underset{\_}{e}}_j,$(9) , (10(10) ${\underset{\_}{\sigma }}_j=S{ }^j{\underset{\_}{ϵ}}_j,$(10) ) and (11(11) $\delta \underset{0.3em\$\$}{{\underset{\_}{ϵ}}_j}=J_{{\eta }_j}\left({\underset{\_}{e}}_j\right)\delta {\underset{\_}{e}}_j\text{where}{J}_{{\eta }_j}{\underset{\_}{e}}_j=\frac{\partial {\underset{0.3em\$\$}{\underset{\_}{\eta }}}_j}{\partial {\underset{0.3em\$\$}{\underset{\_}{e}}}_j}.$(11) ), the virtual work of the element’s elastic forces can be written as

(12) $\delta W_{f_j}^{\left(el\right)}=-\underset{{\underset{\_}{\sigma }}_j^T}{\delta }{\underset{\_}{ϵ}}_j=Q_{e{l}_j}^T\delta {\underset{\_}{e}}_j\text{where}{\underset{\_}{Q}}_{e{l}_j}=-J_{{\eta }_j}^T{\underset{\_}{e}}_{j}Sj{\underset{\_}{\eta }}_j{\underset{\_}{e}}_j.$(12)

Similarly, the virtual work of the element’s damping forces can be written as $\delta W_{f_j}^{(d)}={\underset{\_}{Q}}_{d_j}^T\delta {\underset{\_}{e}}_j$, where the generalized damping forces ${\underset{\_}{Q}}_{d_j}$ can be modelled using various approaches, as discussed by Yoo et al. [73].

Next, we define the virtual work for the external forces. For this, it is assumed that a random point $P_j^i$ of element $j$, located at the beam’s elastic line and at a distance $x_j^i$ from the nodal point $A_j$, is subjected to a point force ${ }^j{\underset{\_}{F}}_i^F(t,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)$ and a point moment ${ }^j{\underset{\_}{M}}_i^F(t,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)$, shown in Figure 4a.

Figure 4. Free body diagrams for DDM and RFEM methods.

Then, the virtual work resulting from these forces and moments can be written as

$\delta W_{f_j}^{(ex)}={(}^j{\underset{\_}{F}}_i^F{)}^T\delta {\underset{\_}{r}}_{O{P}_j^i/F}^F+{(}^j{\underset{\_}{M}}_i^F{)}^T\delta {\underset{\_}{\theta }}_j^i.$

Using equations (4)(4) $\delta {\underset{\_}{r}}_{O{P}_j/F}^F=\frac{\mathrm{\partial }{\underset{\_}{r}}_{O{P}_j/F}^F}{\mathrm{\partial }{\underset{\_}{e}}_j}\delta {\underset{\_}{e}}_j={\mathit{Z}}_j(x_j,{\underset{\_}{e}}_j)\delta {\underset{\_}{e}}_j,$(4) and (6(6) ${\underset{\_}{\theta }}_j={\mathit{S}}_{\theta }(x_j){\underset{\_}{e}}_j\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{S}}_{\theta }(x_j)=\left[\begin{array}{c}{\underset{\_}{s}}_{\bar{\varphi }}(x_j){\mathit{L}}_{\bar{\varphi }}\\ {\underset{\_}{s}}_{\bar{\theta }}(x_j){\mathit{L}}_{\bar{\theta }}\\ {\underset{\_}{s}}_{\bar{\psi }}(x_j){\mathit{L}}_{\bar{\psi }}\end{array}\right].$(6) ), the above expression takes the form

(13) $\delta W_{f_j}^{(ex)}={(}^j{\underset{\_}{F}}_i^F{)}^T{\mathit{Z}}_j^i\delta {\underset{\_}{e}}_j+{(}^j{\underset{\_}{M}}_i^F{)}^T{\mathit{S}}_{\theta }^i\delta {\underset{\_}{e}}_j,$(13)

where ${\mathit{Z}}_j^i={\mathit{Z}}_j(x_j^i,{\underset{\_}{e}}_j)$ and ${\mathit{S}}_{\theta }^i={\mathit{S}}_{\theta }(x_j^i)$. The virtual work of the external generalized forces presented in equation (13)(13) $\delta W_{f_j}^{(ex)}={(}^j{\underset{\_}{F}}_i^F{)}^T{\mathit{Z}}_j^i\delta {\underset{\_}{e}}_j+{(}^j{\underset{\_}{M}}_i^F{)}^T{\mathit{S}}_{\theta }^i\delta {\underset{\_}{e}}_j,$(13) can be extended to accommodate for $n_f$ point forces and $n_m$ point moments. In this case, it can be proven that the total virtual work of the external generalized forces is given as

(14) $\delta W_{f_j}^{(ex)}={\underset{\_}{Q}}_{f_j}^T\delta {\underset{\_}{e}}_j\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\underset{\_}{Q}}_{f_j}=\sum _{i=1}^{n_f}({\mathit{Z}}_j^i{)}^T{(}^j{\underset{\_}{F}}_i^F)+\sum _{k=1}^{n_m}({\mathit{S}}_{\theta }^k{)}^T{(}^j{\underset{\_}{M}}_k^F).$(14)

Finally, the virtual work of the external generalized constraint forces for element $j$ can be calculated as $\delta W_{f_j}^{(c)}={\underset{\_}{Q}}_{c_j}^T\delta \underset{\_}{\bar{q}}$.

2.1.3. Dynamic equilibrium

To derive the flexible needle’s equations of motion we apply the dynamic equilibrium expression based on D’Alembert’s form of the principle of virtual work, which can be written as

$\delta W_f^{(in)}=\delta W_f^{(el)}+\delta W_f^{(d)}+\delta W_f^{(ex)}+\delta W_f^{(c)}.$

The virtual works of the above expression can be calculated considering the contributions from the individual elements as

$\sum _{j=1}^{n_e}\delta W_{f_j}^{(in)}=\sum _{j=1}^{n_e}\left(\delta W_{f_j}^{(el)}+\delta W_{f_j}^{(d)}+\delta W_{f_j}^{(ex)}+\delta W_{f_j}^{(c)}\right).$

Substituting the derived virtual work expressions, the above equation becomes

$\sum _{j=1}^{n_e}{\left({\mathit{M}}_{f_j}{\underset{\_}{\ddot{e}}}_j-{\underset{\_}{f}}_{v_{f_j}}-{\underset{\_}{Q}}_{e{l}_j}-{\underset{\_}{Q}}_{d_j}-{\underset{\_}{Q}}_{f_j}-{\underset{\_}{Q}}_{c_j}\right)}^T\delta {\underset{\_}{e}}_j=0.$

If we introduce the mapping from the local generalized coordinates ${\underset{\_}{e}}_j$ of element $j$ to the flexible beam global generalized coordinates $\underset{\_}{\bar{q}}$, such that ${\underset{\_}{e}}_j={\mathit{L}}_j\underset{\_}{\bar{q}}$, then the above expression can be written as

(15) $({\mathit{M}}_f\overline{\underset{\_}{\ddot{q}}}-{\underset{\_}{f}}_{v_f}-{\underset{\_}{Q}}_{el}-{\underset{\_}{Q}}_d-{\underset{\_}{Q}}_f-{\underset{\_}{Q}}_c{)}^T\delta \underset{\_}{\bar{q}}=0,$(15)

with

(16) ${\mathit{M}}_f=\sum _{j=1}^{n_e}{\mathit{L}}_j^T{\mathit{M}}_{f_j}{\mathit{L}}_j,{\underset{\_}{f}}_{v_f}=\sum _{j=1}^{n_e}{\mathit{L}}_j^T{\underset{\_}{f}}_{v_{f_j}}\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{Q}}_{\omega }=\sum _{j=1}^{n_e}{\mathit{L}}_j^T{\underset{\_}{Q}}_{{\omega }_j}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\omega =\{el,d,f,c\}.$(16)

Due to the explicit integration of the constraint forces ${\underset{\_}{Q}}_c$ in the above expression, the generalized coordinates of the flexible body can be treated as independent and, thus, its equations of motion can be expressed as

(17) ${\mathit{M}}_f\overline{\underset{\_}{\ddot{q}}}={\underset{\_}{f}}_{v_f}+{\underset{\_}{Q}}_{el}+{\underset{\_}{Q}}_d+{\underset{\_}{Q}}_f+{\underset{\_}{Q}}_c.$(17)

2.2. Equations of motion: RFEM approach

2.2.1. Kinematic description

The main idea of the RFEM method is replacing flexible bodies with a system of rigid finite elements (rfes) that are interconnected via spatial spring-damping elements (sdes) [35, pp. 35]. The first step towards describing the kinematics of the proposed method is the division of the flexible body into rfes and sdes. This division can be performed in two stages [35, pp. 37–40]. First, the flexible body is divided into $n$ elements, each of which has a length $\mathrm{\Delta }l=L/n$, where $L$ is the needle’s total length. The sdes are added at the centre point of the initial elements, as shown in Figure 5. In the secondary division, the initial elements are divided again, so that their edge points are interconnected by the sdes. As shown in Figure 5, intermediate rfes have now length equal to $\mathrm{\Delta }L$, while the first and last rfe have a length of $\mathrm{\Delta }L/2$. The configuration of the flexible body can now be described with the help of the resulting $n_e$ rfes. For this, we consider the rfe $j\in [1,n_e]$ shown in Figure 6. Due to the rigidity assumption, the configuration of the rfe $j$ can be completely characterized with the help of the body frame $f_j$, which is rigidly attached to the element’s centre of mass $C_j$ shown in Figure 6. More specifically, the inertial position of an arbitrary point $P_j$ of element $j$ is given as

Figure 5. Needle division to rfes and sdes.

Figure 6. Spatial RFEM element.

(18) ${\underset{\_}{r}}_{O{P}_j/F}^F={\underset{\_}{r}}_{O{C}_j/F}^F+{\mathit{R}}_{f_j}^F{\underset{\_}{r}}_{C_j{P}_j/f_j}^{f_j},$(18)

where ${\underset{\_}{r}}_{O{C}_j/F}^F$ is the origin of frame $f_j$ measured in the inertial frame $F$ and $R_{f_j}^F(\underset{\_}{{\theta }_j})$ the rotation matrix that describes the orientation of frame $f_j$ with respect to $F$. As before, the rotation matrix is parameterized using the set of Euler angles ${\underset{\_}{\theta }}_j$ that describe the sequence of intrinsic rotations $z_0$-$y{ }^{\mathrm{\prime }}$-$x{ }^{\mathrm{\prime }}{ }^{\mathrm{\prime }}$. It should be noted that, due to the rigidity assumption and contrary to the DDM method, the distance ${\underset{\_}{r}}_{C_j{P}_j/f_j}^{f_j}$ is constant. Thus, a complete characterization of the rfe $j$ configuration can be obtained using the generalized coordinates

(19) ${\underset{\_}{e}}_j={\left[\begin{array}{cc}{({\underset{\_}{r}}_{O{C}_j/F}^F)}^T& {\underset{\_}{\theta }}_j^T\end{array}\right]}^T.$(19)

Based on expression (18), the inertial velocity of point $P_j$ is given as

(20) ${\underset{\_}{\dot{r}}}_{O{P}_j/F}^F={\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j){\underset{\_}{\dot{e}}}_j\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j)=\left[\begin{array}{cc}{\mathit{I}}_3& -{\mathit{R}}_{f_j}^F\mathit{S}({\underset{\_}{x}}_j){\mathit{G}}_j\end{array}\right]\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{x}}_j={\underset{\_}{r}}_{C_j{P}_j/f_j}^{f_j}.$(20)

The matrix ${\mathit{G}}_j=\mathit{G}({\underset{\_}{\theta }}_j)$ maps the derivatives of the Euler angles ${\underset{\_}{\dot{\theta }}}_j$ to the rotational velocity of frame $f_j$, such that

${\underset{\_}{\omega }}_{f_j/F}^{f_j}={\mathit{G}}_j({\underset{\_}{\theta }}_j){\underset{\_}{\dot{\theta }}}_j.$

Similarly, the inertial acceleration of point $P_j$ can be defined as

(21) ${\underset{\_}{\ddot{r}}}_{O{P}_j/F}^F={\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j){\underset{\_}{\ddot{e}}}_j+{\underset{\_}{a}}_{v_{f_j}}({\underset{\_}{x}}_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j),$(21)

where

${\underset{\_}{a}}_{v_{f_j}}({\underset{\_}{x}}_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)={\mathit{R}}_{f_j}^F\left({\mathit{S}}^2({\underset{\_}{\omega }}_{f_j/F}^F){\underset{\_}{x}}_j-\mathit{S}({\underset{\_}{x}}_j){\dot{\mathit{G}}}_j{\underset{\_}{\dot{\theta }}}_j\right).$

Furthermore, based on equations (18)(18) ${\underset{\_}{r}}_{O{P}_j/F}^F={\underset{\_}{r}}_{O{C}_j/F}^F+{\mathit{R}}_{f_j}^F{\underset{\_}{r}}_{C_j{P}_j/f_j}^{f_j},$(18) and (20(20) ${\underset{\_}{\dot{r}}}_{O{P}_j/F}^F={\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j){\underset{\_}{\dot{e}}}_j\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j)=\left[\begin{array}{cc}{\mathit{I}}_3& -{\mathit{R}}_{f_j}^F\mathit{S}({\underset{\_}{x}}_j){\mathit{G}}_j\end{array}\right]\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{x}}_j={\underset{\_}{r}}_{C_j{P}_j/f_j}^{f_j}.$(20) ), the inertial virtual displacement of point $P_j$ can be written as

(22) $\delta {\underset{\_}{r}}_{O{P}_j/F}^F=\frac{\mathrm{\partial }{\underset{\_}{r}}_{O{P}_j/F}^F}{\mathrm{\partial }{\underset{\_}{e}}_j}\delta {\underset{\_}{e}}_j={\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j)\delta {\underset{\_}{e}}_j.$(22)

2.2.2. Virtual work

Similarly to Section 2.1.2, the virtual work of the inertial forces can be defined as

$\delta W_{f_j}^{(in)}={\int }_{m_{f_j}}({\underset{\_}{\ddot{r}}}_{O{P}_j/F}^F{)}^T\delta {\underset{\_}{r}}_{O{P}_j/F}^F\mathrm{d}{m}_{f_j}$

which, using equations (21)(21) ${\underset{\_}{\ddot{r}}}_{O{P}_j/F}^F={\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j){\underset{\_}{\ddot{e}}}_j+{\underset{\_}{a}}_{v_{f_j}}({\underset{\_}{x}}_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j),$(21) and (22(22) $\delta {\underset{\_}{r}}_{O{P}_j/F}^F=\frac{\mathrm{\partial }{\underset{\_}{r}}_{O{P}_j/F}^F}{\mathrm{\partial }{\underset{\_}{e}}_j}\delta {\underset{\_}{e}}_j={\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j)\delta {\underset{\_}{e}}_j.$(22) ), leads to

(23) $\delta W_{f_j}^{(in)}=({\mathit{M}}_{f_j}{\underset{\_}{\ddot{e}}}_j-{\underset{\_}{f}}_{v_{f_j}}{)}^T\delta {\underset{\_}{e}}_j,$(23)

with

${\mathit{M}}_{f_j}({\underset{\_}{e}}_j)={\int }_{m_{f_j}}{\mathit{Z}}_j^T({\underset{\_}{x}}_j,{\underset{\_}{e}}_j){\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j)\mathrm{d}{m}_{f_j}$

and

${\underset{\_}{f}}_{v_{f_j}}({\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)=-{\int }_{m_{f_j}}{\mathit{Z}}_j^T({\underset{\_}{x}}_j,{\underset{\_}{e}}_j){\underset{\_}{a}}_{v_{f_j}}({\underset{\_}{x}}_j,{\underset{\_}{e}}_j,{\underset{\_}{\dot{e}}}_j)\mathrm{d}{m}_{f_j}.$

The key component of the RFEM is the definition of the elastic forces that act on each rfe. As shown in Figure 7, the rfe $j$ is subject to the elastic forces ${\underset{\_}{F}}_{e_j},{\underset{\_}{F}}_{e_{j+1}}$ and the elastic moments ${\underset{\_}{M}}_{e_j},{\underset{\_}{M}}_{e_{j+1}}$ due to the deformation of the associated sdes $s_j$ and $s_{j+1}$. Based on Figure 7, these forces and moments can be modelled as [56]

Figure 7. Elastic forces in RFEM.

(24) ${\underset{\_}{V}}_{e_j}^F={\underset{\_}{V}}_{e{l}_j}^l+{\underset{\_}{V}}_{d_j}^l\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{V}}_{e_{j+1}}^F={\underset{\_}{V}}_{e{l}_j}^r+{\underset{\_}{V}}_{d_j}^r\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}V=\{F,M\}.$(24)

The indices $el$ and $d$ refer to spring and damping forces/moments, respectively. These can be modelled as

$\left[\begin{array}{c}{\underset{\_}{F}}_{e{l}_j}^{\omega }\\ {\underset{\_}{M}}_{e{l}_j}^{\omega }\\ {\underset{\_}{F}}_{d_j}^{\omega }\\ {\underset{\_}{M}}_{d_j}^{\omega }\end{array}\right]=\left[\begin{array}{cccc}{\mathit{K}}_s& {\mathit{O}}_3& {\mathit{O}}_3& {\mathit{O}}_3\\ {\mathit{O}}_3& {\mathit{K}}_t& {\mathit{O}}_3& {\mathit{O}}_3\\ {\mathit{O}}_3& {\mathit{O}}_3& {\mathit{C}}_s& {\mathit{O}}_3\\ {\mathit{O}}_3& {\mathit{O}}_3& {\mathit{O}}_3& {\mathit{C}}_t\end{array}\right]\left[\begin{array}{c}{\underset{\_}{r}}_{{\omega }_j}\\ {\underset{\_}{\theta }}_{{\omega }_j}\\ {\underset{\_}{\dot{r}}}_{{\omega }_j}\\ {\underset{\_}{\dot{\theta }}}_{{\omega }_j}\end{array}\right]$

with $\omega =\{l,r\}$ and ${\mathit{O}}_3$ the $3\times 3$ null matrix. Expressions for the constant stiffness and damping matrices ${\mathit{K}}_s,{\mathit{C}}_s,{\mathit{K}}_t$ and ${\mathit{C}}_t$ can be found in [35, pp. 69–82]. In the above expression, vectors ${\underset{\_}{r}}_{r_j}$ and ${\underset{\_}{r}}_{l_j}$ denote the translational deformation of the sdes $s_j$ and $s_{j+1}$ respectively, and can be written as

${\underset{\_}{r}}_{l_j}={\underset{\_}{r}}_{O{A}_j/F}^F-{\underset{\_}{r}}_{O{B}_{j-1}/F}^F\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{r}}_{r_j}={\underset{\_}{r}}_{O{A}_{j+1}/F}^F-{\underset{\_}{r}}_{O{B}_j/F}^F.$

Similarly, the elastic moments in equation (24)(24) ${\underset{\_}{V}}_{e_j}^F={\underset{\_}{V}}_{e{l}_j}^l+{\underset{\_}{V}}_{d_j}^l\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{V}}_{e_{j+1}}^F={\underset{\_}{V}}_{e{l}_j}^r+{\underset{\_}{V}}_{d_j}^r\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}V=\{F,M\}.$(24) are modelled based on the rotational deformations of the associated sdes $s_j$ and $s_{j+1}$, denoted as ${\underset{\_}{\theta }}_{l_j}$ and ${\underset{\_}{\theta }}_{r_j}$. Due to the non-commutative property of Euler angles, these can be extracted from the rotation matrices ${\mathit{R}}_{f_j}^{f_{j-1}}({\underset{\_}{\theta }}_{l_j})$ and ${\mathit{R}}_{f_{j+1}}^{f_j}({\underset{\_}{\theta }}_{r_j})$ accordingly. From the expression above, the element’s generalized elastic forces are nonlinear functions of the system’s generalized coordinates. As a result, if these expressions are treated as forces external to the element $j$, then the virtual work caused by the generalized elastic and damping forces can be written as

(25) $\delta W_{f_j}^{(el)}={\underset{\_}{Q}}_{e{l}_j}^T\delta {\underset{\_}{e}}_j\mathrm{a}\mathrm{n}\mathrm{d}\delta W_{f_j}^{(d)}={\underset{\_}{Q}}_{d_j}^T\delta {\underset{\_}{e}}_j,$(25)

where

${\underset{\_}{Q}}_{{\omega }_j}^T={\mathit{Z}}_{B_j}^T{\underset{\_}{F}}_{{\omega }_j}^r-{\mathit{Z}}_{A_j}^T{\underset{\_}{F}}_{{\omega }_j}^l+{\left({\mathit{R}}_{f_j}^F{\mathit{G}}_j{\mathit{L}}_{\theta }\right)}^T\left({\underset{\_}{M}}_{{\omega }_j}^r-{\underset{\_}{M}}_{{\omega }_j}^l\right)\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\omega =\{el,d\}$

and

${\mathit{Z}}_{D_j}=\left[\begin{array}{cc}{\mathit{I}}_3& -{\mathit{R}}_{f_j}^F\mathit{S}({\underset{\_}{r}}_{C_j{D}_j/f_j}^{f_j}){\mathit{G}}_j\end{array}\right]\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}D=\{A,B\}.$

The locator matrix ${\mathit{L}}_{\theta }$ in the above expression maps the element’s Euler angles to the element’s generalized coordinates such that ${\underset{\_}{\theta }}_j={\mathit{L}}_{\theta }{\underset{\_}{e}}_j$.

Next, the virtual work of the external forces is defined. For this, we consider the point $P_j^i$ of element $j$, located at the position defined by the vector ${\underset{\_}{x}}_j^i$ with respect to frame $f_j$, which is subject to the point force and point moment pair illustrated in Figure 4b. Then, as illustrated in Section 2.1.2, the virtual work resulting from these forces and moments can be written as

$\delta W_{f_j}^{(ex)}={(}^j{\underset{\_}{F}}_i^F{)}^T\delta {\underset{\_}{r}}_{O{P}_j^i/F}^F+{(}^j{\underset{\_}{M}}_i^F{)}^T{\mathit{R}}_{f_j}^F{\mathit{G}}_j\delta {\underset{\_}{\theta }}_j^i.$

Using equation (22)(22) $\delta {\underset{\_}{r}}_{O{P}_j/F}^F=\frac{\mathrm{\partial }{\underset{\_}{r}}_{O{P}_j/F}^F}{\mathrm{\partial }{\underset{\_}{e}}_j}\delta {\underset{\_}{e}}_j={\mathit{Z}}_j({\underset{\_}{x}}_j,{\underset{\_}{e}}_j)\delta {\underset{\_}{e}}_j.$(22) and the matrix ${\mathit{L}}_{\theta }$, the above expression can be written as

(26) $\delta W_{f_j}^{(ex)}=\left[{{(}^j{\underset{\_}{F}}_i^F)}^T{\mathit{Z}}_j^i+{{(}^j{\underset{\_}{M}}_i^F)}^T{\mathit{R}}_{f_j}^F{\mathit{G}}_j{\mathit{L}}_{\theta }\right]\delta {\underset{\_}{e}}_j,$(26)

where ${\mathit{Z}}_j^i=\mathit{Z}({\underset{\_}{x}}_j^i,{\underset{\_}{e}}_j)$. If we consider that the element $j$ is subject to $n_f$ point forces and $n_m$ point moments, it can be shown that the total virtual work of the external generalized forces is given as

(27) $\delta W_{f_j}^{(ex)}={\underset{\_}{Q}}_{f_j}^T\delta {\underset{\_}{e}}_j\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\underset{\_}{Q}}_{f_j}=\sum _{i=1}^{n_f}({\mathit{Z}}_j^i{)}^T{(}^j{\underset{\_}{F}}_i^F)+{\mathit{L}}_{\theta }^T{\mathit{G}}_j{\mathit{R}}_F^{f_j}\sum _{k=1}^{n_m}{(}^j{\underset{\_}{M}}_k^F).$(27)

2.2.3. Dynamic equilibrium

The procedure for deriving the equations of motion for the flexible beam is identical to the one presented in Section 2.1.3. First, we define the mapping from element’s $j$ local generalized coordinates ${\underset{\_}{e}}_j$ to the beam’s global generalized coordinates $\underset{\_}{\bar{q}}$ as ${\underset{\_}{e}}_j=L_j\underset{\_}{\bar{q}}$. Then, following the steps of Section 2.1.3, it can be shown that the equations of motion of the flexible beam using the RFEM method take the form

(28) ${\mathit{M}}_f\underset{\_}{\ddot{\bar{q}}}={\underset{\_}{f}}_{v_f}+{\underset{\_}{Q}}_{el}+{\underset{\_}{Q}}_d+{\underset{\_}{Q}}_f+{\underset{\_}{Q}}_c,$(28)

where the mass matrix and the generalized force vectors can be calculated as in (16).

2.3. Simulation algorithm

As illustrated in the previous sections, the two proposed approaches lead to the same set of equations, given in (17) and (28). This allows the development of a common solution strategy, independent of the method used, as well as the implementation of direct comparisons between them. To provide a more compact form of the above equations we define the vector

$\underset{\_}{g}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\underset{\_}{f}}_{v_f}+{\underset{\_}{Q}}_{el}+{\underset{\_}{Q}}_d+{\underset{\_}{Q}}_f+{\underset{\_}{Q}}_c.$

Then, equations (17)(17) ${\mathit{M}}_f\overline{\underset{\_}{\ddot{q}}}={\underset{\_}{f}}_{v_f}+{\underset{\_}{Q}}_{el}+{\underset{\_}{Q}}_d+{\underset{\_}{Q}}_f+{\underset{\_}{Q}}_c.$(17) and (28(28) ${\mathit{M}}_f\underset{\_}{\ddot{\bar{q}}}={\underset{\_}{f}}_{v_f}+{\underset{\_}{Q}}_{el}+{\underset{\_}{Q}}_d+{\underset{\_}{Q}}_f+{\underset{\_}{Q}}_c,$(28) ) can be rewritten as

(29) $\underset{\_}{\ddot{\bar{q}}}=\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\mathit{M}}_f^{-1}\underset{\_}{g}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}}).$(29)

As stated in Section 1, the aim of the solution algorithm is to provide an accurate and computationally efficient description of the needle’s three-dimensional dynamics, as the rigid base’s centre of mass $C$ follows an arbitrary trajectory imposed by the surgeon’s hand or a robotic gripper (driving constraints). Furthermore, the solution algorithm should describe the profiles of the reaction forces ${\underset{\_}{F}}_C$ and moments ${\underset{\_}{M}}_C$, illustrated in Figure 8, that result from the combined effect of the structure’s rigid body motion and vibration. The following sections provide a detailed description of the process for deriving the required solution.

Figure 8. Rigid base and flexible needle at connection point $A$.

2.3.1. Boundary conditions

Before the proposed algorithm is formulated, some important remarks for equation (29)(29) $\underset{\_}{\ddot{\bar{q}}}=\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\mathit{M}}_f^{-1}\underset{\_}{g}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}}).$(29) must be made. More specifically, it should be noted that the state vector $\underset{\_}{\bar{q}}$ contains the boundary conditions, due to the rigid connection between the flexible needle and the rigid base, as shown in Figure 8. To eliminate the boundary conditions, the state vector $\underset{\_}{\bar{q}}$ can be partitioned in the constrained ${\underset{\_}{q}}_c\in {\mathbb{R}}^{n_c}$ and active ${\underset{\_}{q}}_a\in {\mathbb{R}}^{n_a}$ coordinates such that

(30) $\underset{\_}{\bar{q}}={\left[\begin{array}{cc}{\underset{\_}{q}}_c^T& {\underset{\_}{q}}_a^T\end{array}\right]}^T.$(30)

The constrained vector ${\underset{\_}{q}}_c$ incorporates the prescribed motion of the associated coordinates (boundary conditions), while the active vector ${\underset{\_}{q}}_a$ describes the free (independent) coordinates that define the beam’s motion. Depending on the chosen FMD technique and according to the different sets of generalized coordinates defined in (1) and (19), the constrained vector ${\underset{\_}{q}}_c$ can take one of the following forms:

${\underset{\_}{q}}_c^{DDM}=\left[\begin{array}{c}{\underset{\_}{r}}_{A_1}\\ {\underset{\_}{\theta }}_{A_1}\end{array}\right]\mathrm{o}\mathrm{r}{\underset{\_}{q}}_c^{RFEM}=\left[\begin{array}{c}{\underset{\_}{r}}_{O{C}_1/F}^F\\ {\underset{\_}{\theta }}_1\end{array}\right].$

If ${\underset{\_}{r}}_{OC/F}^F(t)$ and $\underset{\_}{\theta }(t)$ are known functions, with continuous second derivatives, that describe the position and orientation of frame $f_c$ (pose of surgeon’s hand or robotic gripper) with respect to the inertial frame of reference $F$, then, from the above expression and according to Figure 8, it follows that the constrained vector ${\underset{\_}{q}}_c$ is also a known vector function. Furthermore, it follows that the linear and rotational velocity and acceleration of the constrained coordinates can be obtained.

2.3.2. Recursive Newton–Euler method

A widely used approach for solving the equations (29)(29) $\underset{\_}{\ddot{\bar{q}}}=\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\mathit{M}}_f^{-1}\underset{\_}{g}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}}).$(29) while simultaneously obtaining the generalized reaction forces ${\underset{\_}{F}}_C$ and moments ${\underset{\_}{M}}_C$ profiles is the explicit incorporation of the system’s kinematic constraints through the use of the LaGrange multipliers [74, pp. 323–349]. While this technique has been widely employed in the FMD literature for solving systems of interconnected bodies, it leads to a higher-dimensional system of coupled differential and algebraic equations, the solution of which requires special numerical treatment that increases the computational complexity.

Taking advantage of the specific properties of the needle insertion application, our work presents an alternative approach that can accelerate the solution of (29), while obtaining the required generalized reaction force profiles. More specifically, as the examined system is only under two sets of kinematic constraints, i.e. [a)] the connection of the flexible needle to the rigid base at point $A$ (Figure 8) and the driving constraints at the base’s centre of mass C, it allows the efficient implementation of the generalized recursive Newton–Euler method [75], which can lead to accelerated solutions for a low number of kinematic constraints.

To implement the proposed algorithm, the equations of motion (29), based on the partitioning of (30), can be rewritten as

(31) $\left[\begin{array}{cc}{\mathit{M}}_{f_{cc}}& {\mathit{M}}_{f_{ca}}\\ {\mathit{M}}_{f_{ac}}& {\mathit{M}}_{f_{aa}}\end{array}\right]\left[\begin{array}{c}{\underset{\_}{\ddot{q}}}_c\\ {\underset{\_}{\ddot{q}}}_a\end{array}\right]=\left[\begin{array}{c}{\underset{\_}{g}}_c(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\\ {\underset{\_}{g}}_a(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\end{array}\right],$(31)

where

$\left[\begin{array}{c}{\underset{\_}{g}}_c\\ {\underset{\_}{g}}_a\end{array}\right]=\left[\begin{array}{c}{\underset{\_}{f}}_{v_{f_c}}\\ {\underset{\_}{f}}_{v_{f_a}}\end{array}\right]+\left[\begin{array}{c}{\underset{\_}{Q}}_{e{l}_c}\\ {\underset{\_}{Q}}_{e{l}_a}\end{array}\right]+\left[\begin{array}{c}{\underset{\_}{Q}}_{d_c}\\ {\underset{\_}{Q}}_{d_a}\end{array}\right]+\left[\begin{array}{c}{\underset{\_}{Q}}_{f_c}\\ {\underset{\_}{Q}}_{f_a}\end{array}\right]+\left[\begin{array}{c}{\underset{\_}{Q}}_{c_c}\\ {\underset{\_}{Q}}_{c_a}\end{array}\right].$

Note, the components of the vector $\underset{\_}{g}=[{\underset{\_}{g}}_c^T,{\underset{\_}{g}}_a^T{]}^T$ are all known, except for the generalized constrained forces ${\underset{\_}{Q}}_c$. The key idea for obtaining the system’s solution is that the components of the constraint forces corresponding to the active system’s coordinates are zero, ${\underset{\_}{Q}}_{c_a}={\underset{\_}{0}}_{n_a}$. On the other hand, the components ${\underset{\_}{Q}}_{c_c}$, which result from the reaction forces due to the boundary conditions at point $A$, are finite and incorporate information about the connection forces ${\underset{\_}{F}}_A$ and moments ${\underset{\_}{M}}_A$ (Figure 8). Based on this observation, equation (31)(31) $\left[\begin{array}{cc}{\mathit{M}}_{f_{cc}}& {\mathit{M}}_{f_{ca}}\\ {\mathit{M}}_{f_{ac}}& {\mathit{M}}_{f_{aa}}\end{array}\right]\left[\begin{array}{c}{\underset{\_}{\ddot{q}}}_c\\ {\underset{\_}{\ddot{q}}}_a\end{array}\right]=\left[\begin{array}{c}{\underset{\_}{g}}_c(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\\ {\underset{\_}{g}}_a(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\end{array}\right],$(31) can be written as

${\mathit{M}}_{f_{aa}}{\underset{\_}{\ddot{q}}}_a={\underset{\_}{f}}_{v_{f_a}}+{\underset{\_}{Q}}_{e{l}_a}+{\underset{\_}{Q}}_{d_a}+{\underset{\_}{Q}}_{f_a}-{\mathit{M}}_{f_{ac}}{\underset{\_}{\ddot{q}}}_c.$

As the right-hand side of the above equation is known, the equation can be integrated numerically to identify the time evolution of the generalized active coordinates ${\underset{\_}{q}}_a$. Given that both the constrained and active coordinates are now known, the constrained components of the reaction forces are calculated for the current time step as

(32) ${\underset{\_}{Q}}_{c_c}={\mathit{M}}_{f_{cc}}{\underset{\_}{\ddot{q}}}_c+{\mathit{M}}_{f_{ca}}{\underset{\_}{\ddot{q}}}_a-{\underset{\_}{f}}_{f_{v_c}}-{\underset{\_}{Q}}_{e{l}_c}-{\underset{\_}{Q}}_{d_c}-{\underset{\_}{Q}}_{f_c}.$(32)

Based on the above equation, the reaction forces ${\underset{\_}{F}}_C^f$ and moments ${\underset{\_}{M}}_C^f$ can be calculated with the help of the equations of motion of the rigid base. These can be written as

$\sum {\underset{\_}{F}}_C^F=m_r{\underset{\_}{\ddot{r}}}_{OC/F}^F\mathrm{a}\mathrm{n}\mathrm{d}\sum {\underset{\_}{M}}_C^{f_c}={\mathit{I}}_c^{f_c}{\underset{\_}{\alpha }}_{f/F}^f-\mathit{S}({\underset{\_}{\omega }}_{f/F}^f){\mathit{I}}_c^{f_c}{\underset{\_}{\omega }}_{f/F}^f.$

where $m_r$ is the rigid body’s mass, ${\mathit{I}}_c^{f_c}$ its inertial tensor defined in the frame $f_c$ and vector ${\underset{\_}{\alpha }}_{f/F}^f$ the rotational acceleration of frame $f$ or $f_c$, with respect to the inertial frame of reference. From the above equations and Figure 8, the generalized reaction forces that act on the needle’s base are

${\underset{\_}{F}}_C^f=({\mathit{R}}_f^F{)}^T(m_r{\underset{\_}{\ddot{r}}}_{OC/F}^F-{\underset{\_}{W}}_C^F-{\underset{\_}{F}}_A^F),$

and

${\underset{\_}{M}}_C^f={\mathit{I}}_c^{f_c}{\underset{\_}{\alpha }}_{f/F}^f-{\underset{\_}{M}}_A^f-\mathit{S}({\underset{\_}{r}}_{AC/f}^f){\underset{\_}{F}}_A^f-\mathit{S}({\underset{\_}{\omega }}_{f/F}^f){\mathit{I}}_c^{f_c}{\underset{\_}{\omega }}_{f/F}^f.$

The connection forces ${\underset{\_}{F}}_A$ and moments ${\underset{\_}{M}}_A$ in the above expression are derived based on generalized constraint forces of (32).

2.4. Needle–tissue interaction forces

As discussed in Section 1.4, needle–tissue interaction is modelled with either constrained methods, by imposing kinematic constraints between the needle and its surrounding medium, or with unconstrained methods, through the development of needle–tissue interaction force models. Our work is independent of the needle–tissue modelling approach and can be combined with either the constrained or the unconstrained method. However, as discussed in Section 1.4, the constrained approach is highly dependent on the tissue modelling and discretization. Thus, to maintain the general nature of the proposed needle model, we will only study the unconstrained approach. Furthermore, we will only focus on the RFEM formulation, but a similar approach can be used for the DDM method.

First, we consider the last RFEM element $n_e$, illustrated in Figure 9, which incorporates the needle’s bevel tip. Based on the same figure, we define the needle’s tip point $S_t$ as

Figure 9. Last RFEM element and needle tip definition.

${\underset{\_}{r}}_{C_{n_e}{S}_t/f_{n_e}}^{f_{n_e}}=\left[\begin{array}{c}{l}_{n_e}/2\\ 0\\ r_n\end{array}\right],$

with $l_{n_e}$ the element’s length and $r_n$ the needle’s radius. Similarly, the position of the needle’s tip surface centroid $N_t$ can be defined with the help of similar triangles as

${\underset{\_}{r}}_{C_{n_e}{N}_t/f_{ne}}^{f_{n}e}={\underset{\_}{r}}_{B_{n_e}{N}_t/f_{ne}}^{f_{n}e}-\left[\begin{array}{c}{l}_{n_t}\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}{l}_{n_e}/2\\ 0\\ 0\end{array}\right]-\left[\begin{array}{c}{l}_{n_c}/2\\ 0\\ 0\end{array}\right]\mathrm{o}\mathrm{r}{\underset{\_}{r}}_{C_{n_e}{N}_t/f_{ne}}^{f_{n}e}=\frac{1}{2}\left[\begin{array}{c}{l}_{n_e}-l_{n_c}\\ 0\\ 0\end{array}\right].$

As discussed in Section 1.4, the axial forces acting on the needle can be decomposed as a summation of a cutting/stiffness force acting on the needle’s tip and a friction force distributed along the needle’s shaft. The stiffness forces are due to the elastic properties of the tissue and occur right before the tissue/capsule puncture. On the other hand, friction and cutting forces occur after the puncture event.

In our work, both cutting and stiffness forces are modelled as point loads acting on the needle’s tip. We assume that the stiffness force is applied at the needle’s tip point $S_t$ and it remains parallel to the element’s local axis $x_e$, as shown in Figure 10a. The cutting force is assumed to be normal to the needle’s bevel surface and is applied at its centroid $N_t$ [68], as illustrated in Figure 10b.

Figure 10. Stiffness and cutting forces on the last RFEM element.

From the above figures, we can define the stiffness ${\underset{\_}{F}}_s$ and cutting ${\underset{\_}{F}}_c$ forces expressed in the local element frame $f_e$ as

(33) ${{\underset{\_}{F}}_s}^{f_{n_e}}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=\left[\begin{array}{c}-F_s(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\\ 0\\ 0\end{array}\right]\mathrm{a}\mathrm{n}\mathrm{d}{{\underset{\_}{F}}_c}^{f_{n_e}}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=\left[\begin{array}{c}-F_c(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})sin(a)\\ 0\\ F_c(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})cos(a)\end{array}\right].$(33)

Friction forces are distributed along the needle’s shaft as it slides inside the tissue post-puncture. During the needle insertion, the needle’s shaft is also under the influence of distributed lateral forces due to the tissue compression. Both friction and lateral forces, illustrated in Figure. 11, can be considered as traction vectors acting on the needle’s surface and can be analysed with the help the needle’s tangent plane.

Figure 11. Illustration of post-puncture needle–tissue interaction forces.

More specifically, we first consider the rfe $j$, shown in Figure 12. Furthermore, we consider an arbitrary point $P_j$ that belongs to the surface of the rfe $j$. The position of the point with respect to the element’s local frame $f_j$ can be expressed with the help of cylindrical coordinates as

Figure 12. Infinitesimal area and tangent plane vectors for rfe $j$.

(34) ${\underset{\_}{\hat{r}}}_j(x_j,{\alpha }_j)={\underset{\_}{r}}_{C_j{P}_j/f_j}^{f_j}=\left[\begin{array}{c}{x}_j\\ rcos({\alpha }_j)\\ rsin({\alpha }_j)\end{array}\right],$(34)

where $x_j$ is the distance of the point $P_j$ along the $x_{f_j}$ axis of element $j$, illustrated in Figure 12, $r$ is the radius of the needle (or that of the element rfe $j$) and $a_j$ is the complementary polar angle illustrated in the same Figure.

Equation (34)(34) ${\underset{\_}{\hat{r}}}_j(x_j,{\alpha }_j)={\underset{\_}{r}}_{C_j{P}_j/f_j}^{f_j}=\left[\begin{array}{c}{x}_j\\ rcos({\alpha }_j)\\ rsin({\alpha }_j)\end{array}\right],$(34) is the parametric equation of the surface $S_j$ of the element $j$. The surface is a function of two parameters, namely, the longitudinal distance $x_j$ and the angle ${\alpha }_j$. Using this parametric equation we can define the tangents and normals of the surface. More specifically, the unit vectors ${\underset{\_}{n}}_{p_j}$ and ${\underset{\_}{t}}_{p_j}$ at the point $P_j$ shown in Figure 12 are defined as:

${\underset{\_}{n}}_{p_j}^{f_j}=-\frac{\left(\frac{\mathrm{\partial }{\underset{\_}{\hat{r}}}_j}{\mathrm{\partial }{x}_j}\times \frac{\mathrm{\partial }{\underset{\_}{\hat{r}}}_j}{\mathrm{\partial }{\alpha }_j}\right)}{∥\frac{\mathrm{\partial }{\underset{\_}{\hat{r}}}_j}{\mathrm{\partial }{x}_j}\times \frac{\mathrm{\partial }{\underset{\_}{\hat{r}}}_j}{\mathrm{\partial }{\alpha }_j}∥}=\left[\begin{array}{c}0\\ cos({\alpha }_j)\\ sin({\alpha }_j)\end{array}\right]\mathrm{a}\mathrm{n}\mathrm{d}{\underset{\_}{t}}_{p_j}^{f_j}=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right].$

The unit vector ${\underset{\_}{b}}_{p_j}^{f_j}$ at the point $P_j$ is defined as

${\underset{\_}{b}}_{p_j}^{f_j}={\underset{\_}{t}}_{p_j}^{f_j}\times {\underset{\_}{n}}_{p_j}^{f_j}=\left[\begin{array}{c}0\\ -sin({\alpha }_j)\\ cos({\alpha }_j)\end{array}\right].$

Note that vectors ${\underset{\_}{t}}_{p_j}^{f_j}$, ${\underset{\_}{n}}_{p_j}^{f_j}$ and ${\underset{\_}{b}}_{p_j}^{f_j}$ are defined with respect to the element $j$ frame $f_j$. Finally, observing Figure 12, we can define the infinitesimal surface $d{S}_j$ as

$d{S}_j=rd{\alpha }_{j}d{x}_j.$

To model the traction acting on the needle’s surface, we examine the rfe $j$ shown in Figure 13. First, we assume that a mean traction vector ${\underset{\_}{\tilde{F}}}_{q_j}$ acts on the point $P_j$ and it is transmitted across an infinitesimal area $d{S}_j$. As shown in Figure 13, this traction vector can be divided into three components, each of which points in the directions defined by the tangent and normal vectors ${\underset{\_}{t}}_{p_j}^{f_j}$, ${\underset{\_}{n}}_{p_j}^{f_j}$ and ${\underset{\_}{b}}_{p_j}^{f_j}$.

Figure 13. Traction on an infinitesimal area of rfe $j$.

These components represent the distributed contact forces between the needle and the tissue as a result of their relative motion. More specifically, the traction ${\tilde{F}}_{f_j}$ represents the axial friction due to the longitudinal motion of the element, the traction ${\tilde{F}}_{b_j}$ represents the rotational friction due to the needle’s axial rotation and the traction ${\tilde{F}}_{n_j}$ represents the lateral resistance due to tissue displacement/compression. Based on Figure 13 and the definition of the tangent and normal surface vectors of element $j$, we can define the mean traction vector ${\underset{\_}{F}}_{q_j}$ transmitted across the infinitesimal area $d{S}_j$ as

(35) $\begin{array}{c}{\underset{\_}{\tilde{F}}}_{q_j}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j)=-{\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{t}}_{p_j}}^{f_j}(x_j,{\alpha }_j)-{\tilde{F}}_{n_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{n}}_{p_j}}^{f_j}(x_j,{\alpha }_j)\\ -{\tilde{F}}_{b_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{b}}_{p_j}}^{f_j}(x_j,{\alpha }_j)\end{array}$(35)

The definition of the mean traction vector transmitted across the infinitesimal area $d{S}_j$ allows the calculation of the total surface force acting on the element $j$. This is given by integrating over the entire surface $S_j$ of the element $j$

(36) ${{\underset{\_}{F}}_{q_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\int }_{S_j}{\underset{\_}{\tilde{F}}}_{q_j}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j)d{S}_j={\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\underset{\_}{\tilde{F}}}_{q_j}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j)rd{\alpha }_{j}d{x}_j$(36)

where $x_{l_j}$ and $x_{u_j}$ are the lower and upper integration limits for the longitudinal distance of the element $j$. These values depend on the percentage of the element that is in contact with the tissue. Combining equations (35)(35) $\begin{array}{c}{\underset{\_}{\tilde{F}}}_{q_j}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j)=-{\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{t}}_{p_j}}^{f_j}(x_j,{\alpha }_j)-{\tilde{F}}_{n_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{n}}_{p_j}}^{f_j}(x_j,{\alpha }_j)\\ -{\tilde{F}}_{b_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{b}}_{p_j}}^{f_j}(x_j,{\alpha }_j)\end{array}$(35) and (36(36) ${{\underset{\_}{F}}_{q_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\int }_{S_j}{\underset{\_}{\tilde{F}}}_{q_j}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j)d{S}_j={\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\underset{\_}{\tilde{F}}}_{q_j}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j)rd{\alpha }_{j}d{x}_j$(36) ) the total surface force on element $j$ is given as

${{\underset{\_}{F}}_{q_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={{\underset{\_}{F}}_{f_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})+{{\underset{\_}{F}}_{n_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})+{{\underset{\_}{F}}_{b_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})$

with

(37) ${{\underset{\_}{F}}_{f_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{t}}_{p_j}^i}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j$(37)

(38) ${{\underset{\_}{F}}_{n_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{n_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{n}}_{p_j}}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j$(38)

and

(39) ${{\underset{\_}{F}}_{b_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{b_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{b}}_{p_j}}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j.$(39)

Equations (33)(33) ${{\underset{\_}{F}}_s}^{f_{n_e}}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=\left[\begin{array}{c}-F_s(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\\ 0\\ 0\end{array}\right]\mathrm{a}\mathrm{n}\mathrm{d}{{\underset{\_}{F}}_c}^{f_{n_e}}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=\left[\begin{array}{c}-F_c(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})sin(a)\\ 0\\ F_c(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})cos(a)\end{array}\right].$(33) , (37(37) ${{\underset{\_}{F}}_{f_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{t}}_{p_j}^i}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j$(37) ), (38(38) ${{\underset{\_}{F}}_{n_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{n_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{n}}_{p_j}}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j$(38) ) and (39(39) ${{\underset{\_}{F}}_{b_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{b_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{b}}_{p_j}}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j.$(39) ) can provide a full characterization of the needle–tissue interaction. These can be included in the RFEM needle model, using equation (27)(27) $\delta W_{f_j}^{(ex)}={\underset{\_}{Q}}_{f_j}^T\delta {\underset{\_}{e}}_j\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{\underset{\_}{Q}}_{f_j}=\sum _{i=1}^{n_f}({\mathit{Z}}_j^i{)}^T{(}^j{\underset{\_}{F}}_i^F)+{\mathit{L}}_{\theta }^T{\mathit{G}}_j{\mathit{R}}_F^{f_j}\sum _{k=1}^{n_m}{(}^j{\underset{\_}{M}}_k^F).$(27) . To provide a complete characterization, however, it is essential to define expressions for the functions of stiffness, cutting, axial friction, rotational friction and lateral compression forces.

As discussed in Section 1.4, there is no unique modelling approach for characterizing these forces. For example, in the study by Okamura et al. [1], cutting forces were modelled as constant, while Elgezua et al. [76] proposed a quadratic cutting force model. Similarly, there are various formulations for the distributed loads. For example, friction forces in the study by Okamura et al. [1] were assumed to be uniformly defined over the needle’s length (${\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j)={\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})$) with a magnitude characterized by the modified Karnopp model. Other models, such as the LuGre friction model, have also been proposed [66]. Applying a uniformly distributed Karnopp model for the friction forces of the element $j$ in (37) leads to

${{\underset{\_}{F}}_{f_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]d{\alpha }_{j}d{x}_j=2\pi r(x_{u_j}-x_{l_j})\left[\begin{array}{c}{\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\\ 0\\ 0\end{array}\right],$

with ${\tilde{F}}_{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})$ defined by Okamura et al. [1]. A similar expression can be derived for the rotational friction and lateral tissue resistance of equations (39)(39) ${{\underset{\_}{F}}_{b_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{b_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{b}}_{p_j}}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j.$(39) and (38(38) ${{\underset{\_}{F}}_{n_j}}^{f_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})=-r{\int }_{x_{l_j}}^{x_{u_j}}{\int }_0^{2\pi }{\tilde{F}}_{n_j}(\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}},x_j,{\alpha }_j){{\underset{\_}{n}}_{p_j}}^{f_j}(x_j,{\alpha }_j)d{\alpha }_{j}d{x}_j$(38) ). For complex distributed load expressions, the surface traction integrals can be computed numerically (e.g. through Gaussian quadrature) and applied as equivalent point loads on the RFEM elements.

This Section provided some guidelines regarding the integration of the proposed needle model with unconstrained needle–tissue interaction forces. However, characterizing the needle–tissue interaction forces is beyond the scope of this work. This preserves the general nature of our work and allows the reader to develop and integrate their own needle–tissue interaction models.

2.5. Real-time considerations

As discussed in Section 1, one of the key aims of our work is to provide both highly accurate and computationally efficient mathematical models that allow real-time execution. In this regard, except for the solution strategy presented above, which has been tailored to the requirements of the application to maximize computational efficiency, our work considers three more modelling and implementation aspects that contribute to the acceleration of the overall system’s performance.

The first strategy takes advantage of the structure’s geometry to decrease the computational complexity of the solution. More specifically, our work only considers beam finite elements for modelling the dynamics of flexible medical needles. The reason for this is twofold. First, beam elements have been specifically designed to tackle the particularities of beam-like structures, such as a biopsy needle and, thus, are expected to have superior performance when compared to generic finite elements, such as tetrahedrals or hexahedrals [8]. Second, while 3D beam elements and tetrahedral elements might have the same number of nodal coordinates, tetrahedral-based meshes require the use of a significantly higher number of elements for covering the needle’s geometry. Furthermore, the requirement for mesh symmetry leads to a further increase in the number of required tetrahedral elements, while the employed beam elements are inherently symmetric [8]. Thus, by considering only beam elements, our work aims to minimize the required computational complexity and accelerate the model’s solution.

The second strategy aims to tackle the numerical particularities of the problem. More specifically, due to the thinness of the examined structure, the needle’s axial and torsional stiffness is significantly higher than that in the transversal directions. This leads to a high ratio in the eigenvalues of the Jacobian $J_f=\mathrm{\partial }\underset{\_}{f}/\mathrm{\partial }\underset{\_}{\bar{q}}$ of the vector function $\underset{\_}{f}$ in equation (29)(29) $\underset{\_}{\ddot{\bar{q}}}=\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\mathit{M}}_f^{-1}\underset{\_}{g}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}}).$(29) and, subsequently, to a significant deviation between the response time scales of the $\underset{\_}{\bar{q}}$ vector components. This phenomenon, known as differential equation stiffness, creates significant numerical problems and needs special treatment.

In our work, different approaches have been considered for tackling the stiffness of the differential equations (29(29) $\underset{\_}{\ddot{\bar{q}}}=\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\mathit{M}}_f^{-1}\underset{\_}{g}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}}).$(29) ), while allowing real-time execution. This includes the analysis of explicit, implicit, semi-implicit and multistep numerical integration approaches. Explicit integration techniques are not suitable for stiff differential equations as they lead to solution instability, unless an extremely small integration step size is chosen (which significantly increases the computational time) [77]. On the other hand, implicit integration techniques can guarantee stability without the need for extremely small step sizes, at the cost of solving a set of nonlinear algebraic equations at each time step. To cover the numerical requirements of all of the proposed methods, our work employs diagonally implicit Runge-Kutta methods with an inner Newton-Raphson algorithm parallelized in multiple processing cores for improved performance. This technique leads to stable solutions and real-time execution. As will be discussed in Section 4, the numerical properties of the RFEM method also allowed the use of Rosenbrock-Type methods (or linearly implicit Runge-Kutta methods) [77] leading to accuracy, stability and even further computational efficiency.

The third strategy stems from the concurrency that characterizes the proposed algorithms. More specifically, both the DDM and the RFEM methods include concurrent components that can be executed in parallel allowing for a significant improvement in the overall computational efficiency. Examples of these components are the generalized forces and mass properties that can be calculated independently among the different elements. By exposing the concurrency of the proposed algorithms and employing multithreading execution techniques, the integration time of Eq. (29)(29) $\underset{\_}{\ddot{\bar{q}}}=\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\underset{\_}{f}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}})={\mathit{M}}_f^{-1}\underset{\_}{g}(t,\underset{\_}{\bar{q}},\underset{\_}{\dot{\bar{q}}}).$(29) was significantly accelerated.

3. Results

To demonstrate the performance of the proposed algorithms, three numerical examples are presented in this section. These were designed to give insight into the accuracy and computational efficiency of the proposed methods for both the static and dynamic analysis of the system, while aspects of linear and nonlinear deflection were also considered.

For the validation of the proposed methods, the simulation scenarios were also solved with the help of commercial simulation software. More specifically, the Ansys software was used for the static analysis of the needle’s deflection, while dynamic simulations were performed with the help of MSC Adams. The simulations were implemented on a 4-core Intel® Core™ i7-7700HQ CPU running at 2.80 GHz, using the C++ Armadillo library [78] with the Intel Math Kernel Library (MKL) integration. Aspects of algorithmic parallelization were implemented with the help of the OpenMP API [79].

The needle’s base was treated as a solid cylinder with length $l_h$ = 0.153 m, radius $r_h$ = 0.017 m and a mass of $m_h$ = 0.09 kg. The flexible needle was modelled as a homogeneous stainless steel beam with a circular cross-section of radius $r=0.635\mathrm{m}\mathrm{m}$ and a length of $L=0.2623\mathrm{m}$. Both the mechanical and geometrical properties of the needle were assumed to remain constant along its span. The properties of the system were based on commercially available bevelled tip LATP needles. Next, the three numerical examples are presented.

3.1. Steady-state analysis: tip force loading

First, we examine the steady-state deflection of the flexible needle, using the two proposed FMD approaches. The needle’s deflection results from a constant point force loading applied at its tip (point B, Figure 2). To analyse the performance of the proposed techniques, we consider three spatial loading conditions:

(40) ${ }^{\omega }{\underset{\_}{F}}_{tip}^F={\left[\begin{array}{ccc}0& d_{\omega }& d_{\omega }\end{array}\right]}^T\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\omega =\{a,b,c\}\mathrm{a}\mathrm{n}\mathrm{d}\left\{\begin{array}{c}{d}_a=0.2\mathrm{N}\\ d_b=0.6\mathrm{N}\\ d_c=1.0\mathrm{N}\end{array}\right..$(40)

The gradual increment in the loading allows the analysis of the effect of the large needle deflection on the accuracy of the proposed methods. To examine the accuracy of the methods, the problem is also solved using Ansys (both linear and nonlinear solver). Figure 14 shows the transversal $y$ and $z$ deflection for the three tip force loading conditions $a$, $b$ and $c$ of (40). Figure 15 shows the error measure

Figure 14. Needle response under the three tip force loading conditions of (40).

Figure 15. Error $\hat{e}$ under the three tip force loading conditions of (40).

(41) $\hat{e}(x)=\sqrt{{(y-y_v)}^2+{(z-z_v)}^2},$(41)

where $y_v$ and $z_v$ correspond to the deflections calculated from the Ansys nonlinear solver.

3.2. Transient analysis: step force loading and sinusoidal base movement

Next, we examine the dynamic behaviour of the flexible needle. For this simulation scenario, the needle’s rigid base undergoes a sinusoidal motion, while a step force loading is applied at its tip. More specifically, the input driving constraint has the form

(42) ${\underset{\_}{r}}_{OC/F}^F={\left[\begin{array}{ccc}0& 0& a(t)\end{array}\right]}^T\mathrm{a}\mathrm{n}\mathrm{d}\underset{\_}{\theta }={\left[\begin{array}{ccc}0& a(t)& 0\end{array}\right]}^T,$(42)

where $a(t)=\hat{a}sin(2\pi f_{d}t)$. The excitation magnitude $\hat{a}$ and frequency $f_d$ were arbitrarily chosen as $\hat{a}=0.2\mathrm{m}$ and $f_d=0.5\mathrm{H}\mathrm{z}$. The external force vector, applied at the needle’s tip B, is chosen as

(43) ${\underset{\_}{F}}_B^F={\left[\begin{array}{ccc}0& 0& F_{B}tanh(t)\mathcal{H}(t-2.5)\end{array}\right]}^T.$(43)

The magnitude $F_B$ has a value of $0.6\mathrm{N}$, while $\mathcal{H}(t)$ is the Heaviside function. For the simulation results shown in Figure. 16(a-d), a spatial discretization of four elements is chosen for both methods. The system’s solution using the MSC Adams software acts as a reference point for testing the accuracy of the proposed approaches.

Figure 16. (a) Response of tip’s position in $x$ direction. (b) Response of tip’s position in $z$ direction. (c) Reaction force at needle’s handle in $z$ direction. (d) Reaction moment at needle’s handle in $y$ direction.

To analyse the convergence, accuracy and computational efficiency of the two methods, the following error metrics and ratios are introduced. First, the error metrics that define the time-averaged deviation of the tip position ${\hat{e}}_t$ and the reaction force ${\hat{e}}_f$ from the MSC Adams results are defined as

(44) ${\hat{e}}_t=\frac{1}{n_t}\sum _{i=1}^{n_t}{\parallel }^{\omega }{\underset{\_}{r}}_{OB/F}^F(t_i){-}^a{\underset{\_}{r}}_{OB/F}^F(t_i)\parallel$(44)

and

(45) ${\hat{e}}_f=\frac{1}{n_t}\sum _{i=1}^{n_t}{\parallel }^{\omega }{\underset{\_}{F}}_C^F(t_i){-}^a{\underset{\_}{F}}_C^F(t_i)\parallel ,$(45)

where $n_t$ is the number of points used for the temporal discretization, $\omega =\{RFEM,DDM\}$ and $a=\{Adams\}$. Next, the computational efficiency of the algorithm is determined using the ratio of the algorithm’s execution time to the simulated time. Note that a real-time simulation requires that this ratio is less than one (real-time execution line Figure. 17). The proposed error metrics and the execution-time-to-real-time ratio, with respect to the number of elements used for the beam’s spatial discretization, are illustrated in Figure. 17.

Figure 17. (a)/(b) Mean error ${\hat{e}}_t$ for the tip position and execution time to real-time ratio with respect to the number of elements for RFEM/DDM. (c)/(d) Mean error ${\hat{e}}_f$ for the handle force and execution time to real-time ratio with respect to the number of elements for RFEM/DDM.

3.3. Needle–tissue interaction: Winkler foundation model

The final simulation scenario is illustrated in Figure 18. The transversal needle–tissue interactions are modelled with the help of spring-damper elements along the needle’s length. This model, which originates from the Winkler foundation modelling approach [80], has been extensively used in the needle insertion literature to provide a simplified model of the laterally distributed forces that are present during needle insertion procedures [81–83].

Figure 18. Needle tissue interaction model based on Winkler foundation.

Let $T_i$ be an arbitrary point at the flexible needle, which for the undeformed configuration is written as $T_i^0$. As illustrated in Figure 18, the displacement of the point due to the needle’s deflection can be defined as

(46) ${\underset{\_}{r}}_{\mathrm{\Delta }{T}_i}={\underset{\_}{r}}_{O{T}_i/F}^F-{\underset{\_}{r}}_{O{T}_i^0/F}^F.$(46)

If we assume that the sping-damper forces are only applied in the inertial $z$ direction, as shown in Figure 18, the total force exerted at point $T_i$ can be modelled as

(47) ${\underset{\_}{F}}_{T_i}^F={\left[\begin{array}{ccc}0& 0& -k_t{r}_{\mathrm{\Delta }{T}_i}^z-c_t{\dot{r}}_{\mathrm{\Delta }{T}_i}^z\end{array}\right]}^T,$(47)

where $r_{\mathrm{\Delta }{T}_i}^z$ is the projection of the displacement vector (46) in the inertial $z$ direction, ${\dot{r}}_{\mathrm{\Delta }{T}_i}^z$ is the associated velocity component, while $k_t$ and $c_t$ are the constants of the spring and damper elements, respectively.

The proposed Winkler foundation model was solved using both the DDM and the RFEM approaches, as well as MSC Adams. The stiffness constant was chosen as $k_t=300.0\mathrm{N}/\mathrm{m}$, based on the values reported by Glozman and Shoham [81], while the damping ratio was chosen as $c_t=2.0\mathrm{N}\mathrm{s}/\mathrm{m}$, based on the estimated stiffness-to-damping ratio for soft tissues reported by Khadem et al. [14]. It should be noted that these values are used for providing a qualitative approximation of the types of forces that are present during needle insertion procedures and do not aim to provide an accurate tissue model. The distributed loading was approximated with the help of evenly distributed points that spanned half of the needle’s length. For this simulation, the needle’s rigid base was constrained to follow an arbitrary trajectory

(48) ${\underset{\_}{r}}_{OC/F}^F={\left[\begin{array}{ccc}0& 0& a(t)\end{array}\right]}^T\mathrm{a}\mathrm{n}\mathrm{d}\underset{\_}{\theta }={\left[\begin{array}{ccc}0& 0& 0\end{array}\right]}^T,$(48)

where $a(t)=\hat{a}sin(2\pi f_{d}t)$, with $\hat{a}=0.1\mathrm{m}$ and $f_d=0.5\mathrm{H}\mathrm{z}$. The simulation times are indicative of the time required for a single sample acquisition in LATP procedures. The simulation results are illustrated in Fig. 19.

Figure 19. (a) Response of tip’s position in $x$ direction. (b) Response of tip’s position in z direction. (c) Reaction force at needle’s handle in z direction. (d) Reaction moment at needle’s handle in y direction. (e) Error $e_t$ for the tip position. (f) Error $e_f$ for the handle force.

4. Discussion

The simulation results presented in Section 3 allow for a direct accuracy and efficiency comparison between the different modelling approaches. As shown in Figure 14, linear modelling techniques pose significant limitations for capturing the deflection profile of the flexible needle in the case of large deformations. More specifically, the linear Ansys solver remains valid only for small deformations (Fig. 14a and 15a), while in the case of large needle deflections (Fig. 14b, 15b, 14c and 15c), the linear approach leads to unrealistic and non-volume preserving deflection profiles. On the other hand, as illustrated in the same figures, due to their inherent geometric nonlinearity, both the DDM and the RFEM methods lead to highly accurate deflection profiles even for large needle deformations. This is also illustrated in Figure 15, in which the DMM and the RFEM methods lead to a near-zero value of the performance measure of equation (41)(41) $\hat{e}(x)=\sqrt{{(y-y_v)}^2+{(z-z_v)}^2},$(41) , while the linear approach leads to significant errors especially for large needle deflections. Due to their poor accuracy, linear solutions were not included in the subsequent investigations.

The second simulation scenario was designed to provide insight into the dynamic behaviour of the proposed modelling approaches. As illustrated in Fig. 16 and 17, both the DDM and RFEM approaches capture the vibrational behaviour of the needle with high accuracy and quickly converge to the desired solution, leading to a close to zero deviation from the MSC Adams simulation results. It should be noted that both the equations of motion and the associated algorithms were structured in an identical way for the two methods, which allowed the implementation of direct efficiency comparisons between them. As shown in Fig. 17, while both methods can be used for acquiring real-time solutions, the DDM method leads to significantly higher execution times compared to the RFEM (for the same accuracy level). More specifically, comparison of the mean error ${\hat{e}}_t$ in Fig. 17(a) and 17(b) reveals that for the examined simulation the RFEM approach allows real-time execution using a discretization of up to six elements. On the contrary, a real-time simulation using the DDM approach can only be achieved with no more than three elements, resulting in a comparatively lower accuracy of the method. The same observations can be made for the mean error ${\hat{e}}_f$ illustrated in Fig. 17(c) and 17(d). The performance difference between the two methods can be partly attributed to the fact that DDM requires twice the number of DOFs per element than the RFEM. Additionally, it was observed that, during the implicit integration of the equations of motion, the DDM method generally required more iterations of the inner Newton–Raphson for the solution’s convergence, especially in the cases where a rigid body motion was imposed on the system. Further analysis of the system revealed that the source of the numerical problems, which led to more iterations, was inherent to the method as a large condition number was observed on the generalized strain Jacobian defined in equation (11)(11) $\delta \underset{0.3em\$\$}{{\underset{\_}{ϵ}}_j}=J_{{\eta }_j}\left({\underset{\_}{e}}_j\right)\delta {\underset{\_}{e}}_j\text{where}{J}_{{\eta }_j}{\underset{\_}{e}}_j=\frac{\partial {\underset{0.3em\$\$}{\underset{\_}{\eta }}}_j}{\partial {\underset{0.3em\$\$}{\underset{\_}{e}}}_j}.$(11) . On the other hand, the RFEM method showed superior numerical behaviour allowing fast integration and stable solutions even with the use of semi-implicit integration schemes that can further improve the overall system’s performance. Furthermore, the sparse structure of the system’s matrices and the property of the RFEM algorithm for destructuring and concurrent execution allowed further optimizations on the system’s numerical efficiency.

The third simulation scenario was designed to analyse the effect of complex loading conditions and driving constraints that resemble those encountered during actual needle insertion procedures. As illustrated in Fig. 19, both methods captured the needle dynamics with high accuracy, leading to negligible errors for both the tip position and the reaction force at the surgeon’s hand. It should be noted that the performance plots for this simulation scenario followed a similar structure to the ones presented in Fig. 17. More specifically, both methods allowed real-time execution, with the RFEM outperforming the DDM method in terms of computational efficiency. In both dynamic simulations, the execution-time-to-real-time ratio reported by MSC Adams was close to one. However, as the software does not give the user the ability to extract the exact details of the numerical solution, these results have not been reported in Fig. 17.

A systematic investigation of the different FMD modelling approaches reveals that the RFEM method constitutes one of the most optimal candidates for providing robust and accurate solutions for the examined problem. The properties of the RFEM algorithm, combined with the proposed solution strategy (optimized to the particular properties of the problem) allow high computational efficiency and real-time execution capabilities without the need for modelling assumptions that compromise the overall system’s accuracy. Furthermore, even though the parametric nature of our models allows its application to a variety of needle insertion procedures, our focus remains on long, slender and moderately flexible needles, such as the ones used in LATP and prostate brachytherapy procedures.

It should be noted that the computational times reported in the previous examples do not take into account the deformation of the surrounding tissue, but they solely refer to the computational effort for capturing the needle’s deflection under a generalized force field that quantitatively emulates the one encountered during needle insertion procedures. Ensuring numerical stability and accuracy when modelling soft tissues undergoing large localized strain induced by the needle is the remaining part of the challenge, as tissue modelling is a highly challenging modelling and computational endeavour. The identification of the tissue properties, the proper modelling of fracture and contact dynamics and the treatment of spatial discontinuities during needle insertion constitute open modelling challenges.

However, this work studies the needle’s dynamic model and its vibrational behaviour independent of the tissue and their interaction. It does not attempt to solve the whole problem of needle insertion, but it attempts to provide high fidelity models for the needle’s dynamics. Such needle models, are invaluable for identifying the properties of needle insertion as the needle can be seen as a filter through which the generalized needle–tissue interaction forces can be measured.

5. Conclusions

In this paper, we presented two novel three-dimensional dynamic rigid/flexible models of brachytherapy and LATP needles under a general 3D force field. The proposed models were chosen based on a thorough investigation of different FMD methods and their ability to tackle the numerical particularities of thin flexible medical needles. By employing the theory of FMD and incorporating methods that account for geometric nonlinearities, the proposed approaches minimized the need for introducing modelling assumptions that could compromise their performance. They provided a complete characterization of the system’s rigid motion and vibrational behaviour for both small and large deformations and allowed the correlation of the system’s state with input parameters, such as insertion velocity and needle orientation.

Furthermore, this work presented a novel simulation algorithm, based on the generalized recursive Newton–Euler formulation that characterized and quantified the correlation between the system’s inertial and external forces with the generalized reaction forces acting on its base. The proposed algorithm, tailored to the requirements of the specific application, was combined with parallel architecture implementation strategies to allow for computationally efficient solutions and real-time execution.

The models were assessed in terms of accuracy and computational efficiency and their performance was validated with the help of commercial simulation software. It was shown that linear approaches lead to significant limitations and are not suitable for providing accurate and realistic solutions. On the other hand, both of the employed nonlinear approaches presented high levels of accuracy for both static and dynamic system responses. The RFEM method is shown to be more robust and numerically efficient than the DMM method, and is recommended as the proposed technique for further investigations on medical needle dynamics.

Future work will further analyse the RFEM model in order to validate its performance in experimental studies. The model will also be used in conjunction with computationally efficient and accurate tissue models and will act as the basis for modelling the dynamics of virtual surgical instruments of a fully featured medical simulation solution. Future extensions of our work will also allow the implementation of the proposed model with the help of the graphics processing unit, aiming at further improvement of the model’s computational efficiency.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Additional information

Funding

This work was supported by EPSRC under Grant EP/S021795/1.