# Continuum concentric push-pull robots: a Cosserat rod model

Journal Title
XX(X):1–27

DOI:

<figure><img src="image_1.png" /></figure>

Matthias Tummers<fnref n="1" />, Frédéric Boyer<fnref n="2" />, Vincent Lebastard<fnref n="2" />, Alexis Offermann<fnref n="1" />, Jocelyne Troccaz<fnref n="1" />, Benoît Rosa<fnref n="3" />, and M. Taha Chikhaoui<fnref n="1" />

## Abstract

Various approaches and structures emerged recently to design continuum robots. One of the most promising designs regards a new concept of continuum concentric push-pull robots (CPPRs) that have the characteristic of combining several key advantages of tendon actuated, multi-backbone, and concentric tube ones (direct curvature actuation, small outer/inner diameter ratio, free lumen, etc.). Geometrically-exact models of such recently introduced robots are yet to be developed to gain leverage of their full potential. This article extends beyond usual definitions of Cosserat rod theory in order to take into account this new type of continuum robots, constituted by sliding rods, in a shape of tubes whose cross-sections are neither uniform nor symmetrical along their entire length. The introduced model is capable of considering versatile design options, external loads, 3D deformations, an arbitrary number of tubes and profiles of the centroid lines, as well as a new actuation method consisting of an input rotation. Numerical simulations and experiments on CPPR prototypes validate our model.

## Keywords

Continuum robots, Cosserat rod modeling, Kinematics, Medical Robots and Systems, Flexible Robotics, Mechanics

## 1 Introduction

Continuum robots are thin, slender flexible structures that can conform to curvilinear paths and are inherently compliant. Such robots have been developed over the last decades, especially in applications that require dexterity in confined, difficult to access spaces. Many designs have been proposed for applications ranging from nuclear decontamination to medical applications (Burgner-Kahrs et al. 2015; Dupont et al. 2022). Concentric agonist-antagonist robots (CPPRs), were introduced recently by Oliver-Butler et al. (2017). Similarly to concentric tube continuum robots (CTRs), CPPRs are composed of nested, concentric tubes actuated at their base. Thus, they benefit from the same interesting robot outer diameter/open lumen diameter ratio, and can be deployed in a telescopic manner. However, unlike CTRs, the tubes are asymmetrically notched and attached together at one of their tips (see Figure 1). Notching the tubes shifts the line of the centers of mass of their cross-sections, or centroid line. This offset causes a bending moment when one of the tubes is pushed or pulled, causing it to slide relative to the other. Moreover, the same tubes can be rotated at their base to produce a torsion along the robot, which, coupled with the bending, enhances its deformation possibilities. Thanks to this new actuation principle, CPPRs do not undergo elastic instabilities known to occur with concentric tube robots (Peyron et al. 2019). Moreover, CPPRs can bend over large ranges of curvatures and have the potential of replacing tendon actuated continuum robots (TACRs) in many applications, owing to their simplicity of assembly. In addition, similarly to CTRs, CPPRs can be miniaturized to very small scales. Swaney et al. (2016) show that it is possible to manufacture notched tubes of diameters

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Photo of an archetypal two-tube CPPR. From left to right: The outer and inner notched tubes separately. The CPPR with its tubes nested together, at rest, and bent under the effect of a pulling force applied to the tip of the inner tube. Observe how the tubes slide coaxially in relation to each other.</figcaption></figure>

below 0.5 mm. Finally, as the actuation does not rely on tendons but on rods, as in multi-backbone continuum robots (MBCRs) (Dupont et al. 2022), these robots have

**Corresponding author:**

M. Taha Chikhaoui, TIMC Laboratory, Pavillon Taillefer, Rond Point de la Croix de Vie, 38700 Grenoble, France

Email: taha.chikhaoui@univ-grenoble-alpes.fr

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France</footnote>

<footnote n="2">²LS2N Laboratory, Institut Mines Telecom Atlantique, 44307 Nantes, France</footnote>

<footnote n="3">³Université de Strasbourg, CNRS, INSERM, ICube UMR7357, Strasbourg, France</footnote>
bidirectional (push-pull) actuation capabilities and a higher payload than their tendon-driven counterparts (Oliver-Butler et al. 2022). These numerous advantages have led to CPPRs being the topic of multiple preliminary studies for the treatment of various diseases (Rox et al. 2018; Qin et al. 2023). Among these, urology is one of the most promising application field of CPPRs, as they are envisioned to be inserted through the long working channels of flexible endoscopes (Childs and Rucker 2023).

In order to deploy CPPRs and control their shape or tip position, a precise model linking the input variables (i.e. the tube actuation forces or displacements) to the output (i.e. the robot shape) is required. Indeed, modeling CPPRs will ultimately allow to leverage intervention and/or patient-specific design optimization procedures, to provide a solid basis for planning surgical intervention using realistic path planning and environmental interaction. Furthermore, such model can be a key component to advance recent research regarding force sensing, model-based motion control, and stiffness control. Such advances are compulsory to develop safe, reliable, and efficient medical devices.

A constant curvature model of CPPRs was proposed by Oliver-Butler et al. (2017). The model was extended in (Oliver-Butler et al. 2022) by considering variable notch designs using a piecewise constant curvature model. This last model enables approximating planar robot shapes with variable curvature. The model was subsequently used for designing CPPRs able to achieve specific shapes. The piecewise constant curvature assumption can be embedded in mechanical models to account for environment interaction (Rone and Ben-Tzvi 2012; Yuan et al. 2019; Chikhaoui and Rosa 2022). Nevertheless, the model in (Oliver-Butler et al. 2022) exploits the piecewise constant curvature assumption in a pure geometric way, therefore fundamentally limiting their potential to account for environment interactions.

A promising possibility for developing a full mechanical model for CPPRs would be to use Cosserat rod theory (Cosserat and Cosserat 1909). In recent years, Cosserat rod models have been derived for different types of standard continuum robots that emerged, including CTRs, TACRs, MBCRs, and parallel continuum robots (Dupont et al. 2010; Rucker et al. 2010; Rucker and Webster 2011; Till and Rucker 2017; Chen et al. 2021). Exploiting the continuous formulations that are based on this model, static and dynamic simulations have been carried out, first with the shooting algorithm (Rucker and Webster 2011; Till and Rucker 2017), and more recently with a Rayleigh-Ritz reduction method applied to Cosserat strains (Renda et al. 2020; Boyer et al. 2021; Santina and Rus 2020; Li et al. 2023; Tummers et al. 2023). Despite all this success, the Cosserat model cannot be applied to CPPRs in the usual way. This is essentially due to the fact that unlike other systems such as CTRs, the cross-sections of the tubes in a CPPR are neither symmetrical nor uniform. In contrast to the classic theory, the tubes must therefore be modeled as rods whose centroid line differs from their centerline and may even be curved according to variations in the depth and orientation of the notches. In addition, due to the non-uniformity of the tubes, the cross-sectional parameters of the CPPR (including stiffness) depend on the relative position of the tubes as they slide

coaxially relative to each other (the sliding is visible in Figure 1).

These difficulties were addressed for the first time in the recently published article (Childs and Rucker 2023), where Cosserat's model was applied to the case of multi-tube CPPRs whose centroid lines are coplanar with the centerlines. In this approach, the tube cross-sections are defined perpendicular to the centerlines, and Kirchhoff rod theory is used. Based on this model, the coaxiality conditions, which impose the coincidence of the tube centerlines, are "embedded" in a reduced, but not minimal, parameterization of the CPPR. This parameterization consists of the two fields of curvature along the centerline of the CPPR in its Bishop frame, plus the tube curvature fields around the offset vectors of their centroids (which support the radial distance between the centerline and the line of centroids). On the basis of this continuous reduced parameterization of the CPPR, the remaining fields are projected onto a Rayleigh-Ritz basis of polynomials, and the potential energy principle is applied to the whole system. This approach reduces the kinetostatic solution to a finite-dimensional optimization problem, where, due to the non-minimal nature of the continuous parameterization, the distal boundary conditions appear as constraints in a constrained optimization problem. This model, which can take into account plane bending in any direction and the effects of out-of-plane loading, has been successfully tested on a set of CPPRs manufactured to very high standards (laser-notched Nitinol tubes with sandblasted surfaces).

Elaborated in parallel but independently, the present paper proposes an alternative paradigm for kinetostatic modeling of CPPRs. As in (Childs and Rucker 2023), the tubes are initially modeled as Cosserat rods. However, thanks to a relabeling process recently proposed for the dynamics of sliding rods (Boyer et al. 2022), the coaxiality constraints are here used to express all the inner tube configuration variables in terms of those of the outer tube and a single additional field of angle, capable of capturing the combined effects of torsion and bending in three-dimensional (3D) deformations. In this parameterization, which is minimal and purely material (it does not use Bishop's reference frame), the outer tube becomes the leader that guides the kinematics of the inner tubes, considered as followers. Applying the Rayleigh-Ritz method to this reduced kinematics finally provides a set of static equilibrium equations in conventional unconstrained Lagrangian form, which can be solved using a standard root solver.

For pedagogical reasons, this modeling approach is first detailed on an archetypal two-tube planar CPPR, and then extended to the wider context of multi-tube systems capable of 3D deformation. This includes the case of tubes whose centroid and centerlines are not in the same plane (for example with helical centroid lines that wrap around the centerline of the CPPR), as well as the case of CPPRs subjected to the combined effects of push/pull forces and torques applied to the base of their follower tubes. The proposed model is validated in 2D through an extensive set of experiments using various 3D-printed 2-tube CPPR designs, with and without external forces applied at the robot tip. The number of experimental conditions for experiment to simulation comparison reaches
28, 44, and 20 for, respectively, various regularly noched CPPRs, various CPPRs with variable curvature, and CPPR subject to external loading. The experiments show that the proposed model is capable of capturing the behavior of a CPPR subject to external loads, while the tip errors stay below 6.33% of the robot length in the unloaded cases, and 6.57% for the loaded cases. In the 3D context, the capabilities of the approach are illustrated by a comparison with the 3D experimental results of (Childs and Rucker 2023), i.e. for a 2-tube CPPR subjected to a load out of the plane of bending, and a 3-tube system capable of bending in any plane containing its central axis. Beyond these restricted conditions, first simulations and experiments illustrating more complex situations combining push/pull forces and torques applied to the tubes, are reported. To ease prototyping of CPPR applications, the associated code can be found at https://github.com/TIMClab-CAMI/Cosserat-Rod-Modeling-of-Continuum-Concentric-PushPull-Robots/.

The remainder of this paper is organized as follows. Section 2 introduces the modeling approach on the case of 2-tube CPPRs in planar deformation. In Section 3, this modeling approach is extended to the case of multi-tube CPPRs able to deform in the three dimensions. In Section 4, the model is first validated against experiments conducted under 2D conditions, and then applied to simulations and experiments of 3D scenarios that go beyond the state of the art in this field. The obtained results are discussed and put in perspective in Section 5 that concludes the paper. For the sake of research reproducibility and to ease prototyping of CPPR applications, the associated code will be released upon the publication.

## 2 Kinetostatic Modeling of an Archetypical 2-Tube Planar CPPR

In this section, the kinetostatic modeling approach of the paper is introduced through its application to the archetypal case of a planar 2-tube CPPR, as shown in Figure 2. This example will allow to present in a pedagogical way all the key ingredients of the approach that go from design to the static simulation of these new continuum robots, and to explain how the Cosserat rod model can be applied to them. Although illustrated here in a restricted context, all the concepts proposed in this section are formulated within the general geometrical framework of SE(3), which will be exploited in the next section to extend the approach to the wider context of multi-tube CPPRs capable of deforming in three dimensions.

### 2.1 Working Principle

As illustrated in Figure 2, a 2-tube CPPR is composed of an inner tube nested within an outer tube. The two tubes are coaxial, i.e. they share a common centerline, and are attached together at their distal end. Each tube is notched along its length with asymmetric notches that offset its centroid line with respect to its centerline. By varying the geometry of the notches along tubes, it is possible to control these *centroid offsets*. By analogy with the well-known TACRs, the centroid and centerlines here play the role of the tendons and the central backbone respectively, the rigidity of the backbone being transferred to the two nested tubes. In contrast to

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Composition, working principle and frames of a CPPR. Outer (a) and inner (b) tubes at rest with their centroid and centerlines. The assembled tubes, forming the CPPR in a resting configuration, with the cross-sectional frames attached to the tubes (c). Actuating the CPPR with $\Delta a > 0$, bends it in a deformed configuration and makes the frames move (d).</figcaption></figure>

TACRs, actuation stresses are transmitted through tubes (not tendons), that can be both pushed and pulled (see also Extension 1). Without loss of generality, the outer tube is clamped in the baseplate, while the inner tube can be pushed or pulled axially through the same plate of a controlled length $\Delta a$ (see Figure 2(d)). This axial movement lengthens or shortens the part of the inner tube contained in the outer tube, which bends the CPPR due to the differences in length and position of the tube centroid lines. Remarkably, all the features of more standard continuum robots are here integrated and distributed in a very compact design (Oliver-Butler et al. 2017).

### 2.2 Design Parameters

Throughout the article, *i* denotes the index of tubes, with *i* = 1 for the outer tube and *i* = 2 for the inner tube. Quantities that are related to tube *i* are labeled with the subscript *i*. The geometry of our 2-tube CPPR is defined in the resting (straight) configuration of Figure 2(c), in a frame $\mathcal{F}_s = (O, e_X, e_Y, e_Z)$ located on the baseplate, with $e_X$ supporting the centerlines of tubes 1 and 2, denoted $C_1$ and $C_2$ respectively. This coincidence of tube centerlines $C_1 = C_2$ is called the *coaxiality constraint*. It is required by the operating principle of CPPRs and will play a key role in its kinetostatic model when it deforms. Considered separately, each tube consists of alternating closed rings, or *teeth*, and cut-out portions, called *notches* (see Figure 2(a) and (b)). The notches¹ are cut in an initially axi-symmetrical tube according to a longitudinal and radial geometry described in Figure 3. The longitudinal geometry of a tube of index *i*, is defined by the length (along $C_i$) of its teeth and notches, which are constant along the tube and denoted by $k_i$ and $c_i$ respectively (see Figure 3(left)). Designed for planar deformations, the radial geometry is here only defined by three independent design parameters $\rho_i$, $\bar{\rho}_i$, $\alpha_i$ (see Figure 3(right)). $\rho_i$ and $\bar{\rho}_i$ are respectively the inner and outer radii of tube *i*, while $\alpha_i$ governs the depth $\gamma_i = \bar{\rho}_i + \alpha_i$ of its notches. In the present section, only the depth
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Longitudinal and radial design parameters of notches of tubes 1 and 2, displayed on a partial front (left) and sectional view (right) of a CPPR. Remind that $\gamma_i$ and $\alpha_i$ are related by $\gamma_i = \bar{\rho}_i + \alpha_i$.</figcaption></figure>

parameter $\alpha_i$ is adjusted from one notch to the next, in order to design the stiffness and centroid offset of the tube along its length. In Section 3, we will add to $(\underline{\rho}_i, \bar{\rho}_i, \alpha_i)$ an angle of rotation of the notches in the plane of their cross-sections. This additional design parameter will allow to describe, e.g., tubes whose centroid line wraps helically around their centerline. For the time being, we will keep to the above restricted design, and in the next subsection present a model which approximates the two tubes as Cosserat rods whose notch depth varies continuously along their length.

## 2.3 Modeling the Tubes as Equivalent Rods

Firstly, it is observed that when loaded separately, each of the tubes of a CPPR deforms along its length, with smooth shapes typical of those of a beam, i.e. dominated by bending and torsion. Therefore, they are approximated as equivalent Cosserat rods whose cross-sectional geometry is asymmetric and varies continuously along their length. Such approximation is justified by the following functional analysis and is proved to be valid thanks to the experimental results of Section 4. Unlike other continuum robot designs (e.g. TACR, CTR), the asymmetry and non-uniformity of cross-sections is necessary to capture the effects of notches and their possible depth variations along the tubes. Secondly, applying a functional analysis to the design of the tubes shows that the function of the teeth is to ensure coaxial guidance of one tube into the other, without compromising the function of the notches, which is to control centroid offsets and the rigidity of rod-like deformations. When this is the case (such as for the prototypes in this article), the longitudinal stress forces transmitted between two notches by a tooth are channeled along its lower part, i.e. they do not penetrate deeply into the tooth 3D annular domain, which can be considered rigid. Based on this observation, which was confirmed by a preliminary finite element analysis on the Ansys software (see also FEM simulations in (Childs and Rucker 2023)), the deformation along each tooth can be approximated by that of a rod segment whose cross-section has the shape of those of the adjacent notches, and whose dimensions are an interpolation of them. To implement this idea, let us consider the tubes assembled in the resting configuration of a CPPR, with abscissa X along the CPPR axis supported by $e_X$ (see Figure 2(c)). The cross-sections of the equivalent rod modeling tube i have the U-shape of its notches (see Figure 4), with constant internal and external

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Rod with continuously variable sections (dark red), equivalent to a tube with variable notches (light and dark red). Top: in red: longitudinal section (or profile) of the tube; in blue: interpolation of the tube depth profile with equation (1). Bottom: for any value of the continuous label X, the blue function can be evaluated, yielding the continuous U-shaped cross-sections along the equivalent rod. See also Extension 2.</figcaption></figure>

radii $\underline{\rho}_i$, $\bar{\rho}_i$, and notch depths, defined for all X, by the cubic spline interpolation² of the piece-wise constant profile $(\vec{X}_i, \vec{\alpha}_i)$ of $\alpha_i$

$$ \alpha_i(X) = S(X, (\vec{X}_i, \vec{\alpha}_i)), \quad (1) $$

where $\vec{X}_i$ is the vector of the ordered values of X at the middle of the notches of the tube i, while $\vec{\alpha}_i = \alpha_i(\vec{X}_i)$ gathers the corresponding discrete values of $\alpha_i$.

## 2.4 Cosserat Rod Model of the Tubes at Rest

Through the above approximation process, the two tubes are now considered as Cosserat rods whose asymmetric cross-sections vary continuously along their length. According to this model, a frame is attached to each of the rigid cross-sections of the rod. In the following subsections, these cross-sectional frames are defined (drawn) on each tube i considered separately and in its resting configuration. In Section 2.5, the tubes will be assembled and the same frames will move as the CPPRs deform.

### 2.4.1 Specification of the Cross-Sectional Frames
Recalling the foundations of the Cosserat rod theory, a rod is considered to be a continuous stack of rigid material cross-sections along a material reference line drawn along the largest dimension of its resting (unstressed) configuration (Antman 2005). The definition of the reference line and cross-sections is a choice left to the modeler. In the case of a CPPR, it is convenient to take as reference lines of the rods, their centroid lines or neutral bending lines, noted $G_i$, while their cross-sections are chosen to be orthogonal to the centerlines of the tubes $C_i$ (see Extension 2). From now on, we denote by $X_i$ the arc-length along $G_i$ of total length $l_i$ oriented positively from the base to the tip of the CPPR for tube 1, and from the tip (distal end) of the CPPR to the base for tube 2, as shown in Figure 2(d).

**Remark 1:** It is important to note that in contrast to common practice in continuum robotics, the reference (centroid) lines $G_i$ do not coincide with the centerlines $C_i$ of the tubes. In particular, the former may not be parallel to the latter, or may even be curved, while the latter are always straight in a CPPR at rest. As a consequence, the cross-sections of a rod are in general not orthogonal to its reference line when
<figure><img src="image_6.png" /><figcaption>Figure 5. (a) : Lateral view of a piece of a 2-tube planar CPPR at rest, with non-parallel curved centroid lines. Referring to Section 2.5, $\beta_i$ is a design shearing parameter measuring the angle between $b_{iX}$ and $t_i = (dr_i/dX_i)$ (see also Extension 2). (b) : Geometry of the same CPPR in its resting configuration with curved centroid lines, in the baseplate frame $(O, e_X, e_Y, e_Z)$. $X_1$ and $X_2$ are the arc lengths of the continuous green and red lines, counted from the base and tip of the CPPR respectively, with their "straight" approximations on both sides as in Remark 4.</figcaption></figure>

the rod is at rest (see Figure 5). This is in contrast with the usual theory. In the next section, the Cosserat rod model will naturally accommodate for this choice, by declaring a shear field (consequence of design) along each tube. In Section 3, similar reasons will drive us to declare another design torsion field capturing how $G_i$ wraps around $C_i$ in the rest configuration. •

We are now in a position to define the field of cross-sectional frames that will be used in the model of the tubes. To each cross-section labeled by $X_i$, or $X_i$-cross-section, we rigidly attach a frame $\mathcal{F}_i(X_i) = (G_i, b_{iX}, b_{iY}, b_{iZ})(X_i)$. The origin of this frame $G_i(X_i) \in \mathcal{G}_i$ is the centroid of the $X_i$-cross-section, while $b_{iX}$ is its unit normal vector, oriented towards increasing $X_i$, and $b_{iY}$, $b_{iZ}$ are such that $e_Z = b_{iZ} = b_{iX} \times b_{iY}$, as indicated in Figure 2(c).

### 2.4.2 Cross-Sectional Inertia Parameters of the Tubes
In the cross-sectional frame $\mathcal{F}_i(X_i)$, the 6 × 6 matrix of the inertia tensor of the $X_i$-cross-section of tube $i$, reads

$$ \mathcal{I}_i(X_i) = \begin{pmatrix} I_i & 0_{3 \times 3} \\ 0_{3 \times 3} & A_i 1_{3 \times 3} \end{pmatrix}, \quad (2) $$

where $A_i(X_i) \in \mathbb{R}$ is the area of the cross-section, while the non-diagonal components are zero, since they are proportional to the position vector of the section centroid $G_i(X_i)$ in $\mathcal{F}_i(X_i)$, i.e. in a frame located on $G_i(X_i)$. Throughout this manuscript, the notations $0_{n \times m}$ and $1_{n \times m}$ stand for $n \times m$ zero and identity matrices, respectively. As regards the angular matrix component of (2), since the $\mathcal{F}_i$-frame respects the symmetries of the U-shaped $X_i$-cross-section, we have $I_i(X_i) = \text{diag}(I_{iX}, I_{iY}, I_{iZ})$, which from Huygens transport theorem are given by

$$ I_{iX} = J_{iX} - A_i D_i^2, \quad I_{iY} = J_{iY}, \quad I_{iZ} = J_{iZ} - A_i D_i^2. \quad (3) $$

where $J_{iX} = J_{iY} + J_{iZ}$, $J_{iY}$ and $J_{iZ}$, are the usual angular inertia momenta around the $(C_i, b_{iX})$, $(C_i, b_{iY})$ and $(C_i, b_{iZ})$-axes, respectively, with $C_i(X_i)$ the center point of

<figure><img src="image_7.png" /><figcaption>Figure 6. Varying the depth of the notches along the tubes varies the position of their offset centroid line. (a,b) Tubes with convergent (divergent) centroid lines. The notch depth linearly decreases (increases) towards the tip of the CPPR. (c) Cross-sectional frames and centroid offsets of the inner (top) and outer (bottom) tubes.</figcaption></figure>

the $X_i$-cross-section. In (3), $D_i(X_i)$ represents the only non-zero component of the centroid offset vector

$$ \overrightarrow{C_i G_i}(X_i) = D_i(X_i) b_{iY}(X_i), \quad (4) $$

which we will call the centroid offset function of tube $i$ at $X_i$ (see Figure 6(c)). Finally, the detailed expressions of the four independent parameters ($A_i, D_i, J_{iY}, J_{iZ}$) are given in Appendix B for rectangular notches. As expected, they depend on the cross-sectional parameters $(\rho_i, \bar{\rho}_i, \alpha_i)$ of the equivalent rods of Section 2.3, and are functions of $X_i$, through the only design parameter which is considered to vary along the length of our 2-tube CPPR, namely $\alpha_i$.

### 2.4.3 Cross-Sectional Stiffness Parameters of the Tubes
The cross-sectional matrices $\mathcal{I}_i$ are fully diagonal, and assuming that both tubes are made of a same homogeneous material, their stiffness can be modeled by cross-sectional 6 × 6 Hooke matrices in usual diagonal form

$$ \mathcal{H}_i(X_i) = \text{diag}(GI_{iX}, EI_{iY}, EI_{iZ}, EA_i, GA_i, GA_i), \quad (5) $$

with $E$ being Young's modulus, $G = E(1 + \nu)/2$ the shear modulus, and $\nu$ Poisson's ratio of the material.

**Remark 2:** Although (3) to (5) are valid for CPPRs moving in 3D, in the rest of this section, a kinetostatic model of a 2-tube CPPR able to deform in the $(O, e_X, e_Y)$-plane will be proposed. As a consequence of this restriction, this model will only depend on the centroid offset functions $D_1$ and $D_2$ of the two tubes 1 and 2, their cross-section area functions $A_1, A_2$ (gravity is a priori considered), and their second inertia moments $I_{1Z}$ and $I_{2Z}$. All these parameters defining functions of $X_1$ and $X_2$ respectively •

### 2.4.4 Illustration on a Specific Design
In Figure 6(a) and Figure 6(b), we schematized the profile of $D_i$ along the centerline of a 2-tube CPPR as the depth of the notches increases (a), and decreases (b) linearly along its length from the base to its tip. Referring to the position of the centroid lines, in the first case, the design of notches is said to be convergent, while it is divergent in the second. In the same
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Varying inertia parameters ($A_i, D_i, I_{iZ}$)($X_i$), $i = 1, 2$ along a 2-tube CPPR (design $R_3$ of Table 1), with convergent tube centroid lines.</figcaption></figure>

way, the design is said to be regularly notched when the depth of notches is uniform. To illustrate the developments of this section, Figure 7 displays the profile of ($A_i, D_i, I_{iZ}$)($X_i$) of the two tubes for the convergent design ($R_3$ of Table 1). These parameter functions are obtained with (3) and the expressions of ($A_i, D_i, J_{iZ}$)($X_i$) in Appendix B. They are all based on the interpolation of (1).

**Remark 3:** It should be noted that, despite its conceptual simplicity, this calculation process is subject to unavoidable modeling approximations, and could be improved by more sophisticated models of asymmetric rods (Young et al. 2012), or replaced by a direct identification of the tube stiffness coefficients. This is the choice adopted by Childs and Rucker (2023), where the bending and torsion stiffness, as well as the centroid offsets, are directly deduced from a preliminary finite element analysis of the isolated tubes •

## 2.5 Tube Deformation Model

Now that the Cosserat rod model of tubes at rest has been defined in Sections 2.3 and 2.4, let us consider deformations. Using $\mathcal{F}_s = (O, e_X, e_Y, e_Z)$ of Figure 2(d) as an inertial reference frame, for all $X_i \in [0, l_i]$, the (inertial) pose of $\mathcal{F}_i(X_i)$ with respect to $\mathcal{F}_s$ is parameterized by the homogeneous transformation $g_i(X_i) \in SE(3)$

$$g_i(X_i) = \begin{pmatrix} R_i(X_i) & r_i(X_i) \\ 0_{1 \times 3} & 1 \end{pmatrix}, \qquad (6)$$

where $R_i(X_i) \in SO(3)$ is the orientation matrix of the basis of $\mathcal{F}_i(X_i)$ in $\mathcal{F}_s$ and $r_i(X_i) \in \mathbb{R}^3$ is the position vector of its origin in the same frame. The deformation of each tube is defined by the field of (space-rate) twists in $se(3) \cong \mathbb{R}^6$

$$\xi_i(X_i) = \left( g_i^{-1} \frac{\partial g_i}{\partial X_i} \right)^{\vee}, \qquad (7)$$

where $\xi_i = (K_i^T, \Gamma_i^T)^T$, and $K_i = (K_{iX}, K_{iY}, K_{iZ})^T$ and $\Gamma_i = (\Gamma_{iX}, \Gamma_{iY}, \Gamma_{iZ})^T$ are respectively the vector of twist ($K_{iX}$) and curvatures ($K_{iY}, K_{iZ}$), and that of stretch ($\Gamma_{iX}$) and transverse shears ($\Gamma_{iY}, \Gamma_{iZ}$). The Cosserat rod model can capture several rod sub-models by allowing only a sub-vector of $\xi_i$ to be free (Boyer et al. 2021). In the present case, the two tubes are considered as Kirchhoff rods, straight at rest, and moving in the $(O, e_X, e_Y)$ plane of

Figure 2(d). Therefore, the free (strain) fields are the two curvatures $K_{1Z}(X_1)$ and $K_{2Z}(X_2)$, along $e_Z = b_{iZ}$ (normal to the motion plane), while all other components of $\xi_1$ and $\xi_2$ are some constants fixed as follows. First, since the deformation holds in a plane normal to $e_Z$, $\Gamma_{iZ} = K_{iX} = K_{iY} = 0$. Second, as Kirchhoff rods are inextensible, $\Gamma_{iX} = 1$. Finally, with the definition of the material cross-sections and reference lines of Section 2.4.1, $\Gamma_{iY}(X_i)$ fixes the time-independent small tilt angle between the unit tangent to $G_i$ at $X_i$, and the unit normal to the $X_i$ cross-section, respectively defined by $(dr_i/dX_i)(X_i)$ and $b_{iX}(X_i)$. Referring to the Cosserat rod model, $\Gamma_{iY}$ can be interpreted as offset shear fields capturing the non-uniformity of the asymmetric cross-sections along the tubes. Note here that in accordance with Kirchhoff’s theory, these fields are not degrees of freedom of the CPPR (and do not store internal energy), but are instead fixed by design through the relations

$$\beta_i(X_i) = \frac{dD_i}{dX_i}(X_i), \qquad (8)$$

which are deduced from simple geometric considerations applied to Figure 5(a), and define the fields of shearing angles $\beta_i$. Finally, under these modeling assumptions the twists of (7), take the detailed form

$$\xi_i(X_i) = (0, 0, K_{iZ}, 1, \beta_i, 0)^T(X_i), \qquad (9)$$

where the shearing (tilt) angles $\beta_i$ are fixed by design, while the curvature fields $K_{iZ}$ are related to the deformation of the tubes of a CPPR.

## 2.6 Kinematics of a 2-Tube Planar CPPR

So far the two tubes have been considered either separately or assembled in the resting configuration of the CPPR of Figure 2(c). In the assembled CPPR design, both tubes can either be pulled or pushed relative to each other in an agonist-antagonistic way. Without loss of generality, in the remainder of this article, the CPPR is actuated by pulling or pushing the tip of tube 2 through the baseplate. Referring to Figure 2(c) and (d) the oriented length of the piece of tube 2 below the baseplate is $a$, with $a^o$ its nominal value when the CPPR is in its resting configuration, and $\Delta a = a - a^o$ its actuated variation.

### 2.6.1 Geometric Assumptions of the Model
A CPPR is a slender (composite) body. Thus, by noting $d = 2\bar{\rho}_1$ its diameter and $l = l_1$ its length (quantities related to the outer tube), the aspect ratio $|d/l| \simeq \varepsilon$ defines a small positive number characteristic of all small quantities in the CPPR model. According to the geometrically-exact model of Cosserat rods (Simo and Vu-Quoc 1988), the tubes are assumed to endure finite deformations and small strains. This means first of all that $|K_{iZ}d|$ is a small quantity compared to the unit. Second, since the tilt $\beta_i(X_i)$ of the centroid lines along the CPPR centerline is captured by the shear strain of the model, these tilts must be small as well. Third, the length variation of the portion of tube 2 below the baseplate, noted $\Delta a$ in Figure 2(d), characterizes the retraction or elongation of tube 2 above the baseplate. This is also a small quantity compared to $l$ (typically the displacement actuation $|\Delta a|$ is of the order of the diameter $d$ of the CPPR). All these
assumptions that are compatible with the operating principle of a CPPR, can be summarized as follows:

$$|K_{iZ}d| \approx \varepsilon, \quad \left|\frac{dD_i}{dX_i}\right| = |\beta_i| \approx \varepsilon, \quad \left|\frac{\Delta a}{l}\right| \approx \varepsilon. \quad (10)$$

In the rest of this section, we derive a model for a 2-tube CPPR consistent at the first order with respect to $\varepsilon$, i.e. a model which preserves the geometrically-exact character of the Cosserat rod theory. In Section 3, we will extend this model to multi-tube CPPRs capable of twisting and bending. In this broader context, the approximations (10) will be supplemented by others consistent with small torsional strains. Finally, as in (10), any first order approximation will be indicated by the symbol $\simeq$.

**Remark 4:** As a first illustration of such an approximation, let us consider our 2-tube CPPR in its resting (straight) configuration of Figure 2(c). Using the coordinate $X$ along $e_X$, geometrical considerations emphasized in Figure 5(b), show that we have $dX = \pm \cos(\beta_i)dX_i \simeq \pm(1 - \beta_i^2/2)dX_i$ (+ if $i = 1$, - if $i = 2$), which gives after integration, and to the first order in $\varepsilon$

$$X \approx X_1 \approx l_1 - X_2, \quad \forall (X_1, X_2) \in [0, l_1] \times [0, l_2 - a^o], \quad (11)$$

and $l_1 \simeq l_2 - a^o$. As will be explained in the next section, the correspondence (11) between $X_1$ and $X_2$ no longer holds when the CPPR deforms. Indeed, as soon as the CPPR bends, the two tubes slide relative to each other, and a non-negligible (first order) deformation-dependent shift between the initially coincident cross-sections $X_1$ and $X_2 = l_1 - X_1$ appears (see Figure 8(a) and (b)). To model these effects we need to introduce new concepts detailed in the next section •

**2.6.2 Face-to-Face Function and Relabeling Process** When tube 2 is pulled ($\Delta a > 0$) or pushed ($\Delta a < 0$) across the baseplate (see Figure 2(d) and Figure 8(a) and (b)), it slides into tube 1. To capture the effects of these sliding motions in the CPPR model, it is convenient to define a function $\tilde{X}_2 : X_1 \in [0, l_1] \to [0, l_1 - a]$, named face-to-face function. This function gives in any (deformed) configuration of a 2-tube CPPR, the label (arc length) $X_2$ along tube 2, noted $\tilde{X}_2(X_1)$, of the point where the $X_1$-cross-section plane intersects the centroid line $G_2$ of tube 2. In other words, it permits to look at the cross-section of tube 2 that faces a given section of tube 1. This face-to-face function allows any function $f_2$ of $X_2 \in [0, l_2 - a]$, to be changed into a function $\tilde{f}_2 = f_2 \circ \tilde{X}_2$ of $X_1 \in [0, l_1]$ by the pullback operation

$$\tilde{f}_2(X_1) = f_2(\tilde{X}_2(X_1)), \quad (12)$$

which stands for a relabeling of a field, here $f_2$, defined along tube 2, with the $X_1$-label of tube 1. Since the Cosserat kinematics handle first order derivatives with respect to rod labels (see (7) and (8)), we will also use the following key relationship, tangent to (12)

$$\frac{d\tilde{f}_2}{dX_1} = \left(\frac{d\tilde{X}_2}{dX_1}\right) \left(\frac{df_2}{dX_2}\right)^{\sim} = h(X_1) \left(\frac{df_2}{dX_2}\right)^{\sim}, \quad (13)$$

which defines $h$, and allows transforming the $X_2$-derivative of any field $f_2$, originally defined as a function of $X_2 \in [0, l_1 - a]$, into the $X_1$-derivative of a function of $X_1$ on $[0, l_1]$.

**Remark 5:** Note that the pullback process of (12) and (13) rationalizes the change of variable on functions: $f_2(X_2(X_1))$, and their differentials: $(df_2/dX_1) = (dX_2/dX_1)(df_2/dX_2)$. In the rest of the article, the exhaustive notations of (12) and (13) will be systematically used to avoid any confusion. In short while $X_1$ and $X_2$ are variables (labels along the centroid lines of tubes), $\tilde{X}_2$ and any tilde-quantity, originally related to tube 2 (and indexed 2), denotes a function of the label $X_1$ along tube 1. In Section 2.8, the detailed expression of this pullback operation will be derived, in order to relabel the Cosserat rod model of tube 2, along tube 1 •

**2.6.3 Kinematic Sliding Constraint** When tube 2 is pulled or pushed, it slides into tube 1 as well as into the tubular guide below the baseplate. These contacts between the two tubes and between tube 2 and the guide, are all considered rigid and friction-less, and can therefore be fully modeled by the following kinematic constraints. First, the coaxiality constraint, imposes the coaxiality of the tube centerlines $C_1 = C_2$. Second, the guiding constraint guides the tube 2 inside the baseplate. Extension 3 presents a visual representation of these contact constraints. For practical calculation purposes, $(E_X, E_Y, E_Z)$ will denote the canonical (numerical) base of $\mathbb{R}^3$, i.e. $E_X = (1, 0, 0)^T$, $E_Y = (0, 1, 0)^T$, and $E_Z = (0, 0, 1)^T$.

**2.6.4 Coaxiality Constraint Between Tubes** Coaxiality can be split in two sub-conditions, which are the coplanarity and the concentricity of all cross-sections that face each other along the tubes. Coplanarity means that any $X_1$-cross-section of tube 1 shares a common normal straight line with the $\tilde{X}_2(X_1)$-cross-section of tube 2. Equivalently, introducing the relative transformation $g_{1,2}$ between these two cross-sections, such that

$$\forall X_1 \in [0, l_1] : g_2(\tilde{X}_2(X_1)) = g_1(X_1)g_{1,2}(X_1), \quad (14)$$

coplanarity means that $g_{1,2}$ must leave any line supported by $b_{1X}(X_1)$ unchanged. This condition requires the rotational part of $g_{1,2}$ to be of the form $R_{1,2} = \exp(\theta \hat{E}_Z) \triangleq R_Z(\theta)$, with $\theta = 0$ or $\pi$. With our choice of cross-sectional frames (see Figure 8(c)) and because of the opposite orientation of the centroid lines $G_i$, $R_{1,2}$ must flip $b_{1X}$, and we have $R_{1,2} = R_Z(\pi)$. With such a rotation, $g_{1,2}$ changes a $X_1$-cross-section on tube 1 into a $\tilde{X}_2(X_1)$-cross-section on tube 2 lying in the same plane. Now these two cross-sections need to share a common center, i.e. $g_{1,2}$ must also fulfill the concentricity sub-condition

$$\forall X_1 \in [0, l_1] : \overrightarrow{C_1\tilde{C}_2} = \overrightarrow{C_1G_1} + \overrightarrow{G_1\tilde{G}_2} + \overrightarrow{\tilde{G}_2C_2} = \vec{0}, \quad (15)$$

with $C_1(X_1) \in C_1$ the center of the $X_1$-cross-section of tube 1, and $\tilde{C}_2(X_1) = C_2((\tilde{X}_2(X_1)) \in C_2$ the center of the cross-section of tube 2 facing the $X_1$-cross-section of tube 1 (see Figure 8(c)). Expressing (15) in the $X_1$-cross-sectional frame, with the definition of centroid offsets $\overrightarrow{C_iG_i} = D_i b_{iY}$ of Section 2.4.2, leads to the positional component $r_{1,2}$, of $g_{1,2}$

$$r_{1,2}(X_1) = -D(X_1)E_Y, \quad (16)$$
where we have introduced the offset between the centroid lines of the two tubes (see Figure 8(c))

$$D(X_1) = D_1(X_1) + \tilde{D}_2(X_1), \quad (17)$$

and remind that $\tilde{D}_2(X_1) = D_2(\tilde{X}_2(X_1))$. As a result, in (14), $g_{1,2}$ must take the detailed form

$$g_{1,2} = \begin{pmatrix} R_Z(\pi) & -DE_Y \\ 0_{1 \times 3} & 1 \end{pmatrix}. \quad (18)$$

Finally, the coaxiality of tubes can be entirely modeled by the pointwise constraint (14) with $g_{1,2}$ given by (18). In Section 3, the form (14) will be used to formulate constraints capable of capturing coaxiality of a multi-tube CPPR subject to bending and twisting in 3D.

### 2.6.5 Guiding Constraint of Tube 2 Below the Baseplate

Similarly, the effect of the rigid tubular guide (see Figure 2(d) and Extension 4) on the part of tube 2 below the baseplate can be modeled by another set of point-wise geometric kinematic constraints

$$E_Y^T r_2(X_2) = E_Z^T r_2(X_2) = 0, \quad R_2(X_2) = 1_{3 \times 3}. \quad (19)$$

which hold $\forall X_2 \in [l_2 - a, l_2]$, and prevent the tube to rotate or translate laterally inside the guide.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a) and (b): Cross-section label shift due to the sliding between deformed tubes for a regularly notched CPPR. (c): Two coincident cross-sections of tube 1 and 2. Remind that for any quantity $f_2$ related to tube 2, we have $\tilde{f}_2(X_1) = f_2(\tilde{X}_2(X_1))$ with $\tilde{X}_2(X_1)$ the label of the cross-section of tube 2 that coincides with the $X_1$-cross-section of tube 1. The distance $D(X_1)$ is the offset between centroid lines at $X_1$. See also Extension 4.</figcaption></figure>

### 2.6.6 Computation of the Face-to-Face Function Applying definitions (12) and (13) to the coaxiality constraint (14), provides after calculations detailed in Appendix C, the face-to-face function as the solution of the backward integration (i.e. from $X_1 = l_1$ to $X_1 = 0$), of the ordinary differential equation (ODE) initialized by $\tilde{X}_2(l_1) = 0$

$$\frac{d\tilde{X}_2}{dX_1}(X_1) = h(X_1) = -(1 + K_{1Z}D), \quad (20)$$

where we recognize the offset $D$ between the centroid lines of the two tubes of (17). Finally, note that the term $K_{1Z}D$ in (20) models the effect of the axial shift of the cross-section labels (or "sliding"), due to deformation (see Extension 4). Referring to Section 2.6.1, this term is of the order of $\varepsilon$, and cannot be ignored in a geometrically-exact model of the CPPR.

## 2.7 Kinetostatic Model of a 2-Tube Planar CPPR

This section presents a closed formulation that models the static balance between internal elastic forces and external forces exerted on each of the tubes considered in isolation. These external forces are induced by contacts of tube 2 with tube 1, and of tube 2 with the guide in the baseplate, by the gravity field; the actuation push-pull force $-T_+e_Z$ applied at the tip $X_2 = l_2$ of tube 2 (that generates $\Delta a$); and a possible external force $f_+$, applied at the tip of the CPPR. To introduce this closed formulation, we must define the model of stress. Remind that in the Cosserat rod theory, stresses are in the dual space of $X$-rates $\xi$. Therefore, the stress state of tube $i$ is modeled by a field of wrench $\Lambda_i \in se(3)^* \cong \mathbb{R}^6$, along the tube length $[0, l_i]$. In detail, $\Lambda_i(X_i) = (M_i^T, N_i^T)^T(X_i)$ is the internal wrench exerted by the piece of rod $[X_i, l_i]$ onto the piece $[0, X_i]$ across the $X_i$-cross-section, with $M_i = (M_{iX}, M_{iY}, M_{iZ})^T$, and $N_i = (N_{iX}, N_{iY}, N_{iZ})^T$, named the *moment-stress* and *force-stress* respectively, all being expressed in the $X_i$-cross-sectional frame. Note that in our case, the CPPR moving in the plane $(O, e_X, e_Y)$, we have $M_{iX} = M_{iY} = N_{iZ} = 0$. Moreover, since the rods are Kirchhoff, $M_{iZ}$ captures a field of elastic moments, while $N_{iY}$ and $N_{iZ}$ model some internal reaction forces preventing the shearing and stretching along tubes. Based on this model of stress, in the (quasi)static regime, the deformations of tubes 1 and 2 are governed by the set of the Cosserat rod ODEs (Tummers et al. 2023), which hold for $i = 1, 2$

$$\frac{dg_i}{dX_i} = g_i \hat{\xi}_i, \quad \frac{d\Lambda_i}{dX_i} = ad_{\xi_i}^T \Lambda_i - \bar{F}_i. \quad (21)$$

From left to right, the first ODEs of (21) stand for a continuous geometric model of the tubes, while the second represent their static wrench-stress balances, with $\bar{F}_i \in se(3)^* \cong \mathbb{R}^6$, the density of external wrench exerted along tube $i$ (including contact forces), expressed in the cross-sectional frames and detailed hereafter. These ODEs are supplemented with boundary and connectivity conditions. By denoting $g_{1,-}$ the inertial pose of the first cross-section of tube 1, one has the boundary condition (BC) at $X_1 = 0$

$$g_1(0) = g_{1,-}. \quad (22)$$

For tube 2, the BC at $X_2 = l_2$ is

$$\Lambda_2(l_2) = F_{+,2} = (0, 0, 0, T_+, 0, 0)^T, \quad (23)$$

where $T_+$ represents the actuation push-pull force applied to the tip of tube 2 (i.e. $X_2 = l_2$). At $X_1 = l_1$ and $X_2 = 0$ the tubes are clamped in a same rigid cross-section at the tip of the CPPR, which imposes the connectivity condition (CC) on poses

$$g_2(0) = g_1(l_1)g_{1,2}(l_1), \quad (24)$$

where $g_{1,2}(l_1)$ is given by (18) in which $D(l) = D_1(l) + D_2(0)$. The connection of tubes at the tip of the CPPR, also imposes the condition on wrenches

$$\Lambda_1(l_1) = Ad_{g_{1,2}(l_1)}^{-T} \Lambda_2(0) + Ad_{g_{1,s}}^{-T} F_+, \quad (25)$$

where $F_+ = (0_{1 \times 3}, f_+^T)^T$ is the wrench of a possible external concentrated force $f_+$, expressed in the inertial
frame and imposed at the tip of the CPPR, while $g_{1,s}^{-1}$ is the pose of the tip cross-sectional frame of tube 1 in a frame defined by the inertial basis located at the application point of the force. For the CPPRs here studied, the external wrench density $\bar{F}_i$ of (21) is detailed as

$$ \bar{F}_i = \bar{F}_{\text{ext},i} + \bar{F}_{c,i}. \qquad (26) $$

Where $\bar{F}_{\text{ext},i}$ models the effect of all distributed external forces applied to the CPPRs except those due to contacts. In this paper, the distributed external forces are restricted to gravitational forces which take the detailed form

$$ \bar{F}_{\text{ext},i}(X_i) = (0_{1 \times 3}, (\mu_i A_i)(X_i)(R_i^T(X_i)a_g)^T)^T, \quad (27) $$

with $\mu_i A_i$, the mass density per unit of length of tube $i$, and $a_g$, the uniform gravitational acceleration field expressed in the inertial frame $\mathcal{F}_s$. The second component $\bar{F}_{c,i}$ of (26) stands for a field of contact wrenches whose components play the role of Lagrange multipliers in charge of forcing the coaxiality and guiding constraints (14) and (19)$^3$. Finally, one needs to consider the constitutive relation along the tubes, which consistently with the Cosserat small strain assumption, takes the linear form

$$ \Lambda_i = \mathcal{H}_i(\xi_i - \xi_i^o), \qquad (28) $$

where $\xi_i^o$ is the rest value of $\xi_i$. Then, using (9) in (28) with (5), we have for a 2-tube planar CPPR

$$ M_{iZ} = EI_{iZ}K_{iZ}, \qquad (29) $$

where $EI_{iZ}$ is the flexural stiffness, i.e. the only component of Hooke's tensor (5) that holds for a Kirchhoff rod moving in the $(O, e_X, e_Y)$ plane. Finally, the equations (21) to (27) and (29) with the constraints (14) and (19) define a closed kinetostatic formulation for a 2-tube planar CPPR.

**Remark 6:** Note that introducing $K_{iZ}$ into (9) and the result into the geometric models of (21), allows integrating (numerically) these two ODEs from $X_i = 0$ to $X_i = l_i$, and to reconstruct the pose fields $g_1$ and $g_2$ along the tubes. As a consequence, the two fields $K_{1Z}$ and $K_{2Z}$ entirely parameterize the configuration of a 2-tube planar CPPR, and the above formulation is no more than its kinetostatic model on the configuration space $\mathcal{C} = \mathcal{C}_1 \times \mathcal{C}_2$, with: $\mathcal{C}_i = \{K_{iZ} : X_i \in [0, l_i] \mapsto K_{iZ}(X_i) \in \mathbb{R}\} \bullet$

## 2.8 Lagrangian Reduction

The space $\mathcal{C} = \mathcal{C}_1 \times \mathcal{C}_2$ is a non-minimal configuration space (with more dimensions than degrees of freedom) of the 2-tube CPPR. This redundancy is reflected by the presence of constraints (14) and (19) and their associated Lagrange multipliers $\bar{F}_{c,i}$ in (26). For the purpose of simulation, it is convenient to remove the constraints (and Lagrange multipliers) from the formulation and to approximate the residual independent fields on a basis of Ritz functions (Boyer et al. 2021). This process, which will be detailed hereafter, is called reduction. It is here performed in two stages. The first stage is completed at a continuous level. It consists in isolating the subsystem above the baseplate, and re-parameterizing its configurations by

the curvature field $K_{1Z}$ of rod 1 only. In short, tube 1 is considered as a leader for tube 2, whose curvature $K_{2Z}$ follows $K_{1Z}$, in order to fulfill the coaxiality constraint along tube 1. As a consequence of this re-parameterization, the configuration of our 2-tube CPPR can now be entirely reconstructed from $K_{1Z} \in \mathcal{C}_1$, which becomes the new reduced configuration space of the CPPR. Moreover, the constraints (14) and (19) and their Lagrangian multipliers of (26), become useless and can be removed from the previous redundant formulation which becomes "minimal". Based on this minimal formulation, in a second stage, $K_{1Z}$ is projected on a Ritz basis of (modal) functions, according to the strain based reduction process of Boyer et al. (2021). At the end, the static model of the CPPR will take the usual Lagrangian matrix form

$$ Q_{\text{int}}(q) + Q_{\text{ext}}(q) = 0_{n \times 1}, \qquad (30) $$

where $q \in \mathbb{R}^n$ is a vector of generalized coordinates, while $Q_{\text{int}}$ and $Q_{\text{ext}}$ denote generalized internal and external force vectors, respectively. Friction being neglected, $Q_{\text{int}}$ (respectively $Q_{\text{ext}}$) will consist only of the elastic restoring forces (respectively the forces due to gravity, the push and pull forces of the actuation, and any tip loads).

### 2.8.1 First Reduction Stage: Continuous Reduction
From the above redundant formulation (14), (19), (21) to (27), and (29), the first reduction stage proceeds in two steps. The first step removes the guidance constraint (19). The second step removes the coaxiality constraint (14). Note that these two steps require all the equations related to tube 2 to be relabeled with the label $X_1$ of tube 1, as indicated in Section 2.6.2.

(1) First Step: Removing the Guiding Constraint: The relabeling of the guiding constraint (19) is only possible if we ignore the piece of tube 2 below the baseplate. Or in other words, if the CPPR is modeled by isolating the material system contained at each instant in a control volume above the baseplate. In fact, as it has been recently demonstrated in the broader context of sliding Kirchhoff rod dynamics by (Boyer et al. 2022), this simplification is fully justified under the quasi-static conditions studied here. Under these conditions, the relabeled version of the guiding constraint (19) degenerates into a single BC imposed on the tube 2, at the guide opening

$$ \Lambda_2(\tilde{X}_2(0)) = \tilde{\Lambda}_2(0) = F_{+,2} = (0, 0, 0, T_+, 0, 0)^T, \quad (31) $$

which replaces the BC (23) of the redundant formulation.

(2) Second Step: Removing the Coaxiality Constraint: We are looking for a relationship between $K_{2Z}$ and $K_{1Z}$, which is compatible with the coaxiality constraint (14). To derive such an expression, one needs to consider the $X_1$-derivative of the coaxiality constraint on poses (14)

$$ \frac{d\tilde{g}_2}{dX_1} = \tilde{g}_2 (Ad_{g_{2,1}}\xi_1 + \xi_{2/1})^\wedge, \qquad (32) $$

were $\xi_{2/1} = (g_{1,2}^{-1}(dg_{1,2}/dX_1))^\vee$ is the field of relative twists between the two tubes. Using (13), (17), and (18) applied to $f_2 = D_2$, and the definition of shearing angles (8), $\xi_{2/1}$ reads

$$ \xi_{2/1} = \begin{pmatrix} 0_{3 \times 1} \\ (dD/dX_1)E_Y \end{pmatrix} = \begin{pmatrix} 0_{3 \times 1} \\ (\beta_1 + h\tilde{\beta}_2)E_Y \end{pmatrix}, \quad (33) $$
which only depends on the deformation through h. Moreover, using (20) and invoking (10), one can use the consistent approximation

$$ \xi_{2/1}(X_1) \approx \begin{pmatrix} 0_{3 \times 1} \\ (\beta_1 - \tilde{\beta}_2)E_Y \end{pmatrix}. \quad (34) $$

Note that (32) being exclusively labeled by $X_1$, it defines an ODE along tube 1, with $\xi_1$ and $\xi_{2/1}$ defined by (7) and (34). Integrating this ODE from $X_1 = 0$ to $X_1 = l_1$, allows reconstructing the configuration of tube 2 (the $g_2$ field), from the knowledge of $K_{1Z}$ only, while preserving the coaxiality of both tubes (see Extension 4).

**Remark 7:** As well as being used to reconstruct the $g_2$ field, (7), (32), and (34) can be used to express $K_{2Z}$ directly in terms of $K_{1Z}$. Indeed, using (13) with $f_2 = g_2$ in the definition (7) of $\xi_2$, and introducing (32) in the result, gives

$$ \tilde{\xi}_2 = h^{-1}(Ad_{g_{2,1}}\xi_1 + \xi_{2/1}). \quad (35) $$

Now, using (34) and

$$ Ad_{g_{2,1}} = \begin{pmatrix} R_Z(\pi) & 0_{3 \times 3} \\ -D\hat{E}_Y R_Z(\pi) & R_Z(\pi) \end{pmatrix}, \quad (36) $$

in (35), and pre-multiplying the result by $(E_Z^T, 0_{1 \times 3})$, provides the relation between the curvature fields along the two tubes

$$ \tilde{K}_{2Z} = h^{-1}K_{1Z} \Rightarrow \tilde{K}_{2Z} = \frac{-K_{1Z}}{1 + \tilde{K}_{1Z}D}, \quad (37) $$

which, testifies of the minimal parameterization of the 2-tube planar CPPR configuration through $K_{1Z}$ •

**Remark 8:** When the two centroid lines are parallel, i.e. $(dD/dX_1) = 0$, the relation (37) can be directly deduced from planar geometric consideration based on $R_1 = R_2 + D$, with $R_1$ and $R_2$ the radii of the osculating circles tangent to the centroid lines $G_1$ and $G_2$ respectively •

**Remark 9:** The three relations (32), (35), and (37) are all consequences of the coaxiality. In fact, they indicate how the kinematic variables of tube 2, $g_2$, $\xi_2$, $K_{2Z}$, adapt to those of tube 1 in order to preserve this condition along deformation •

At this stage of the process, introducing (32), (35), and (37) along with (31) in the redundant formulation (14), (19), (21) to (27), and (29), allows removing the coaxiality and guiding constraints (14) and (19), and their Lagrange multipliers $\tilde{F}_{c,i}$ in (26).

To complete the process, one needs to relabel all the other equations related to tube 2 with respect to $X_1$, in which the above minimal parameterization is introduced. This is done as follows. In the balance of stress of (21), we replace $\xi_2$ by $\tilde{\xi}_2$ given by (35), and apply (13) to $f_2 = \Lambda_2$. This allows rewriting the stress balance of tube 2 into

$$ \frac{d\tilde{\Lambda}_2}{dX_1} - ad_{(Ad_{g_{2,1}}\xi_1+\xi_{2/1})}^T\tilde{\Lambda}_2 + h\tilde{F}_{\text{ext},2} = 0_{6 \times 1}, \quad (38) $$

where $\tilde{F}_{\text{ext},2}(X_1)$ now denotes the density of external (gravitational) wrenches applied on tube 2 per unit of $X_2$,

at $X_1$ along tube 1. Reconsidering (27), one finds

$$ \tilde{\tilde{F}}_{\text{ext},2} = (0_{1 \times 3}, \mu_2 \tilde{\Lambda}_2 (\tilde{R}_2^T a_g)^T)^T. \quad (39) $$

Using (37), the constitutive relation (29) of tube 2 reads

$$ \tilde{M}_{2Z} = E \tilde{I}_{2Z} \tilde{K}_{2Z} = - \frac{E \tilde{I}_{2Z} K_{1Z}}{1 + K_{1Z} D}. \quad (40) $$

The BC (22) and CC (24) do not change, while in the CC (25), replacing $\Lambda_2(0)$ by $\tilde{\Lambda}_2(l_1)$, gives

$$ \Lambda_1(l_1) = Ad_{g_{1,2}(l_1)}^{-T} \tilde{\Lambda}_2(l_1) + Ad_{g_{1,s}}^{-T} F_+. \quad (41) $$

**Remark 10:** At the end of this first reduction stage, the continuous redundant formulation (14), (19), (21) to (27), and (29) is replaced by a minimal continuous formulation defined for tube 2 by (32) and (38) to (40) with BCs (22), (24), and (31). For tube 1, the equations are unchanged, except that the coaxiality constraints and the associated Lagrange multipliers $\tilde{F}_{c,1}$ have been removed. This continuous reduced formulation represents the kinetostatic model of our 2-tube CPPR on the $\mathcal{C}_1$ configuration space of its reference tube •

### 2.8.2 Second Reduction Stage: Ritz Approximation
The second reduction stage consists in parameterizing $\mathcal{C}_1$ on a functional truncated basis of strain functions of tube 1. As proposed by Boyer et al. (2021), this is formally achieved by decomposing the allowed strain fields $K_{1Z}$ on a truncated basis of n Ritz functions $\Phi = (\Phi_1, \Phi_2, \dots, \Phi_n)$ as

$$ K_{1Z}(X_1) = \Phi(X_1)q, \quad (42) $$

where $q = (q_1, q_2, \dots, q_n)^T$ stands for the set of generalized coordinates introduced in (30). In this set of coordinates, the model takes the final form of (30), and is canonically deduced from the principle of virtual work (Tummers et al. 2023)

$$ \delta W_{\text{int}} + \delta W_{\text{ext}} = 0, \quad (43) $$

which holds for any variation of the configuration (virtual displacement field) compatible with the kinematics of the CPPR. In this balance, $\delta W_{\text{int}}$ and $\delta W_{\text{ext}}$, denote the virtual work of internal and external forces, respectively, which both take the generic form

$$ \delta W = \int_{0}^{l_1} \delta \xi_1^T \Lambda_1 dX_1 + \int_{0}^{l_2-a} \delta \xi_2^T \Lambda_2 dX_2, \quad (44) $$

where $\Lambda_i$ stands for the opposite of restoring wrench stress of (28) in the case of $\delta W_{\text{int}}$, and for the wrench stress that balance the external gravitational, actuation push-pull, and tip forces, in the case of $\delta W_{\text{ext}}$. In (43), the variations $\delta \xi_1 = (0, 0, \delta K_{1Z}, 0_{1 \times 3})^T$ and $\delta \xi_2 = (0, 0, \delta K_{2Z}, 0_{1 \times 3})^T$ are compatible with the coaxiality constraint (14), i.e. using the variation of (35) with (34), they are related by

$$ \delta \tilde{\xi}_2 = h^{-1}(Ad_{g_{2,1}}\delta \xi_1 - \delta h \tilde{\xi}_2), \quad (45) $$

whose the only non-zero component is

$$ \delta \tilde{K}_{2Z} = \frac{-\delta K_{1Z}}{(1 + K_{1Z} D)^2}, \quad (46) $$
and can be directly deduced from the variation of (37). Finally, let us change the integral along tube 2 of (44) into an integral along tube 1 through the relabeling process of Section 2.6.2. Then, using (45) with the Ritz reduction (42), (44) can be reformulated in the form

$$ \delta W = \delta q^T Q(q), \quad (47) $$

with

$$ Q = \int_{0}^{l_1} \Phi^T (E_Z^T, 0_{1 \times 3}) \left[ \Lambda_1 - \left( Ad_{g_{2,1}}^T - \left[ \frac{\partial h}{\partial \xi_1} \right] \tilde{\xi}_2^T \right) \tilde{\Lambda}_2 \right] dX_1, \quad (48) $$

where we noted $(\partial h / \partial \xi_1) = (0, 0, \partial h / \partial K_{1Z}, 0, 0, 0)^T$. The formula (48) allows converting any pair of wrench-stress fields $(\Lambda_1, \tilde{\Lambda}_2)$ transmitted across the two tubes, into the corresponding generalized force vector $Q$ applied to the CPPR. In particular, for the generalized vector of internal forces $Q_{int}$, using (29), (40) and (42), the non-zero components of $\tilde{\Lambda}_1 = -\mathcal{H}_1(\xi_1 - \xi_1^o)$ and $\tilde{\Lambda}_2 = -\tilde{\mathcal{H}}_2(\tilde{\xi}_2 - \tilde{\xi}_2^o)$, are

$$ M_{1Z} = -EI_{1Z}\Phi q, \quad \tilde{M}_{2Z} = \frac{E\tilde{I}_{2Z}\Phi q}{(1 + D\Phi q)}, \quad (49) $$

which, once introduced in (48), provides the detailed expression for $Q_{int}$

$$ Q_{\text{int}}(q) = - \int_{0}^{l_1} \Phi^T \left( EI_{1Z} + \frac{E \tilde{I}_{2Z}}{(1 + D\Phi q)^2} \right) \Phi q \, dX_1. \quad (50) $$

Ultimately, to practically compute $Q_{ext}$, a computed-torque style algorithm is applied. This algorithm, which also computes (50), is detailed in the next section.

## 2.9 Simulation Algorithm of a 2-Tube Planar CPPR

### 2.9.1 Algorithm for Computing $Q_{int}$ and $Q_{ext}$
While $Q_{int}(q)$ can be directly calculated for any $q$ with the explicit expression (50), the calculation of $Q_{ext}(q)$ requires more attention. This essential difference is explained by the fact that $Q_{int}$ is directly generated by a field of stress wrench (governed by the constitutive relation (28)), while in contrast, the external gravity, push-pull, and tip forces are primarily defined as wrenches of forces⁴ $\tilde{F}_{ext,i}$, $F_{+,2}$, and $F_+$ respectively. A natural way to calculate $Q_{ext}(q)$ is to first convert $\tilde{F}_{ext,i}$, $F_{+,2}$ and $F_+$ into wrench stress fields $\tilde{\Lambda}_1$ and $\tilde{\Lambda}_2$, and then to introduce these wrench fields into (48). This conversion is achieved with an inverse (computed-torque) algorithm, proposed by Boyer et al. (2021), adapted to the statics of CPPRs. This algorithm is structured in two passes. The first one is a forward kinematic pass that computes, for a given $q$, the pose field of the two tubes from the base $X_1 = 0$ of tube 1, to the intersection of tube 2 with the baseplate $\tilde{X}_2(0) = l_2 - a$, via the rigid end-effector connection. The second pass is a backward pass that calculates the stress generated by the external forces along the structure, starting from the intersection of tube 2 with the baseplate, and ending at the base of tube 1, going through the rigid end-effector

connection. In a last step, these two stress fields are used in (48) to compute $Q_{ext}$. Finally, this algorithm also computes $Q_{int}$ through (50), which requires the numerical integration of the face-to-face function (20). It is summarized as follows.

#### ■ Inputs of the algorithm:

$$ \xi_1 = (0, 0, \Phi q, 1, \beta_1, 0)^T. \quad (51) $$

#### ■ Forward kinematics (1st pass):

$$ \text{Initial conditions: } g_1(0) = g_{1,-}. \quad (52) $$

Integrate from $X_1 = 0$ to $X_1 = l_1$

$$ \frac{dg_1}{dX_1} = g_1 \hat{\xi}_1. \quad (53) $$

Initial conditions

$$ \tilde{X}_2(l_1) = 0, \quad \tilde{g}_2(l_1) = g_1(l_1)g_{1,2}(l_1). \quad (54) $$

Integrate from $X_1 = l_1$ to $X_1 = 0$:

$$ \frac{d\tilde{X}_2}{dX_1} = h = -(1 + D\Phi q), \quad (55) $$

$$ \tilde{\xi}_2 = h^{-1} (Ad_{g_{2,1}} \xi_1 + \xi_{2/1}), \quad (56) $$

$$ \frac{d\tilde{g}_2}{dX_1} = \tilde{g}_2(h\tilde{\xi}_2)^{\wedge}. \quad (57) $$

$$ \text{Store } (g_1, \tilde{X}_2, \tilde{g}_2, \tilde{\xi}_2)(\cdot). \quad (58) $$

#### ■ Backward statics (2nd pass):

Initial conditions:

$$ \tilde{\Lambda}_2(0) = (0_{1 \times 3}, T_+, 0, 0)^T. \quad (59) $$

Integrate from $X_1 = 0$ to $X_1 = l_1$

$$ \frac{d\tilde{\Lambda}_2}{dX_1} = h ad_{\tilde{\xi}_2}^T \tilde{\Lambda}_2 - h \tilde{F}_{\text{ext},2}. \quad (60) $$

Initial conditions:

$$ \Lambda_1(l_1) = Ad_{g_{1,2}(l_1)}^{-T} \tilde{\Lambda}_2(l_1) + Ad_{g_{1,s}}^{-T} F_+. \quad (61) $$

Integrate from $X_1 = l_1$ to $X_1 = 0$

$$ \frac{d\Lambda_1}{dX_1} = ad_{\tilde{\xi}_1}^T \Lambda_1 - \tilde{F}_{\text{ext},1}. \quad (62) $$

$$ \text{Store } (\tilde{\Lambda}_2, \Lambda_1)(\cdot). \quad (63) $$

#### ■ Calculation of the outputs of the algorithm:

Compute:

$$ Q_{\text{ext}} = \int_{0}^{l_1} \Phi^T \left( E_Z^T \right)^T \left[ \Lambda_1 - \left( Ad_{g_{2,1}}^T - \left[ \frac{\partial h}{\partial \xi_1} \right] \tilde{\xi}_2^T \right) \tilde{\Lambda}_2 \right] dX_1 \quad (64) $$

$$ Q_{\text{int}} = - \int_{0}^{l_1} \Phi^T \left( EI_{1Z} + \frac{E_2 \tilde{I}_{2Z}}{(1 + D\Phi q)^2} \right) \Phi \, dX_1 q. \quad (65) $$

The space integration of $Q_{int}$ is performed by first computing $\tilde{X}_2$ in each $X_1$ from (55), and then by computing the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Schematic of the two passes of the algorithm. For both passes, the integration starts at the circled bullet with the corresponding initial condition and follows through the arrows from one tube to the other. The dotted lines show the associated differential equations. The light gray lines link the quantities that are necessary for the implementation of the algorithm to the corresponding equation numbers of the paper. For a step-by-step version of the graphical algorithm, see Extension 5.</figcaption></figure>

corresponding values of $I_{2Z}$, i.e. $\tilde{I}_{2Z}(X_1) = I_{2Z}(\tilde{X}_2(X_1))$, with (3), (104), and (107), applied to the interpolated cross-sectional geometry of Section 2.3. A similar process holds for the computation of $\tilde{D}_2$, $\tilde{\beta}_2$ and $\tilde{A}_2$ in (55) and (57); and $\tilde{F}_{\text{ext},2}$ of (60). This algorithm is summarized in Figure 9 where the equation numbers of the corresponding definitions are given (see also Extension 5). Ultimately, this algorithm defines a numerical input-output map that produces $Q_{\text{ext}}(q)$ and $Q_{\text{int}}(q)$ from $q$, $T_+$, and $f_+$.

### 2.9.2 Quasi-Static Simulation of a 2-Tube Planar CPPR

The simulation of a given CPPR begins by calculating the cross-section parameters (3) to (5) of the rods equivalent to its tubes, which are themselves deduced from the interpolation process of Section 2.3. Using these parameters, the static balance of (30), where $Q_{\text{int}}(q)$ and $Q_{\text{ext}}(q)$ are computed by the above algorithm defines a set of nonlinear

algebraic equations, which can be set into the usual form

$$ \mathfrak{R}\mathfrak{e}\mathfrak{s}(q) = Q_{\text{int}}(q) + Q_{\text{ext}}(q) = 0_{n \times 1}, \quad (66) $$

with $\mathfrak{R}\mathfrak{e}\mathfrak{s}(q)$, the residual vector of the system, and $Q_{\text{int}}(q)$, $Q_{\text{ext}}(q)$ computed by the above algorithm. The system (66) can be solved with respect to $q$ at each step of load ($T_+, F_+$) with standard root-finding algorithms (e.g. Newton-Raphson, Levenberg-Marquardt). In the next section, this approach is applied to different robot designs. The iterative search for the roots of (66) is performed with Matlab fsolve function with the Levenberg-Marquardt algorithm with all parameters set to defaults. The Jacobian matrix of the residual vector $(\partial \mathfrak{R}\mathfrak{e}\mathfrak{s}/\partial q)(q)$ is numerically approximated by finite differences. The spatial integrations of the above computed-torque algorithm are performed with a spectral method (Trefethen 2000).

**Remark 11:** Note that so far we have only considered an external tip wrench $F_+$, the detailed model of which will be required by the experimental processing in Section 4.7. However, the above algorithm can take into account loads concentrated at any arbitrary point along the length of the robot. In this case, one can (virtually) cut the outer tube into two pieces connected at the point where the load is applied to. Then, the backward pass of the stress calculation algorithm is stopped at this junction point, and the concentrated load is added to the final conditions of the first piece to feed the initial conditions of the second. This approach has been implemented on TACRs in (Tummers et al. 2023), where discontinuities in tendon paths introduce concentrated stresses at junctions between segments •

## 3 Extension to Three-Dimensional (3D) Multi-Tube CPPRs

So far, we have considered the case of a 2-tube CPPR capable of 2D bending. We will now see how we can apply the above modeling approach to the more general case of a multi-tube CPPR capable of both twisting and bending in three dimensions. In line with our approach, this means that we now have to tackle a tree-like system of Cosserat rods, in which the proximal end of tube 1 is attached to the base, while its other end supports a rigid moving platform to which $N-1$ other tubes, numbered $k=2, 3, \dots, N$, are rigidly connected. All these tubes can be bent and twisted in 3D while sharing a common centerline, or more simply they must respect the coaxiality condition. In the most general case, the centroid lines of all the tubes are not necessarily regular at rest, and even more so may spiral around the central axis of the CPPR, which is assumed to be always straight in the configuration at rest. As in the planar case, tube 1 of this CPPR will be considered as the reference (or leader) for the others (the followers) in a two-stage kinematic reduction, the first of which is based on exploiting the coaxiality condition. Another difference with the planar case, is that now the follower tubes can be subject both to a push-pull force $T_+$ along $e_X$, and a torque $\Upsilon_+$ around $e_X$, at the level of the baseplate. In the rest of this section, we will follow the 2D approach step by step. This begins with the 3D parameterization of a CPPR at rest, and the statement of its coaxiality constraints, which will then be
used to deduce a minimal continuous parameterization of the deformations. We then continue with the Ritz reduction, and explain how the associated generalized forces entering into the static equilibrium can be calculated in this wider context. In order not to overload the presentation, we will not go into as much detail as for the 2D modeling. In particular, while we will insist on the kinematic aspects of the first stage of the reduction (Section 2.8.1 before Remark 9), we will not detail the reduction of the stress balance of the tubes (21), which under the assumption of frictionless contacts, can be performed with the same arguments as those used in Section 2.8.1 (after Remark 9). Despite these shortcuts in the presentation, all the key relationships required by the quasi-static simulation of a multi-tube CPPR in 3D are shown below and can be used on their own.

## 3.1 3D Parameterization of the CPPR at Rest

As in the planar case, we must first define a parameterization of the tubes assembled in the resting configuration of the CPPR. This parameterization fixes the design of the CPPR in the frame of its baseplate (O, e_X, e_Y, e_Z). The CPPR is straight at rest, with its centerline supported by e_X and denoted C^o, where from now on, a right upper index o, indicates that a quantity or a geometric object is related to the reference configuration. The cross-sectional frames (G^o_i, b^o_{iX}, b^o_{iY}, b^o_{iZ})(X_i) of tube i are fixed consistently with the choices of Section 2.4, i.e. G^o_i is the centroid of the X_i-cross section, and

$$b^o_{iX} = e_X, \quad b^o_{iY} = \overline{G^o_i C^o_i} / D_i, \quad b^o_{iZ} = b^o_{iX} \times b^o_{iY}, \quad (67)$$

where C^o_i(X_i) is the orthogonal projection of G^o_i(X_i) onto C^o, while D_i(X_i) = ||G^o_i C^o_i|| still denotes the offset between the centerlines and the centroid lines as represented for i = 1, in Figure 10(a).

The centroid lines at rest G^o_i are now 3D curves that can be spiraled around C^o_i. Referring to Section 2.5, this more general context requires the introduction of other design offsets than the shear angles of the planar case, in the definition of the tube space-twist fields ξ_i, for each tube i. These offsets are naturally defined as the components of ξ^o_i, which describes how the frames of a tube i vary as we move along its centroid line G^o_i in the direction of increasing X_i. As we shall see in the next section, thanks to the coaxiality condition C^o_i = C^o, i = 1, 2 .. N, all the ξ^o_k, k = 2, 3 .. N, can be defined from the knowledge of

$$\xi_1^o = \left( (g_1^o)^{-1} \frac{dg_1^o}{dX_1} \right)^\vee, \quad (68)$$

and the value at rest of a further new angle that will describe how the followers rotate inside the leader, when the CPPR deforms. To calculate ξ^o_1, we use (68) in which g^o_1 = (R^o_1, r^o_1) is the pose of the frames along tube 1 defined by the conventions (67), with r^o_1(X_1) the inertial position of G^o_1(X_1) along G^o_1, and R^o_1 = (b^o_{1X}, b^o_{1Y}, b^o_{1Z}) the orientation of the X_1-cross-sectional frame in the inertial frame (O, e_X, e_Y, e_Z). After calculation, one finds

$$\begin{align}
\xi_1^o &= \begin{pmatrix} K_1^o \\ \Gamma_1^o \end{pmatrix} = \begin{pmatrix} (R_1^{oT} R_1^{o'})^\vee \\ R_1^{oT} r_1^{o'} \end{pmatrix} \tag{69} \\
&= \begin{pmatrix} K_{1X}^o E_X \\ E_X + D_1^T E_Y + D_1 K_{1X}^o E_Z \end{pmatrix}, \nonumber
\end{align}$$

<figure><img src="image_11.png" /><figcaption>Figure 10. (a) Perspective view of one tube (number 1) in the resting configuration of the CPPR. In red are indicated the offset twists due to design. (b): Section view of the leader tube (numbered 1) and a follower (numbered k) in a CPPR able to deform in 3D. Tube k can rotate of an angle θ_k, with respect to tube 1 around the centerline of the CPPR. Note that this figure generalizes the planar context of Figure 8(c) where we had θ_k = 0.</figcaption></figure>

where, from now on, the derivation d./dX_1 will be denoted by a'. In (69), K^o_{1X} models how the frames of tube 1 are twisted along the CPPR, while as in the planar case, D'_1 is a shear offset which models the way in which the centroid line moves towards or away from the CPPR centerline. Finally, it should be noted that although we started from the definition of the cross-sectional frames (67) to deduce the field of offset twist (69), in practice, it is easier to first define D_1 and K^o_{1X}. Then, by introducing these 2 functions into (69), one can integrate g^o_1' = g^o_1 \hat{\xi}_1^o, from g^o_1(0) = g_{1,-} (with g_{1,-} inertial pose of the cross-section of tube 1 at the baseplate), and obtain the frame field (67) of tube 1. This second procedure has the advantage of avoiding the singularity D_1 = 0 of (67). For the other tubes, k = 2, 3 .. N (the followers), their cross-sectional frames will be positioned in relation to those of tube 1 (the leader) both in the resting and the deformed configuration of the CPPR. This is made possible by the coaxiality condition and an additional torsion angle, as described in the next section.

## 3.2 Coaxiality Constraints in 3D

The reduced model of a planar CPPR with 2 tubes can naturally be extended to model the 3D deformations of a CPPR with an arbitrary number of tubes. To this end, we will continue the ideas developed in the planar case, and express, as far as possible, the deformation of tubes k = 2, 3 .. N (the followers), as a function of the strain fields of tube 1 (the leader). To achieve this, we reapply the coaxiality condition (14) to all the follower tubes k = 2, 3 .. N

$$\forall X_1 \in [0, l_1] : g_k(\tilde{X}_k(X_1)) = g_1(X_1)g_{1,k}(X_1), \quad (70)$$

In this expression, g_{1,k} is the relative pose of the coplanar cross-sections of the two tubes, which is now defined by (see Figure 10(b))

$$g_{1,k}(X_1) = \begin{pmatrix} R_Z(\pi)R_X(\tilde{\theta}_k) & R_Z(\pi)(D_11_{3\times3} + \tilde{D}_k R_X(\tilde{\theta}_k))E_Y \\ 0_{1\times3} & 1 \end{pmatrix}, \quad (71)$$
where $\tilde{\theta}_k(X_1)$ is the angle of a planar rotation $R_X(\tilde{\theta}_k) = \exp(\tilde{\theta}_k \tilde{E}_X)$ of the $\tilde{X}_k(X_1)$-section of tube $k$, with respect to the $X_1$-cross section of tube 1, around the (possibly deformed) centerline $C$ of the CPPR. Since such an angle exists for every cross-section along tube 1, this defines a torsion angle field along the CPPR which, in the most general case, takes the form

$$ \tilde{\theta}_k(X_1) = \tilde{\theta}_k^o(X_1) + \Delta\tilde{\theta}_k(X_1), \quad (72) $$

where $\tilde{\theta}_k^o$ is a design torsion angle field that models how the centroid line of the tube $k$ wraps around the central axis of the CPPR in its resting configuration. On the other hand, $\Delta\tilde{\theta}_k$ defines an additional torsion angle field due to deformation. Along $[0, l_1]$, $\Delta\tilde{\theta}_k$ defines an independent kinematic field of the CPPR deformation.

**Remark 12:** Although the torsional fields of tube $k$ were naturally introduced in (71), as functions of $X_1$, as the tilde above them indicates, they are first naturally defined as functions of $X_k$. As a result, $\tilde{\theta}_k^o$ and $\Delta\tilde{\theta}_k$ of (72) are defined by a relabeling process similar to that used in the planar case (12), but now operating between each of the follower tubes and the leader. In particular, the pullback formula on the derivatives of functions (13) now holds along any tube $k$ with $h_k(X_1) = (d\tilde{X}_k/dX_1)(X_1)$, and we have for any function $f_k$ of $X_k$

$$ \left( \frac{d\tilde{f}_k}{dX_1} \right) = h_k \left( \frac{df_k}{dX_k} \right)^{\sim} \bullet \quad (73) $$

**Remark 13:** Note that the angular field $\theta_k^o$ is a design parameter that must be added to the 3 independent parameters $(\rho_i, \bar{\rho}_i, \alpha_i)$ of the planar design from Section 2.2. Like the notch depth field $\alpha_i$, this design torsion angle field can be deduced by interpolating its values notch by notch along the robot, as in (1). As regards tube 1, a similar design parameter must be defined either by an angle field $\theta_1^o$ or through the torsion rate $\theta_1^{o'} = K_{1X}^o$ of Section 3.1 •

Differentiating (70) with respect to $X_1$, with (73), provides the coaxiality constraints on space-twists

$$ \tilde{\xi}_k = h_k^{-1} \left( \tilde{g}_k^{-1} \frac{d\tilde{g}_k}{dX_1} \right)^{\vee} = h_k^{-1} \left( Ad_{g_{k,1}} \xi_1 + \xi_{k/1} \right) . \quad (74) $$

This relationship has the same form as (35) but with $Ad_{g_{k,1}}$ and $\xi_{k/1} = (g_{1,k}^{-1}(dg_{1,k}/dX_1))^{\vee}$ now calculated with (71). Performing these calculations, we obtain

$$ \xi_{k/1} = \Delta\xi_{k/1} + \xi_{k/1}^o, \quad (75) $$

with

$$ \Delta\xi_{k/1} = \begin{pmatrix} E_X \\ \tilde{D}_k E_Z \end{pmatrix} \Delta\tilde{\theta}_k', \quad (76) $$

$$ \xi_{k/1}^o = \begin{pmatrix} \tilde{\theta}_k^{o'} E_X \\ (ck D'_1 + \tilde{D}'_k)E_Y - (sk D'_1 - \tilde{D}_k \tilde{\theta}_k^{o'})E_Z \end{pmatrix}, \quad (77) $$

where we used the notations $sk = \sin(\tilde{\theta}_k)$, $ck = \cos(\tilde{\theta}_k)$. Note that imposing $\theta_k = 0$ in (77), one obtains again the expression (33) of the planar case.

### 3.3 Minimal Parameterization of a CPPR in 3D

Unlike the 2D case, the expression (75) to (77) of $\xi_{k/1}$ has a non-zero shear component which depends on the deformation via the torsion $\Delta\theta_k$ of tube $k$. On the basis of these considerations, we use the inextensible but shearable rod model for all tubes. With this choice, the 3D configurations of the CPPR can be parameterized by the $\Delta K_1 = K_1 - K_1^o$ vector field of the leader, which enters in its twist field $\xi_1$, as a deformation component added to the offset field (69)

$$ \xi_1 = \xi_1^o + \begin{pmatrix} \Delta K_1 \\ D_1(E_X^T \Delta K_1) E_Z \end{pmatrix}, \quad (78) $$

and by the $N-1$ fields of torsion angles $(\Delta\theta_2, \Delta\theta_3, \dots, \Delta\theta_k)$ of the followers. Introducing these angles into the coaxiality constraints (74), with (75) to (77), then provides the $\tilde{\xi}_k$ for each follower tube $k$.

**Remark 14:** Note that although $D_1(E_X^T \Delta K_1) = D_1 \Delta K_{1X}$ appears in (78) as a shear component, it only depends on the torsion deformation $\Delta K_{1X}$. This is also the case for follower tubes with the term $\tilde{D}_k \Delta\tilde{\theta}_k'$ in (76). These apparent shears are a kinematic consequence of torsion and eccentricity of centroids. Geometrically, they are due to the fact that the coaxiality condition requires $\xi_1$ and all the vectors $\xi_k$ to be referred to the center (and not the centroid) of the sections •

As in the planar case, using the inextensibility of all tubes, the projection of the linear component of the constraints (74) onto $E_X$ provides the detailed expression of the $X_1$-rate of the face-to-face function between the leader tube and all the followers $k = 2, 3, \dots, N$

$$ \begin{aligned} h_k(X_1) &= \frac{d\tilde{X}_k}{dX_1}(X_1) \\ &= - \left( 1 + D_1 K_{1Z} + \tilde{D}_k (ck K_{1Z} + sk K_{1Y}) \right). \end{aligned} \quad (79) $$

Backward integrating (79) along the leader, defines the face-to-face function $\tilde{X}_k(X_1)$ of the relabeling process along each of the followers. Finally, note that by imposing $\theta_k = 0$ in (79), we again obtain expression (20) of the planar case.

**Remark 15:** Similarly to the 2D case, $\Delta K_1$ and all the fields $\Delta\tilde{\theta}_k$ can be used to reconstruct the configuration of the CPPR, i.e. the inertial pose fields of all its constitutive tubes. For the leader, this is performed by integrating, from $X_1 = 0$ to $X_1 = l_1$, the ODE

$$ g'_1 = g_1 \tilde{\xi}_1, \quad (80) $$

initialized by $g_1(0) = g_{1,-}$, and with $\xi_1$ given by (78). For the follower tubes, one has to $X_1$-integrate, from $X_1 = l_1$ to $X_1 = 0$, the ODEs

$$ \tilde{g}'_k = \tilde{g}_k (h_k \tilde{\xi}_k)^{\wedge}, \quad (81) $$

initialized by $\tilde{g}_k(l_1) = g_1(l_1)g_{1,k}(\theta_k^o(l_1))$, with $\tilde{\xi}_k$ given by $\xi_1$ and (74) to (77). Note that $\tilde{g}_k$ can be equivalently reconstructed by integrating the ODE (32) with $k$ replacing $2$ •
**Remark 16:** Note that in the above model, (75) to (77) require to compute $\tilde{D}'_k$, $\Delta\tilde{\theta}'_k$, and $\tilde{\theta}'_k$, i.e. the $X_1$-derivatives of functions that originally depend on $X_k$. To practically compute them, we use the pullback formula (73)

$$ \tilde{D}'_k = h_k(dD_k/dX_k)\sim, \quad \Delta\tilde{\theta}'_k = h_k(d\Delta\theta_k/dX_k)\sim, \\ \tilde{\theta}'_k = h_k(d\theta_k^o/dX_k)\sim. \quad (82) $$

Moreover, as it has been done for the shear offsets $\tilde{D}'_k$ that were already present in the planar case, one can apply the further approximations: $|dD_k/dX_k| \simeq |d\theta_k^o/dX_k| \simeq |d\Delta\theta_k/dX_k| \simeq \varepsilon$, which simply means that in any configuration (deformed or not), the winding of the centroid line of tube $k$ around the central axis of the CPPR remains moderate. Using (79) in (82), and the notation $(df_k/dX_k)\sim = (f'_k)\sim$, we have at first order w.r.t. $\varepsilon$: $\tilde{D}'_k \simeq -1(D'_k)\sim$, $\Delta\tilde{\theta}'_k \simeq -1(\Delta\theta'_k)\sim$, $\tilde{\theta}'_k \simeq -1(\theta'_k)\sim$, and (75) to (77) are changed into

$$ \xi_k/1 \simeq \nu_k(\Delta\theta'_k)\sim + \mu_k, \quad (83) $$

with

$$ \nu_k = -(1, 0, 0, 0, 0, \tilde{D}_k)^T \quad (84) $$

and

$$ \mu_k = (-(\theta_k^{o'})\sim, 0, 0, 0, ckD'_1 - (D'_k)\sim, \\ -skD'_1 - \tilde{D}_k(\theta_k^{o'})\sim)^T. \quad (85) $$

Finally, let us note that in contrast to $\tilde{D}'_k$, $\tilde{\theta}'_k$ and $\Delta\tilde{\theta}'_k$ of (76) and (77), the derivatives $(D'_k)\sim$, $(\theta'_k)\sim$, and $(\Delta\theta'_k)\sim$ in the above expressions are practically (and numerically) easy to compute from the knowledge of the original functions $\theta_k^o$, $\Delta\theta_k$, and $D_k$ of $X_k$. In fact, we just have to differentiate these functions with respect to $X_k$, change their sign, and then apply the relabeling from $X_k$ to $X_1$ •

Finally, (74) and (83) to (85) are the relations that the space-twists of the tubes must respect in order to preserve the coaxiality of the CPPR when it deforms in 3D.

## 3.4 Kinetostatic Model and Simulation in 3D

The kinematic model (74), (78), and (83) to (85) corresponds to the first (continuous) stage of the kinematic reduction in the 3D case. To achieve the second stage, all the independent kinematic variables of this reduced continuous model are approximated on truncated Ritz functional basis as

$$ \Delta K_1(X_1) = \Phi_{\epsilon}(X_1)q_{\epsilon}, \quad \Delta \tilde{\theta}_k(X_1) = \Phi_{\theta_k}(X_1)q_{\theta_k}, \quad (86) $$

where $\Phi_{\epsilon}$ and $\Phi_{\theta_k}$ are $3 \times n_1$ and $1 \times n_k$ matrices of functions (typically polynomials). Note here that since the follower tubes are fixed to the tip rigid section of the CPPR ($X_1 = l_1$), we have $\Delta\tilde{\theta}_k(l_1) = 0$, and we need to take all the polynomial sets $\Phi_{\theta_k}$, such that

$$ \Phi_{\theta_k}(l_1) = 0. \quad (87) $$

As in the planar case, (86) defines the reduced space of generalized coordinates, here detailed as $q = (q_{\epsilon}^T, q_{\theta_2}^T, q_{\theta_3}^T, \dots, q_{\theta_N}^T)^T \in \mathbb{R}^{n_1} \times \mathbb{R}^{n_2} \times \dots \times \mathbb{R}^{n_N}$. On this configuration space, the static balance of a 3D CPPR

takes the Lagrangian form:

$$ \begin{pmatrix} Q_{\epsilon, \text{int}} \\ Q_{\theta, \text{int}} \end{pmatrix} + \begin{pmatrix} Q_{\epsilon, \text{ext}} \\ Q_{\theta, \text{ext}} \end{pmatrix} = \begin{pmatrix} 0_{n_1 \times 1} \\ 0_{(n_2+n_3+\dots+n_N) \times 1} \end{pmatrix}, \quad (88) $$

where $(Q_{\epsilon, \text{int}}^T, Q_{\theta, \text{int}}^T)^T$ is again the vector of the generalized restoring elastic forces inside the system and $(Q_{\epsilon, \text{ext}}^T, Q_{\theta, \text{ext}}^T)^T$ is that of the external forces induced by gravity, external forces applied along the robot, and the control forces $T_{k,+}$ and torques $\Upsilon_{k,+}$, exerted on the follower tubes at the level of the baseplate. As in the planar case, the static balance (88) defines a system of nonlinear algebraic equations that can be solved with respect to $q$ with the Levenberg-Marquardt algorithm at each step of load. Once $q$ is calculated, it is used to calculate $\xi_1$ and all the $\tilde{\xi}_k$, and finally to reconstruct the pose fields $g_1$ and $\tilde{g}_k$ as indicated in Remark 15.

To calculate the generalized forces of (88), we do as in the planar case, i.e. we first compute $\Lambda_1$ and all the $\Lambda_k$ using the computed-torque algorithm of Section 2.9, now applied to a tree like-structure, and project them onto the reduced space of the modal generalized coordinates $q = (q_{\epsilon}^T, q_{\theta}^T)^T$. To obtain the expression of these projection operators, the virtual strain twists $\delta\xi_1, \delta\xi_2, \dots, \delta\xi_N$, along which the stress wrenches $\Lambda_1, \Lambda_2, \dots, \Lambda_N$ work, are first expressed in terms of $\delta K_1, \delta\theta_2, \delta\theta_3, \dots, \delta\theta_N$, using the variation of (74), (78), and (79), and finally as functions of $\delta q = (\delta q_{\epsilon}^T, \delta q_{\theta}^T)^T$ with the variation of (86). Then using the identity of virtual works

$$ \begin{aligned} \delta W &= \int_0^{l_1} \delta \xi_1^T \Lambda_1 dX_1 + \sum_{k=2}^{N} \int_0^{l_k-a_k} \delta \xi_k^T \Lambda_k dX_k \quad (89) \\ &= \delta q_{\epsilon}^T Q_{\epsilon} + \delta q_{\theta}^T Q_{\theta}, \end{aligned} $$

provides by simple identification, the generalized forces $Q_{\epsilon}$ and $Q_{\theta}$ produced by any set of $(\Lambda_1, \Lambda_2, \dots, \Lambda_N)$. This dualization process gives after a few algebraic operations, the $n_1 \times 1$ vector

$$ Q_{\epsilon} = \int_{0}^{l_1} \Phi_{\epsilon}^{T} (1_{3 \times 3} \quad D_1 E_X E_Z^{T}) \begin{bmatrix} \Lambda_1 \\ -\sum_{k=1}^{n_k} \left( Ad_{g_{k,1}}^{T} - \left[ \frac{\partial h_k}{\partial \xi_1} \right] \tilde{\xi}_k^{T} \right) \tilde{\Lambda}_k \end{bmatrix} dX_1, \quad (90) $$

as well as $Q_{\theta} = (Q_{\theta_1}^T, Q_{\theta_2}^T, \dots, Q_{\theta_N}^T)^T$, with the $n_k \times 1$ vectors

$$ \begin{aligned} Q_{\theta_k} = & -\int_0^{l_1} \Phi_{\theta_k}^{T} \left[ \nu_k^T \tilde{\Lambda}_k \right] dX_1 \\ & - \int_0^{l_1} \Phi_{\theta_k}^{T} \left[ \frac{\partial \mu_k}{\partial \theta_k} + \nu_k^T a d_{(Ad_{g_{k,1}} \xi_1)}^{T} \right. \\ & \left. - \left( \frac{\partial h_k}{\partial \theta_k} \right) \tilde{\xi}_k^T \right] \tilde{\Lambda}_k dX_1, \end{aligned} \quad (91) $$

where we noted $(\partial h_k / \partial \xi_1) = ((\partial h_k / \partial K_1)^T, 0_{1 \times 3})^T$, with $(\partial h_k / \partial K_1) \in \mathbb{R}^3$ the gradient vector of $h_k$ with respect to the components of $K_1$; $(\partial h_k / \partial \theta_k) \in \mathbb{R}$, is the usual derivative of $h_k$ w.r.t. $\theta_k$, and $(\partial \mu_k / \partial \theta_k) \in \mathbb{R}^6$ is the vector whose components are the $\theta_k$-derivatives of those of $\mu_k$ given by (85), i.e. $(\partial \mu_k / \partial \theta_k) = D'_1(0, 0, 0, 0, -sk, -ck)^T$. Note that (90) is the 3D generalization of (48), whereas (91) was absent from the 2D case.
Finally, introducing in (90) and (91), the stress fields $(\Lambda_1, \Lambda_2, \dots, \Lambda_N)$ produced by the external forces and the internal (restoring) ones, provides the expected generalized external and internal forces of the Lagrangian balance (88). To numerically compute these sets of stress fields, we proceed as follows:

* As regards the external forces, the corresponding stress $(\Lambda_1, \Lambda_2, \dots, \Lambda_N)$ are computed with the two passes algorithm of Section 2.9, where a loop on the index $k$ replaces each of the calculations related to tube 2, while the backward ODEs on stress are now initialized by $F_{k,+} = (\Upsilon_{k,+}, 0, 0, T_{k,+}, 0, 0)$.

* As regards the internal (restoring forces), we use the full Hookean constitutive law of Cosserat rods (28), with $H_i$ given by (5), to get for $k = 2, 3, \dots, N$

$$ \Lambda_1 = H_1(\xi_1 - \xi_1^o), \quad \tilde{\Lambda}_k = \tilde{H}_k(\tilde{\xi}_k - \tilde{\xi}_k^o), \quad (92) $$

where $\xi_1^o$ and $\tilde{\xi}_k^o$ are the values of the twist fields $\xi_1$ and $\tilde{\xi}_k$ in the resting configuration of the CPPR. In detail, $\xi_1 - \xi_1^o$ is given by (78), whereas $\tilde{\xi}_k^o$ is obtained by applying (74) with (83) to (85) in the rest configuration, i.e. with $(\Delta\theta'_k) \sim 0$ in (83), and $\Delta\theta_k = 0$ in (85).

# 4 Experimental Validation

This section starts with a set of experiments with various 2-tube CPPR designs (see Table 1 and Extension 6) in order to validate our modeling approach in the planar case (i.e. the reduced model of Section 2.8). Designs with different dimensions and centroid lines are first validated in free-space, then externally loaded. Following, we focus on the 3D cases where out-of-plane loading is highlighted on a 2-tube CPPR and the 3-tube CPPR design is validated. The latter two scenarios were introduced experimentally in (Childs and Rucker 2023). Finally, extending beyond the state-of-the-art, more complex 3D cases exhibiting combined bending and torsion are introduced, namely CPPRs with helicoidal centroid lines actuated in translation, and a novel actuation approach consisting in combined translation and rotation at the base of each component tube of a CPPR.

## 4.1 Materials

Each of the tested prototypes is identified with a symbol $R_k$, with $k = 1..4$. The prototypes are subjected to different push-pull forces $T_+$, while preliminary tests with tip loads are presented. The basis of the two tubes of the CPPRs are attached to the upper two carriages of an actuation unit that translates the tubes (see Figure 11). The carriages are translated with EC-max22 brushless 12W Maxon motors via a PTGSG-10X2-01-R IGUS 2mm pitch lead screw. To ensure position and orientation precision, the carriages are guided with two µm precision linear guides. The motors are equipped with MR 2048 CPT Maxon encoders and controlled with a DMC-4040 Galil control board. In the experiments presented in this paper, tube 1 remained clamped in the baseplate at all times and only one stage of the actuation unit was used to translate tube 2 vertically.

The setup is equipped with two 6 degrees of freedom electromagnetic (EM) sensors (Aurora, NDI) in order to

<figure><img src="image_12.png" /><figcaption>Figure 11. Multi-stage CPPR actuation unit with a mounted prototype used for experimental validation. In the experiments presented in this paper, only one stage was used to translate one CPPR tube with respect to the other. (1) CPPR; (2) tip electromagnetic sensor; (3) camera; (4) nylon thread; (5) idler deflection pulleys; (6) reference electromagnetic sensor; (7) baseplate; (8) translation stages; (9) calibration weight.</figcaption></figure>

measure the tip position in the inertial frame. One reference 6 DoF sensor is attached at the basis of the actuation unit, while the other is attached to the tip of the prototypes (in a purposely designed slot on a rigid extension of the tip). The setup is also equipped with a camera which is used for qualitative and quantitative shape estimation. For each experiment with a planar CPPR, the camera intrinsic and extrinsic (3D pose of the camera w.r.t. the CPPR base frame) are calibrated. The images are then re-projected onto the robot bending plane (which is also the case for the images shown in Figures 12 to 14 and 16) and processed using standard image processing techniques in order to extract the metric robot shape. It is then compared to the model-predicted robot shape to compute a shape RMS error.

All the tubes are 3D printed with PLA (polylactic acid) using a Raise3D Pro2 printer in exactly the same position and orientation on the printing bed to limit printing variability. The height $k$ of all notches and the distance $c$ between all notches is 5 mm (see Figure 3).

## 4.2 Experimental Protocol

A CPPR is said to be *loaded* when it is subjected to a non-zero external tip force $f_+ \neq 0$, and *unloaded* otherwise,
<table><caption>**Table 1.** Geometrical and mechanical parameters of the prototypes.</caption><thead><tr><th rowspan="2">Robot design</th><th>ρ̄<sub>1</sub></th><th>ρ̄<sub>2</sub></th><th>ρ̄<sub>3</sub></th><th>ρ̄<sub>4</sub></th><th>k</th><th>c</th><th>l</th><th>γ<sub>1</sub></th><th>γ<sub>2</sub></th><th>f<sub>D</sub></th><th>E</th><th>Δa<sub>max</sub></th></tr><tr><th>[mm]</th><th>[mm]</th><th>[mm]</th><th>[mm]</th><th>[mm]</th><th>[mm]</th><th>[mm]</th><th>[mm]</th><th>[mm]</th><th>[-]</th><th>[GPa]</th><th>[mm]</th></tr></thead><tbody><tr><td>R<sub>1</sub> (regular)</td><td>3.0</td><td>2.0</td><td>1.5</td><td>0.5</td><td>5</td><td>5</td><td>150</td><td>5.0</td><td>2.0</td><td>1.27</td><td>5.39</td><td>7</td></tr><tr><td>R<sub>2</sub> (regular)</td><td>5.5</td><td>4.5</td><td>4.0</td><td>3.0</td><td>5</td><td>5</td><td>150</td><td>10.0</td><td>7.0</td><td>1.27</td><td>5.39</td><td>7</td></tr><tr><td>R<sub>3</sub> (convergent)</td><td>4.5</td><td>3.7</td><td>3.2</td><td>2.4</td><td>5</td><td>5</td><td>150</td><td>8.20 - 7.20*</td><td>5.60 - 4.95*</td><td>1.32</td><td>5.39</td><td>11</td></tr><tr><td>R<sub>4</sub> (divergent)</td><td>4.5</td><td>3.7</td><td>3.2</td><td>2.4</td><td>5</td><td>5</td><td>150</td><td>7.20 - 8.20*</td><td>4.95 - 5.60*</td><td>1.06</td><td>5.39</td><td>11</td></tr></tbody></table>

i.e. $f_+ = 0$. In the unloaded experiments, all prototypes defined in Table 1 (see. Sections 4.5 and 4.6 for more details) were actuated by alternating incremental 1 mm steps $\Delta a \in \{0, \pm 1 \text{ mm}, \pm 2 \text{ mm}, \dots, \pm \Delta a_{\text{max}}\}$. $\Delta a_{\text{max}}$ was determined once for, respectively, the regularly and variably notched CPPRs such as to cover a representative workspace, i.e. bending angle of approximately 90°. In the case of the regularly notched CPPRs ($R_1$ and $R_2$ detailed in Section 4.5), $\Delta a_{\text{max,reg}} = 7 \text{ mm}$ was determined with the CPPR of smallest diameter. In the case of the variably notched CPPRs $R_3$ and $R_4$ detailed in Section 4.6, as discussed in Section 5, we observed that the highest actuation values of the experiments translated to an internal load in simulations that goes beyond the theoretical estimation of critical buckling load $T_{\text{crit}}$. For this reason, validation of the model is limited to actuation values up to $\Delta a_{\text{max, var}} = 11 \text{ mm}$. In all tests, while alternating positive and negative actuation, the resting configuration (home position, $\Delta a = 0$) is systematically measured and compared to its initial reference value in order to detect potential hysteresis (due to plastic deformations or dry frictions), that might occur during bending. Extension 6 proposes a visual representation of the experimental protocol.

The loaded experiments were carried out over 5 fixed actuation values ($\Delta a \in \{0, \pm 3, \pm 6\} \text{ mm}$). Calibration weights of 20 g were attached to the tip and routed through idler deflection pulleys located symmetrically in the bending plane at both sides of the robot, as shown in Figure 11.

All tests (each robot design, actuation value, and loading condition) were repeated three times with duplicate prototypes. A total of 276 measurements (600 if counting the zero measurements) were performed for experimental validation. The 95% confidence interval for the tip measurement repeatability across the three repetitions with identical but distinct prototypes is found in the range [0, 2.36] mm, with a median at 0.6 mm, suggesting a very good repeatability between experiments with the same conditions. The 95% confidence interval for accuracy of the Aurora measurement apparatus is specified by the manufacturer to be [0, 0.88] mm. Two outliers were however found in one set of experiments compared to the two other ones of the same case (design $R_1$ with $f_+ = 0$, and $\Delta a \in \{-7, -6\} \text{ mm}$). After analysis, it was found that a motor stall was the source of the problem, and these two data points were removed from the subsequent data analysis. In the following results, for each test, the mean tip position across repetitions is used.

## 4.3 Use of the Simulator for Comparisons

For the purpose of comparisons with experiments, the simulator of Section 2.9 is used as follows. Firstly, note that the prototypes are actuated with $\Delta a$ and not with $T_+$. To address this difference, we use a linear PI controller that allows any prescribed $\Delta a_d$ to be first converted into a tension $T_+ = k_p e_{\Delta a} + k_I \int e_{\Delta a} dt$, with $e_{\Delta a} = \Delta a_d - \Delta a$, a positional error, and $k_p, k_I$, some proportional and integral gains. Once $e_{\Delta a} = 0$ (in practice $e_{\Delta a} \approx 10^{-6}$), $T_+ = k_I \int e_{\Delta a} dt$ defines the simulation actuation input of the current trial. It should be noted that this approach exploits the most general form of the kinetostatic model, where actuation is by force, while allowing the benefits of actuation by displacement to be exploited, which eliminates the influence of friction in particular. Secondly, the Ritz reduction method of Boyer et al. (2021) is here applied to $K_{1Z}$ with the first 5 Legendre polynomials, which proved to be sufficient over all cases for the solution to have converged. Finally, the weight of tubes has been observed to have a limited impact on the deformation and $\bar{F}_{\text{ext},i} = 0_{6 \times 1}$ everywhere in the algorithm.

## 4.4 Model Calibration

3D printing allows fabricating a variety of CPPR designs. There is, however, inherent variability due to the printing process used, which may impact the tube diameters, notch depth, etc. It is therefore necessary to calibrate the model before comparing simulation results with experiments. A quick analysis of the above model shows that the behavior of a CPPR is essentially governed by the balance between the elastic restoring moments, and those produced by the actuation push and pull force $T_+$ and external tip load $f_+$. Based on this observation, we chose to calibrate two model parameters: the distance function between the centroid lines of (17), and the Young modulus $E$ of the flexural stiffness $EI_{iZ}$. The offset $D$ was calibrated by seeking the value of a unique variable $f_D$ that minimizes the sum of the squares of the tip errors between experiments and simulations performed with $D_{\text{sim}} = f_D \times D$. Formally, $f_D$ is a constant scaling coefficient of $D$, solution of

$$f_D = \arg \min_{f_D} \sum_{k,j} \|r_{l,\text{sim}}(f_D, \mathcal{R}_k, \Delta a_j) - r_{l,\text{exp}}(\mathcal{R}_k, \Delta a_j)\|^2 . \quad (93)$$

The sum over the j-index means that this optimization is performed for all the $\Delta a_j$ in the actuation spaces described in Section 4.2. Regarding the other sum index k, the

<footnote>* Values of the proximal (γ₀) and distal (γ₁) depth of the notches varying linearly with X along the CPPRs (cf. Section 4.6).</footnote>
optimization problem is solved once for the two regularly notched prototypes ($R_1$ and $R_2$) and once for each of the two variable notch models ($R_3$ and $R_4$).

**Remark 17:** Note that the bending stiffness of the tubes depend on the centroid offset $D$ via Huygens theorem (see (3)). However, this stiffness being calibrated in a second stage through the Young modulus $E$, this dependency is ignored in this first calibration and absorbed by the second ●

**Remark 18:** When tube 2 is pushed ($\Delta a < 0$), the abrupt step in the inertia of the inner tube at the level of baseplate $\tilde{X}_2(0) = l_2 - a$ (see Figure 7) was found to slow down the convergence of the root-finding algorithms. Preliminary simulations allowed to identify a smaller value for $I_{2Z}(X_2 > l)$ that is yet sufficiently large to have a negligible impact on the results ●

Once the offset function $D$ rescaled, the Young modulus $E$ is calibrated in a similar way, but now with the loaded experiments (i.e. $f_+ \neq 0$), by solving the further problem

$$E = \arg \min_{E} \sum_{k,j,f_+} \|r_{l,\text{sim}}(E, \mathcal{R}_k, \Delta a_j, f_+) - r_{l,\text{exp}}(\mathcal{R}_k, \Delta a_j, f_+)\|^2, \quad (94)$$

where $k \in \{1, 2\}$, and $f_+$ stands for the tip loads experimented in the test of Section 4.7. The calibrated values of $f_D$ and $E$ deduced from this two stage process are given in Table 1. Note that the value of $E$ is identified with $\mathcal{R}_1$ and $\mathcal{R}_2$, and then used for all robots. It falls within the typical range for PLA (i.e. [2.7, 16] GPa).

**Remark 19:** Since the CPPRs are actuated by imposing $\Delta a$, and not $T_+$ (see Section 4.3), a variation of $E$ modifies $T_+$, but not their shape in the unloaded cases. In other words, the unloaded cases cannot be used to identify $E$. This explains why $E$ must be identified with several non-zero values of $f_+$ ●

## 4.5 Regularly Notched CPPRs

In these tests, two different regularly notched CPPR designs ($\mathcal{R}_1$ and $\mathcal{R}_2$) with outer diameters of respectively 6 and 11 mm, and centroid offsets such that $D_1 < D_2$, were studied. The geometrical parameters are given in Table 1, referring to Figure 3 for the notations. Reconsidering the face-to-face function (20) with $X_1 = 0$, and remarking that $\tilde{X}_2(0) = l_2 - a \simeq l - a$, shows after some algebraic manipulations, that we have $|K_{1Z}D| \simeq |\Delta a/l|$. This means that for a fixed $\Delta a$, the curvature of a CPPR increases as $D$ decreases. This tendency is confirmed by the fact that $\mathcal{R}_1$ has a greater range of curvature than $\mathcal{R}_2$ when each is subjected to the same actuation $\Delta a$. We used here the 14 actuation values $\Delta a \in \{0, \pm 1, \pm 2, \dots, \pm 7\}$ mm, which yields $2 \times 14 = 28$ simulations for $3 \times 28 = 84$ experiments. The tip errors between the simulation and experimental results for both regularly notched robots are reported in Figure 12. As shown by the plots of Figure 12(bottom), the tip errors do not exceed 3.4% of the total length. The shape RMS errors have an average of 1.7 mm (1.1% of the robot length) for

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Top, snapshots of the experiments and simulation results for the regularly notched prototypes, here with prototype $\mathcal{R}_1$. The colored lines represent for each case the centroid line of tube 1 and 2, respectively, with corresponding colors. The blue marks are the tip EM sensor data. The purple marks are the simulated robot tips. Bottom, tip error as function of the actuation $\Delta a$ for both regularly notched prototypes. See also Extension 6.</figcaption></figure>

both $\mathcal{R}_1$ and $\mathcal{R}_2$, with a maximum value of 6.8 mm (4.5%) and 2.4 mm (1.6%), respectively. This confirms the good agreement of our model with the experimental robot shape. These good results are confirmed by Figure 12(top), which displays some snapshots obtained with the robot $\mathcal{R}_1$, over its entire workspace, both in experiments and simulation. It should be noted that, as in the rest of the article, these snapshots are representative of all the others, which have not been reported here for the sake of visualization.

## 4.6 CPPRs with Variable Centroid Lines

With the design shear of (8), the model should be able to reproduce the deformations of a CPPR when the centroid lines of the tubes are not parallel. To validate this capability of the model, we constructed two additional prototypes, called $\mathcal{R}_3$ and $\mathcal{R}_4$, with convergent and divergent centroid lines respectively, as illustrated in Figure 6.

In such CPPRs, the depth of the notches decreases or increases from the base to the tip of the robot, in the resting configuration, according to the linear function of $X$ along $e_X$ (see Figure 2(c))

$$\gamma(X) = \gamma_0 + \left( \frac{\gamma_l - \gamma_0}{l} \right) X, \quad (95)$$

where $\gamma_0$ and $\gamma_l$ are, respectively, the proximal and distal notch depths given by the value pairs for $\gamma$ in Table 1. Remind that the interpolation (1) provides some $\alpha_i$ such that $\alpha_i(X) = \gamma_i(X) + \bar{\rho}_i$ with $\gamma_i$ defined by
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Top, snapshots of the experiments and simulation results for the convergent centroid line prototypes ($R_3$). The snapshots are split (odd/even translation actuation values) over two plots to prevent overlapping. The colored lines represent for each case the centroid line of tube 1 and 2, respectively, with corresponding colors. The blue marks are the tip EM sensor data. The purple marks are the simulated robot tips. Bottom, tip error as function of the actuation $\Delta a$. See also Extension 6.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Top, snapshots of the experiments and simulation results for the divergent centroid line prototypes ($R_4$). The snapshots are split (odd/even translation actuation values) over two plots to prevent overlapping. The colored lines represent for each case the centroid line of tube 1 and 2, respectively, with corresponding colors. The blue marks are the tip EM sensor data. The purple marks are the simulated robot tips. Bottom, tip error as function of the actuation $\Delta a$. See also Extension 6.</figcaption></figure>

(95). The two robots are subjected to the 11 actuation values $\Delta a \in \{0, \pm 1, \pm 2, \dots, \pm 11\}$ mm, that is to say a total of 44 simulations, for $3 \times 44 = 132$ experiments. The comparisons between the experimental and simulation results are presented in Figures 13 and 14 for prototypes $R_3$ and $R_4$ respectively. As shown by the plots of Figure 13(bottom) and Figure 14(bottom), the tip errors do not exceed 2.3% and 6.5% of the total length, respectively. Once again, examination of the shape of the robots as shown in the few representative snapshots of Figure 13(top) and Figure 14(top) confirms this trend. Such results are also confirmed by the shape RMS errors, which have an average value of 1.9 mm (1.2% of the robot length) and 1.8 mm (1.2%) for $R_3$ and $R_4$, respectively. The maximum shape RMS errors are respectively 5.0 mm (3.4%) and 7.7 mm (5.1%).

## 4.7 Externally Loaded CPPRs

This last series of experiments aims at assessing the ability of the model to capture the effect of an external load $f_+$ applied at the tip of a CPPR. To this end, as described in Section 4.2, a nylon thread is attached at the tip of the outer tube and routed around the left (respectively right) return pulley to a calibrated weight of mass $m$. The thread thus transmits a force of known intensity to the robot applied on its left (respectively right) side. Such a tip load $f_+$ is applied to the two robots $R_1$ and $R_2$, initially deformed under the effect of an imposed $\Delta a$. Although the intensity of $f_+$ is fixed once for all by taking $m = 20$ g, its direction is a priori unknown and depends on the interaction of the combined effects of $\Delta a$ and $f_+$ on the CPPR structure. To feed the simulator, $f_+$ is calculated from the experimental data as follows. It is assumed that $f_+$ is applied at the point of attachment of the nylon thread with the outer tube of the CPPR. Using
<figure><img src="image_16.png" /><figcaption>Figure 15. Schematic of the pulley system for the loading tests.</figcaption></figure>

notations of Section 2.1 and the context of Figure 15, the inertial position of this point is

$$r_{f+} = r_1(l_1) + (\bar{\rho}_1 \mp D_1(l_1)) b_{1Y}(l_1), \quad (96)$$

where $b_{1Y} = R_1 E_Y$, while the -/+ signs indicate whether the tip force is applied to the left (-) or right (+) of the robot. An examination of Figure 15 shows also that $f_+$ can be expressed as

$$f_+ = m a_g (\sin(\vartheta), \pm \cos(\vartheta), 0)^T, \quad (97)$$

where again the +/- signs indicate a left/right tip force, while $\vartheta$ is deduced from trigonometric considerations on triangles

$$\vartheta = \cos^{-1} \left( \frac{e_Y^T (r_{ply} - r_{f+})}{\|r_{ply} - r_{f+}\|} \right) + \sin^{-1} \left( \frac{\rho_{ply}}{\|r_{ply} - r_{f+}\|} \right), \quad (98)$$

with $r_{ply}$ and $\rho_{ply}$ the known inertial position and radius of the pulleys. Once $f_+$ has been computed in this way, it is introduced into the wrench $F_+$ of (25), that feds the backward pass of the algorithm of Section 2.9.1, when the stress wrench is transferred from tube 2 to tube 1.

For both robots $R_1$ and $R_2$, $f_+$ is applied successively on the left and on the right, with $\Delta a \in \{0, \pm 3, \pm 6\}$ mm, which gives a total of $2 \times 2 \times 5 = 20$ comparisons between the simulations and the experiments, and thus $3 \times 20 = 60$ experimental trials. The tip errors between the simulation and experimental results for all loaded experiments are reported in Figure 16. As shown in Figure 16(bottom), these errors do not exceed 6.3% of the total CPPR length. Once again, shape RMS errors confirm these results, with an average of 2.4 mm (1.6% of the robot length) and a maximum of 5.9 mm (4.0%) for those experiments.

Figure 16(top and middle) show the 10 snapshots obtained with $f_+$ applied on the right (top) and left (middle) of robot $R_1$ for a given repetition. By comparing, two by two and from left to right, the shaded shapes of the background with those of the foreground, one can see how the prototypes are deviated when forces are applied.

## 4.8 A CPPR Loaded out of its Bending Plane

Unlike the previous validation, which concerned the planar case, in this and the following (see also Extension 7), the

<figure><img src="image_17.png" /><figcaption>Figure 16. Top, snapshots of the experiments and simulation results for loaded tests with the regularly notched prototype $R_1$. The snapshots are split (right/left tip forces, identified by the blue arrow) over two plots to prevent overlapping. The shaded shapes on the background correspond to the free-space shapes at identical actuation (match left to right by pairs). The colored lines represent for each case the centroid line of tube 1 and 2, respectively, with corresponding colors. The blue marks are the tip EM sensor data. The purple marks are the simulated robot tips. Bottom, tip error as function of the actuation $\Delta a$ for all loaded tests with both regularly notched prototypes. See also Extension 6.</figcaption></figure>

performance of the 3D model is evaluated for CPPRs with 3D deformations, under different conditions and designs. In all these tests, the curvature fields are approximated with the first 3 Legendre polynomials, and the torsion fields, with the first two Legendre polynomials after they have been adapted to ensure (87). In this section we start by a test taken from (Childs and Rucker 2023) where a 2-tube CPPR 33 mm long is first bent in the horizontal plane with a push and pull force, then subjected to a weight of mass $m = 50$g attached to the end-effector. Note that in (Childs and Rucker 2023), the CPPRs have much smaller notches than our 3D printed mock-ups. This scenario can be simulated with a 2-tube version of the full 3D algorithm of Section 3, where the vertical load is modeled by a wrench $F_+ = (0_{1 \times 3}, N_+^T)^T$,
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Boxplots of the tip position error between experimental results and simulations for all tests. The first three groups are unloaded cases ($f_+ = 0$). The ‖ symbol is related to the number of experiments.</figcaption></figure>

with $N_+ = 9.87 m R_1(l_1)^T e_Z$. Based on the results of the FEM analysis of Childs and Rucker (2023), the 3D model is fed with the following flexural and torsional stiffness. For tube 1, ($EI_{1Y}, EI_{1Z}, GI_{1X}$) = (17.83, 1.25, 1.08) Ncm², and $D_1$ = 1.09 mm; for tube 2: ($EI_{2Y}, EI_{2Z}, GI_{2X}$) = (20.7, 1.85, 1.71) Ncm², $D_2$ = 0.94 mm. Figures 18 and 19 show top, frontal, and perspective views of the centerline of the loaded CPPR for $\Delta a_d$ varying from -1 to 1 mm every 0.25 mm. As shown by the frontal view, such an out-of-plane load generates a slight 3D deformation.

<figure><img src="image_19.png" /><figcaption>**Figure 18.** 2-tube CPPR loaded out of the plane of bending. Left: Top view, in the X-Y (horizontal) plane. Right: front view in the Y-Z (vertical) plane.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** (top left) X-Y view, (top right) Y-Z view and (bottom) perspective view of the 3-tube CPPR of (Childs and Rucker 2023)</figcaption></figure>

## 4.9 A CPPR with 3 Regularly Notched Tubes

This second test is based on (Childs and Rucker 2023). It consists of a 30 mm-length CPPR with 3 tubes regularly notched and symmetrically arranged every $2\pi/3$ around the central axis. In these conditions, the tube centroid

lines are straight and parallel to the central axis, and pushing and pulling each of the tubes generates a planar bending in any direction. The full 3D model of Section 3 is here used. As in (Childs and Rucker 2023), tube 1 is clamped into the base and stands for the leader tube, while the two others are followers that are pushed and pulled at the level of the baseplate. Based on the results of the FEM analysis of Childs and Rucker (2023), the 3D model is fed with the following flexural and torsional stiffness. For tube 1, ($EI_{1Y}, EI_{1Z}, GI_{1X}$) = (1.96, 1.78, 2.73) Ncm², and $D_1$ = 1.12 mm; for tube 2: ($EI_{2Y}, EI_{2Z}, GI_{2X}$) = (1.81, 1.78, 2.78) Ncm², $D_2$ = 0.983 mm, and for tube 3: ($EI_{3Y}, EI_{3Z}, GI_{3X}$) = (1.83, 1.89, 2.91) Ncm² and $D_3$ = 0.764 mm. Moreover, the 3D model of Section 3 also requires the shear stiffness of the two followers. Practically, we here used ($G_1A_1, G_2A_2, G_3A_3$) = (1.13, 1.45, 1.23)$10^4$ Ncm² deduced from ($G_1, G_2, G_3$) = (21.15, 27.61, 25.74) GPa and the dimensions of the tube cross-sections (ignoring the notches) given in (Childs and Rucker 2023). While tube 1 is maintained fixed, the two others are subject to 3 sequences of translation steps of 0.5 mm each, such that ($\Delta a_{d2}, \Delta a_{d3}$) ∈ [-1, 1] × {0} (mm) (sequence 1), ($\Delta a_{d2}, \Delta a_{d3}$) ∈ {0} × [-1, 1] (mm) (sequence 2), ($\Delta a_{d2}, \Delta a_{d3}$) ∈ [-1, 1]² (mm) (sequence 3). In Figure 20, the configurations of the CPPR's centerline are displayed in a lateral plane (a), a top plane (b) and a perspective view (c), under these actuation conditions. The parallel straight lines of the centroids of the tubes in this CPPR allow the fields of curvature of the tubes to be decoupled. This is illustrated in Figure 20(top-right) where it can be seen that the 3 sequences of tube translations move the CPPR in a pure bending deformation along the 3 almost straight branches of a star. Although difficult to quantify, the perspective view of Figure 20(bottom) is close to that of (Childs and Rucker 2023).

## 4.10 A CPPR Actuated with a Torque

This test, like the two following illustrates some of the capabilities of the model that go beyond the current state of the art. To this end, we consider a 2-tube CPPR with parallel opposite centroid lines. The geometric properties and Young modulus of this CPPR are those of prototype $R_1$ (see Table 1). The shearing modulus is $G = 1$ GPa for the two tubes. The follower tube is subject to a torque that deforms the CPPR in torsion only. Like for translation, a PI controller is implemented to ensure $\Delta\theta_2(0) = \Delta\theta_d$, with $\Delta\theta_d$ a prescribed angle variation, here such that $\Delta\theta_d \in \{-\pi, -\pi/2, 0, \pi/2, \pi\}$ rad. In Figure 21, the leader tube is in green, the follower is in red, and can rotate at the level of the baseplate along the red circle. The black solid line indicates the material line linking the tip of 2 tubes at the level of the distal platform. The perspective view of Figure 21 shows how the centroid lines are twisted around their common centerline. The follower tube is the most twisted because of its lower torsional stiffness. As expected, the top view shows each of the centroid lines spiralling around a cylinder centered on the centerline of the CPPR.
<figure><img src="image_21.png" /><figcaption>**Figure 20.** (top left) X-Z (lateral) view, (top right) Y-Z (top) view, and (bottom) perspective view of the 3-tube CPPR of (Childs and Rucker 2023).</figcaption></figure>

base with prescribed translations $\Delta a_d$ varying from -10 mm to 10 mm by a step of 1 mm.

<figure><img src="image_22.png" /><figcaption>**Figure 22.** 2-tube CPPR with helicoidal centroid lines actuated in translation. (left) Configuration at rest. (right) Snapshots of a deformation due to a translation of tube 2 from -10 mm to 10 mm. See also Extension 7.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 21.** 2-tube CPPR actuated in rotation only. Perspective view (left). Top view (right). The red arc of circle represents the path taken by the centroid of tube 2 at the base of the CPPR, while the green arc of circle represents the path taken by the centroid of tube 1 at the top of the CPPR. The solid black line represents the rigid segment connecting the centroids of the two tubes at the tip of the CPPR. Their intersection lies on the central axis of the CPPR. See also Extension 7.</figcaption></figure>

## 4.11 CPPRs with Helicoidal Centroid Lines

In this test, we consider the same 2-tube CPPR as in the previous test, except that the centroid line of the leader at rest now has a non-zero torsional offset $K_{1X}^o = 2\pi/l_1$. This forces the two opposing centroid lines to wrap helically around the centerline of the CPPR at rest, as shown on Figure 22(left). Figure 22(right) shows the snapshots of a simulation where the follower tube is pushed and pulled at its

## 4.12 A CPPR Actuated with Push-Pull Forces and Torques

In this test, we carried out an initial series of comparisons between the 3D simulator and the experiment. To do this, the 2-tube prototype $R_1$ is used as in the two previous sections. However, the shear moduli of the two tubes are considered independent of each other and calibrated with respect to the experiments in the pure torsion test of Figure 23. Note that this rough estimate leads to $(G_1I_{1X})/(G_2I_{2X}) \approx 0.65$ (an absolute estimate of each stiffness would require a measurement of the actuating torque), i.e. a leader tube twice as torsionally flexible as the follower tube. Based on this calibration, the comparison involves a scenario in which the follower tube is first pushed and pulled to its base with prescribed translations $\Delta a_d = \pm 3$ mm, then rotated at the base with $\Delta\theta_d \in \{-\pi/2, 0, \pi/2\}$ rad. As illustrated in Figure 24, these loading conditions couple bending and torsion to generate a slightly spiraled shape in 3D and a finite tip deflection. Such actuation capability is validated experimentally using end-effector pose measurements, similar to the planar cases. Additionally, we performed shape measurements using a mechanical digitizer FAROArm (FARO Technologies, Florida, USA) equipped with a 3D laser scanner. Afterwards, the 3D point clouds are processed and registered in the robot's base frame. In Figure 24, these experimental data are represented by the tip frames, and two scanned configurations (one at rest, the second deformed).
<figure><img src="image_24.png" /><figcaption>**Figure 23.** Top-view of the 2-tube CPPR ($R_1$) actuated in rotation only with $\Delta\theta_d \in \{-\pi, -\pi/2, 0, \pi/2, \pi\}$ rad. The blue arrows represent one of the base vectors of the frame of reference attached to the tip of the CPPR and measured by the Aurora sensing device. The values of $G_1$ and $G_2$ are calibrated and the corresponding vectors of the model (black solid lines) are shown. See also Extension 7.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 24.** 3D shape measurements of a 2-tube CPPR (prototype $R_1$) actuated in translation and rotation. The frames at the tip of the snapshots are defined by pose measurements. The thick red and green lines are the simulated centroid lines of the tubes. See also Extension 7.</figcaption></figure>

## 5 Discussion and Conclusion

In this paper, a mechanics-based model for CPPRs based on Cosserat rod theory is proposed as an alternative to (Childs and Rucker 2023). Extensive experimental validation was first carried out on four 3D-printed 2-tube CPPR models (different dimensions and notch geometries), with external forces applied to the tip of two of them. Consistent with (Childs and Rucker 2023) and previous studies applied to other continuum robots (Rucker and Webster 2011; Rucker et al. 2010; Black et al. 2018; Orekhov et al. 2017) the results obtained show that the model performs well, with mean and median tip errors of 1.47% and 1.06% of the robot length, respectively, over the set of experiments. Figure 17 provides a graphical summary of the tip error statistics and number of experiments for each category. Overall, in 94.6% of the experiments, the tip error was below 5% of the robot length, including experiments under external loading. In the planar cases, the RMS shape errors show a similar trend, with an average error of 1.9 mm (1.2% of the robot length), and all experiments yielding an error below 5.1% of the robot length. 92% of the experiments yielded a shape RMS error below 3% of the robot length.

As in the case in most of continuum robots, model errors tend to increase further from the base, and are generally at their maximum at the end-effector. The same applies to configurations at the boundary of the workspace of such robots, where the actuation (should it be the translation or rotation of rods/tendons/tubes) amplitude increases. In fact, adding to fabrication imperfections, backlash, interplay, and other assembly artifacts, the unmodeled phenomena, such as friction (and tendon elongation in TACR) tend to reveal when soliciting such deformable systems. Recalling most of continuum robots are fixed and actuated at their base, such effects accumulate towards the end-effector.

Systematic comparison of the resting configuration of the CPPRs with that obtained by relaxing $\Delta a$ to zero after each measured deformation, showed us that the CPPRs endured limited hysteresis (e.g. plastic deformation), with tip displacements lower than 2 mm in all cases.

Cases with actuation values $\Delta a > 11$ mm were discarded from the results as mentioned in Section 4.2. Indeed, the application of a pushing force at the tip of the inner tube is countered by a reaction force at its opposite end, from the outer tube. This reaction force is opposed to the applied pushing force and loads the inner tube into compression. Conversely, the application of a pulling force on the inner tube exerts a pair of opposing compressing forces at the ends of the outer tube. Therefore, due to their operating principle, CPPRs can buckle beyond a certain load which can be approximated by calculating Euler's critical load for a pin ended rod

$$T_{\text{crit}} \approx \frac{\pi^2 EI_Z}{l^2}, \qquad (99)$$

where $I_Z = \max(I_{1Z} + I_{2Z})$. Applying this approximation to the robots $R_3$ and $R_4$, the critical load is plotted in Figure 25 (dotted line), with the simulated values of the push and pull force $T_+$ as a function of $\Delta a$. Based on this analysis, note first that the model accuracy decreases when approaching the theoretical critical load, i.e. when $\Delta a$ is close to $\pm 9$ mm for the convergent $R_3$ and divergent $R_4$ designs. Second, above this critical value, we observed that
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Axial load in the variable centroid line prototypes. The absolute value of the load increases with increasing actuation translation. Beyond $|Δa| > 11$ mm the axial load was found to exceed Euler's critical load of the rods composing the CPPRs.</figcaption></figure>

the equilibrium configurations calculated with the model, jump abruptly to another family of solutions. Furthermore, in some cases, the root-finding algorithm may not converge or converge very slowly. This explains why these pathological cases cannot be included in the simulation analysis. It should be noted that this buckling phenomenon is not an intrinsic limitation of CPPRs, but is rather due to our experimental conditions, especially those related to external loading.

Although encouraging, the above results call for further investigation. The proposed model has the potential of addressing the full 3D deformations of a CPPR subject to any combination of bending and torsion with arbitrary variable notch designs (possibly spiraled around its central axis). However, over and above its theoretical consistency and the encouraging preliminary results obtained in simulation, the approach will need to be validated by experiments under these extended conditions and further extended to other possibilities for these systems, such as S-shaped CPPRs, for example, or multi-segment CPPRs. A few physical phenomenons remain un-modeled in our approach and their inclusion could further improve the accuracy of the model. The first phenomenon is the slight non-tangency between the tubes. As with many models of robots with concentric tubes (Dupont et al. 2010; Rucker and Webster 2008), our model assumes that the free space between the tubes is ideally sufficient to allow their relative motion without friction, while keeping their tangents coincident. As shown by Ha et al. (2019), the free space modeling allows incorporating the friction between the tubes, which improves the accuracy of the model. However, the integration of these phenomenons into a Lagrangian formulation based on the Cosserat rod model in the context of CPPRs remains an open problem.

Another element worth discussing is the model calibration. 3D printing enables easy manufacturing of CPPRs with various designs, as evidenced in the experimental validation with e.g. convergent centroid line designs. 3D printing is however subject to imperfections, which can lead to large errors in key dimensions. For example, in simulation, increasing the notch depth by 100 µm for design $R_1$ leads to a tip error up to 4.1 mm in the workspace with respect to the nominal model results. In order to compensate for such effects, we chose to perform a calibration of the parameter

$D$ (see Section 4.4). However, it is worth mentioning that by using more accurate measuring equipment and more sophisticated manufacturing (laser notching of Nitinol tubes, as in York et al. (2015)), such a calibration procedure can be avoided as this has been done in the 3D simulations of Sections 4.8 and 4.9, which are based on the much more mature fabrication of Childs and Rucker (2023). Added to such aspects, further CPPR development would benefit from other factors introduced in previous research, such as the critical load in (Oliver-Butler et al. 2022), maximum curvature and notch closure conditions introduced in (Childs and Rucker 2023). These can be essential factors for designers and modelers in order to comply to the target application requirements and enhance integration possibilities.

Encouragingly, past years have seen the technical development and assessment of CPPRs in medical applications. The concept was first introduced by Oliver-Butler et al. (2017) and patented by Riojas et al. (2019). In 2018, the concept of its integration through the working channel of a bronchoscope was shown, as a preliminary step towards medical applications (Rox et al. 2018). They proved the feasibility of central airway tumor resection during *ex vivo* sheep trials. Further, Qin et al. (2023) investigated the design of multi-section CPPRs integrated in a colonoscope for treatment of colorectal cancer. The research eventually led to the development of CPPRs within the company “EndoTheia, Inc, Nashville, TN, USA” as pre-clinical systems for the treatment of biliary disease, kidney stones, and middle ear disease. The mentioned advancements highlight the potential of CPPRs in translating into surgical robots. This article provides ground for addressing challenges related to their modeling within more realistic scenarios including complex 3D deformations.

Finally, the modeling and simulation approach developed in this article is the natural extension of our previous works on TACRs and continuum parallel robots, and the dynamics of continuum robots in general (Tummers et al. 2023; Briot and Boyer 2023; Boyer et al. 2021, 2024). Exploiting this common framework, the results of the present paper can be extended to the dynamics of CPPRs following (Boyer et al. 2021, 2024), or exploited for stability analysis, as it has been done in (Briot and Boyer 2023) for continuum parallel robots. In addition, the structure of CPPRs, which can be considered as tree-like systems of Cosserat rods organized in coaxial bundles, has necessitated the development of new concepts to meet the challenges of sliding rods with asymmetrical cross-sections, the reduction of coaxiality constraints, and the calculation, using stress/strain duality, of the generalized forces of a kinetostatic model in Lagrangian form. The workflow of this paper can thus be adapted to other structures with comparable requirements, such as MBCRs. All in all, the model proposed in this paper is a ground for further structural and functional analyses, application specific design optimization, intrinsic shape/force sensing, and motion/force control.

## Funding

This work was supported by grants ANR-11-LABX-0004-01, ANR-19-P3IA-0003, ANR-20-CE33-0001, ANR-10-IAHU-0002, ANR-18-CE19-0012, and ANR-21-ESRE-0015.
## Notes

1. In this paper we consider rectangular notches, however, the derivation of the proposed model is independent of this design consideration. Moreover, a similar procedure to the one described in Appendix B can be applied to any kind of asymmetric notches. Such notch profiles can have any front profile shape as well as any cross-section shape, including square or hexagonal tubes as in (Oliver-Butler et al. 2017).

2. Any regular parameterization could be used here. However, being piecewise polynomial functions composed of several cubic polynomials, cubic splines provide a smooth approximation to the data, which is highly desirable in our case (Poirier 1973).

3. Note that only $\bar{F}_{c,2} \neq 0_{6 \times 1}$ below the baseplate, while inside the CPPR, we have $\bar{F}_{c,1} = -\bar{F}_{c,2}$ in virtue of the action-reaction principle.

4. Note that, throughout this paper, the difference in nature between stresses and forces is indicated by the notations $\Lambda$ and $F$.

## References

Antman SS (2005) *Nonlinear Problems of Elasticity*. 2nd edition. New York, NY, USA: Springer-Verlag.

Black CB, Till J and Rucker DC (2018) Parallel continuum robots: Modeling, analysis, and actuation-based force sensing. *IEEE Transactions on Robotics* 34(1): 29–47.

Boyer F, Gotelli A, Tempel P, Lebastard V, Renda F and Briot S (2024) Implicit time integration simulation of robots with rigid bodies and cosserat rods based on a newton-euler recursive algorithm. *IEEE Transactions on Robotics, In press*, (10.1109/TRO.2023.3334647). hal-04291835.

Boyer F, Lebastard V, Candelier F and Renda F (2021) Dynamics of continuum and soft robots: A strain parameterization based approach. *IEEE Transactions on Robotics* 37(3): 847–863.

Boyer F, Lebastard V, Candelier F and Renda F (2022) Extended hamilton’s principle applied to geometrically exact kirchhoff sliding rods. *Journal of Sound Vibration* 516: 116511. DOI: 10.1016/j.jsv.2021.116511.

Briot S and Boyer F (2023) A geometrically-exact assumed strain modes approach for the geometrico- and kinemato-static modellings of continuum parallel robots. *IEEE Transactions on Robotics* 39(2): 1527–1543. DOI:10.1109/TRO.2022.3219777.

Burgner-Kahrs J, Rucker DC and Choset H (2015) Continuum robots for medical applications: A survey. *IEEE Transactions on Robotics* 31(6): 1261–1280.

Chen Y, Wu B, Jin J and Xu K (2021) A variable curvature model for multi-backbone continuum robots to account for inter-segment coupling and external disturbance. *IEEE Robotics and Automation Letters* 6(2): 1590–1597.

Chikhaoui MT and Rosa B (2022) Chapter 8 - modeling and control strategies for flexible devices. In: *Endorobotics*. New York, NY, USA: Academic. ISBN 978-0-12-821750-4, pp. 187–213.

Childs JA and Rucker C (2023) A kinostatic model for concentric push-pull robots. *IEEE Transactions on Robotics*: 1–19DOI: 10.1109/TRO.2023.3327811.

Cosserat E and Cosserat F (1909) *Théorie des corps déformables*. Paris, France: Hermann & Cie.

Dupont PE, Lock J, Itkowitz B and Butler E (2010) Design and control of concentric-tube robots. *IEEE Transactions on Robotics* 26(2): 209–225.

Dupont PE, Simaan N, Choset H and Rucker DC (2022) Continuum robots for medical interventions. *Proceedings of the IEEE* 110(7): 847–870.

Ha J, Fagogenis G and Dupont PE (2019) Modeling tube clearance and bounding the effect of friction in concentric tube robot kinematics. *IEEE Transactions on Robotics* 35(2): 353–370.

Li H, Xun L and Zheng G (2023) Piecewise linear strain cosserat model for soft slender manipulator. *IEEE Transactions on Robotics* 39(3): 2342–2359.

Oliver-Butler K, Childs JA, Daniel A and Rucker DC (2022) Concentric push-pull robots: Planar modeling and design. *IEEE Transactions on Robotics* 38(2): 1186–1200.

Oliver-Butler K, Epps ZH and Rucker DC (2017) Concentric agonist-antagonist robots for minimally invasive surgeries. In: *Proceedings of SPIE Medical Imaging*. p. Art. no. 1013511.

Orekhov AL, Aloi VA and Rucker DC (2017) Modeling parallel continuum robots with general intermediate constraints. In: *Proceedings of IEEE International Conference on Robotics and Automation*. pp. 6142–6149.

Peyron Q, Rabenorosoa K, Andreff N and Renaud P (2019) A numerical framework for the stability and cardinality analysis of concentric tube robots: Introduction and application to the follow-the-leader deployment. *Mechanism and Machine Theory* 132: 176–192.

Poirier DJ (1973) Piecewise regression using cubic splines. *Journal of the American Statistical Association* 68(343): 515–524.

Qin T, Connor P, Dang KT, Alterovitz R, Webster RJ and Rucker DC (2023) Computational analysis of design parameters for a bimanual concentric push-pull robot. In: *Proceedings of the Hamlyn Symposium on Medical Robotics*. pp. 19–20.

Renda F, Armanini C, Lebastard V, Candelier F and Boyer F (2020) A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation. *IEEE Robotics and Automation Letters* 5(3): 4006–4013.

Riojas KE, Webster RJ, Rucker DC, Butler KO and Ponten R (2019) Surgical device tip with deflectable joint. Patent: US 2019/0133705 A1.

Rone WS and Ben-Tzvi P (2012) Continuum manipulator statics based on the principle of virtual work. In: *ASME International Mechanical Engineering Congress and Exposition, volume Volume 4: Dynamics, Control and Uncertainty, Parts A and B*. pp. 321–328.

Rox M, Riojas KE, Emerson M, Oliver-Butler K, Rucker DC and Webster RJ (2018) Luminal robots small enough to fit through endoscope ports: Initial tumor resection experiments in the airways. In: *Proceedings of the Hamlyn Symposium on Medical Robotics*. pp. 63–64.

Rucker DC, Jones BA and Webster RJ (2010) A geometrically exact model for externally loaded concentric-tube continuum robots. *IEEE Transactions on Robotics* 26(5): 769–780.

Rucker DC and Webster RJ (2008) Mechanics-based modeling of bending and torsion in active cannulas. In: *Proceedings of IEEE/RAS-EMBS International Conference on Biomedical Robotics and Biomechatronics*. pp. 704–709.

Rucker DC and Webster RJ (2011) Statics and dynamics of continuum robots with general tendon routing and external loading. *IEEE Transactions on Robotics* 27(6): 1033–1044.
Santina CD and Rus D (2020) Control oriented modeling of soft robots: The polynomial curvature case. *IEEE Robotics and Automation Letters* 5(2): 290–298.

Simo JC and Vu-Quoc L (1988) On the dynamics in space of rods undergoing large motions - a geometrically exact approach. *Computer Methods in Applied Mechanics and Engineering* 66(2): 125–161.

Swaney PJ, York PA, Gilbert HB, Burgner-Kahrs J and Webster RJ (2016) Design, fabrication, and testing of a needle-sized wrist for surgical instruments. *Journal of Medical Devices* 11(1): 014501–1–xx.

Till J and Rucker DC (2017) Elastic stability of Cosserat rods and parallel continuum robots. *IEEE Transactions on Robotics* 33(3): 718–733.

Trefethen LN (2000) *Spectral methods in MATLAB*. Philadelphia, PA, USA: SIAM.

Tummers M, Lebastard V, Boyer F, Troccaz J, Rosa B and Chikhaoui MT (2023) Cosserat rod modeling of continuum robots from newtonian and lagrangian perspectives. *IEEE Transactions on Robotics* 39(3): 2360–2378.

York PA, Swaney PJ, Gilbert HB and Webster RJ (2015) A wrist for needle-sized surgical robots. In: *Proceedings of IEEE International Conference on Robotics and Automation*. pp. 1776–1781.

Young WC, Budynas RG and Sadegh AM (2012) *Roark's formulas for stress and strain*. McGraw-Hill Education.

Yuan H, Zhou L and Xu W (2019) A comprehensive static model of cable-driven multi-section continuum robots considering friction effect. *Mechanism and Machine Theory* 135: 130–149. DOI:10.1016/j.mechmachtheory.2019.02.005.

# Appendix

## A Index to Multimedia Extensions

<table><thead><tr><td>Extension</td><td>Media Type</td><td>Description</td></tr></thead><tbody><tr><td>1</td><td>Video</td><td>Working principle of CPPR</td></tr><tr><td>2</td><td>Video</td><td>Modeling the notched tubes as Cosserat rods</td></tr><tr><td>3</td><td>Video</td><td>Modeling the contacts</td></tr><tr><td>4</td><td>Video</td><td>Model reduction</td></tr><tr><td>5</td><td>Video</td><td>Solving the model</td></tr><tr><td>6</td><td>Video</td><td>Experimental validation</td></tr><tr><td>7</td><td>Video</td><td>CPPRs with 3D deformations</td></tr></tbody></table>

## B Cross-Sectional Inertia Parameters of Notched Tubes

As the expressions below are valid for both tubes, the subscript *i* is omitted for brevity. This appendix gives the expressions of the parameters (*A*, *D*, *J*<sub>Y</sub>, *J*<sub>Z</sub>)(*X*) of Section 2.4.2. To that end, we use a frame not centered on the centroid *G*(*X*), but on the center *C*(*X*) of a cross-section. Thus, (*Y*, *Z*) below designates the coordinates of a point on the cross-section *X*, expressed in (*C*, *b*<sub>X</sub>, *b*<sub>Y</sub>, *b*<sub>Z</sub>)(*X*). Remind that the basis (*b*<sub>Y</sub>, *b*<sub>Z</sub>) respects the symmetries of the U-shaped sections.

### B.1 Centroids and Cross-Sectional Area

Following the methodology of Swaney et al. (2016); York et al. (2015), we first compute the areas $\bar{A}$ and $\underline{A}$ of the pieces of full disk of radius $\bar{\rho}$ and $\underline{\rho}$ for $|Y| > \alpha$ (see Figure 3), as well as the offsets of their centroids $\bar{D}$ and $\underline{D}$, with respect to the CPPR centerline

$$ \bar{A} = \frac{\bar{\rho}^2(\bar{\phi} - \sin(\bar{\phi}))}{2}, \quad \underline{A} = \frac{\underline{\rho}^2(\underline{\phi} - \sin(\underline{\phi}))}{2}, \quad (100) $$

$$ \bar{D} = \frac{4\bar{\rho}\sin^3(\frac{1}{2}\bar{\phi})}{3(\bar{\phi} - \sin(\bar{\phi}))}, \quad \underline{D} = \frac{4\underline{\rho}\sin^3(\frac{1}{2}\underline{\phi})}{3(\underline{\phi} - \sin(\underline{\phi}))}, \quad (101) $$

where we introduced the notations

$$ \bar{\phi} = 2 \cos^{-1}(\alpha/\bar{\rho}), \quad \underline{\phi} = 2 \cos^{-1}(\alpha/\underline{\rho}), \quad (102) $$

with $\alpha = \gamma - \bar{\rho}$ and $\gamma$ the depth of the notches as shown in Figure 3. Then, using two such pieces of full disk, one with positive, the other with negative area, the cross-sectional area of the notched tube is given by their sum

$$ A = \iint dY dZ = \bar{A} - \underline{A}, \quad (103) $$

while the offset of its centroid *G* is

$$ D = \frac{1}{\bar{A}} \iint Y \, dY \, dZ = \frac{\bar{D} \bar{A} - \underline{D} \underline{A}}{\bar{A} - \underline{A}}. \quad (104) $$

In this implementation, extra attention must be paid to check whether the notches reach deeper than the opposing inner wall of the tube ($\gamma > \underline{\rho} + \bar{\rho}$). In this case, $\underline{A} = 0$ and the offset is simply $D = \bar{D}$.

### B.2 Cross-Sectional Inertia

Following a similar additive process of positive and negative area, provides the expression of the second inertia momenta

$$ J_X = J_Y + J_Z, \quad (105) $$

$$ J_Y = \iint Z^2 \, dY \, dZ = \frac{\bar{\rho}^4}{48} \Psi_Y(\bar{\rho}) - \frac{\rho^4}{48} \Psi_Y(\underline{\rho}), \quad (106) $$

$$ J_Z = \iint Y^2 \, dY \, dZ = \frac{\bar{\rho}^4}{16} \Psi_Z(\bar{\rho}) - \frac{\rho^4}{16} \Psi_Z(\underline{\rho}), \quad (107) $$

with

$$ \Psi_Y(\rho) = 6\pi - 12\psi - 8\sin(2\psi) - \sin(4\psi), \quad (108) $$

$$ \Psi_Z(\rho) = 2\pi - 4\psi + \sin(4\psi), \quad (109) $$

and the further notation $\psi(\rho) = \sin^{-1}(\alpha/\rho)$. Finally, it should be noted that thin-walled open sections can be subject to phenomena of combined deformation modes thus virtually affecting the neutral line position and the shear modulus (Young et al. 2012). For the prototypes simulated in this paper, it was verified that the used formulas remained within their valid range of operation.
C Derivation of the ODE (20) that Governs the Face-to-Face Function of a Planar 2-Tube CPPR

To derive the expression of the face to face function, we first split the coaxiality condition on poses given by (14) and (18), into its linear component

$$ \forall X_1 \in [0, l_1] : \tilde{r}_2(X_1) = r_1(X_1) - D(X_1)R_1(X_1)E_Y, \quad (110) $$

and its angular one

$$ \forall X_1 \in [0, l_1] : \tilde{R}_2(X_1) = R_1(X_1)R_Z(\pi). \quad (111) $$

Then, differentiating (110) with respect to $X_1$ and using (13) with $f_2 = r_2$, gives

$$ h \left( \frac{dr_2}{dX_2} \right)^{\sim} = \frac{dr_1}{dX_1} - R_1 \left( \frac{dD}{dX_1} E_Y + DK_1 \times E_Y \right). \quad (112) $$

Further, using (111) in (112) and rearranging yields

$$ h R_Z(\pi) \tilde{R}_2^T \left( \frac{dr_2}{dX_2} \right)^{\sim} = R_1^T \frac{dr_1}{dX_1} - \left( \frac{dD}{dX_1} E_Y - DK_1 Z E_X \right), \quad (113) $$

in which we recognize the translational component of (7)

$$ \Gamma_i = R_i^T \left( \frac{dr_i}{dX_i} \right). \quad (114) $$

Hence, one can rewrite (113) as

$$ h R_Z(\pi) \tilde{\Gamma}_2 = \Gamma_1 - \left( \frac{dD}{dX_1} E_Y - DK_1 Z E_X \right). \quad (115) $$

But, with the model of deformations of Section 2.5, $\Gamma_1$ and $\Gamma_2$ can be decomposed as

$$ \Gamma_i = (1, \beta_i(X_i), 0)^T = E_X + \beta_i(X_i) E_Y. \quad (116) $$

Therefore, introducing the detailed expressions (116) into (115) and projecting the result onto $E_X$ gives the expression of $h$ that once equalized with its definition (13) provides (20), where the minus sign comes from $R_Z(\pi)$.