Open Advanced Search
Get 20M+ Full-Text Papers For Less Than $1.50/day.
Start a 14-Day Trial for You or Your Team.
Learn More →
Modelling and real-time dynamic simulation of flexible needles for prostate biopsy and brachytherapy
Modelling and real-time dynamic simulation of flexible needles for prostate biopsy and brachytherapy
Martsopoulos, Athanasios; Hill, Thomas L.; Persad, Rajendra; Bolomytis, Stefanos; Tzemanaki, Antonia
MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 2023, VOL. 29, NO. 1, 1–40 https://doi.org/10.1080/13873954.2022.2158875 RESEARCH ARTICLE Modelling and real-time dynamic simulation of flexible needles for prostate biopsy and brachytherapy a,b a c c Athanasios Martsopoulos , Thomas L. Hill , Rajendra Persad , Stefanos Bolomytis a,b and Antonia Tzemanaki a b School of Civil, Aerospace and Mechanical Engineering, University of Bristol, Bristol, UK; Department of Medical Robotics, Bristol Robotics Laboratory, Bristol, UK; Bristol Urological Institute, Southmead Hospital, Bristol, UK ABSTRACT ARTICLE HISTORY Received 29 October 2021 Percutaneous needle insertion constitutes a widely adopted techni- Accepted 09 December 2022 que for performing minimally invasive operations. Robot-assisted needle placement and virtual surgical training platforms have the KEYWORDS potential to significantly improve the accuracy of these operations. Surgery simulation; needle For this, the development of mathematical models that provide insertion; flexible multibody a complete characterization of the underlying dynamics of medical dynamics; rigid finite needles is considered of paramount importance. In this paper, we element method 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 pre- vious 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 effi - ciency 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 effi - cient solution for capturing the dynamics of flexible medical needles. 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 . 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 . 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 . 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 CONTACT Athanasios Martsopoulos firstname.lastname@example.org School of Civil, Aerospace and Mechanical Engineering, University of Bristol, Bristol, UK © 2023 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 2 A. MARTSOPOULOS ET AL. in anatomical structures among patients make navigation through tissue particularly complicated and decrease the procedure’s overall accuracy . 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 : [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 . 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 experimen- tally 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. , 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 beha- viour 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 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 3 compared with experimental studies . 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 mod- els, 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  and Goksel et al.  developed quasi-static methods based on nonlinear 2D and 3D Finite Element Method (FEM). Similarly, Adagolodjo et al.  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 assum- ing 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 . 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 . Due to these limitations, it is necessary to develop modelling methodologies that incorporate the underlying needle dynamics. Khadem et al.  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. , 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 4 A. MARTSOPOULOS ET AL. and Hashtrudi-Zaad  extended the static angular spring FEM approach of Goksel et al. : 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. , 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 con- straints 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.  and Haddadi and Hashtrudi-Zaa , 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. . 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. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 5 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  a detailed review of available techni- ques 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 . 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) . As discussed by Shabana , 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 non- linearities 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 geo- metric nonlinearity and a computationally efficient solution, it suffers from the numer- ical problem of element locking, as described by Reddy . 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 . 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 . 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 . 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 sim- plifies the associated equations of motion and accelerates computations . 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 6 A. MARTSOPOULOS ET AL. 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) . To address this, different approaches have been proposed in the literature, such as the development of higher-order thin elements , the uncoupling between the shear deformation and antisymmetric bending  and the redefinition of the system’s degrees of freedom (DOFs) . While some of the proposed techniques have the potential to alleviate the problem of Poisson locking, they introduce additional modelling complexity  or remove the initial benefits of the ANCF, such as the formulations of Gruber et al.  and von Dombrowski  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 . A comparative analysis between DDM and a modified (locking-free) version of ANCF in the study by Schwab and Meijaard  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 deforma- tion 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 , 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. 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 deflec- tions, 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 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 7 Figure 1. Large deflections and large rigid body motion in prostate brachytherapy and biopsy. 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 , during needle insertion the deformed needle shaft compresses the surrounding tissue which, in turn, applies forces to the former influen- cing its deflection. As a result, needle deflection and tissue deformation are coupled effects that influence each other . 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 interac- tions 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 . The computational cost of this remeshing makes the application of this technique unsuitable for interactive applications. One possible solu- tion is using a meshless discretization, as in the study by Xu et al. , or using the technique of local constraints presented by Wang and Hirai . 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 8 A. MARTSOPOULOS ET AL. plausible or applicable in real-time applications or in vivo tests . 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 depen- dent 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 . One of the most widely used unconstrained models is the one proposed by Okamura et al. . Experimental studies reported in the study Okamura et al.  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. , significant research has been devoted to the definition of mathema- tical models that can accurately capture these force profiles. For example, whilst Okamura et al.  proposed a modified Karnopp friction model, different friction modelling approaches were then adopted by Asadian et al.  and Yang et al. . 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 formula- tions [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. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 9 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 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 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 r which describes the position of a vector AP with respect AP=f to frame f expressed in frame f [71, pp. 13–16]. The following sections present the 1 2 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 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 jP½1; n � can be acquired based on the inertially measured cross-sectional position and orientation at its nodes A and B , while intermediate points are defined with j j the help of interpolation functions. The element’s vector of generalized coordinates can be defined as 10 A. MARTSOPOULOS ET AL. Figure 3. Spatial DDM element at reference and deformed state. h i T T T T e ¼ r θ r θ ; (1) j A A B B j j j j F F where r ¼ r and r ¼ r are the inertial positions of the nodal cross-sections’ A OA =F B OB =F j j j j midpoints A and B , while θ and θ are the Euler angles, based on the sequence of j j A B j j 0 00 intrinsic rotations z -y -x , that define the cross-sectional orientation at nodes A and B . 0 j j These are defined as h i � θ ψ θ ¼ D D with D ¼ fA; Bg: j j D D j With the help of the Euler angles θ and θ , the set of triad unit vectors, shown in A B j j Figure 3, are defined as F F e ¼ R ðθ Þ e with r ¼ fx; y; zg and D ¼ fA; Bg: r f D D D D j j j F F The rotation matrices R and R describe the orientation of the local frames f and f A B f f j j A B j j f f f f f f A A A B B B j j j j j j with respect to the inertial frame of reference F. The e ; e ; e and e ; e ; e are x y z x y z A A A B B B j j j j j j constant unit vectors pointing in the associated directions x; y and z. The inertial position of point P of element j, illustrated in Figure 3, can be obtained with the help of the generalized element coordinates and interpolation functions �ðxÞ as i j F F F r ¼ � ðxÞ r þ � ðxÞ e þ � ðxÞ r þ � ðxÞ e : (2) 1 j 2 j 3 j 4 j OP =F A x B x j j j A B j j Expressions for the �ðxÞ; i ¼ f1; . . . 4g, can be chosen as in the study by Jonker and i j Meijaard , while the scalar x is the point’s distance along the element’s elastic line measured using the co-rotational frame introduced by Jonker and Meijaard . Based on equation (2), the inertial velocity of point P can be expressed as _ _ r ¼ Z ðx ; eÞ e ; (3) j j j j OP =F j MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 11 where h i f f A B j j F F Z ðx ; eÞ ¼ � ðxÞI � ðxÞR Sðe ÞG � ðxÞI � ðxÞR Sðe ÞG : j j j 1 j 3 2 j x A 3 j 3 4 j x B f A j f B j A j B j j j The matrices G ¼ Gðθ Þ and G ¼ Gðθ Þ in the above expression map the derivatives A B j A j B j j _ _ of the Euler angles θ and θ to the respective frame rotational velocities, such that A B j j f f A B j j _ _ ω ¼ G θ and ω ¼ G θ : A B j A j B f =F j f =F j A B j j Based on equation (3) the inertial virtual displacement of point P is given as @r OP =F δr ¼ δe ¼ Z ðx ; eÞδe ; (4) j j OP =F j j j @e while the acceleration of point P takes the form € r ¼ Z ðx ; eÞ€ e þ a ðx ; e ; e _ Þ: (5) j j j j j v j j OP =F j f In the above expression, the Coriolis/centrifugal component of acceleration is given as A B j j a ðx ; e ; e _ Þ ¼ � ðxÞ a þ � ðxÞ a ; v v j 2 j 4 j v j j f f where � � f f f D D D D j j j j F 2 _ _ a ¼ R S ðω Þe Sðe ÞG θ with D ¼ fA; Bg: v x x D f f D D j D D f =F j j j Similarly, the orientation of the cross-sections along the element’s length can be writ- ten as h i 2 3 > � 2 3 > � ¼ � ; � � A B > j j s �ðxÞ � j > � � j j < h i 6 7 4 5 � θ ¼ j ¼ s ðxÞ θ with ; 4 � 5 θ ¼ θ ; θ j j A B θ j j j j ψ h i s ðxÞ ψ > T j j ψ � ψ ¼ ψ ; ψ A B j j where expressions for the interpolation vector functions s ðxÞ, s�ðxÞ and s ðxÞ can be � j j j � ψ � found in the study by Schwab and Meijaard . If the mapping from the nodal components of the Euler angles to the generalized element’s coordinates is defined as � r ¼ L e , with r ¼ f�; θ; ψg, then the above expressions take the form j � r 2 3 s ðxÞL � � � � 4 5 s�ðxÞL� θ ¼ S ðxÞ e where S ðxÞ ¼ j : (6) θ j θ j θ θ j j s ðxÞL � j ψ 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 12 A. MARTSOPOULOS ET AL. thickness. Based on this assumption, the virtual work of the inertial forces of element j can be defined as ð ð ðinÞ T F F δW ¼ ð€ r Þ δr dm þ ρI � δ� dx ; f p j j j OP =F OP =F j j f j j m 0 where m is the mass, l the length, I the polar moment of area and ρ the density of f ej p j j element j. The second term of the above equations employs the first component of the sequence of Euler angles � to describe the inertial forces due to the element’s torsional dynamics. Using equations (4), (5) and (6), the above expression can be rewritten as ðinÞ δW ¼ ðM € e f Þ δe ; (7) j j j t r t r where M ¼ M þ M . The mass matrices M and M are defined as j f f f f j j j j ð ð t T r T T M ðeÞ ¼ Z ðx ; eÞZ ðx ; eÞ dm and M ¼ L ρI s s dx L �: j j j f � p � � j j j j � � f j j f � j � j j m 0 The generalized Coriolis/centrifugal vector takes the form f ðe ; e _ Þ ¼ Z ðx ; eÞ a ðx ; e ; e _ Þ dm : j j f j j j j v j j j v 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 . The generalized strains, which remain invariant under arbitrary rigid body motions and orthogonal coordinate transformations , take the form 2 3 l l j e � � 6 7 T T 2 3 F F F F 6 7 ðe Þ e ðe Þ e =2 � z y y z � A B A B 1 6 j j j j 7 6 7 6 7 j � � F 2 6 7 l ðeÞ e 6 7 j i z 6 A 7 6 7 � � 6 7 6 7 T � �¼ ¼ ; (8) 6 F 7 6 7 l ðeÞ e � e 4 6 j i z 7 6 7 j 6 7 4 5 � � 6 j 7 j F 6 l ðeÞ e 7 j i y � j 6 6 7 4 5 l ðeÞ e j i y � � where l ¼k r r k and e ¼ ðr r Þ=l . In equation (8), the generalized strain j j B A l B A j j j j � � ; is used to describe the beam’s axial elongation, � � is for torsion, while � � ; . . . ; � � 1 2 3 6 j j j 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 � ¼ � �ðeÞþ ω ~ð� �ðeÞÞ ¼ ηðeÞ; (9) j j j j j j MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 13 where expressions for the function ω ~ð� �ðeÞÞ, which incorporates the second-order terms, and the function ηðeÞ, which maps the generalized coordinates to generalized strain measures, are given by Jonker and Meijaard . 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  σ ¼ S �; (10) j j where S is the element’s constant symmetric stiffness matrix, defined by Schwab and Meijaard . Based on equation (9), the infinitesimal change of the generalized strains can be calculated as � � � � @η δ� ¼ J e δe where J e ¼ : (11) η η j j j j j j @e With the help of equations (9), (10) and (11), the virtual work of the element’s elastic forces can be written as � � � � ðelÞ T T T j δW ¼ σ δ� ¼ Q δe where Q ¼ J e S η e : (12) j j j η j j j el el j j j j Similarly, the virtual work of the element’s damping forces can be written as ðdÞ δW ¼ Q δe , where the generalized damping forces Q can be modelled using various j d d j j approaches, as discussed by Yoo et al. . Next, we define the virtual work for the external forces. For this, it is assumed that i i a random point P of element j, located at the beam’s elastic line and at a distance x from j j the nodal point A , is subjected to a point force F ðt; e ; e _ Þ and a point moment i j j j F M ðt; e ; e _ Þ, shown in Figure 4a. i j j Then, the virtual work resulting from these forces and moments can be written as ðexÞ j T j T i F F F δW ¼ ð F Þ δr þð M Þ δθ : f i OP =F i j Using equations (4) and (6), the above expression takes the form ðexÞ j T j T F i F i δW ¼ ð F Þ Z δe þð M Þ S δe ; (13) f i j j i θ j i i i i where Z ¼ Z ðx ; eÞ and S ¼ S ðxÞ. The virtual work of the external generalized forces j θ j j j θ j presented in equation (13) can be extended to accommodate for n point forces and n f m point moments. In this case, it can be proven that the total virtual work of the external generalized forces is given as f n X X ðexÞ T j T j T i F k F δW ¼ Q δe with Q ¼ ðZ Þ ð F Þþ ðS Þ ð M Þ: (14) j j i θ k f f f j j i¼1 k¼1 Finally, the virtual work of the external generalized constraint forces for element j can be ðcÞ calculated as δW ¼ Q δ� q. j c j 14 A. MARTSOPOULOS ET AL. (a) Free body diagram of element j for the (b) Free body diagram of element j for the DDM method. RFEM method. Figure 4. Free body diagrams for DDM and RFEM methods. 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 ðinÞ ðelÞ ðdÞ ðexÞ ðcÞ δW ¼ δW þ δW þ δW þ δW : f f f f f The virtual works of the above expression can be calculated considering the contributions from the individual elements as n n e e � � X X ðinÞ ðelÞ ðdÞ ðexÞ ðcÞ δW ¼ δW þ δW þ δW þ δW : f f f f f j j j j j j¼1 j¼1 Substituting the derived virtual work expressions, the above equation becomes � � n T M € e f Q Q Q Q δe ¼ 0: j j j el d f c v j j j j j¼1 If we introduce the mapping from the local generalized coordinates e of element j to the flexible beam global generalized coordinates � q, such that e ¼ L � q, then the above expression can be written as ðM € q f Q Q Q Q Þ δ� q ¼ 0; (15) el d f c with n n n e e e X X X T T T M ¼ L M L ; f ¼ L f and Q ¼ L Q with ω ¼ fel; d; f; cg: f f j j j j j ω ω v v f f j¼1 j¼1 j¼1 (16) Due to the explicit integration of the constraint forces Q 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 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 15 M € q ¼ f þ Q þ Q þ Q þ Q : (17) el d f c 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 �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 �L, while the first and last rfe have a length of �L=2. The configuration of the flexible body can now be described with the help of the resulting n rfes. For this, we consider the rfe jP½1; n � e 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 , which is rigidly attached to the element’s centre of mass C shown in Figure 6. More specifically, the inertial position of an arbitrary point P of element j is given as F F F r ¼ r þ R r ; (18) OP =F OC =F f C P =f j j j j j j F F where r is the origin of frame f measured in the inertial frame F and R ðθÞ the j j OC =F f j j rotation matrix that describes the orientation of frame f with respect to F. As before, the rotation matrix is parameterized using the set of Euler angles θ that describe the 0 00 sequence of intrinsic rotations z -y -x . It should be noted that, due to the rigidity assumption and contrary to the DDM method, the distance r is constant. Thus, C P =f j j j a complete characterization of the rfe j configuration can be obtained using the general- ized coordinates h i F T e ¼ ðr Þ θ : (19) j j OC =F Based on expression (18), the inertial velocity of point P is given as h i F F j I R SðxÞG r _ ¼ Z ðx ; eÞ e _ with Z ðx ; eÞ ¼ 3 j and x ¼ r : (20) j j j j j j j j f j OP =F j C P =f j j j j The matrix G ¼ GðθÞ maps the derivatives of the Euler angles θ to the rotational j j velocity of frame f , such that ω ¼ GðθÞ θ : j j f =F Similarly, the inertial acceleration of point P can be defined as € r ¼ Z ðx ; eÞ € e þ a ðx ; e ; e _ Þ; (21) OP =F j j j v j j j j f j 16 A. MARTSOPOULOS ET AL. Figure 5. Needle division to rfes and sdes. Figure 6. Spatial RFEM element. where � � F 2 F a ðx ; e ; e _ Þ ¼ R S ðω Þ x SðxÞG θ : v j j j f f =F j j j f j j Furthermore, based on equations (18) and (20), the inertial virtual displacement of point P can be written as @r OP =F δr ¼ δe ¼ Z ðx ; eÞδe : (22) j j j j OP =F @e 2.2.2. Virtual work Similarly to Section 2.1.2, the virtual work of the inertial forces can be defined as MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 17 ðinÞ T F F δW ¼ ð€ r Þ δr dm f OP =F OP =F j j j j which, using equations (21) and (22), leads to ðinÞ δW ¼ ðM € e f Þ δe ; (23) j j f j with M ðeÞ ¼ Z ðx ; eÞZ ðx ; eÞ dm f j f j j j j j j j and f ðe ; e _ Þ ¼ Z ðx ; eÞa ðx ; e ; e _ Þ dm : j j j j j v j j j 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 F ; F and the elastic e e j jþ1 moments M ; M due to the deformation of the associated sdes s and s . Based on j jþ1 e e j jþ1 Figure 7, these forces and moments can be modelled as  F l l F r r V ¼ V þ V and V ¼ V þ V with V ¼ fF; Mg: (24) e el d e el d j j j jþ1 j j The indices el and d refer to spring and damping forces/moments, respectively. These can be modelled as 2 3 2 3 2 3 el ω K O O O j s 3 3 3 j 6 7 6 7 M 6 7 O K O O ω 6 el 7 6 j 7 j 3 t 3 3 6 7 6 7 ¼ 6 7 4 5 r _ O O C O ω 4 5 4 j 5 d 3 3 s 3 M O O O C θ 3 3 3 t d ω j j with ω ¼ fl; rg and O the 3� 3 null matrix. Expressions for the constant stiffness and damping matrices K ; C ; K and C can be found in [35, pp. 69–82]. In the above s s t t expression, vectors r and r denote the translational deformation of the sdes s and r l j j s respectively, and can be written as jþ1 F F F F r ¼ r r and r ¼ r r : l OA =F OB =F r OA =F OB =F j j j 1 j jþ1 j Similarly, the elastic moments in equation (24) are modelled based on the rotational deformations of the associated sdes s and s , denoted as θ and θ . Due to the non- j jþ1 l r j j commutative property of Euler angles, these can be extracted from the rotation matrices f f j 1 j R ðθ Þ and R ðθ Þ accordingly. From the expression above, the element’s generalized l r f j f j j jþ1 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 ðelÞ ðdÞ T T δW ¼ Q δe and δW ¼ Q δe ; (25) f j f j j el j d j j 18 A. MARTSOPOULOS ET AL. Figure 7. Elastic forces in RFEM. where � � � � T T r T l F r l Q ¼ Z F Z F þ R G L M M with ω ¼ fel; dg j θ B ω A ω f ω ω ω j j j j j j j and h i Z ¼ I R Sðr ÞG with D ¼ fA; Bg: D 3 j j f j C D =f j j j The locator matrix L in the above expression maps the element’s Euler angles to the element’s generalized coordinates such that θ ¼ L e . j j Next, the virtual work of the external forces is defined. For this, we consider the point i i P of element j, located at the position defined by the vector x with respect to frame f , j 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 ðexÞ j F T F j F T F i δW ¼ ð F Þ δr i þð M Þ R G δθ : i i f j f OP =F j j Using equation (22) and the matrix L , the above expression can be written as h i T T ðexÞ j j F i F F δW ¼ ð F Þ Z þð M Þ R G L δe ; (26) j θ i j i f j j j i i where Z ¼ Zðx ; eÞ. If we consider that the element j is subject to n point forces and n f m j j j point moments, it can be shown that the total virtual work of the external generalized forces is given as f n X X ðexÞ T j j j T i F T F δW ¼ Q δe with Q ¼ ðZ Þ ð F Þþ L G R ð M Þ: (27) f j j i θ F k j f f j j i¼1 k¼1 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 e to the beam’s global generalized coordinates � q as e ¼ L � q. j j 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 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 19 M � q ¼ f þ Q þ Q þ Q þ Q ; (28) el d f c 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 gðt; � q; � qÞ ¼ f þ Q þ Q þ Q þ Q : el d f c Then, equations (17) and (28) can be rewritten as € _ _ _ � q ¼ fðt; � q; � qÞ where fðt; � q; � qÞ ¼ M gðt; � q; � 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 F and moments M , illustrated in C C 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. 2.3.1. Boundary conditions Before the proposed algorithm is formulated, some important remarks for equation (29) must be made. More specifically, it should be noted that the state vector � 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 � q n n c a can be partitioned in the constrained q PR and active q PR coordinates such that c a h i T T q q � q ¼ : (30) c a The constrained vector q incorporates the prescribed motion of the associated coordi- nates (boundary conditions), while the active vector q 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 q can take one of the following forms: � � � � DDM A RFEM 1 OC =F q ¼ or q ¼ : c c 1 1 If r ðtÞ and θðtÞ are known functions, with continuous second derivatives, that OC=F describe the position and orientation of frame f (pose of surgeon’s hand or robotic gripper) with respect to the inertial frame of reference F, then, from the above expression 20 A. MARTSOPOULOS ET AL. Figure 8. Rigid base and flexible needle at connection point A. and according to Figure 8, it follows that the constrained vector q 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) while simultaneously obtaining the generalized reaction forces F and moments M profiles is the explicit incorporation C C 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 , 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 " # � �� � € g ðt; � q; � qÞ M M f f cc ca ¼ ; (31) € q _ M M g ðt; � q; � qÞ f f ac aa where " # " # � � � � � � � � g Q Q Q f f c el d c c c c c c ¼ þ þ þ þ : g f Q Q Q Q el d f c a v a a a a fa T T Note, the components of the vector g ¼ ½g ; g � are all known, except for the general- c a ized constrained forces Q . 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, Q ¼ 0 . On the other hand, the components Q , which result from the reaction c a c a c MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 21 forces due to the boundary conditions at point A, are finite and incorporate information about the connection forces F and moments M (Figure 8). Based on this observation, A A equation (31) can be written as € € M q ¼ f þ Q þ Q þ Q M q : f f aa ac el d f a v a a a 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 q . 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 Q ¼ M € q þ M € q f Q Q Q : (32) f f cc ca c el d f c c a f c c c f f Based on the above equation, the reaction forces F and moments M can be calculated C C with the help of the equations of motion of the rigid base. These can be written as X X f f f f F F c f f c c F ¼ m € r and M ¼ I α Sðω ÞI ω : C OC=F C c f=F f=F c f=F where m is the rigid body’s mass, I its inertial tensor defined in the frame f and vector r c α the rotational acceleration of frame f or f , with respect to the inertial frame of f=F reference. From the above equations and Figure 8, the generalized reaction forces that act on the needle’s base are f T F F F F F ¼ ðR Þ ðm € r W F Þ; C f OC=F C A and f f f f f f f f f c c M ¼ I α M Sðr ÞF Sðω ÞI ω : C c A A c f=F AC=f f=F f=F The connection forces F and moments M in the above expression are derived based on A A 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 , illustrated in Figure 9, which incorporates the needle’s bevel tip. Based on the same figure, we define the needle’s tip point S as t 22 A. MARTSOPOULOS ET AL. 2 3 l =2 4 5 r ¼ 0 ; C S =f n t n e e with l the element’s length and r the needle’s radius. Similarly, the position of the n n needle’s tip surface centroid N can be defined with the help of similar triangles as 2 3 2 3 2 3 2 3 l l =2 l =2 l l n n n n n t e c e c f e f e f e n n n 4 5 4 5 4 5 4 5 r ¼ r 0 ¼ 0 0 or r ¼ 0 : C N =f B N =f C N =f ne t ne ne t ne ne t ne 0 0 0 0 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 and it remains parallel to the element’s local axis x , 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 , as illustrated in Figure 10b. From the above figures, we can define the stiffness F and cutting F forces expressed s c in the local element frame f as 2 3 2 3 F ð� q; � qÞ sinðaÞ F ð� q; � qÞ f f n n e 4 5 e 4 5 _ _ F ð� q; � qÞ ¼ and F ð� q; � qÞ ¼ 0 : (33) s c F ð� q; � qÞ cosðaÞ 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. More specifically, we first consider the rfe j, shown in Figure 12. Furthermore, we consider an arbitrary point P that belongs to the surface of the rfe j. The position of the point with respect to the element’s local frame f can be expressed with the help of cylindrical coordinates as Figure 9. Last RFEM element and needle tip definition. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 23 (a) Stiﬀness forces. (b) Cutting forces. Figure 10. Stiffness and cutting forces on the last RFEM element. 2 3 4 5 ^ r ðx ; αÞ ¼ r ¼ r cosðαÞ ; (34) j j C P =f j j j r sinðαÞ where x is the distance of the point P along the x axis of element j, illustrated in j j f Figure 12, r is the radius of the needle (or that of the element rfe j) and a is the complementary polar angle illustrated in the same Figure. Equation (34) is the parametric equation of the surface S of the element j. The surface is a function of two parameters, namely, the longitudinal distance x and the angle α . j j Using this parametric equation we can define the tangents and normals of the surface. More specifically, the unit vectors n and t at the point P shown in Figure 12 are p p j j defined as: � � 2 3 2 3 @^ r @^ r j j 0 1 @x @α f j j f j j 4 5 4 5 � � cosðαÞ n ¼ ¼ and t ¼ 0 : p j p j j @^ r @^ r � � j j � � � sinðαÞ 0 @x @α j j The unit vector b at the point P is defined as p j 2 3 f f f j j j 4 5 b ¼ t � n ¼ sinðαÞ : p p p j j j cosðαÞ f f f j j j Note that vectors t , n and b are defined with respect to the element j frame f . Finally, p p p j j j observing Figure 12, we can define the infinitesimal surface dS as dS ¼ rdα dx : j j 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 F acts on the point P and it is q j j 24 A. MARTSOPOULOS ET AL. Figure 11. Illustration of post-puncture needle–tissue interaction forces. Figure 12. Infinitesimal area and tangent plane vectors for rfe j. transmitted across an infinitesimal area dS . As shown in Figure 13, this traction vector can be divided into three components, each of which points in the directions defined by f f f j j j the tangent and normal vectors t , n and b . p p p j j 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 F represents the axial friction due to the longitudinal motion of the element, the traction F represents the rotational friction due to the needle’s axial rotation and the traction F 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 F transmitted across the infinitesimal area dS as j MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 25 f f j j ~ _ ~ _ ~ _ F ð� q; � q; x ; αÞ ¼ F ð� q; � q; x ; αÞt ðx ; αÞ