# Angular Strain Parameterization for Solving Static Cosserat Rods

Radhouane Jilani, Pierre-Frédéric Villard, Erwan Kerrien

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

## To cite this version:

Radhouane Jilani, Pierre-Frédéric Villard, Erwan Kerrien. Angular Strain Parameterization for Solving Static Cosserat Rods. ICRA@40, IEEE, Sep 2024, Rotterdam, Netherlands. hal-04774573

HAL Id: hal-04774573

https://inria.hal.science/hal-04774573v1

Submitted on 8 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Angular Strain Parameterization for Solving Static Cosserat Rods

Radhouane Jilani*, Pierre-Frédéric Villard*, Erwan Kerrien*

*Université de Lorraine, CNRS, Inria, LORIA, F-54000 Nancy, France

**Abstract**— We show that static Cosserat rod clamped-free boundary value problems can be solved by only employing an angular strain parameterization. Consequently, the angular strain is discretized, while the linear strain remains continuous. This is made possible by leveraging the known boundary conditions at both ends of the rod and solving the balance equation of generalized moments. To assess this approach, our method is compared against the well-validated shooting and strain parameterization methods through two benchmarks. Compared to strain parameterization, our approach reduces the size of the residual vector, which considerably improves the computational time. Furthermore, since the linear strain is continuous, it converges more quickly to the solution, making shear and extension deformations more accurate. We have made our code available at https://gitlab.inria.fr/rjilani/ICRAX40.

## I. INTRODUCTION

Cosserat rod theory is currently being used as a replacement for the full 3D finite element method (FEM) in many continuum and soft robot problems [1]. In recent years, there have been two main approaches used to solve Cosserat rod boundary value problems (BVPs). The first is the continuous approach, primarily solved using the shooting method. The second is the reduction-based approach, primarily solved through a strain parameterization (SP) [1], [2]. The SP method has many advantages. Compared to the shooting method, SP is more suited for problems involving large loads because it provides better convergence guarantees [2]. Additionally, complex robots are easier to implement with SP than with the shooting method due to the use of Lagrangian mechanics [2]. Compared to FEM, SP requires fewer degrees of freedom (DOF) for solving BVPs, making it a powerful approach for control [1]. However, the shooting method is still the most used technique in the robotics community. One of the main reasons is that SP method is computationally expensive, primarily because of its high-dimensional residual vector [2]. In this work, we aim to address this issue by proposing a method that decreases the size of the residual vector without compromising accuracy.

## II. BOUNDARY VALUE PROBLEM

First, the static Cosserat rod BVP is presented. In the following, (.)' means $\frac{\partial(r)}{\partial X}$ with $X \in [0, l]$ being the arc-length and $l$ the length of the rod before deformation.

$$r' = R\Gamma, \quad r(0) = 0_{3 \times 1} \qquad (1)$$

$$R' = R\hat{K}, \quad R(0) = 1_{3 \times 3} \qquad (2)$$

$$N' = -\hat{K}N - \hat{N}, \quad N(l) = N_+ \text{ or } N(l) = R(l)^T n_+ \qquad (3)$$

$$C' = -\hat{K}C - \hat{\Gamma}N - \bar{C}, \quad C(l) = C_+ \text{ or } C(l) = R(l)^T c_+ \qquad (4)$$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Left: 45° bending, the tip force is applied in the inertial frame. Right: follower tip force, the force is applied in the cross-sectional frame.</figcaption></figure>

The first two lines represent the kinematics of Cosserat rods in the inertial frame (IF), while the other two represent the static balance equations in the cross-sectional frame (CF). $r$ is the position of the rod centerline, and $R$ is the rotation matrix of the cross-section orientation. $\Gamma$ and $K$ indicate the rates of change of $r$ and $R$ in the CF, respectively. For conciseness, the following components are represented in linear and angular DOF, respectively. $N$ and $C$ denote the stresses, $\bar{N}$ and $\bar{C}$ are the external distributed loads. $N_+$ and $C_+$ refer to the external tip loads in the CF, while $n_+$ and $c_+$ are their IF counterparts. The hat operator $\hat{\cdot}$ denotes the skew-symmetric matrix of a vector. Next, the BVP is supplemented with linear elastic constitutive laws

$$N = \mathcal{H}_l(\Gamma - \Gamma_0) \qquad (5)$$

$$C = \mathcal{H}_a(K - K_0) \qquad (6)$$

where $H_l$ and $H_a$ are the stiffness matrices and $\Gamma_0$ and $K_0$ are the rates of change in the initial configuration.

## III. ANGULAR STRAIN PARAMETERIZATION

The SP method proposed in [1], [2] leverages the following fact: if one knows $\Gamma$ and $K$, and uses the known boundary conditions at the proximal end of the rod ($r(0)$ and $R(0)$), one can forward integrate (from $X = 0$ to $X = l$) the kinematics (1, 2) to determine the rod configuration, $r$ and $R$. In this work, we go a step further by noting that if one knows $K$ and leverages the known boundary conditions at both ends of the rod, the rod configuration can be determined with the following algorithm (Algorithm 1):

1) Forward integrate (2).

2) Backward integrate (3,4), where $\Gamma$ is computed by inverting (5).
<figure><img src="image_4.png" /><figcaption>Fig. 2. Follower tip force: extension stress w.r.t the arc-length X.</figcaption></figure>

3) Forward integrate (1), where $\Gamma$ is computed by inverting (5).

In order to express $K$, we employ an angular strain parameterization (ASP)

$$K(X) = \Phi(X)q + K_0 \quad (7)$$

where $q \in \mathbb{R}^k$ is the generalized angular strain coordinates and $\Phi \in \mathbb{R}^{3 \times k}$ is a block diagonal matrix of a chosen shape function. Here, $k = \sum_{i=1}^3 k_i$, with $k_i$ being the number of modes (i.e. shape functions) for bending and torsion DOFs. The unknown $q$ is computed by solving the static balance of generalized moments:

$$K_{\epsilon\epsilon}q + Q_{\text{ext}} = 0 \quad (8)$$

where $K_{\epsilon\epsilon} \in \mathbb{R}^{k \times k}$ is the matrix of generalized (bending and torsion) stiffness, and $Q_{\text{ext}} \in \mathbb{R}^k$ is a vector of generalized external moments. They are computed as follows: $K_{\epsilon\epsilon} = \int_0^l \Phi^T H_a \Phi dX$ and $Q_{\text{ext}} = - \int_0^l \Phi^T C dX$ [1]. $C$ is computed with the first two lines of Algorithm 1 with $K$ determined using (7).

To capture all Cosserat rod's DOFs, [1] and [2] solve the balance equation of generalized forces and moments with the generalized strain coordinates $\tilde{q} \in \mathbb{R}^k$, where $\tilde{k} = \sum_{j=1}^6 \tilde{k}_j$, and $\tilde{k}_j$ being the number of modes for bending, torsion, shear and extension DOFs. In contrast, our approach implicitly accounts for shear and extension through the backward integration of (3).

## IV. NUMERICAL APPLICATIONS

The objectives of this section are to compare our ASP to the SP of [2] in terms of accuracy and computational efficiency. The accuracy of the two methods was tested w.r.t the shooting method. All methods were implemented in Python, and used Newton-Raphson for solving $q$ and $\tilde{q}$ with a termination tolerance of $10^{-7}$. Integration was done with Runge-Kutta 4(5). For each numerical application, we ran three simulations, each with a different number of modes: $k_i$ and $\tilde{k}_j \in \{2, 4, 8\}$. Legendre polynomials were used as shape functions. For conciseness, the physical parameters of the rods will not be reported. We invite the reader to refer to the referenced papers.

<table><caption>TABLE I
COMPARISON BETWEEN STRAIN PARAMETERIZATION (SP) AND ANGULAR STRAIN PARAMETERIZATION (ASP)</caption><thead><tr><th>Method</th><th>Avg. e<sub>r</sub> (%)</th><th>Avg. e<sub>R</sub> (deg)</th><th>Avg. nb<sub>i</sub></th><th>Avg. nb<sub>f</sub></th><th>Avg. time (s)</th></tr></thead><tbody><tr><td>ASP 8 modes</td><td>5.29e-06</td><td>1.21e-05</td><td>4.48</td><td>112</td><td>11.8</td></tr><tr><td>SP 8 modes</td><td>5.31e-06</td><td>1.21e-05</td><td>4.48</td><td>219</td><td>27.6</td></tr><tr><td>ASP 4 modes</td><td>5.22e-03</td><td>1.24e-04</td><td>4.48</td><td>58.2</td><td>8.18</td></tr><tr><td>SP 4 modes</td><td>5.24e-03</td><td>1.24e-04</td><td>4.48</td><td>112</td><td>17.8</td></tr><tr><td>ASP 2 modes</td><td>2.1</td><td>4.41e-2</td><td>4.33</td><td>30.3</td><td>4.9</td></tr><tr><td>SP 2 modes</td><td>2.1</td><td>4.41e-2</td><td>4.33</td><td>56.3</td><td>10.5</td></tr></tbody></table>

### A. 45° bending

In the initial configuration, the rod is 1/8 of a circular arc, then an out-of-plane tip force ($n_+$) in the x-direction is applied incrementally in 20 steps, from 0 to 3000N (cf. Fig. 1) [3]. Table I shows the average error, across all loading steps, of tip position and tip rotation. The errors are computed as follows: $e_r = \frac{|r(l) - r_s(l)|}{l} \times 100$ for positions and $e_R = \cos^{-1}\left(\frac{\text{trace}(R_s(l)R^T(l))-1}{2}\right)$ for rotations, where subscript $s$ means that the component belongs to the shooting method. As expected, by increasing the number of modes, the errors decrease. Furthermore, our method is slightly more accurate due to continuous linear strain. Table I also shows the average, across all loading steps, of the number of iterations executed by Newton-Raphson ($nb_i$), the number of evaluations of the residual function ($nb_f$) (8), and the execution time. It reveals that $nb_i$ is identical for the two methods. However, ASP's $nb_f$ is approximately half that of SP's, due to the residual vector size ratio $k/\tilde{k} = 0.5$. Consequently, the execution time is also approximately halved. Furthermore, ASP solves Cosserat BVPs with fewer DOFs than SP, which makes it better for control.

### B. Follower tip force

In the initial configuration, the rod is straight, then a follower tip force ($N_+$) in the x-direction is applied incrementally in 20 steps, from 0 to 130kN (cf. Fig. 1) [1]. Fig. 2 plots the extension stress of the last loading step w.r.t $X$. We can conclude that our method converges more quickly to the solution compared to the SP method. This is due to the continuous linear strain model.

## ACKNOWLEDGMENT

This work was supported by the French Agence Nationale de la Recherche under Grant ANR-20-CE45-0011.

## REFERENCES

[1] F. Boyer, V. Lebastard, F. Candelier, and F. Renda, “Dynamics of Continuum and Soft Robots: A Strain Parameterization Based Approach,” *IEEE Transactions on Robotics*, vol. 37, pp. 847–863, June 2021.

[2] M. Tummers, V. Lebastard, F. Boyer, J. Troccaz, B. Rosa, and M. T. Chikhaoui, “Cosserat Rod Modeling of Continuum Robots from Newtonian and Lagrangian Perspectives,” *IEEE Transactions on Robotics*, vol. 39, pp. 2360–2378, June 2023.

[3] R. Jilani, P.-F. Villard, and E. Kerrien, “An Orthogonal Collocation Method for Static and Dynamic Cosserat Rods,” in *2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*, pp. 4328–4333, Oct. 2023.