Finite Arc Method: Fast-Solution Extended Piecewise Constant Curvature Model of Soft Robots with Large Variable Curvature Deformations

Abstract

Accurate deformation modeling of soft flexural robots is of high practical importance, especially for high-risk tasks such as surgery. In this study, a new mechanistic model, that is, finite arc method (FAM), for soft robots, for example, tendon-drive, was proposed and validated. First, the catheter was modeled as a finite number of arcs, each with a constant bending curvature, hence the name FAM. Afterward, using a validated Bezier shape approximation, the deformation was parameterized, and the governing equations of the robot were derived. Also, a fast and recursive algorithm was proposed and implemented for the mechanical solution of the robot's deformation. To validate the proposed method, two validation studies were performed. In Study I, the FAM's predicted deformations for eight load cases in each two-dimensional and three-dimensional space on a 40 mm long flexure were compared with the nonlinear finite element method (FEM). In Study II, a representative set of lateral forces on a cardiac catheter (obtained in our previous study) was used to find its FAM-based deformation and was compared with the experimental reference. The error between FAM and FEM deformations was 0.23 ± 0.89 mm with computation times of 3 mseconds (FAM) versus 1244 mseconds (FEM). Also, the error of FAM compared with ground-truth in Study II was 1.41 ± 1.47 mm with a computation time of 7 mseconds. The proposed method showed acceptable performance for the accurate prediction of highly complex large deformations in real time.

Introduction

Background

Flexural robots are a class of soft robots that have been used in interventional robotics frequently. For example, flexural robots have been used in catheter ablation therapy in the left atrium.^1,2 This treatment has become the gold standard for the treatment of cardiac atrial fibrillation, the most prevalent cause of cardiac arrhythmia and the leading cause of hospitalization in cardiac electro-physiology.³

Figure 1 depicts a schematic view of the RFA catheter inside the atrium during the ablation procedure. This procedure is minimally invasive; thus, the surgeon does not have a direct vision of the trajectory of the shape of the catheter. Instead, they obtain such information through real-time X-ray imaging, that is, fluoroscopy.

FIG. 1.

The schematic shape of steerable ablation catheters in the left atrium during ablation interventions.

To avoid the health risks associated with X-ray exposure, remote leader–follower robotic platforms with tendon-driven catheters have been proposed recently.⁴ Such systems have enabled surgeons to not only move away from the X-ray exposure risk zone but also do intervention from miles away.^5,6 Similar to other novel medical technologies, simulation-based training is a clinical need for shortening surgeons' learning curves and increasing the adoption rate of flexural medical robots.

To this end, having a fast and accurate computational model that can reliably predict the deformed shape of the flexural robot, for example, catheter, in real time is a requirement. In this study, we have proposed and validated a fast and accurate model for flexural robots that has superior performance in computational time compared with the most accurate computational method in the literature, for example, the finite element method (FEM).

Related studies

In the field of soft robotics, the mechanical modeling of these systems has evolved through various phases, each marked by distinct methodologies and theoretical advances. Initially, the focus was on discretization strategies to manage the infinite degrees of freedom inherent in soft robotic systems.^7,8 However, while simplifying the equations of motion, this approach often led to a loss of critical information and inaccurate deformation predictions.

As a countermeasure, researchers began incorporating continuous shape interpolation functions, such as Bezier splines, for better deformation modeling and integration into kinematic balance and constitutive equations.^9–11

The modeling of large deformations in flexural robots presented unique challenges, inspiring a variety of methods.^9,12 In this regard, Cardoso and Furuie¹³ modeled a guidewire using a series of elastically connected vertices. The series of vertices would model the backbone of a flexible catheter. Their model was fairly accurate for small deformations, that is, deflection less than 15°, but would become erroneous at large deformations.

Several methods for modeling large deflections of flexures have been proposed, for example, the elliptical integral approach,¹⁴ numerical integration approach with iterative shooting techniques,^15,16 incremental finite element or finite difference method with Newton–Raphson iteration techniques,^17,18 mechanics-based static and dynamic models,^19,20 Cosserat rod model approach,²¹ and data-driven method,^22–24 These approaches have been adopted in various studies on flexural robots.

For example, Khoshnam et al.¹¹ used the beam theory to simulate the distal portion of a catheter. On the other hand, traditional beam theories are not immediately applicable to large deflections. With deflection angles larger than 10°–15°, the beam theory becomes erroneous. Further, Cotin et al.²⁵ constructed an incremental finite element model based on the beam theory, assuming the catheter is made up of wire-like segments. Each of these approaches addressed specific aspects of large deformation modeling but also came with limitations such as computational inefficiency or limited scalability.²⁶

Two particular approaches that gained attention were the Piecewise Constant Curvature (PCC) and variable curvature models.^27,28 The PCC method simplified the robot into segments with constant curvature, easing control and planning but potentially overlooking nuances in behavior under significant deformation. Variable curvature models, offering a more dynamic representation, provided greater accuracy in complex environments but at the cost of increased computational demand and complexity in control strategies.²⁵

For a tendon-driven soft manipulator whose shape is described by the PCC, Camarillo et al.²⁹ presented a comprehensive derivation for both forward and inverse kinematics utilizing parallel algorithms. Sadati et al.³⁰ provided a critical analysis of the PCC parametrization model. Their study indicated that the application of this model encounters limitations, particularly in its standard form. On further examination, it is revealed that the relaxation of the PCC constraints yields more favorable outcomes. This observation underscores the need for a nuanced approach in the implementation of the PCC parametrization, suggesting that its rigid application may not be universally effective across various scenarios.

In recent years, the Cosserat rod theory has emerged as a significant advancement, presenting a unified set of equations for both large and small deformations of flexural robots. While offering a sophisticated modeling approach, Cosserat rod models often lack analytical or quick-solving solutions, posing challenges for real-time applications.^31–33

The continual evolution in soft robot modeling reflects a broader trend within robotics, where the complexity of theoretical tools has exponentially grown, often intersecting with other scientific disciplines such as continuum mechanics, computational mechanics, and computer graphics. These advancements highlight the shift from simplistic models to more complex, accurate representations, acknowledging the trade-off between accuracy and numerical usability for robotic purposes.³⁴

Soft robots, unlike traditional rigid robots, are continuum systems characterized by highly nonlinear partial differential equations that are generally not analytically integrable. This nature of soft robots necessitates a more nuanced understanding of continuum mechanics, where the robot's configuration is parameterized through positional fields dependent on time and material coordinates. Modeling these systems, thus, requires a deep understanding of the dynamic principles and constitutive laws that govern the stress-strain relationships within the continuum medium.³³ The mechanical modeling in soft robotics is marked by an ongoing quest to balance computational feasibility with the need for accurate, detailed representations of complex systems.

From initial discretization strategies to sophisticated models such as the Cosserat rod theory, each method has contributed to a richer understanding of these intricate systems, underscoring the dynamic and ever-evolving nature of research in soft robotics. To alleviate the aforementioned limitations for fast and accurate large deformation and variable curvature computation on flexural robots, we have proposed a new method based on Bezier shape approximation and an analytical solution on a finite set of arcs along the flexure. The method is called the finite arc method (FAM), was preliminarily tested for accuracy and speed, and was compared with finite element analysis as the most accurate method in the literature.

Contributions

The main contributions of this study are: 1.

Derivation of large deformation kinematics of soft flexural robots based on Bezier shape approximation.

Analytical solution of deformation of a single two-dimensional (2D) arc robot subjected to generalized forces and moments.

Proposing a universal FAM model by the generalization of the 2D single arc solution to three-dimensional (3D) flexures.

Implementation and validation of the proposed FAM model for representative 2D and 3D robot actuation scenarios.

Experimental validation of the proposed FAM model for a surgical interventional use-case, that is, robotic cardiovascular catheter.

Methodology

In this section, first the derivation of the kinematics, constitutive equation, and moment balance of a single-segment beam with PCC deformation assumption (single arc) is provided. Afterward, the proposed FAM model is constituted and analytically solved by aggregation of multiple single-arcs hence, extending the PCC assumption. By decreasing the discretization resolution, we have shown that the proposed methodology exhibits high computational performance and accuracy in modeling large beam deformations with variable curvature.

Kinematics: Single arc

The kinematics derivation of the flexural robot was based on the discretization of the flexure into several finite arcs, each represented with a third-order Bezier curve that satisfied a single-plane constant bending radius condition. Single Bezier arc kinematics for a flexural robot has been previously derived and validated in Ref.⁹ Figure 2 depicts the deformed shape of a flexural robot as a slender beam. In a previous study, the author showed that given the tip position of a circular arc in bending plane $Γ$ , that is, $(x_{\circ}, y_{\circ})$ , the shape of the flexure's centerline $c (s)$ can be accurately approximated with:

FIG. 2.

Schematic deformed shape of the tendon-driven continuum arm Bezier spline centerline in Γ-plane.⁹

c (s) = (\begin{matrix} x (s) \\ y (s) \\ z (s) \end{matrix}) = (\begin{matrix} p_{0} & p_{1} & p_{2} & p_{3} \end{matrix}) (\begin{matrix} {(1 - s)}^{3} \\ 3 s {(1 - s)}^{2} \\ 3 s^{2} (1 - s) \\ s^{3} \end{matrix}),

(1)

where $s \in [0, 1]$ is the normalized length parameter, $p_{i} \in Γ$ are the control points determining the shape of $c (s)$ , and $Γ$ -plane is the bending plane. Thus, assuming length-preservation, the radius of curvature of the arc, r, and bending angle $φ$ were: $r = \frac{x_{\circ}^{2} + y_{\circ}^{2}}{2 y_{\circ}}, φ = \frac{L}{r} .$ (2)

Also, the control points in the $Γ$ -plane would be of the form:

with $x_{\circ}$ and $y_{\circ}$ the horizontal and vertical coordinates of the tip in the bending plane, $Γ$ , and $λ \in [0, 1]$ a constant controlling the curvature throughout shape:⁹

which is only a function of $φ$ . From the geometry in Figure 2, the tip point was parameterized in terms of $φ$ as: $x_{\circ} = \frac{L sin φ}{φ}, y_{\circ} = \frac{L (1 - cos φ)}{φ} .$ (5)

For the described bending, the only kinematic variable was the geodesic bending curvature $κ_{g}$ :¹⁰ $κ_{g} (s, φ) = - \frac{x' (s, φ) y'' (s, φ) - y' (s, φ) x'' (s, φ)}{{({x'}^{2} (s, φ) + {y'}^{2} (s, φ))}^{\frac{3}{2}}},$ (6)

with $(.)' = \frac{d}{d s} (.)$ as the derivation operator. For demonstration, Figure 2 depicts circular arcs, their Bezier fittings, and their corresponding $κ_{g}$ -s. As shown, the proposed shape interpolation resulted in fair curvature approximations. However, for bendings with $| φ | > π$ , adoption number of control points higher than 3 would be necessary for accurate curvature approximation.⁹

The Bezier curve approximation used in the FAM represents a significant simplification when compared with the PCC model. In the FAM, a single parameter, the bending angle, is employed to define the arc of the soft body. This approach contrasts with the PCC model, where the radius of curvature or bending angle for each segment is considered independently.

The Bezier curve provides a more intuitive and computationally efficient means of capturing the continuous nature of bending in soft robots. Unlike the PCC model, which requires handling each segment's curvature separately, the Bezier curve approach allows for a smoother and more integrated representation of the entire structure, facilitating faster computations without compromising the accuracy required for applications such as surgical robotics. In addition, using a Bezier shape approximation allows for choosing various classes of such hodograph curves, for example, twist-free, minimal twist, length-preserving, etc.,³⁵ which would allow for different reduced order modeling of soft robots.

Constitutive equation

The Euler-Bernoulli constitutive relationship between the bending deformation, that is, $κ_{g} (s) - κ_{g}^{⋆} (s)$ , with $κ_{g}^{⋆} (s)$ as the initial curvature, and the internal moment, $m (s)$ along the continuum arm was: $m (s) = E I (\begin{matrix} 0 \\ 0 \\ κ_{g} (s) - κ_{g}^{⋆} (s) \end{matrix}),$ (7)

with $E I$ as the flexural rigidity of the flexure. The constitutive equation implies that the internal bending moment has only an out-of-plane component perpendicular to $Γ$ (in this case, XY-plane). Also, by assuming the initial shape of the flexure as a straight line, $κ_{g}^{⋆} (s) = 0$ was assumed.

Moment balance: Single arc

For the sake of generality, it was assumed that flexure bears an external force and an external moment at the tip. This loading is valid for tendon-driven catheters as a representative example. Figure 4 shows the free-body diagram (FBD) of the flexure at equilibrium. The flexure is bent as a result of an external tip force $f^{T i p}$ and an external tip bending moment $m^{T i p}$ .

In the case of tendon-driven catheters, the tip force and moment are the results of the pulling in the catheter's steering tendons. Quasi-static equilibrium of the FBD depicted in Figure 4 shows that the force and moment at the base of the flexure applied to it from its anchorage are: $f_{0} = - f^{T i p},$ (8) $m_{0} = - m^{T i p} - (\begin{matrix} x_{\circ} \\ y_{\circ} \\ 0 \end{matrix}) \times f^{T i p} .$ (9)

Also, the moment balance equation at any given point $p (s, φ) = {(\begin{matrix} x (s, φ) & y (s, φ) & 0 \end{matrix})}^{T}$ throughout the flexure is: $m (s, φ) = m^{T i p} + (\begin{matrix} x_{\circ} - x (s, φ) \\ y_{\circ} - y (s, φ) \\ 0 \end{matrix}) \times f^{T i p} .$ (10)

Substituting the constitutive definition of $m (s)$ from Equation (7), the balance equation, which relates the shape of the flexure to the tip force and moment, was obtained as: $m^{T i p} + (\begin{matrix} x_{\circ} - x (s, φ) \\ y_{\circ} - y (s, φ) \\ 0 \end{matrix}) \times f^{T i p} = E I (\begin{matrix} 0 \\ 0 \\ κ_{g} (s, φ) \end{matrix}) .$ (11)

This equation is a vectorial equation with two trivial scalar equations ( $0 = 0$ ) from the first and second rows and a nonlinear bi-variable equation in terms of s and $φ$ from the third row: $m^{T i p} + f_{y} (x_{\circ} (φ) - x (s, φ)) - f_{x} (y_{\circ} (φ) - y (s, φ)) = E I κ_{g} (s, φ),$ (12)

with $m^{T i p}$ as the external tip moment applied through pulling the steering tendons, f_x, and f_y as the horizontal and vertical components of the tip force in the bending angle, respectively. Also, the tip bending moment was known in real time by measuring the torque of the driving motor. In this equation, the only unknown is the bending angle $φ$ .

Deformation solution

As stated earlier, Equation (12) is under-determined with two unknowns $φ$ , s and one equation. However, the balance equation must be satisfied for all the points along the flexure's length, that is, $s \in [0, 1]$ . Therefore, the balance equation was formulated as a minimization problem over $s \in [0, 1]$ domain. It is noteworthy that to find the deformation along the flexure, the solution starts from the last segment and then propagates to the next segments by transferring the internal moments and internal forces by knowing the deformation of the next segments.

To the best of the authors' knowledge, this approach is unprecedented in the literature. To find the $φ$ in a given segment, functional $Π (φ)$ was defined as:

Functional $Π (φ)$ represents the summation of the residual momentum balance equation along the length of the catheter. Minimization of $Π (φ)$ , in fact, enforces the balance of forces and moments [Eq. (12)] along the length of the flexure, leading to an accurate representation of the deformation. The functional was further simplified to Equation (17) by introducing: $ξ_{1} (s, φ) = f_{y} (x_{\circ} (φ) - x (s, φ)),$ (14) $ξ_{2} (s, φ) = - f_{x} (y_{\circ} (φ) - y (s, φ)),$ (15)

ξ_{3} = m^{T i p} - E I κ_{g} (s, φ),

(16)

Afterward, the functions $ξ_{1, 2, 3} (s, φ)$ were replaced with their Taylor expansion around $s = \frac{1}{2}$ (the mid-point of the integral domain on s) with $N = 15$ terms. Without loss of generality, $N = 15$ was selected empirically. This expansion led to polynomial representation of the integrand in Equation (13) in terms of s.

where Given the polynomial representation of $ξ_{i}$ , the integrand was polynomial of degree $2 N = 30$ ( $2 N$ because the integrand was a complete square). The integral was carried out analytically over $s \in [0, 1]$ . The analytical solution was of the form:

with $α_{k} (φ)$ nonlinear co-efficients of the resulting polynomial. Afterward, by reversing the order of summation over k and integration over s, Equation (19) is simplified as:

Imposing the necessary condition for minimizing $Π (φ)$ to find the unknown bending angle $φ$ , Equation (20) is simplified as:

To solve Equation (21), the analytical form of ${(\sum_{k = 0}^{2 N} \frac{1}{k + 1} \frac{\partial α_{k} (φ)}{\partial φ})}^{2}$ was minimized using a fast graphics processing unit (GPU)-based particle swarm optimizer in C++ programming language that led to high computational efficiency, especially for real-time applications.

Finite arcs: 3D generalization

In general, the deformation occurring along the flexure is not planar if forces and moments are 3D. Nevertheless, we hypothesized that the general shape of the flexure could be obtained by the aggregation of finite deformations of discretized arcs. We assumed that each arc would undergo a planar deformation and a torsional deformation.

Thus, each arc would necessarily bend in a plane perpendicular to the total bending moment vector acting on it. Thus, in the following, we describe finding the moment plane ( $Γ$ -plane), calculating the deformation of each arc in that plane, and then mapping it to the global coordinate system (GCS).

As shown in Figure 3, the $Γ$ -plane with normal $\hat{n}$ creates a new coordinate system, denoted as bending coordinate system (BCS), in which the deformation of the flexure occurs. Normal vector $\hat{n}$ has been obtained in Equation (23). The normal vector of the bending plane is derived from the direction of the total bending moment. The total bending moment is defined as the cumulative effect of all bending moments along the soft robot's structure.

FIG. 3.

A representative continuum arc before bending (straight shape) with external loading. The arc's deformation occurs in the bending plane Γ denoted with normal n in response to the total external force f_t and bending moment m_t.

This approach assumes that the normal to the bending plane is parallel to the total bending moment vector. We chose this assumption based on the observation that, in most practical applications, the dominant bending moments align with the plane's normal. Equation (19) reflects this by excluding a component in the x-direction. This exclusion is justified as, typically, bending moments inducing significant deformation in soft robots, such as catheters primarily occur in the plane perpendicular to the robot's primary axis.

This approach differs from the PCC model, where an additional parameter defines the bending plane, offering a more streamlined and computationally efficient method while still capturing the essential dynamics of bending. As a result, the first step was to obtain the bending plane and, by using the normal of this plane, define a new coordinate system in which the deformation is calculated.

The normal vector of bending plane $Γ$ can be found from the total bending moment vector m _t . $m_{t} = (\begin{matrix} 0 \\ m_{y} \\ m_{z} \end{matrix}),$ (22)

while the x-axis is selected along the flexure, and y and z-axes were selected arbitrarily. Therefore, the normal vector $\hat{n}$ of the bending plane $Γ$ was: $\hat{n} = \frac{m_{t}}{\sqrt{m_{t} \cdot m_{t}}},$ (23)

with $\cdot$ representing the dot product of two vectors. By selecting the x-axis as the second local axis ( $\hat{u}$ ) of a new coordinate system, the third axis was calculated as given next: $\hat{v} = \hat{n} \times \hat{u},$ (24)

with $\times$ representing the cross-product of two vectors. Equation (21) requires the in-plane component of the total force at the tip of the flexure, that is, $f^{T i p}$ as well as the total tip bending moment $m^{T i p}$ that using force and moment projection onto the bending plane were obtained as: $f^{T i p} = f_{Γ} = (\begin{matrix} (f - f_{n}) \cdot \hat{u} \\ (f - f_{n}) \cdot \hat{v} \\ 0 \end{matrix}) .$ (25) $m^{T i p} = m_{t} \cdot \hat{n} .$ (26)

Using the vectors $\hat{u}, \hat{v}, \hat{n}$ , the transformation matrix from BCS to GCS can be determined as: $^{g} R_{b} = (\begin{matrix} \hat{u} & \hat{v} & \hat{n} \end{matrix}),$ (27)

which was used to transfer the deformation from BCS to GCS. The rotation matrix can also be obtained using Rodrigues' rotation formula, which leads to the same result.

Torsion

An introducer sheath is used to introduce a catheter into the vasculature, generating friction that hinders mobility. This friction demonstrates the minimum force required to move the catheter during an interventional operation.³⁶ During interventional operation, when a steerable catheter has one or two degrees of freedom, the tip cannot bend in all directions at the same time. Therefore, a twisting moment (torque) must be applied to overcome vascular friction during the operation and rotate the catheter. The torsional deformation, despite its nature, results in a planar effect due to the specific geometry and constraints of the soft robotic structure, particularly in the context of flexural robots such as catheters.

To include the torsion in the solution schema, the continuum robot was assumed to have a uniform cross-section along its length. In this case: $θ = \frac{L}{G J} T,$ (28)

where $θ$ is the angle of twist in radians, T is the applied torsion moment (torque), L is the segment length, G is the Modulus of rigidity (shear modulus) of the material, and J is the torsional constant. The torsional rigidity (GJ) is selected as a function of flexural rigidity (EI), which is a property of the continuum robot and is generally constant as given next: $G J = \frac{E I}{1 + V} .$ (29)

Therefore, the torsional rotation matrix R _t around local axis $\hat{u}$ can be defined as: $R_{t} = (\begin{matrix} 1 & 0 & 0 \\ 0 & cos θ & - sin θ \\ 0 & sin θ & cos θ \end{matrix}) .$ (30)

Hence, the total rotation matrix to transfer the deformation, including moment and torsion, can be defined as: $R = {R_{t}}^{g} R_{b} .$ (31)

Deformation solution

Based on the solution method provided in Equation (21), it was hypothesized to discretize the flexure into N segments that each segment's deformation would be a single bending radius flexure. This way, the total deformation of the flexure would be obtained by aggregating all the deformations from the base $s = 0$ to the tip $s = 1$ of the flexure. Since the flexure would undergo large deformation, the internal moment and force at the boundaries of each segment would be obtained by translation and rotation of the forces on the segments next to it to the endpoints of the segment, and then the deformation of the segment would be solved using Equation (21). Thus, the deformation solution starts from the last segment and then propagates to the segments before it by transferring the internal moments and internal forces by knowing the deformation of the next segments.

Figure 4 depicts the proposed segmentation schematically. To the best of the author's knowledge, this approach is unprecedented in the literature. This angle for the next section $(φ_{i - 1})$ can be determined with the same process by transferring the resultant force and moment to section $i - 1$ . For this transformation, first, the rotation matrix of i to $i - 1$ and translation matrix of endpoint in coordinate $i - 1$ are as given next: $^{i - 1} R_{i} = 2 (\begin{matrix} {q_{0}}^{2} + {q_{1}}^{2} - \frac{1}{2} & q_{1} q_{2} - q_{0} q_{3} & q_{1} q_{3} + q_{0} q_{2} \\ q_{1} q_{2} + q_{0} q_{3} & {q_{0}}^{2} + {q_{2}}^{2} - \frac{1}{2} & q_{2} q_{3} - q_{0} q_{1} \\ q_{1} q_{3} - q_{0} q_{2} & q_{2} q_{3} + q_{0} q_{1} & {q_{0}}^{2} + {q_{3}}^{2} - \frac{1}{2} \end{matrix}),$ (32)

FIG. 4.

The FBD of a deformed flexural robot with multiple circular arc segments. FBD, free-body diagram.

where: $(\begin{matrix} q_{0} \\ q_{1} \\ q_{2} \\ q_{3} \end{matrix}) = (\begin{matrix} cos \frac{φ_{i}}{2} \\ sin \frac{φ_{i}}{2} {\hat{n}}_{x} \\ sin \frac{φ_{i}}{2} {\hat{n}}_{y} \\ sin \frac{φ_{i}}{2} {\hat{n}}_{z} \end{matrix}) .$ (33) $^{i - 1} p_{i} =^{i - 1} R_{i} (\begin{matrix} x_{i} \\ y_{i} \\ z_{i} \end{matrix}) =^{i - 1} R_{i} (\begin{matrix} r_{i} sin φ_{i} \\ r_{i} (1 - cos φ_{i}) \\ 0 \end{matrix}) .$ (34)

As a result, the amount of force and moment exerted on the endpoint of section $i - 1$ is obtained from the resultant force and moment from section i in addition to the external force and moment applied directly on section $i - 1$ , which is formulated in Equations (35) and (36). $f_{i - 1} =^{i - 1} R_{i} f_{i} + f_{i - 1}^{e x t} .$ (35)

m_{i - 1} =^{i - 1} p_{i} \times f_{i} + Σ_{k = i}^{N} m_{k} + Σ_{k = i}^{N} m_{k}^{e x t} .

(36)

This process continues in a similar way to transmit force and momentum to all sections from the beam's free end to the base. In the end, the final deformed shape of the flexure was constructed by augmenting the deformed shape of segment $i + 1$ at the end of segment i. Eventually, the deformed position of the endpoint of the n-th segment, P _n in the GCS was obtained as: $^{0} p_{n} = \frac{L_{n}}{φ_{n}} {\prod_{i = 1}^{n}}^{i - 1} T_{i} (φ_{i}) (\begin{matrix} sin φ_{n} \\ (1 - cos φ_{n}) \\ 0 \end{matrix}),$ (37)

where L $_{i}$ -s are the length of the i-th segment. Also, $T_{i} (φ_{j})$ represents the transformation matrix for each segment, which includes both rotation and translation components structured as given next:

Our method of initiating the solution from the last segment and proceeding backward is highly effective for analyzing problems with pose-independent external forces, such as actuation forces and moments in tendon-driven catheters and soft robots. We acknowledge this limitation when dealing with pose-dependent forces, such as gravity or contact forces whose component resolution in the “deformed” coordinate systems are unknown when solving Equation (21).

Nevertheless, for many intraluminal procedures such as those in cardiovascular applications, the catheter's weight is negligible compared with the actuation forces; thus, the proposed method may not lead to significant errors. We have shown the validity of this speculation in the presented validation studies of this work. To alleviate this limitation in the future, an implicit solution schema of a system of Equation (21)-s for all the segments could be investigated.

Validation Studies

Study I: Simulation-based validation

Protocol

Figure 5 shows eight planar and eight spatial load-cases solved by FAM and compared with FEM. The developed FE models were of 40 × 6 mm size. The flexural rigidity of the robots was 750 Nmm². The FE simulations were performed in Abaqus CAE 2020—Student Edition.

FIG. 5.

Eight load cases on a representative flexure for Study I in (a) 2D space and (b) 3D space.

Derivation of the analytical form of Equation (21) was performed in Maple 2019 (Maplesoft; Cybernet Systems Co., Ltd.), and FAM was developed in C++ programming language. For linear algebraic operations and numerical calculations, GPU-accelerated Armadillo library was used. The FE model of the robots was meshed with quadratic beam elements of 50 elements.

The FE results were not sensitive to mesh refinement for the number of elements used in the simulation. Also, the number of segments for the FAM solution was $N = 12$ . The rationale for this selection and the sensitivity of the results of FAM to the segmentation resolution are discussed in the Results section.

Results

To compare the accuracy and computational performance of the FAM method with FEM, two series of loadings that represent various extrinsic loadings on soft robots, for example, tendon-driven robots, were performed, that is, planar deformation (where loadings were co-planar) and spatial (where loadings were not co-planar).

Planar deformation

Figure 6a shows the deformed shape of the catheters after applying force and moment in eight different cases (Fig. 5a). In each load case, each segment was first $λ$ obtained using the proposed shape interpolation method [Eq. (4)]. Table 1 summarizes the results of the deformation simulation.

FIG. 6.

(a) comparison of FAM and FEM results in deflection estimations for eight load cases on soft flexures, (b) distribution of error of deformation estimation with FEM as the reference along the x-direction and (c) y-direction. FAM, finite arc method; FEM, finite element method.

Table 1.

Comparison of Finite Arc Method and Finite Element Method Solutions for Eight Load Cases in Two-Dimensional Simulation

Case	Tip displacement x^Tip (mm)			Tip displacement y^Tip (mm)
Case	FEM	FAM	Error	FEM	FAM	Error
1	11.35	11.53	−0.18	−31.55	−31.61	+0.06
2	33.51	33.16	+0.35	7.65	7.09	+0.56
3	28.17	29.67	−1.50	−15.83	−15.90	+0.07
4	−5.12	−6.38	+1.27	−7.79	−8.76	+0.97
5	39.53	39.26	+0.27	5.30	6.91	−1.60
6	36.84	37.42	−0.58	−14.21	−13.35	−0.86
7	33.96	34.11	−0.15	17.88	17.66	+0.22
8	32.53	33.81	−1.28	−18.98	−17.56	−1.41
Mean ± SD	—	—	−0.27 ± 0.90	—	—	−0.25 ± 0.94

FAM, finite arc method; FEM, finite element method; SD, standard deviation.

The comparison of the estimated tip position with the reference showed that the tip position estimation along the x-axis had an average error of −0.27 ± 0.90 mm and along the y-axis had an average error of −0.25 ± 0.94 mm. The mean-absolute-error (MAE) of the deflection estimation with respect to the reference was 0.29 ± 0.48 along x direction and 0.39 ± 0.58 along y direction with root-mean-square of 0.4865 and 0.5732 mm respectively (Fig. 6b, c).

Spatial deformation

Figure 5b shows eight different load-cases in the 3D simulation. Eight simulations were performed and compared with the FEM (Fig. 7 and Table 2). The catheter model was similar to before, with flexural rigidity of 1500 Nmm². The comparison of the estimated tip position with the reference showed that the tip position estimation along the x-axis had an average error of +0.36 ± 0.98 mm and along the y-axis had an average error of −0.18 ± 0.56 mm, and along the z-axis had an average error of −0.13 ± 0.37 mm.

FIG. 7.

Comparison of FAM and FEM results in deflection estimations for eight load cases on soft flexures in 3D simulation.

Table 2.

Comparison of Finite Arc Method and Finite Element Method Solutions for Eight Load Cases in Three-Dimensional Simulation

Case	Tip displacement x^Tip (mm)			Tip displacement y^Tip (mm)			Tip displacement z^Tip (mm)
Case	FEM	FAM	Error	FEM	FAM	Error	FEM	FAM	Error
1	24.38	24.37	+0.01	15.47	14.48	−0.01	−21.92	−21.90	−0.02
2	28.49	28.52	−0.03	7.36	7.49	−0.13	24.00	23.97	+0.03
3	20.66	19.07	+1.59	21.20	22.12	−0.92	21.61	21.91	−0.30
4	28.24	29.65	−1.41	25.68	24.23	+1.45	−0.12	0.31	−0.43
5	21.55	19.03	+2.52	20.70	21.69	−0.99	20.79	21.69	−0.99
6	31.77	31.74	+0.03	−4.61	−4.56	−0.05	20.77	20.89	−0.12
7	20.15	20.17	−0.02	19.65	19.64	+0.01	−19.65	−19.64	−0.01
8	37.16	36.98	+0.18	6.19	6.96	−0.77	−6.19	−6.96	+0.77
Mean ± SD	—	—	+0.36 ± 0.98	—	—	−0.18 ± 0.56	—	—	−0.13 ± 0.37

The MAE of the deflection estimation with respect to the reference was 0.44 ± 0.28 mm along the x direction, 0.52 ± 0.20 mm along the y direction, and 0.51 ± 0.14 mm along the z direction with root-mean-square of 0.5222, 0.5590, and 0.5284 mm respectively.

Segmentation effect

To investigate the segmentation effect, load case 4 was selected as the worst-case scenario and solved with a series of segmentation numbers. Figure 8a shows a representative set of final shapes of the deformed catheter with various segmentations. The segmentation was changed between 1 and 5000, and the absolute error of the flexure's tip position was studied (Fig. 8b). It was observed that by increasing the number of segments in the model from 1, the error dropped drastically and hit a minimum of 12 segments.

FIG. 8.

Effect of segmentation, (a) FAM deflection solution for a different number of segments and (b) absolute tip displacement error for a different number of segments.

The absolute error value was less than 3 mm (7.5% of length) between N = 12 and 500. For N > 500, the error steadily increased. The reason might be because each arc's deformation is obtained through minimization. At extreme segmentation (N > 500), the arcs tend to straight lines; thus, the residual error of the minimization increases as it cannot find a finite number for the curvature ( $φ \to 0$ ).

Torsional warping effect

Load case 1 in the 3D simulation was selected to investigate the effect of torsion and solved by increasing the torsion moment from 0 to 50 mN.m (Fig. 9). Table 3 summarizes the tip position error for different torque values applied to the continuum robot's tip. Also, by defining normalized torsion $\hat{τ}$ as:

FIG. 9.

(a) comparison of FAM versus NL-FEM results with tip torques and (b) variation of tip displacement error with tip torques. NL-FEM, non-linear FEM.

Table 3.

Comparison of Finite Arc Method and Finite Element Method Solutions for Different Values of Torsion Applied to the Tip of the Continuum Robot

T (mNm)	e_x (mm)	e_y (mm)	e_z (mm)	\|\|e\|\| (mm)
0	+0.01	−0.01	−0.00	0.0141
5	+0.21	−0.22	−0.01	0.3043
10	+0.42	−0.41	−0.07	0.5911
15	+0.66	−0.64	−0.15	0.9315
20	+0.89	−0.81	−0.26	1.2312
25	+1.12	−0.97	−0.39	1.5321
30	+1.37	−1.11	−0.53	1.8412
35	+1.63	−1.22	−0.70	2.1530
40	+1.88	−1.31	−0.84	2.4405
45	+2.14	−1.39	−1.02	2.7481
50	+2.41	−1.45	−1.16	3.0424

\hat{τ} = \frac{T}{m g r_{g}},

(39)

with m, g, and r_g as the segment's mass, gravity, and nominal radius of gyration before deformation, the applied torque was normalized. Similarly, the tip position error was normalized with respect to the robot's length. It has been observed that by increasing the twist angle, the error increased because of exiting the assumption of decoupling solution.

Study II: Experimental validation

Setup and protocol

The proposed method was integrated with a representative robot-assisted intervention (RCI) system previously developed by the authors in Refs.^37,38 to test the feasibility of integrating the proposed method with RCI systems (Fig. 10a). The user would transmit the measured rotation and insertion length to the patient module. As shown in Figure 10a, the CataLyst robot (Thermo Fisher, WA, USA) was used to align the entry angle of the femoral vein of the vascular model (SAM plus; Lake Forest Anatomical, Inc., IL, USA).

FIG. 10.

Setup used in validation study-II (a) patient module, (b) deflected catheter in the vascular model with the external forces, and (c) comparison of FAM and ground truth deflections in Study II.

The forces and their location were extracted using an image-based method developed by the authors in Ref.³⁸ The total length of the catheter was 504 mm with a hollow cylinder cross-section with an outer diameter of 2.14 and a thickness of 0.29 mm. Its material is linear elastic with Young's modulus of 1592.3 MPa, and Poisson's ratio of 0.395.³⁸ Three reaction forces were captured in the deformed shape of the catheters depicted in Figure 10b.

The first force f₁ had a magnitude of $f_{1} = 0.463 N$ located at $s = 0.165$ , the second force $f_{2} = - 0.065 N$ was located at $s = 0.577$ , and the third force was $f_{3} = 0.007 N$ located near the tip $s = 0.968$ . The estimated forces from the deformed shape of the catheter were used in the FAM model of the catheter, and the resultant deflections were compared.

The estimation of contact forces from our experimental setups involved a detailed analysis of the deformed shapes of the catheters under applied forces. Utilizing an image-based method developed in our previous work,^38,39 we captured the deformation patterns of the catheters. From these patterns, the location and magnitude of contact forces were inferred. This estimation process involves correlating specific deformations with known force applications, allowing us to reverse-engineer the forces responsible for observed deformations. This method provides a robust approach to understanding and quantifying the interaction forces between the catheter and its environment, which is crucial for the accurate modeling of soft robotic systems in real-world applications.

Results

Figure 10c compares the magnitude of the reference with the estimated deformation of a catheter. As the comparison showed, the FAM results were in fair agreement with the observed deformation. The maximum absolute error of deflection along vertical axis y was 4.95 mm (8.5% of the maximum deflection), and it was at the tip. The FAM method resulted in a larger estimation error at the tip, which might be due to error accumulation through the segments along the flexure.

A similar trend was observed in Study I on the simulated deformation. Nevertheless, the average relative error of the experiments was 1.41 ± 1.47 mm. Another finding was that the average computational time of the proposed force estimation method was less than 3 mseconds with C++ implementation. These findings showed that the proposed method could successfully simulate the deformation of flexures in real time.

Table 4 compares the performance of the proposed method with representative studies. The comparison showed that the proposed method exhibited improvement in both accuracy and computational cost. Also, the validation studies showed acceptable accuracy of the proposed method in solving the deformation under various loading conditions. In future studies, this work will be extended to accommodate nonlinear constitutive equations in the solution of the FAM.

Table 4.

Comparison of Accuracy and Computation Time of Finite Arc Method and Representative Studies

Study	Method	Dimensions	Position error (% of length)	Computation time (mseconds)
Roshanfar et al.⁴¹	Cosserat rod	3D	4.2%	27
Bezawada et al.⁴⁹	Piecewise Pythagorean Hodographs	2D	5.0%	—
Garg et al.⁵⁰	Piecewise Constant Curvature	3D	4.0%	—
Till et al.⁵¹	Cosserat rod	3D	3.0%	200
This study	FAM	2D	0.50%	3
This study	FAM	3D	0.55%	4

2D, two-dimensional; 3D, three-dimensional.

The authors have previously proposed such integration with Cosserat rod models in Refs.^40,41 The proposed method also can be used in sensor-based^42–46 and sensor-free bending-based^10,38 force estimation for haptic rendering applications, for example, for RCI systems.^47,48

Conclusion

In this paper, a new modeling method and its solution for soft flexural robots, such as tendon-driven catheters, with large deformation using the FAM was proposed and validated. The fast solution of the model solution was made possible by using the proposed Bezier spline shape approximation and discretizing the catheter into a finite set of circular arc segments. Also, combined with an RCI setup, the proposed modeling approach exhibited an acceptable performance in terms of accuracy for real-time force estimation on flexible catheters in contact with tortuous vascular models. The proposed method showed significant improvement in accuracy and reduction of computational cost for direct and inverse mechanistic problems pertaining to the deformation of flexural robots.

Footnotes

Author Disclosure Statement

No competing financial interests exist.

Funding Information

This study was supported by McGill University, Research Institute of McGill University Health Center (RI-MUHC) [Grant Number 10910], Fond de recherche du Quebec—nature et technologies (FRQNT) [Grant Number 315939], and Montreal General Hospital (MGH) Foundation [Research Award 9204].

Abbreviations Used

References

Jolaei

, Hooshiar

, Dargahi

, et al. Toward task autonomy in robotic cardiac ablation: Learning-based kinematic control of soft tendon-driven catheters. Soft Robotics, 2020; 8(3):340–351.

Jeon

, Hoshiar

, Kim

, et al. A magnetically controlled soft microrobot steering a guidewire in a three dimensional phantom vascular network. Soft Robotics, 2019; 6(1):54–68.

Verma

, Terricabras

, Verma

. The cutting edge of atrial fibrillation ablation. Arrhythm Electrophysiol Rev, 2021; 10(2):101.

Jolaei

, Hooshiar

, Sayadi

, et al. Sensor-Free Force Control of Tendon-Driven Ablation Catheters Through Position Control and Contact Modeling. In 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). IEEE; 2020; pp. 5248–5251.

Hooshiar

, Najarian

, Dargahi

Haptic telerobotic cardiovascular intervention: A review of approaches, methods, and future perspectives. IEEE Rev Biomed Eng, 2020; 13:32–50, 2020.

Patel

, Shah

, Pancholy

. Long distance tele-robotic-assisted percutaneous coronary intervention: A report of first-in-human experience. EClinicalMedicine, 2019; 14:53–58.

Shabana

AA.

Flexible multibody dynamics: Review of past and recent developments. Multibody Syst Dyn, 1997; 1(2):189–222.

Barrientos-Diez

, Dong

, Axinte

, et al. Real-time kinematics of continuum robots: Modelling and validation. Robot Comput Integr Manuf, 2021; 67:102019.

Hooshiar

, Sayadi

, Jolaei

, et al. Accurate Estimation of Tip Force on Tendon-Driven Catheters Using Inverse Cosserat Rod Model. In: 2020. International Conference on Biomedical Innovations and Applications (BIA). IEEE; 2020; pp. 37–40.

10.

Hooshiar

, Sayadi

, Jolaei

, et al. Analytical Tip Force Estimation on Tendondriven Catheters Through Inverse Solution of Cosserat Rod Model. In: 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE; 2021; pp. 1829–1834.

11.

Khoshnam

, Skanes

, Patel

. Modeling and estimation of tip contact force for steerable ablation catheters. IEEE Trans Biomed Eng, 2015; 62(5):1404–1415.

12.

Valentini

, Pennestri

Modeling elastic beams using dynamic splines. Multibody Syst Dyn, 2011; 25(3):271–284.

13.

Cardoso

, Furuie

. Guidewire path determination for intravascular applications. Comput Methods Biomech Biomed Eng, 2016; 19(6):628–638.

14.

Dado

, Al-Sadder

A new technique for large deflection analysis of non-prismatic cantilever beams. Mech Res Commun, 2005; 32(6):692–703.

15.

Lee

Post-buckling of uniform cantilever column under a combined load. Int J Non-Linear Mech, 2001; 36(5):813–816.

16.

Magnusson

, Ristinmaa

, Ljung

Behaviour of the extensible elastica solution. Int J Solids Struct, 2001; 38(46–47):8441–8457.

17.

Liu

, Poirot

, Franson

, et al. Modeling and validation of the three-dimensional deflection of an MRI-compatible magnetically actuated steerable catheter. IEEE Trans Biomed Eng, 2015; 63(10):2142–2154.

18.

Golley

The solution of open and closed elasticas using intrinsic coordinate finite elements. Comput Methods Appl Mech Eng, 1997; 146(1–2):127–134

19.

Grazioso

, Di Gironimo

, Siciliano

A geometrically exact model for soft continuum robots: The finite element deformation space formulation. Soft Robotics, 2019; 6(6):790–811.

20.

Yuan

, Li

Workspace analysis of cable-driven continuum manipulators based on static model. Robot Comput Integr Manuf, 2018; 49:240–252.

21.

Spillmann

, Teschner

Corde: Cosserat Rod Elements for the Dynamic Simulation of One Dimensional Elastic Objects. In: Proceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium on Computer Animation; 2007; pp. 63–72.

22.

Sayadi

, Nourani

, Jolaei

, et al. Force Estimation on Steerable Catheters Through Learning-From Simulation with Ex-Vivo Validation. In: 2021 International Symposium on Medical Robotics (ISMR). IEEE; 2021; pp. 1–6.

23.

, Wang

, Ge

, et al. Analytical modeling and design of generalized pneunet soft actuators with three-dimensional deformations. Soft Robotics, 2021; 8(4):462–477.

24.

, Fernandez

JdG

, Tan

Probabilistic kinematic model of a robotic catheter for 3D position control. Soft Robotics, 2019; 6(2):184–194.

25.

Cotin

, Duriez

, Lenoir

, et al. New Approaches to Catheter Navigation for Interventional Radiology Simulation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer; 2005; pp. 534–542.

26.

Tunay

Spatial Continuum Models of rods undergoing large deformation and inflation. IEEE Trans Robot, 2013; 29(2):297–307.

27.

Chawla

, Frazelle

, Walker

. A Comparison of Constant Curvature Forward Kinematics for Multisection Continuum

Manipulators

. In: 2018 Second IEEE International Conference on Robotic Computing (IRC). IEEE; 2018; pp. 217–223.

28.

Webster III

, Jones

. Design and kinematic modeling of constant curvature continuum robots: A review. Int J Robot Res, 2010; 29(13):1661–1683.

29.

Camarillo

, Milne

, Carlson

, et al. Mechanics modeling of tendon-driven continuum manipulators. IEEE Trans Robot, 2008; 24(6):1262–1273.

30.

Sadati

, Naghibi

, Shiva

, et al. A geometry deformation model for braided continuum manipulators. Front Robot AI, 2017; 4:22.

31.

Janabi-Sharifi

, Jalali

, Walker

. Cosserat rod-based dynamic modeling of tendon-driven continuum robots: A tutorial. IEEE Access, 2021; 9:68703–68719.

32.

Sadati

, Shiva

, Herzig

, et al. Stiffness imaging with a continuum appendage: Real-time shape and tip force estimation from base load readings. IEEE Robot Autom Lett, 2020; 5(2):2824–2831.

33.

Renda

, Armanini

, Lebastard

, et al. A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation. IEEE Robot Autom Lett, 2020; 5(3):4006–4013.

34.

Armanini

, Boyer

, Mathew

, et al. Soft robots modeling: A structured overview. IEEE Trans Robot, 2023; 39(3):1728–1748.

35.

Farouki

RT.

Pythagorean—Hodograph Curves. Springer: Berlin Heidelberg; 2008; pp. 381–391.

36.

Thakur

, Holdsworth

, Drangova

Characterization of catheter dynamics during percutaneous transluminal catheter procedures. IEEE Trans Biomed Eng, 2009; 56(8):2140–2143.

37.

Hooshiar

, Sayadi

, Dargahi

, et al. Integral-free spatial orientation estimation method and wearable rotation measurement device for robot-assisted catheter intervention. IEEE ASME Trans Mechatron, 2021; 27(2):766–776.

38.

Hooshiar

, Bandari

, Dargahi

Image-Based Estimation of Contact Forces on Catheters for Robot-Assisted Cardiovascular Intervention. In: Hamlyn Symposium on Medical Robotics. London, UK: Imperial College London and the Royal Geographical Society; 2018; pp. 119–120.

39.

Hooshiar

Image-Based Force Estimation and Haptic Rendering for Robot-Assisted Car-diovascular Intervention. PhD Thesis, Concordia University; 2021.

40.

Roshanfar

, Dargahi

, Hooshiar

. Toward Semi-Autonomous Stiffness Adaptation of Pneumatic Soft Robots: Modeling and

Validation

. In: 2021 International Conference on Autnomous Systems. IEEE; 2021; pp. 1–5.

41.

Roshanfar

, Sayadi

, Dargahi

, et al. Stiffness Adaptation of a Hybrid Soft Surgical Robot for Improved Safety in Interventional Surgery. In: 2022 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). IEEE; 2022; pp. 4853–4859.

42.

Torkaman

, Roshanfar

, Dargahi

, et al. Embedded Six-Dof Force-Torque Sensor for Soft Robots with Learning-Based Calibration. IEEE Sensors J, 2023; 23(4):4204–4215.

43.

Hooshiar

, Razban

, Bandari

, et al. Sensing Principle for Realtime Characterization of Viscoelasticity in the Beating Myocardial Tissue. In: 2017 IEEE International Conference on Computational Intelligence and Virtual Environments for Measurement Systems and Applications (CIVEMSA). IEEE; 2017; pp. 72–77.

44.

Bandari

, Hooshair

, Packirisamy

, et al. Optical Fiber Array Sensor for Lateral and Circumferential Force Measurement Suitable for Minimally Invasive Surgery: Design, Modeling and Analysis. In: Specialty Optical Fibers. Optical Society of America; Optica Publishing Group; Washington, DC, USA; 2016; p. JTu4A-44.

45.

Bandari

, Hooshair

, Muthukumaran

, et al. Bending-Based Formulation of Light Intensity Modulation for Miniaturization of Optical Tactile Sensors. In: Optical Sensors. Optical Society of America; Optica Publishing Group; Washington, DC, USA; 2018; p. SeM2E-3.

46.

Bandari

, Dargahi

, Packirisamy

Tactile sensors for minimally invasive surgery: A review of the state-of-the-art, applications, and perspectives. IEEE Access, 2020; 8:7682–7708.

47.

Hooshiar

, Payami

, Dargahi

, et al. Magnetostriction-based force feedback for robot-assisted cardiovascular surgery using smart magnetorheological elastomers. Mech Syst Signal Process, 2021; 161:107918.

48.

Hooshiar

, Alkhalaf

, Dargahi

. Development and assessment of a stiffness display system for minimally invasive surgery based on smart magneto-rheological elastomers. Mater Sci Eng C, 2020; 108:110409.

49.

Bezawada

, Woods

, Vikas

. Shape Estimation of Soft Manipulators Using Piecewise Continuous Pythagorean-Hodograph

Curves

. In: 2022 American Control Conference (ACC). IEEE; 2022; pp. 2905–2910.

50.

Garg

, Good

, Revier

, et al. Kinematic Modeling of Handed Shearing Auxetics Via Piecewise Constant Curvature. In: 2022 IEEE 5th International Conference on Soft Robotics (RoboSoft). IEEE; 2022; pp. 423–430.

51.

Till

, Aloi

, Rucker

Real-time dynamics of soft and continuum robots based on cosserat rod models. Int J Robot Res, 2019; 38(6):723–746.