DeepDyve requires Javascript to function. Please enable Javascript on your browser to continue.
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
2023-12-31 00:00:00
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 [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 CONTACT Athanasios Martsopoulos athanasios.martsopoulos@bristol.ac.uk 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 [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 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. [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 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 [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 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 [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 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 [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 4 A. MARTSOPOULOS ET AL. 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 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. [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. 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 [36] 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 [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 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 [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 sim- plifies 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 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) [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 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 [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. 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 [59], 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 [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 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 [22]. 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. [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 8 A. MARTSOPOULOS ET AL. 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 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 [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 mathema- tical 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 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 [55], 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 [55]. 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 [50]. 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 [72]. The generalized strains, which remain invariant under arbitrary rigid body motions and orthogonal coordinate transformations [72], 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 [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] σ ¼ S �; (10) j j where S is the element’s constant symmetric stiffness matrix, defined by Schwab and Meijaard [50]. 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. [73]. 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 [56] 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 [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 " # � �� � € 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 [68], 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) Stiffness 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 ; αÞ F ð� q; � q; x ; αÞ n ðx ; αÞ j j f j j j j n j j j j q j p j p j j j (35) ~ _ F ð� q; � q; x ; αÞ b ðx ; αÞ b j j j j j p The definition of the mean traction vector transmitted across the infinitesimal area dS allows the calculation of the total surface force acting on the element j. This is given by integrating over the entire surface S of the element j ð ð ð x 2π f f j j _ ~ _ ~ _ F ð� q; � qÞ ¼ F ð� q; � q; x ; αÞ dS ¼ F ð� q; � q; x ; αÞr dα dx (36) q j j j j j j j q q j j j S x 0 j l where x and x are the lower and upper integration limits for the longitudinal distance l u j j of the element j. These values depend on the percentage of the element that is in contact with the tissue. Combining equations (35) and (36) the total surface force on element j is given as f f f f j j j j _ _ _ _ F ð� q; � qÞ ¼ F ð� q; � qÞþ F ð� q; � qÞþ F ð� q; � qÞ q f n b j j j j with ð ð x 2π f i _ _ F ð� q; � qÞ ¼ r F ð� q; � q; x ; αÞt ðx ; αÞ dα dx (37) f j j j j j j f j p j j x 0 ð ð x 2π f f j j _ ~ _ F ð� q; � qÞ ¼ r F ð� q; � q; x ; αÞn ðx ; αÞ dα dx (38) n j j j j j j n j p j j x 0 and ð ð x 2π f f j j _ _ F ð� q; � qÞ ¼ r F ð� q; � q; x ; αÞb ðx ; αÞ dα dx : (39) b j j j j j j b j p j j x 0 Equations (33), (37), (38) and (39) can provide a full characterization of the needle–tissue interaction. These can be included in the RFEM needle model, using equation (27). To Figure 13. Traction on an infinitesimal area of rfe j. 26 A. MARTSOPOULOS ET AL. 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 (F ð� q; � q; x ; αÞ ¼ F ð� q; � qÞ) with a magnitude characterized by the f j j f j j 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 2 3 2 3 ð ð ~ _ x 2π F ð� q; � qÞ j j _ ~ _ 4 5 4 5 F ð� q; � qÞ ¼ r F ð� q; � qÞ 0 dα dx ¼ 2πrðx x Þ ; f j j u l f j j j x 0 0 0 � � with F ðq; qÞ defined by Okamura et al. [1]. A similar expression can be derived for the rotational friction and lateral tissue resistance of equations (39) and (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, character- izing 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 execu- tion. 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 compu- tational 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 tetra- hedrals 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 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 27 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=@� q of the vector function f in equation (29) and, subsequently, to a significant deviation between the response time scales of the � 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), 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 algo- rithms. 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) 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] 28 A. MARTSOPOULOS ET AL. with the Intel Math Kernel Library (MKL) integration. Aspects of algorithmic paralleli- zation were implemented with the help of the OpenMP API [79]. The needle’s base was treated as a solid cylinder with length l = 0.153 m, radius r = h h 0.017 m and a mass of m = 0.09 kg. The flexible needle was modelled as a homogeneous stainless steel beam with a circular cross-section of radius r ¼ 0:635mm and a length of L ¼ 0:2623m. 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: < d ¼ 0:2N ω F F ¼ ½ 0 d d � where ω ¼ fa; b; cg and d ¼ 0:6N : (40) ω ω b tip d ¼ 1:0N 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 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ^ eðxÞ ¼ ðy y Þ þðz z Þ ; (41) v v where y and z correspond to the deflections calculated from the Ansys nonlinear solver. v v 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 T T r ¼ ½ 0 0 aðtÞ� and θ ¼ ½ 0 aðtÞ 0� ; (42) OC=F where aðtÞ ¼ ^ a sinð2πf tÞ. The excitation magnitude ^ a and frequency f were arbitrarily d d chosen as ^ a ¼ 0:2m and f ¼ 0:5Hz. The external force vector, applied at the needle’s tip B, is chosen as F ¼ ½ 0 0 F tanhðtÞHðt 2:5Þ� : (43) The magnitude F has a value of 0:6 N, while HðtÞ is the Heaviside function. For the simulation results shown in Figure. 16(a-d), a spatial discretization of four elements is MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 29 (a) Needle deflection in y and z direction under constant tip force loading F . tip (b) Needle deflection in y and z direction under constant tip force loading F . tip (c) Needle deflection in y and z direction under constant tip force loading F . tip Figure 14. Needle response under the three tip force loading conditions of (40). 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. To analyse the convergence, accuracy and computational efficiency of the two meth- ods, the following error metrics and ratios are introduced. First, the error metrics that 30 A. MARTSOPOULOS ET AL. F b (a) Error e ˆ under constant tip force loading F . (b) Error e ˆ under constant tip force loading F . tip tip (c) Error e ˆ under constant tip force loading F . tip Figure 15. Error ^ e under the three tip force loading conditions of (40). ^ ^ define the time-averaged deviation of the tip position e and the reaction force e from the t f MSC Adams results are defined as ω F a F ^ e ¼ k r ðtÞ r ðtÞ k (44) t i i OB=F OB=F i¼1 and ω F a F ^ e ¼ k F ðtÞ F ðtÞ k; (45) i i C C i¼1 where n is the number of points used for the temporal discretization, ω ¼ fRFEM; DDMg and a ¼ fAdamsg. 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 discretiza- tion, are illustrated in Figure. 17. 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 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 31 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. model of the laterally distributed forces that are present during needle insertion proce- dures [81–83]. Let T be an arbitrary point at the flexible needle, which for the undeformed config- uration is written as T . As illustrated in Figure 18, the displacement of the point due to the needle’s deflection can be defined as F F r ¼ r r : (46) �T OT =F OT =F i i 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 can be modelled as � � F z z 0 0 k r c r _ F ¼ ; (47) t t T �T �T i i i z z where r is the projection of the displacement vector (46) in the inertial z direction, r _ �T �T i i is the associated velocity component, while k and c are the constants of the spring and t t 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 ¼ 300:0N=m, based on the values reported by Glozman and Shoham [81], while the damping ratio was chosen as c ¼ 2:0Ns=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 32 A. MARTSOPOULOS ET AL. Figure 17. (a)/(b) Mean error e 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 e for the handle force and execution time to real-time ratio with respect to the number of elements for RFEM/DDM. Figure 18. Needle tissue interaction model based on Winkler foundation. 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 T T r ¼ ½ 0 0 aðtÞ� and θ ¼ ½ 0 0 0� ; (48) OC=F MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 33 where aðtÞ ¼ ^ a sinð2πf tÞ, with ^ a ¼ 0:1m and f ¼ 0:5Hz. The simulation times are d d indicative of the time required for a single sample acquisition in LATP procedures. The simulation results are illustrated in Fig. 19. 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 illu- strated 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), 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 beha- viour 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 ^ e 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 ^ e 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). 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 34 A. MARTSOPOULOS ET AL. 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 for the tip position. (f) Error e for the handle force. t f 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 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 35 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 proper- ties, 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 compro- mise 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 general- ized recursive Newton–Euler formulation that characterized and quantified the correla- tion 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 36 A. MARTSOPOULOS ET AL. 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 relation- ships that could have appeared to influence the work reported in this paper. Funding This work was supported by EPSRC under Grant EP/S021795/1. References [1] A.M. Okamura, C. Simone, and M.D. O’Leary, Force modeling for needle insertion into soft tissue, IEEE Trans. Biomed. Eng. 51 10 (Oct 2004), pp. 1707–1716. doi:10.1109/TBME.2004. [2] G. Ravali and M. Muniyandi, Haptic feedback in needle insertion modeling and simulation: Review. IEEE Reviews in Biomedical Engineering. 10 (2017), pp. 63–77. [3] N. Abolhassani, R. Patel, and M. Moallem, Needle insertion into soft tissue: A survey, Med Eng Phys 29 (4) (2007), pp. 413–431. doi:10.1016/j.medengphy.2006.07.003. [4] H.S. Bassan, R.V. Patel, and M. Moallem, A novel manipulator for percutaneous needle insertion: Design and experimentation, IEEE/ASME Trans J. Mechatron. 14 6 (Dec 2009), pp. 746–761. doi:10.1109/TMECH.2009.2011357 [5] J. Yanof, J. Haaga, P. Klahr et al., CT-integrated robot for interventional procedures: Preliminary experiment and computer-human interfaces, Comput. Aided Surg. 6 (6) (2001), pp. 352–359. doi:10.3109/10929080109146304. [6] J. Hing, A. Brooks, and J. Desai, Reality-based needle insertion simulation for haptic feedback in prostate brachytherapy. In: Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006, Orlando, FL, USA. (2006), pp. 619–624. [7] O. Goksel, K. Sapchuk, and S.E. Salcudean, Haptic simulator for prostate brachytherapy with simulated needle and probe interaction, IEEE Trans Haptics. 4(3), (Jul 2011), pp. 188–198. doi:10.1109/TOH.2011.34 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 37 [8] O. Goksel, E. Dehghan, and S.E. Salcudean, Modeling and simulation of flexible needles, Med Eng Phys 31 (9) (2009), pp. 1069–1078. doi:10.1016/j.medengphy.2009.07.007. [9] S.P. DiMaio and S.E. Salcudean, Needle insertion modeling and simulation, IEEE Trans. Robot. Autom. 19 5 (Oct 2003), pp. 864–875. doi:10.1109/TRA.2003.817044 [10] M. Freutel, H. Schmidt, L. Dürselen et al., Finite element modeling of soft tissues: Material models, tissue interaction and challenges, Clin Biomech 29 (4) (2014), pp. 363–372. doi:10. 1016/j.clinbiomech.2014.01.006. [11] H. Delingette, Toward realistic soft-tissue modeling in medical simulation, Proceedings of the IEEE. 86 (3), (1998), pp. 512–523. [12] S. Jiang, P. Li, Y. Yu et al., Experimental study of needle–tissue interaction forces: Effect of needle geometries, insertion methods and tissue characteristics, J. Biomech. 47 (13) (2014), pp. 3344–3353. doi:10.1016/j.jbiomech.2014.08.007. [13] L. Barbé, B. Bayle, M. De Mathelin et al., Needle insertions modeling: Identifiability and limitations, Biomed. Signal Process. Control 2 (3) (2007), pp. 191–198. doi:10.1016/j.bspc.2007. 06.003. [14] M. Khadem, C. Rossa, N. Usmani et al., A two-body rigid/flexible model of needle steering dynamics in soft tissue, IEEE/ASME Trans J. Mechatron. 21 (5) (Oct 2016), pp. 2352–2364. doi:10.1109/TMECH.2016.2549505 [15] T. Lehmann, C. Rossa, N. Usmani et al., A real-time estimator for needle deflection during insertion into soft tissue based on adaptive modeling of needle–tissue interactions, IEEE/ASME Trans J. Mechatron. 21 (6) (Dec 2016), pp. 2601–2612. doi:10.1109/TMECH.2016.2598701 [16] R.J. Webster I, J.S. Kim, N.J. Cowan et al., Nonholonomic modeling of needle steering, Int J Rob Res 25 (5–6) (2006), pp. 509–525. doi:10.1177/0278364906065388. [17] G.J. Vrooijink, M. Abayazid, S. Patil et al. Needle path planning and steering in a three-dimensional non-static environment using two-dimensional ultrasound images, Int J Rob Res.3310 (20149), pp. 1361–1374 10.1177/0278364914526627 [18] Y. Zhang, Z. Ju, H. Zhang et al., 3-d path planning using improved rrt* algorithm for robot- assisted flexible needle insertion in multilayer tissues, IEEE Can. J. Electr. Comput. Eng. 12 (2021), pp. 1–13. [19] A. Asadian, M.R. Kermani, and R.V. Patel, An analytical model for deflection of flexible needles during needle insertion. In: 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems, San Francisco, CA, USA. (Sep 2011), pp. 2551–2556. [20] N. Abolhassani, R.V. Patel, and F. Ayazi, Minimization of needle deflection in robot-assisted percutaneous therapy, Int. J. Med. Robot. Comput. Assist Surg 3 (2) (2007), pp. 140–148. doi:10.1002/rcs.136. [21] K.E. Bisshopp and D.C. Drucker, Quarterly of Applied Mathematics, Large deflection of cantilever beams. 3 (1945), pp. 272–275. [22] S.P. DiMaio and S.E. Salcudean, Interactive simulation of needle insertion models, IEEE Trans. Biomed. Eng. 52 7 (Jul 2005), pp. 1167–1179. doi:10.1109/TBME.2005.847548 [23] Y. Adagolodjo, L. Goffin, M. De Mathelin et al. Robotic insertion of flexible needle in deformable structures using inverse finite-element simulation, IEEE Trans Rob. 353 (20196), pp. 697–708 10.1109/TRO.2019.2897858 [24] A. Haddadi and K. Hashtrudi-Zaad, Development of a dynamic model for bevel-tip flexible needle insertion into soft tissues. Conference proceedings: Annual International Conference of the IEEE Engineering in Medicine and Biology Society IEEE Engineering in Medicine and Biology Society Conference, Boston, Massachusetts, USA. 2011 (2011), pp. 7478–7482. [25] T. Podder, D. Clark, D. Fuller et al., Effects of velocity modulation during surgical needle insertion. In: 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference, Shanghai, China. (2005), pp. 5766–5770. [26] T. Podder, D. Clark, J. Sherman et al., In vivo motion and force measurement of surgical needle intervention during prostate brachytherapy, Med. Phys. 33 (8) (2006), pp. 2915–2922. doi:10.1118/1.2218061. [27] L. Tan, X. Qin, Q. Zhang et al., Effect of vibration frequency on biopsy needle insertion force, Med Eng Phys 43 (2017), pp. 71–76. doi:10.1016/j.medengphy.2017.02.011. 38 A. MARTSOPOULOS ET AL. [28] C. McGill, J. Schwartz, J. Moore et al., Effects of insertion speed and trocar stiffness on the accuracy of needle position for brachytherapy, Med. Phys. 39 (4) (2012), pp. 1811–1817. doi:10.1118/1.3689812. [29] M. Abayazid, M. Kemp, and S. Misra, 3D flexible needle steering in soft-tissue phantoms using Fiber Bragg Grating sensors. In: 2013 IEEE International Conference on Robotics and Automation, Karlsruhe, Germany. (May 2013), pp. 5843–5849. [30] N. Abolhassani, R. Patel, and F. Ayazi, Effects of different insertion methods on reducing needle deflection. In: 2007 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Lyon, France. (Aug 2007), pp. 491–494. [31] R. Tsumura, Y. Takishita, and H. Iwata, Needle insertion control method for minimizing both deflection and tissue damage. Journal of Medical Robotics Research. 4 (1), (2019), pp. [32] N.A. Wood, K. Shahrour, M.C. Ost et al., Needle steering system using duty-cycled rotation for percutaneous kidney access. In: 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology, Buenos Aires, Argentina. 8 (2010), pp. 5432–5435. [33] S. Badaan, D. Petrisor, C. Kim, Does needle rotation improve lesion targeting?, The International Journal of Medical Robotics + Computer Assisted Surgery: MRCAS. 7 (2), (2011), pp.138–147. doi:10.1002/rcs.381. [34] A. Martsopoulos, P. Rajendra, B. Stefanos et al., Spatial rigid/flexible dynamic model of biopsy and brachytherapy needles under a general force field. In: 2020 IEEE International Conference on Intelligent Robots and Systems (IROS), Las Vegas, NV, USA. (Oct 2020). [35] E. Wittbrodt, I. Adamiec-Wójcik, and S. Wojciech, Dynamics of Flexible Multibody Systems: Rigid Finite Element Method. Springer, Foundations of Engineering Mechanics, Berlin Heidelberg, (2007). [36] A.A. Shabana, Flexible multibody dynamics: Review of past and recent developments, Multibody Syst Dyn. 1 2 (Jun 1997), pp. 189–222. doi:10.1023/A:1009773505418 [37] J.M. Mayo, D. García-Vallejo, and J. Domínguez, Study of the geometric stiffening effect: Comparison of different formulations, Multibody Syst Dyn. 11 4 (May 2004), pp. 321–341. doi:10.1023/B:MUBO.0000040799.63053.d9 [38] A. Shabana, Dynamics of Multibody Systems, Cambridge, England: Cambridge University Press, 2003. [39] E.M. Bakr and A.A. Shabana, Geometrically nonlinear analysis of multibody systems, Comput Struct. 23 6 (Jan 1986), pp. 739–751. doi:10.1016/0045-7949(86)90242-7 [40] J. Mayo, J. Domínguez, and A. Shabana, Geometrically nonlinear formulations of beams in flexible multibody dynamics, J Vib Acoust. 117 4 (Oct 1995), pp. 501–509. doi:10.1115/1.2874490 [41] I. Sharf, Geometrically Non-Linear Beam Element for Dynamics Simulation of Multibody Systems, Int J Numer Methods Eng 39 (5) (1996), pp. 763–786. doi:10.1002/(SICI)1097- 0207(19960315)39:5<763::AID-NME879>3.0.CO;2-X. [42] U. Lugrís, M. Naya, J. Pérez et al., Implementation and efficiency of two geometric stiffening approaches, Multibody Syst Dyn. 20 (2) (Sep 2008), pp. 147–161. doi:10.1007/s11044-008-9114-6 [43] A. Nada, B. Hussein, S. Megahed, Use of the floating frame of reference formulation in large deformation analysis: Experimental and numerical validation, Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics. 224 (1), (Mar 2010). pp. 45–58. [44] J.N. Reddy, An Introduction to Nonlinear Finite Element Analysis Second Edition: With Applications to Heat Transfer, Fluid Mechanics, and Solid Mechanics, Oxford, England: OUP Oxford, 2014. [45] A.A. Shabana, Computational Continuum Mechanics, Cambridge, England: Cambridge University Press, 2008. [46] A. Shabana, An Absolute Nodal Coordinate Formulation for the Large Rotation and Deformation Analysis of Flexible Bodies, Chicago, Illinois: Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, 1996. [47] A.A. Shabana and R.Y. Yakoub, Three dimensional absolute nodal coordinate formulation for beam elements: Theory, Mech. Eng. J. 123 4 (May 2000), pp. 606–613. doi:10.1115/1.1410100 MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 39 [48] J. Sopanen and A. Mikkola, Studies on the stiffness properties of the absolute nodal coordinate formulation for three-dimensional beams, European Journal of Ophthalmology 13 (Suppl. 3) (2003), pp. 0. doi:10.5301/EJO.2008.4238. [49] B.A. Hussein, H. Sugiyama, and A.A. Shabana, Coupled deformation modes in the large deformation finite-element analysis: Problem definition, J. Comput. Nonlinear Dyn. 2 2 (Nov 2006), pp. 146–154. 1. doi:10.1115/1.2447353 [50] A. Schwab and J. Meijaard, Comparison of three-dimensional flexible beam elements for dynamic analysis: Classical finite element formulation and absolute nodal coordinate formulation, J. Comput. Nonlinear Dyn. 51 (2010 Jan; 5), 10.1115/1.4000320. [51] J. Sopanen and A. Mikkola, Description of elastic forces in absolute nodal coordinate formulation, Nonlinear Dyn. 34 1/2 (Oct 2003), pp. 53–74. doi:10.1023/B:NODY. 0000014552.68786.bc [52] J. Gerstmayr and A.A. Shabana, Analysis of thin beams and cables using the absolute nodal co-ordinate formulation, Nonlinear Dyn. 45 1 (Jul 2006), pp. 109–130. doi:10.1007/s11071- 006-1856-1 [53] P.G. Gruber, K. Nachbagauer, Y. Vetyukov et al., A novel director-based Bernoulli–Euler beam finite element in absolute nodal coordinate formulation free of geometric singularities, Int. J. Mech. Sci. 42 (2013 Aug; 4), pp. 279–289 10.5194/ms-4-279-2013 [54] S. von Dombrowski, Analysis of large flexible body deformation in multibody systems using absolute coordinates, Multibody Syst Dyn. 8 4 (Nov 2002), pp. 409–432. doi:10.1023/ A:1021158911536 [55] J.B. Jonker and J.P. Meijaard, A geometrically non-linear formulation of a three-dimensional beam element for solving large deflection multibody system problems, Int J Non Linear Mech. 53 (2013), pp. 63–74. doi:10.1016/j.ijnonlinmec.2013.01.012 [56] H. Kaminski and P. Fritzkowski, Application of the rigid finite element method to modelling ropes, Lat. Am. J. Solids Struct. 10 (1) (2013), pp. 91–99. doi:10.1590/S1679- [57] I. Adamiec-Wójcik, J. Awrejcewicz, L. Brzozowska, Modelling of ropes with consideration of large deformations and friction by means of the rigid finite element method. Springer Proceedings in Mathematics & Statistics. 93 (2014), pp. 115–137. [58] M.A. Gorin, A.R. Meyer, M. Zimmerman et al., Transperineal prostate biopsy with cognitive magnetic resonance imaging/biplanar ultrasound fusion: Description of technique and early results, World J Urol. 388 (20208), pp. 1943–1949. doi:10.1007/s00345-019-02992-4. [59] C. Rossa and M. Tavakoli, Issues in closed-loop needle steering, Control Eng. Pract. 62 (2017), pp. 55–69. doi:10.1016/j.conengprac.2017.03.004. [60] N. Chentanez, R. Alterovitz, D. Ritchie et al., ACM Transactions on Computer Systems, Interactive simulation of surgical needle insertion and steering. 28 (3), (2009). pp. 1–10. [61] D. Gao, Y. Lei, B. Lian et al., Modeling and simulation of flexible needle insertion into soft tissue using modified local constraints, J Manuf Sci Eng. 138 (12), (2016), pp. 1167–1175. [62] J. Xu, L. Wang, K.C.L. Wong et al. A meshless framework for bevel-tip flexible needle insertion through soft tissue. In: 2010 3rd IEEE RAS EMBS International Conference on Biomedical Robotics and Biomechatronics, Tokyo, Japan. (Sep 2010). p. 753–758. [63] L. Wang and S. Hirai, A local constraint method for needle insertion modeling and simulation. In: 2011 IEEE International Workshop on Haptic Audio Visual Environments and Games, Qinhuangdao, China. (Oct 2011). p. 39–44. null. [64] A. Asadian, M.R. Kermani, and R.V. Patel, A compact dynamic force model for needle-tissue interaction. In: 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology, Buenos Aires, Argentina. (Aug 2010). p. 2292–2295. [65] K. Yan, T. Podder, D. Xiao et al., An improved needle steering model with online parameter estimator, Int. J. Comput. Assist. Radiol. Surg. 1 (4) (2006), pp. 205–212. doi:10.1007/ s11548-006-0058-0. [66] A. Asadian, R.V. Patel, and M.R. Kermani, A distributed model for needle-tissue friction in percutaneous interventions. In: 2011 IEEE International Conference on Robotics and Automation, Shanghai, China. (May 2011). p. 1896–1901. 40 A. MARTSOPOULOS ET AL. [67] H. Yang, P.X. Liu, and J. Zhang, Modelling of needle insertion forces for surgical simulation. In: IEEE International Conference Mechatronics and Automation, 2005, Niagara Falls, ON, Canada. 2 (Jul 2005). p. 592–595. [68] M. Khadem, C. Rossa, R.S. Sloboda et al., Mechanics of tissue cutting during needle insertion in biological tissue, IEEE Robot. Autom. Lett. 1 (2) (Jul 2016), pp. 800–807. doi:10.1109/ LRA.2016.2528301 [69] P. Li, S. Jiang, Y. Yu et al., Biomaterial characteristics and application of silicone rubber and PVA hydrogels mimicked in organ groups for prostate brachytherapy, J Mech Behav Biomed Mater. 49 (2015), pp. 220–234. [70] B. Zhao, L. Lei, L. Xu et al. Needle deflection modeling and preoperative trajectory planning during insertion into multilayered tissues. IEEE/ASME Transactions on Mechatronics, 262 (2021-04), pp. 943–954. [71] W. Durham, Aircraft Flight Dynamics and Control Aerospace Series, Hoboken, New Jersey: John Wiley & Sons, 2013. [72] J.P. Meijaard, Validation of flexible beam elements in dynamics programs, Nonlinear Dyn. 9 1 (Feb 1996), pp. 21–36. doi:10.1007/BF01833291 [73] W.S. Yoo, O. Dmitrochenko, and D. Pogorelov, Review of finite elements using absolute nodal coordinates for large-deformation problems and matching physical experiments. ASME 2005 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Long Beach, California, USA, 6 (2005), pp. 1285–1294. [74] A. Shabana, Computational Dynamics, Hoboken, New Jersey, U.S.: Wiley, 2001. [75] A.A. Shabana, Dynamics of flexible bodies using generalized Newton-Euler equations, J Dyn Syst Meas Control. 112 (09 1990), pp.496–503. 10.1115/1.2896170 [76] I. Elgezua, Y. Kobayashi, and M.G. Fujie, Estimation of needle tissue interaction based on non-linear elastic modulus and friction force patterns. In: 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, Chicago, IL, USA. (Sep 2014). p. 4315–4320. [77] E. Hairer and G. Wanner, Solving Ordinary Differential Equations Ii: Stiff and differential- algebraic Problems, Springer, Berlin Heidelberg, 2010. (Springer Series in Computational Mathematics). [78] C. Sanderson and R. Curtin, Armadillo: A template-based c++ library for linear algebra, J. Open Source Softw. 1 (2) (2016), pp. 26. doi:10.21105/joss.00026. [79] L. Dagum and R. Menon, OpenMP: An industry-standard API for shared-memory programming, IEEE Comput. Sci. Eng. 5 1 (Jan 1998), pp. 46–55. doi:10.1109/99.660313 [80] C.S. Yim and A.K. Chopra, Earthquake response of structures with partial uplift on Winkler foundation, Earthq. Eng. Struct. Dyn. 12 (2) (1984), pp. 263–281. doi:10.1002/eqe. [81] D. Glozman and M. Shoham, Image-guided robotic flexible needle steering, IEEE Trans Rob. 23 3 (Jun 2007), pp. 459–467. doi:10.1109/TRO.2007.898972 [82] K. Yan, T. Podder, Y. Yu et al. Flexible needle-tissue interaction modeling with depth- varying mean parameter: Preliminary study. ichange1 IEEE Trans. Biomed. Eng. 56 (2009), pp. 255–262. [83] H. Lee and J. Kim, Estimation of flexible needle deflection in layered soft tissues with different elastic moduli, Med. Biol. Eng. Comput. 52 (2014), pp. 729–740.
http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png
Mathematical and Computer Modelling of Dynamical Systems
Taylor & Francis
http://www.deepdyve.com/lp/taylor-francis/modelling-and-real-time-dynamic-simulation-of-flexible-needles-for-aFWVasxNZa
