Discrete Cosserat Approach for Multi-Section Soft Robots Dynamics

TL;DR

This paper introduces a discrete Cosserat model for multi-section soft robots that captures shear and torsional deformations, providing a more comprehensive and unified dynamics framework than traditional curvature-based models.

Contribution

It presents a novel discrete Cosserat approach that accounts for shear and torsion, bridging the gap between continuous models and practical soft robot dynamics.

Findings

Model accurately predicts both in-plane and out-of-plane motions.

Simulation and experiments validate the model's effectiveness.

Provides a unified framework for soft robot dynamics.

Abstract

In spite of recent progress, soft robotics still suffers from a lack of unified modeling framework. Nowadays, the most adopted model for the design and control of soft robots is the piece-wise constant curvature model, with its consolidated benefits and drawbacks. In this work, an alternative model for multisection soft robots dynamics is presented based on a discrete Cosserat approach, which, not only takes into account shear and torsional deformations, essentials to cope with out-of-plane external loads, but also inherits the geometrical and mechanical properties of the continuous Cosserat model, making it the natural soft robotics counterpart of the traditional rigid robotics dynamics model. The soundness of the model is demonstrated through extensive simulation and experimental results for both plane and out-of-plane motions.

Figures

Tables

Table 1. TABLE I: Design Parameters of the Prototype

Parameter	Value	Parameter	Value
$R_{max}$	$15$ $mm$	$d_1$, $d_2$, $d_3$, $d_4$	$9$ $mm$
$R_{min}$	$4$ $mm$	$d_5$, $d_6$, $d_7$, $d_8$	$6$ $mm$
$L_1$	$98$ $mm$	$d_9$, $d_{10}$, $d_{11}$, $d_{12}$	$3$ $mm$
$L_2$	$203$ $mm$	$gr$	$9.81$ $\frac{m}{s^2}$
$L_3$	$311$ $mm$	$C_x$	$0.01$
$L_4$	$418$ $mm$	${\rho }_w$	$1.022$ $\frac{kg}{d{m}^3}$
$E$	$110$ $kPa$	$C_y$	$2.5$
$\mu$	$300$ $Pa\cdot s$	$C_z$	$2.5$
$\nu$	$0.5$	$B_y$	$1.5$
$\rho$	$1.08$ $\frac{kg}{d{m}^3}$	$B_z$	$1.5$

Equations

$$g=\left(\begin{array}[]{cc}R&u\\ 0&1\end{array}\right).$$

$$\widehat{\xi}=\left(\begin{array}[]{cc}\tilde{\mathsf{k}}&\mathsf{q}\\ 0&0\end{array}\right)\in\mathfrak{se}(3)\;\text{,}\;\xi=\left(\mathsf{k}^{T},\mathsf{q}^{T}\right)^{T}\in\mathbb{R}^{6}\;\text{,}$$

$$\widehat{\eta}=\left(\begin{array}[]{cc}\tilde{\mathsf{w}}&\mathsf{v}\\ 0&0\end{array}\right)\in\mathfrak{se}(3)\;\text{,}\;\eta=\left(\mathsf{w}^{T},\mathsf{v}^{T}\right)^{T}\in\mathbb{R}^{6}\;.$$

$$g^{\prime }=g\hat{\xi }.$$

$${\eta }^{\prime }=\dot{\xi }-{\mathrm{ad}}_{\xi }\eta \text{,}$$

$$\mathrm{ad}_{\xi}=\left(\begin{array}[]{cc}\widetilde{\mathsf{k}}&0\\ \widetilde{\mathsf{q}}&\widetilde{\mathsf{k}}\end{array}\right)\;\text{,}\;\mathrm{ad}_{\xi}^{*}=-\mathrm{ad}_{\xi}^{T}=\left(\begin{array}[]{cc}\widetilde{\mathsf{k}}&\widetilde{\mathsf{q}}\\ 0&\widetilde{\mathsf{k}}\end{array}\right)\;.$$

$${\dot{\eta }}^{\prime }=\ddot{\xi }-{\mathrm{ad}}_{\dot{\xi }}\eta -{\mathrm{ad}}_{\xi }\dot{\eta }.$$

$$\mathcal{M}\dot{\eta }+{\mathrm{ad}}_{\eta }^{\ast }\left(\mathcal{M}\eta \right)={\mathcal{F}}_i^{\prime }+{\mathrm{ad}}_{\xi }^{\ast }{\mathcal{F}}_i+{\bar{\mathcal{F}}}_a+{\bar{\mathcal{F}}}_e\text{,}$$

$${\bar{\mathcal{F}}}_a(X,t)=-\left({\mathcal{F}}_a^{\prime }+{\mathrm{ad}}_{\xi }^{\ast }{\mathcal{F}}_a\right)\text{and}{\bar{\mathcal{F}}}_a(X,t)=0\text{,}$$

$${\mathcal{F}}_i(X)=\Sigma \left(\xi -{\xi }^0\right)+{\rm Y}\dot{\xi }\text{,}$$

$$\begin{array}[]{l}\bar{\mathcal{F}}_{e}=\left(1-\rho_{w}/\rho\right)\mathcal{M}\mathrm{Ad}_{g_{r}g(X)}^{-1}\mathcal{G}-\mathcal{D}||\mathsf{v}||\eta+\delta(X-\bar{X})\mathcal{F}_{p}\;\text{,}\\ \mathcal{M}_{a}=\mathcal{M}+\mathcal{A}\end{array}$$

$$\mathrm{Ad}_{g}=\left(\begin{array}[]{cc}R&0\\ \widetilde{\mathsf{u}}R&R\end{array}\right)\;\text{,}\;\mathrm{Ad}_{g}^{*}=\mathrm{Ad}_{g}^{-T}=\left(\begin{array}[]{cc}R&\widetilde{\mathsf{u}}R\\ 0&R\end{array}\right)\;.$$

$$g(X)=g(L_{n-1})e^{\left(X-L_{n-1}\right){\hat{\xi }}_n}.$$

$$\begin{array}{rl}{e}^{\left(X-L_{n-1}\right){\hat{\xi }}_n}=I_4+\left(X-L_{n-1}\right){\hat{\xi }}_n& \\ +\frac{1}{{\theta }_n^2}\left(1-\cos \left(\left(X-L_{n-1}\right){\theta }_n\right)\right){\hat{\xi }}_n^2& \\ +\frac{1}{{\theta }_n^3}\left(\left(X-L_{n-1}\right){\theta }_n-\sin \left(\left(X-L_{n-1}\right){\theta }_n\right)\right){\hat{\xi }}_n^3& =:g_n(X)\text{,}\end{array}$$

$$e^{\left(X-L_{n-1}\right){\hat{\xi }}_n}=I_4+\left(X-L_{n-1}\right){\hat{\xi }}_n\text{,}$$

$$g(X)=g(L_{n-1})g_n(X)\text{,}$$

$$\begin{array}{rl}\eta (X)=& e^{-\left(X-L_{n-1}\right){\mathrm{ad}}_{{\xi }_n}}\\ & \left(\eta (L_{n-1})+{\int }_{L_{n-1}}^X{e}^{\left(s-L_{n-1}\right){\mathrm{ad}}_{{\xi }_n}}ds{\dot{\xi }}_n\right).\end{array}$$

$$\begin{array}{rl}{e}^{x{\mathrm{ad}}_{{\xi }_n}}=I_6+\frac{1}{2{\theta }_n}(3\sin (x{\theta }_n)-x{\theta }_n\cos (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}& \\ +\frac{1}{2{\theta }_n^2}(4-4\cos (x{\theta }_n)-x{\theta }_n\sin (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}^2& \\ +\frac{1}{2{\theta }_n^3}(\sin (x{\theta }_n)-x{\theta }_n\cos (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}^3& \\ +\frac{1}{2{\theta }_n^4}(2-2\cos (x{\theta }_n)-x{\theta }_n\sin (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}^4& \\ ={\mathrm{Ad}}_{g_n(X)}& ,\end{array}$$

$$\eta (X)={\mathrm{Ad}}_{g_n(X)}^{-1}\left(\eta (L_{n-1})+{\mathrm{AD}}_{g_n}(X){\dot{\xi }}_n\right)\text{,}$$

$$\begin{array}{r}{\mathrm{AD}}_{g_n}(X):={\int }_{L_{n-1}}^X{\mathrm{Ad}}_{g_n(s)}ds=\\ x{I}_6+\frac{1}{2{\theta }_n^2}(4-4\cos (x{\theta }_n)-x{\theta }_n\sin (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}\\ +\frac{1}{2{\theta }_n^3}(4x{\theta }_n-5\sin (x{\theta }_n)+x{\theta }_n\cos (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}^2\\ +\frac{1}{2{\theta }_n^4}(2-2\cos (x{\theta }_n)-x{\theta }_n\sin (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}^3\\ +\frac{1}{2{\theta }_n^5}(2x{\theta }_n-3\sin (x{\theta }_n)+x{\theta }_n\cos (x{\theta }_n)){\mathrm{ad}}_{{\xi }_n}^4.\end{array}$$

$$\begin{array}{rl}\dot{\eta }(X)=& e^{-x{\mathrm{ad}}_{{\xi }_n}}\\ & \left(\dot{\eta }(L_{n-1})+{\int }_{L_{n-1}}^X{e}^{x{\mathrm{ad}}_{{\xi }_n}}\left({\ddot{\xi }}_n-{\mathrm{ad}}_{{\dot{\xi }}_n}\eta \right)ds\right).\end{array}$$

$$\begin{array}{rl}& \dot{\eta }(X)={\mathrm{Ad}}_{g_n(X)}^{-1}\\ & \left(\dot{\eta }(L_{n-1})+{\mathrm{AD}}_{g_n}(X){\ddot{\xi }}_n-{\int }_{L_{n-1}}^X{\mathrm{Ad}}_{g_n(s)}{\mathrm{ad}}_{{\dot{\xi }}_n}\eta (s)ds\right)\text{.}\end{array}$$

$${\mathrm{Ad}}_{g_n(s)}{\mathrm{ad}}_{{\dot{\xi }}_n}\eta (s)={\mathrm{ad}}_{{\mathrm{Ad}}_{g_n(s)}{\dot{\xi }}_n}\left(\eta \left(L_{n-1}\right)+{\mathrm{AD}}_{g_n}(s){\dot{\xi }}_n\right).$$

$${\mathrm{Ad}}_{g_n(s)}{\mathrm{ad}}_{{\dot{\xi }}_n}\eta (s)={\mathrm{ad}}_{{\mathrm{Ad}}_{g_n(s)}{\dot{\xi }}_n}\eta \left(L_{n-1}\right)\text{,}$$

$$\dot{\eta }(X)={\mathrm{Ad}}_{g_n}^{-1}\left(\dot{\eta }(L_{n-1})-{\mathrm{ad}}_{{\mathrm{AD}}_{g_n}{\dot{\xi }}_n}\eta \left(L_{n-1}\right)+{\mathrm{AD}}_{g_n}{\ddot{\xi }}_n\right).$$

$$\eta (X)=\sum _{i=1}^n\left(\prod _{j=n}^i{\mathrm{Ad}}_{g_j\left(min\left(L_j,X\right)\right)}^{-1}\right){\mathrm{AD}}_{g_i}\left(min\left(L_i,X\right)\right){\dot{\xi }}_i\text{.}$$

$$\overrightarrow{\xi }={\left[{\xi }_1^T{\xi }_2^T\cdots {\xi }_N^T\right]}^T\in {\mathbb{R}}^{6N}\text{,}$$

$$\eta (X)=J(X)\dot{\overrightarrow{\xi }}\text{,}$$

$$J(X)=\left\{\begin{array}[]{lcl}\left[\left(\mathrm{Ad}_{g_{1}}^{-1}\mathrm{AD}_{g_{1}}\right)(X)\;\;0_{6}\;\;\cdots\right]&\mbox{if}&0<X\leqslant L_{1}\\ \\ \left[\mathrm{Ad}_{g_{2}(X)}^{-1}\left(\mathrm{Ad}_{g_{1}}^{-1}\mathrm{AD}_{g_{1}}\right)(L_{1})\;\;\left(\mathrm{Ad}_{g_{2}}^{-1}\mathrm{AD}_{g_{2}}\right)(X)\;\;0_{6}\;\;\cdots\right]&\mbox{if}&L_{1}<X\leqslant L_{2}\\ \\ \left[\mathrm{Ad}_{g_{3}(X)}^{-1}\mathrm{Ad}_{g_{2}(L_{2})}^{-1}\left(\mathrm{Ad}_{g_{1}}^{-1}\mathrm{AD}_{g_{1}}\right)(L_{1})\;\;\mathrm{Ad}_{g_{3}(X)}^{-1}\left(\mathrm{Ad}_{g_{2}}^{-1}\mathrm{AD}_{g_{2}}\right)(L_{2})\;\;\cdots\;\;0_{6}\right]&\mbox{if}&L_{1}<X\leqslant L_{2}\\ \\ \vdots&\vdots&\vdots\\ \\ \left[\prod\limits_{1}^{j=N}\mathrm{Ad}_{g_{j}\left(min\left(L_{j},X\right)\right)}^{-1}\mathrm{AD}_{g_{1}}(L_{1})\;\;\cdots\;\;\left(\mathrm{Ad}_{g_{N}}^{-1}\mathrm{AD}_{g_{N}}\right)(X)\right]&\mbox{if}&L_{N-1}<X\leqslant L_{N}\end{array}\right.$$

$${}^{s}J(X)={\mathrm{Ad}}_{g(X)}J(X).$$

Full text

I Introduction

Since the beginning of the soft robotics field, many researchers have contributed in the development of mathematical modeling approaches which could be able to describe the kinematics and dynamics of such infinite Degrees of Freedom (DoF) robots, while addressing the challenging requirements imposed by their robotic applications [24], [11]. In order to meet the standards achieved in traditional rigid robotics, a model for soft robotics should be at the same time computational inexpensive and sufficiently accurate. Furthermore, it should be able to shed light on the mathematical sub-models and to encompass them in a unified framework. Such a modeling framework is the necessary condition for developing the physical designs and control architectures of these new soft robots as well as their task-related motions and path planning.

Despite the short history of soft robotics, important results have been already achieved and several complementary modeling approaches have been proposed to date. Those approaches can be divided into three main categories: Piece-wise Constant Curvature (PCC) models, continuum Cosserat models and 3D Finite Elements Models (FEM).

The PCC modeling approach is by far the most adopted in the soft robotics community [28]. It represents the soft robot as a finite collection of circular arcs, which can be described by only three parameters (radius of curvature, angle of the arc and bending plane), a simplification which drastically reduces the number of variables needed. Originally devoted to kinematics modeling [10], this approach has been extended and improved over the years with excellent results as in [9, 8]. In spite of this success, the constant curvature assumption is not always valid, especially when the robot is subject to non-negligible external loads including gravity.

The continuum Cosserat approach is an infinite DoF model where the soft robot is represented by continuously stacking an infinite number of infinitesimal micro-solids. It has been primarily used in the context of hyper-redundant robot [6], and more recently applied to soft robotics locomotion [3, 2] and manipulation [23], in both static [18] and dynamic [19] conditions. This approach has been also extended to shell-like soft robots for underwater locomotion inspired by cephalopods [20], [21]. Despite their accuracy and fidelity to the continuous mechanics of the soft robots, the resulting partial differential equations are computationally demanding and difficult to use for control purposes.

Finally, a FEM based approach has also been explored for modeling and real-time control of soft robots [12]. It is so far limited to quasi-static conditions, and needs linearisation of the structural elasticity that may not apply to many soft robot geometries.

Although it might be impossible, due to physical reasons, to achieve the same elevate standard reached by the mathematical models for rigid robotics, the research outlined above constitute a significant attempt in this direction. In the present paper, we build upon the two main pillars achieved so far to obtain, in the authors opinion, one of the most promising approach towards a unified mathematical framework between traditional and soft robotics. Going further into details, the continuous model developed with the Cosserat approach is discretized in order to implement the PCC idea of reducing the dimension of the configuration space by assuming a piece-wise constant deformation along the soft manipulator. As a consequence, the soft manipulator is completely described by a finite set of strain vectors which plays the same role as that of the joint vector for traditional robotics.

The strains allowed by the Cosserat approach include torsion and shears along with curvature and elongation. Thus, we call this method Piece-wise Constant Strain (PCS) model. With respect to the PCC model, the PCS model not only takes into account shears and torsion, which are both essential to cope with out-of-plane external loads, but also shares a common geometric structure with the equations of motion of their rigid robotics counterpart. As a matter of fact, the PCS model provides a direct forward kinematics between the joint space and the task space without any intermediate map. Furthermore, based on the $SE(3)$ geometry of the Cosserat approach, it guarantees a closer relation with the rigid body geometry of the traditional robotics. Finally, the discrete Cosserat framework allows the adaptation of different actuation solutions and external loads models, including the interaction with a dense medium, without any significant changes in the structure of the model, and is so more independent from the specific applications.

Part of the present work has been presented in the conference paper [16]. Beyond this work, the full multi-section dynamics is addressed here and a recursive algorithm for calculating the coefficients of the dynamics equations is presented. Furthermore, the homogeneity with the standard rigid robotics theory is highlighted and extensive simulations along with experimental results are shown for the multi-section dynamics case. In the following, in section II the continuous Cosserat model is briefly reminded in order to introduce the discretization developed in the subsequent section III. Finally, the model is corroborated through extensive simulation and experimental results in sections IV and V for the general case of a soft manipulator operating in a dense medium like water.

II Continuous Cosserat Model

In the Cosserat theory, the configuration of a micro-solid of a soft body with respect to the inertial frame at a certain time is characterized by a position vector $u$ and an orientation matrix $R$, parameterized by the material abscissa $X\in[0,L]$ along the robot arm. Thus, the configuration space is defined as a curve $g(\cdot):X\mapsto g(X)\in SE(3)$ with

[TABLE]

Then, the strain state of the soft arm is defined by the vector field along the curve $g(\cdot)$ given by $X\mapsto\widehat{\xi}(X)=g^{-1}\partial g/\partial X=g^{-1}g^{\prime}\in\mathfrak{se}(3)$, where the hat is the isomorphism between the twist vector representation and the matrix representation of the Lie algebra $\mathfrak{se}(3)$. The components of this field are specified in the (micro-)body frames as:

[TABLE]

where $\mathsf{q}(X)$ represents the linear strains, and $\mathsf{k}(X)$ the angular strains. The tilde is the isomorphism between three dimensional vectors and skew symmetric matrices.

The time evolution of the configuration curve $g(\cdot)$ is represented by the twist vector field $X\mapsto\eta(X)\in\mathbb{R}^{6}$ defined by $\widehat{\eta}(X)=g^{-1}\partial g/\partial t=g^{-1}\dot{g}$. This field can be detailed in terms of their components in the (micro-)body frames as:

[TABLE]

Where $\mathsf{v}(X)$ and $\mathsf{w}(X)$ are respectively the linear and angular velocity at a given time instant.

II-A Continuous Kinematics

Given the above construction, we can obtain the kinematic equations relating the strains of the robot arm $\xi$ with the position $g$, velocity $\eta$ and acceleration $\dot{\eta}$ for each infinitesimal micro-solid constituting the robot. By definition, the first equation is given by:

[TABLE]

Then, the equality of mixed partial derivatives $(\dot{g})^{\prime}=\dot{(g^{\prime})}$ gives the following compatibility equation between strain and velocity:

[TABLE]

where $\mathrm{ad}$ is the adjoint map defined as (together with the coadjoint map $\mathrm{ad}^{*}$):

[TABLE]

Finally, by taking the derivative of (2) with respect to time, we obtain the continuous model of acceleration:

[TABLE]

II-B Continuous Dynamics

In [4] it is shown that Cosserat beam dynamics can be directly derived from the extension to continuum media of a variational calculus originally introduced by H. Poincaré [14]. In contrast to usual Lagrangian mechanics, this calculus allows deriving the dynamics of a system whose the configuration space definition requires the structure of Lie group. In this context, the dynamics of the Cosserat medium can be entirely derived from a Lagrangian density $\mathfrak{T}(\eta)-\mathfrak{U}(\xi)$, where $\mathfrak{T}$ and $\mathfrak{U}$ are functions of the Lie algebra vectors modelling the densities of kinetic and elastic energy of the Cosserat beam per unit of material length $X$. Applying this variational calculus to such a density leads to the strong form of a Cosserat beam with respect to the micro-solid frames.

[TABLE]

where $\mathcal{F}_{i}(X)=\partial\mathfrak{U}/\partial X$ is the wrench of internal forces, $\bar{\mathcal{F}}_{a}(X,t)$ is the distributed actuation loads, $\bar{\mathcal{F}}_{e}(X)$ is the external wrench of distributed applied forces and $\mathcal{M}(X)$ is the screw inertia matrix. Let us specify the angular and linear components of the internal and external wrenches: $\mathcal{F}_{i}=\left(\mathsf{M}_{i}^{T},\mathsf{N}_{i}^{T}\right)^{T}$, $\bar{\mathcal{F}}_{a}=\left(\mathsf{m}_{a}^{T},\mathsf{n}_{a}^{T}\right)^{T}$, $\bar{\mathcal{F}}_{e}=\left(\mathsf{m}_{e}^{T},\mathsf{n}_{e}^{T}\right)^{T}\in\mathbb{R}^{6}$, where $\mathsf{N}_{i}(X)$ and $\mathsf{M}_{i}(X)$ are the internal force and torque vectors, $\mathsf{n}_{a}(X,t)$ and $\mathsf{m}_{a}(X,t)$ are the actuation force and torque inputs, while $\mathsf{n}_{e}(X)$ and $\mathsf{m}_{e}(X)$ are the external force and torque for unit of $X$. By choosing a local (micro-)body frame oriented as in figure 2, with the $X$ axis pointing toward the tip of the robot arm and the $Y$ and $Z$ axes laying on the plane of the section (considered symmetric), the screw inertia matrix is equal to: $\mathcal{M}=diag(J_{x},J_{y},J_{z},A,A,A)\rho$, where $\rho$ is the body density, $A$ is the section area and $J_{y}$, $J_{z}$, $J_{x}$ are respectively the bending and torsion second moment of inertia of the beam cross-section.

Let us now specify the models of the distributed actuation, external load and internal forces appearing in (4) for the general case of a soft robot arm moving in a dense surrounding medium like water. Considering the two most important actuation systems implemented in soft robotics, the cable driven and the fluidic actuation [24], we have respectively:

[TABLE]

where $\mathcal{F}_{a}$ is the cable wrench acting on the micro-solid given by the cable tension and the cable path from the tip to the base [19], [23], [17]. The model of the fluidic actuator, widely used in soft robotics nowadays [15], condensates the action of the pressure in a concentrated load at the tip of the section (Fig. 1).

Regarding the wrench of internal passive forces, a linear visco-elastic constitutive model, based on the Kelvin Voigt assumptions, is chosen [19].

[TABLE]

where $\Sigma$ and $\Upsilon$ are constant screw stiffness and viscosity matrices, equal to $\Sigma=diag(GJ_{x},EJ_{y},EJ_{z},EA,GA,GA)$, $\Upsilon=diag(J_{x},3J_{y},3J_{z},3A,A,A)\upsilon$, $E$ being the Young modulus, $G$ the shear modulus and $\upsilon$ the shear viscosity modulus; $\xi^{0}=[0\;0\;0\;1\;0\;0]^{T}$ stands for the zeros strain vector in the reference straight configuration. No other assumptions except the constitutive model are needed to describe the elastic behavior of the robot arm.

As for the external loads, we have considered the general case of underwater operation, i.e. distributed loads due to gravity and buoyancy, drag, added mass and a concentrated/point load due to externally applied loads or contacts [19]:

[TABLE]

where $\rho_{w}$ is the water density, $\mathcal{G}=[0\;0\;0\;-9.81\;0\;0]^{T}$ is the gravity twist with respect to the inertial frame (in accordance with the choice of inertial frame given in figure 2), $g_{r}$ is the transformation between the spatial frame and the base frame of the soft manipulator, $\mathcal{D}(X)$ is the drag fluid dynamics coefficient, $\delta(\cdot)$ is the Dirac distribution, $\mathcal{F}_{p}$ is the wrench corresponding to the point load applied at $\bar{X}$ and $\mathcal{A}(X)$ is the added mass fluid dynamics coefficient. Note here that replacing $\mathcal{M}$ by $\mathcal{M}_{a}$ in (4) allows modeling inertial hydrodynamics forces exerted along the arm. Finally, we have introduced the Adjoint representation ($\mathrm{Ad}$) of Lie group $SE(3)$, defined as (together with the coAdjoint map $\mathrm{Ad}^{*}$):

[TABLE]

Finally, when the soft arm is working in a sparse medium like air we will let $\rho_{w}$, and consequently $\mathcal{D}$ and $\mathcal{A}$, be equal to zero.

III Discrete Cosserat Model

Equations (1), (2), (3) and (4) of the continuous Cosserat model are suitable to model the kinematics and dynamics of soft robots expressing a non-constant deformation, as it has been presented in [19] (and [20], [21] for bi-dimensional bodies). In the subsequent development, we unify the constant and non-constant cases under the same mathematical framework. To that end, the continuous model is discretized by an analytic spatial integration. This is allowed by the piece-wise constant strain assumption which provides the condition to analytically integrate the continuum model and leads to the extension of the piece-wise constant curvature model, by including torsion and shears, without any additional effort. Furthermore, a profound and useful parallelism with the rigid manipulators theory can be achieved, which leads to the soft robot counterpart of the Lagrangian model of rigid serial manipulators.

III-A Piece-wise Constant Strain Kinematics

At any instant $t$, considering the strain field $\xi(X)$ constant along each of the $N$ sections of the soft arm, we can replace the continuous field with a finite set of $N$ twist vectors $\xi_{n}$ ($n\in\{1,2,..,N\}$), which play the role of the joint vectors of traditional rigid robotics. Under this assumption, equation (1) becomes an homogeneous, linear, matrix differential equation with constant coefficients, which can be analytically solved at any section $n$ using the matrix exponential method with the appropriate interval of $X$ and initial value [7]. Going further into details, the material abscissa $X\in[0,L]$ is divided into $N$ sections of the form $[0,L_{1})$, $(L_{1},L_{2})\dots(L_{N-1},L_{N}]$ (with $L_{N}=L$) and the initial value for the differential equation of the section $n$ is given by the solution at the right boundary of the previous section ($X=L_{n-1}$). In other words, the solutions are glued together, one on top of the other. With these considerations, the integration of (1) at a certain instant $t$ becomes:

[TABLE]

It turns out that the infinite series of the exponential in (8) can be expressed in a compact way as follows [25]:

[TABLE]

where $\theta_{n}^{2}=\mathsf{k_{n}}^{T}\mathsf{k_{n}}$. For straight configurations of the section, we have $\widehat{\xi}_{n}^{2}=0_{4}$ and hence:

[TABLE]

which allows circumventing the well known singularity of straight arm pose of the PCC models [27], [22]. Equation (9) can be viewed as the $SE(3)$ counterpart of the Rodrigues formula in $SO(3)$. Calling $g_{n}(X)$ the exponential function in (9), equation (8) can be written in the more familiar way:

[TABLE]

which recursively returns the position and orientation of the micro-solid at $X$ knowing the set of strains $\xi_{n}$ only.

Similarly, the velocity of each micro-solid $\eta(X)$ can be obtained by a piece-wise integration of the continuum model (2). Under constant strains condition, at each section $n$ and time $t$, equation (2) is a non-homogeneous, linear, matrix differential equation with constant coefficients (reminds that also $\dot{\xi}_{n}$ is piece-wise constant) which can be analytically solved using the variation of parameters method, with the appropriate initial value [7].

[TABLE]

Again, the exponential function in (11) can be expressed with a finite number of terms [25] (for the sake of presentation, $x=X-L_{n-1}$ holds in the following).

[TABLE]

where for straight configurations we have $\mathrm{ad}_{\xi_{n}}^{2}=0_{6}$ and thus, taking the limit for $\theta_{n}\rightarrow 0$, $e^{x\mathrm{ad}_{\xi_{n}}}=I_{6}+x\mathrm{ad}_{\xi_{n}}$. Thanks to the fact that $e^{\mathrm{ad}_{\xi}}=\mathrm{Ad}_{e^{\widehat{\xi}}}$ ([1] pg. 403 Lemma 7.5.9), we can notice that the exponential function (12) is nothing else but the Adjoint representation of the Lie group transformation $g_{n}(X)$ of (9). With this definition at hand, equation (11) can be rewritten as follows:

[TABLE]

where we have defined:

[TABLE]

Remarkably, equation (13) recursively compute the velocity of any micro-solid at $X$ along the soft arm as a function of the set of strains $\xi_{n}$ and strain rates $\dot{\xi}_{n}$.

Finally, the acceleration of any micro-solid at $X$ ($\dot{\eta}(X)$) can be calculated at any time $t$, by means of a piece-wise integration of the continuous equation (3). Considering constant strains along one section, equation (3) is a non-homogeneous, linear, matrix differential equation with non-constant coefficients (given by the term $\mathrm{ad}_{\dot{\xi}_{n}}\eta$ which is not constant with respect to $X$ due to $\eta(X)$). A direct application of the variation of parameters method with the appropriate initial value gives:

[TABLE]

Then, by virtue of the definitions of $\mathrm{Ad}_{g_{n}}$, $\mathrm{AD}_{g_{n}}$, we obtain:

[TABLE]

Let us focus on the term $\mathrm{Ad}_{g_{n}(s)}\mathrm{ad}_{\dot{\xi}_{n}}\eta(s)$ inside the integral of the right end side. First, by means of equation (13) and the properties of the adjoint map, we can write:

[TABLE]

Then, evoking the linearity and anticommutativity of the adjoint map, and using equations (12) and (14), we obtain the equivalence

[TABLE]

which once substituted in (16) (and using twice the anticommutativity of the adjoint map to make appear $\mathrm{AD}_{g_{n}}$), gives the model of accelerations as follows:

[TABLE]

Again, equation (17) returns the acceleration of any micro-solid at $X$ by means of the set of strains $\xi_{n}$, strain rates $\dot{\xi}_{n}$ and rates of strain rate $\ddot{\xi}_{n}$ only.

In order to develop the discrete Cosserat dynamic model for soft robots a relation between the kinematics quantities $\eta$, $\dot{\eta}$ and a joint vector for soft robotics needs to be established. To do so, we back track to the base the velocity term $\eta(L_{n-1})$ on the right end side of (13), which becomes:

[TABLE]

where $j$ is a descending index, we have considered a fixed base ($\eta(0)=0_{6\times 1}$) and $L_{n-1}<X\leqslant L_{n}$. Introducing the soft robots joint vector

[TABLE]

equation (18) can be expressed as:

[TABLE]

which define the softs robot geometric Jacobian $J(X)\in\mathbb{R}^{6\times 6N}$, shown in (20).

It is important to notice that the Jacobian (20) is calculated directly from the strains $\xi_{n}$ by means of equations (12) and (14). For this reason, in accordance with the rigid manipulators theory, it is referred to as geometric Jacobian, in contrast with the analytic Jacobian. Furthermore, the action of the Jacobian $J(X)$ on the joint vector $\overrightarrow{\xi}$ returns the body velocity $\eta(X)$ which is expressed in the (micro-)body coordinate frame. Accordingly, $J(X)$ is called body Jacobian. The relation with the corresponding spatial Jacobian ${}^{s}J(X)$, which returns the spatial velocity ${}^{s}\eta(X)$ expressed in the fixed spatial frame, is obtained by multiplying both side of (19) with $\mathrm{Ad}_{g(X)}$ and reads:

[TABLE]

Finally, by taking the time derivative of (19) the acceleration vector $\dot{\eta}(X)$ is obtained as:

[TABLE]

where $\dot{J}(X)$ is obtained by a lengthy but straightforward calculation. Defining the $6\times 6$ components of the Jacobian as $J(X)=\left[S_{1}(X)\;\;S_{2}(X)\;\;\cdots\;\;S_{N}(X)\right]$, the time derivative of the Jacobian can be expressed as:

[TABLE]

where $L_{n-1}<X\leqslant L_{n}$ and we have defined $J_{i}(X)$ as the Jacobian containing $0_{6}$ elements except for the $i$-th: $J_{i}(X):=[0_{6}\cdots S_{i}(X)\cdots 0_{6}]\in\mathbb{R}^{6\times 6N}$. Alternatively, equation (21) and the expression of Jacobian derivative (22) can be obtained by back tracking the acceleration $\dot{\eta}(L_{n-1})$ and velocity $\eta(L_{n-1})$ terms on the right side of (17).

Comparison with the PCC Model

The development above led us to three kinematics equations (8), (19) and (21), which give a model to calculate all the kinematic quantities from the knowledge of the joint space of the piece-wise soft arm, in a very similar fashion to traditional rigid manipulators. Compared to the PCC model, the discrete Cosserat approach presented here is able to handle not only constant curvature and elongation, but also shear and torsion, which are fundamental to deal with the strong interactions with the environment characteristic of locomotion and manipulation.

Furthermore, the joint space $\overrightarrow{\xi}$ composed by the $N$ constant strains $\xi_{n}$ is directly related to the configuration kinematics through the equations (8), (19) and (21), while the PCC model needs an additional map between the joint space and the arc parameters space, composed by the length, the curvature and the plane of bending of the section. This allowed us to build a geometric Jacobian instead of an analytic Jacobian, which preserves the natural geometric structure of the motion.

Finally, the intrinsic geometry of the soft robots is reveled. In fact, recognizing (9) as a screw motion in space, we can conclude that each section forms an arc of screw whose parameters are determined by the constant strain $\xi_{n}$ by adapting the formulas normally used for time-twist [17].

III-B Piece-wise Constant Strain Dynamics

In this section we derive the generalized equation of motion of the multi-section piece-wise constant strain model. To that end, we reconsider the continuous dynamics (4), that we restate in the weak form of virtual works, i.e. for any field: $\delta\zeta(\cdot):X\mapsto\delta\zeta(X)\in\mathfrak{se}(3)$:

[TABLE]

Note that the above weak form can be derived from the extended Poincaré variational calculus of [5]. Though being equivalent to the strong form (4), this weak form has the advantage of being directly usable to shift the dynamics from the continuous to our piece-wise discrete approach. In fact, to derive the discrete dynamics corresponding to the discrete kinematics (10), it suffices to introduce the relation: $\delta\zeta(X)=J(X)\overrightarrow{\delta\xi}$ in addition to the kinematics relations (19) and (21). In these conditions, (23) becomes:

[TABLE]

which leads to the following generalized dynamics equation once the external loads (7) and cable driven actuation (5) have been introduced in (24):

[TABLE]

where, when needed, $n$ represents the section corresponding to the running value of $X$ inside the integrals and $|\cdot|_{\mathrm{v}}$ takes the norm of the translational part of the operand according to equation (7).

In the remaining part of the section we will describe the different components of (25), let us start with the internal elastic and actuation load, those loads are traditionally called

[TABLE]

Due to the linearity of the integral, each element $\tau_{n}$ has the form:

[TABLE]

where we note that by definition $S_{n}(X)=0_{6}$ for $X\leqslant L_{n-1}$ (Fig. 3). Each of the integrals in the series except of the first can be directly solved analytically making use of the identity $\mathrm{Ad}_{g}^{*}\left(\mathcal{F}^{\prime}+\mathrm{ad}_{\xi}\mathcal{F}\right)=\left(\mathrm{Ad}_{g}^{*}\mathcal{F}\right)^{\prime}$, while the first one can be analytically solve with an integration by part with the additional use of the identity $\mathrm{AD}_{g}^{T\prime}=\mathrm{Ad}_{g}^{T}=\mathrm{Ad}_{g^{-1}}^{*}$. Applying this operations, we obtain the internal elastic and actuation load for the section $n$ as follows.

[TABLE]

where $l$ is the length of the section equal to $L_{n}-L_{n-1}$ and we have assumed elastic and actuation loads constant along the section, i.e., $\mathcal{F}_{a}(L_{n-1}<X<L_{n})=\mathcal{F}_{a}$ and $\mathcal{F}_{i}(L_{n-1}<X<L_{n})=\mathcal{F}_{i}$ are constants.

In order to calculate the sum in (26), we exploit the boundary condition at each section. For the cable-driven actuation case they are given below.

[TABLE]

where the cables are assumed to run from the point of anchorage to the base of the manipulator. The contribution of the cables attached at $L_{n}$ is indicated with $\mathcal{F}_{an}$ and the constant internal load of the section $n$ with $\mathcal{F}_{in}$. As expected, crossing an anchoring edge $L_{n}$ causes a jump in both the internal elastic and actuation load due respectively to the concentrated load of the cables anchored at that position and the suddenly increase of the number of cable running through the section. Substituting (27) into (26), results in a brutal cancellation of the first term (the sum), which becomes:

[TABLE]

For what concern the fluidic actuation case, the boundary condition are as follows.

[TABLE]

where we have taken into account the load exerted at the bottom of the section (Fig. 1) in the jump from $L_{n}^{+}$ to $L_{n}^{-}$ and the fact that there is no distributed load along the section. Substituting (29) into (26), results in a cancellation of the elastic load in the first term and of the actuation load in the second term, which leads to:

[TABLE]

with $\mathcal{F}_{a(N+1)}=0_{6\times 1}$

The second term on the right end side of equation (25) represents the generalized external concentrated load and is usually referred to as

[TABLE]

where each elements is simply:

[TABLE]

Finally, with those definition at hand and naming the coefficients matrices in squared parenthesis of (25), we obtain the piece-wise constant strain dynamic equation:

[TABLE]

where we recognize the structure of the Lagrangian model of rigid serial manipulators.

Let us now break down each matrix coefficients of the dynamic equation (32). Looking at (25), the mass matrix $M\left(\overrightarrow{\xi}\right)\in\mathbb{R}^{6N\times 6N}$ is a symmetric, positive define matrix and his $6\times 6$ block-element of block-row $n$ and block-column $m$ is calculated as follows.

[TABLE]

where we have exploited the fact that for $i<max(n,m)$ or equivalently $X<L_{max(n,m)-1}$ either $S_{n}(X)$ or $S_{m}(X)$ is equal to $0_{6}$ (Fig. 3). Similarly, for the Coriolis matrices $C_{1}\left(\overrightarrow{\xi},\dot{\overrightarrow{\xi}}\right)\text{, }C_{2}\left(\overrightarrow{\xi},\dot{\overrightarrow{\xi}}\right)\in\mathbb{R}^{6N\times 6N}$ we obtain

[TABLE]

[TABLE]

while, for the drag matrix $D\left(\overrightarrow{\xi},\dot{\overrightarrow{\xi}}\right)\in\mathbb{R}^{6N\times 6N}$, we get:

[TABLE]

With the same reasoning, the block-element of block-row $n$ of the gravitational-buoyancy matrix $N\left(\overrightarrow{\xi}\right)\in\mathbb{R}^{6N\times 6}$ is as follows

[TABLE]

We have now all the ingredients to process the joint dynamic (32) and reconstruct the shape, velocity and acceleration of the soft manipulator with (8), (19) and (21).

IV Simulation Results

In this section the PCS dynamic model (32) is tested through different simulations. First, a plane motion of three sections manipulator is shown, then an out-of-plane motion involving the torsion of all the three sections, which is not possible with the PCC model, is performed. Finally, in order to show how the PCS model copes with non constant external load, the model is compared with a cantilever beam, simulated using the continuous Cosserat model (4). Before that, an efficient recursive algorithm aiming to calculate the coefficient matrices of (32) is presented.

IV-A Recursive Algorithm

The basic idea for the recursive algorithm is that each section of the soft manipulator contributes to a very specific set of block-elements of the coefficient matrices. In particular, the non zero block-elements due to section $n$ of the mass matrix $M$, the first Coriolis matrix $C_{1}$ and the drag matrix $D$ are those located in the square block-matrix of block-rows 1 to $n$ and block-column 1 to $n$, while for the gravitational matrix $N$, they compose the block-rows from 1 to $n$ of the only block-column and finally, for the second Coriolis matrix $C_{2}$, the non zero block-elements form a rectangular block-matrix of block-rows 1 to $n$ and block-column 1 to $n-1$ (Figure 4). This can be seen by splitting the integrals in the coefficient matrices of equation (25) into the $N$ integrals corresponding to each section, then the non zero square block-matrix rise from the varying structure of the Jacobian $J(X)$ (and the modified Jacobian for $C_{2}$) with respect to $X$ as shown in figure 3.

An efficient way to implement this technique, is to benefit from the results of the calculations given by the previous section. Going further into details, at a certain $L_{n-1}<X\leqslant L_{n}$ the Jacobian element $S_{n-1}(X)$ has only one member which actually depends on $X$, all the rest being inherited from the last evaluation of the same quantity in the previous section $\left(S_{n-1}(L_{n-1})\right)$, as it can be visualized by inspecting equation (20). This is used in the calculation of $J(X)$ and the adjoint elements of $C_{2}$ (after multiplication with $\dot{\xi}_{n-1}$). Furthermore, $\eta(L_{n-1})$ and $g(L_{n-1})$ are calculated respectively through equation (13) and (10) to obtain the co-adjoint member in $C_{1}$ and the Adjoint member in $N$.

IV-B Plane & Out-of-Plane Motion

In this section the feasibility of the model to perform highly dynamic motion both in plane and out-of-plane is shown. The simulated soft manipulator is composed by three cylindrical sections of length $l=250$ $mm$, radius equal to $10$ $mm$, Young modulus $E=110$ $kPa$, shear viscosity modulus $\upsilon=0.3$ $kPasec$, Poisson modulus equal to $0.5$ and mass density $\rho=1080$ $kg/m^{3}$. The manipulator lays upside-down as shown in figure 2 and shares the $Z$ axis with the inertial frame, therefore the transformation map between the spatial frame and the base frame is: $g_{r}=diag(-1,-1,1,1)$. The actuation load in both cases is fluidic and imposed over time through a ramp starting from zero, with a 1 second width, and reaching in the planar case:

[TABLE]

and in the out-of-plane case:

[TABLE]

where all the value are in $10^{-3}Nm$. The gravity load has been neglected.

Few snapshots of the plane motion are shown in figure 5, while the out-of-plane motion is shown in figure 6. The screws associated with the last configuration of the out-of-plane motion are also shown in figure 6. Using the terminology of screw [13], the colored arrows represent the axis $\hat{a}_{n}$ of the three screws, around which the sections rotate of an amount equal to the magnitude $m_{n}$, while the black arrows indicate the amount of translation in the direction of the screw, given by $h_{n}m_{n}$, $h_{n}$ being the pitch of the screw. For this particular configuration we obtained:

[TABLE]

IV-C Cantilever Beam Comparison

Even if one could manage to design the actuation of a soft manipulator in order to be constant in each section, non-constant loads due to gravity, external forces and inertial forces are unavoidable in realistic condition. For this reason, it is important to know how the discrete Cosserat model handles a non-constant load and what are the effects of such loads in terms of accuracy of the result. To do so, we notice that each member of the dynamic equation (32) is pre-moltiplied by $J^{T}(X)$ or, block-element-wise, by $S_{n}^{T}(X)$. Now, looking at the Jacobian (20), we see that an element $S_{n}^{T}(X)$ will first map the considered load to the base of the section $n$ through $\mathrm{Ad}^{*}$ and then integrate all the ”re-maps” of this load up to $X$ through

[TABLE]

which gives, by definition of integral, $l$ times the mean of this load on the section $n$. It is worth to highlight here, that this is in essence how the discerete Cosserat model relates the ideal assumption of piece-wise constant strains with the real continuously varying counterpart.

Intuitively, the wider is the interval in which the mean is evaluated the larger is the discrepancy with the real distribution. In order to show this fact and test the model with a non-constant load scenario, the continuous Cosserat model (4) and the discrete Cosserat model (32) are applied in the following to a cantilever beam with vertical tip load. The simulated beam is of cilindrical shape, with length $L=250$ $mm$, radius equal to $10$ $mm$, Young modulus $E=110$ $kPa$, shear viscosity modulus $\upsilon=0.3$ $kPasec$, Poisson modulus equal to [math] and mass density $\rho=2000$ $kg/m^{3}$. The beam lays on the right side of the inertial frame (($e_{1},e_{2},e_{3})$ in figure 2) and share the $Z$ axis with this frame, therefore the map between the spatial frame and the base frame is:

[TABLE]

Finally, the external tip load points in the positive $y$ direction with respect to the fixed base frame and is applied at $\bar{X}=L$, thus it has the form:

[TABLE]

with unit reference of $10^{-3}N$.

In figure 7, on top, is shown the resulting curvature of the continuous cantilever as a function of space and time, followed by the curvatures of the discrete cantilever divided in one, two and three sections (blue lines). We immediately notice that the oscillation frequency for the one section case is much higher than that of the continuous cantilever. This discrepancy is quantified at each time by the tip position error expressed in percentage of the total length $L$ (red markers, in both the $e_{1}$ and $e_{2}$ directions). Intuitively, this can be explained by the fact that the additional constraint of constant strains applied to the Cosserat micro-solids in the discrete model makes the system more rigid. As expected, the oscillation frequency progressively slows down toward the continuous value with the increase of the number of sections while, accordingly, the error gradually decreases.

The additional rigidity due to the constant strain constraint is confirmed by the steady state comparison. As a matter of fact, in all the three cases with one, two and three sections, the steady state tip position is above the real cantilever tip position, or in other words the beam is less deformed. Again, the error decreases in both the $e_{1}$ and $e_{2}$ directions with the increase of sections. The respective steady-state error values are shown below.

[TABLE]

V Experimental Results

In this section the PCS dynamic model performances are compared against experimental data. In order to evaluate the results with respect to the continuous Cosserat model, we have used the same prototype, parameters and experimental data provided in [19], which we refer to for more exhaustive details on the experimental platform and measurement set up.

In short, the prototype is composed of a single conical piece of silicone, with a base radius $R_{max}$ and a tip radius $R_{min}$, actuated by 12 cables embedded inside the robot body. The cables run parallel to the midline at a distance $d_{j}$ ($j\in\{1,2,\dots,12\}$) and are anchored four at a time at three different lengths along the robot arm ($L_{1}$, $L_{2}$, $L_{3}$) and with a relative angle of 90 degrees (Figure 8). During operation, the cable tensions, driven by servomotors, are measured by force sensors while the motion of the arm is recorded with two high speed cameras. The 3-D motion is then reconstructed through a process based on the direct linear transformation (DLT).

The soft manipulator has been tested for three different conditions, a single bending motion produced by cable 11, an in-plane multi-bending produced by a sequence of activation of cables 9, 11, 1, 3 and an out-of-plane multi-bending produce by cables 11, 5, 2. The details of the cable activaion is reported in Figure 9 (top three graphs).

In [19], it has been found that the drag and added mass matrix in this case can be expressed as

[TABLE]

where $\mathsf{D}(X)=diag(1/2\pi RC_{x},RC_{y},RC_{z})\rho_{w}$ and $\mathsf{F}(X)=diag(0,AB_{y},AB_{z})\rho_{w}$, $R(X)$ being the radius of the soft arm and $C_{x}$, $C_{y}$, $C_{z}$, $B_{y}$, $B_{z}$ being fluid dynamics coefficients. The mechanical and geometrical parameters of the arm are summarized in Table I.

V-A Comparison

In order to exploit the dynamics equations developed in III-B, the soft manipulator has been modeled as a stack of four cilindrical constant-strain sections defined by $L_{1}$, $L_{2}$, $L_{3}$, $L_{4}$, with a radius equal to the mean of the prototype radius for each section ($R_{1}$, $R_{2}$, $R_{3}$ and $R_{4}$ in Figure 8). The dynamics and kinematics equations have been solved by implementing the recursive algorithm presented above for the three load conditions of the experiments. The results of the tests are reported in Figure 9 together with the one obtained in [19] with the continuous Cosserat model. The error is calculated as the normalized mean at each time step of the Euclidean distance between the simulated and real markers positioned a the tip of each section.

V-B Discussion

As can be observed from Figure 9, the results of the present discrete model are comparable or even better than the ones obtained with the continuous model, as a matter of fact, the average error in the three cases are respectively 5.1%, 5.2% and 5.4% for the continuous model versus 5.7%, 4.2% and 4.3% for the discrete model. In the authors opinion, the reason for this improvement lays mainly on the higher numerical stability shown by the PCS discrete model, in particular, on the different management of the internal point load exerted where the cables are fastened. In the continuous model an internal point load is modeled with a Dirac function that has to be discretized during the numerical integration, while in the PCS model there is no such approximation and the concentrated load is introduce naturally with the boundary conditions for each section.

In order to further improve the accuracy of the model, the friction of the cables should be included which in turn models the hysteresis behavior of the load-unload cycle. As it is highlighted in [19], the hysteresis behavior is clear from the plane bending experiments, in which the error increase drastically after the relaxation of the cable which is when the load is mainly driven by the friction of the cable against the silicone body. From a geometric point of view, the number of discrete sections could be increased in order to capture non-negligible variation of the strain due to external loads. Furthermore, the model could be able to better take into account the variation with respect to $X$ of the mass $\mathcal{M}(X)$ in the calculus of the mass matrix $M$ as well as the variation of the stiffness and viscosity matrices $\Sigma(X)$, $\Upsilon(X)$ in the calculus of the internal elastic load of $\overrightarrow{\tau}$, which are due to the conical shape of the manipulator.

VI Conclusion

In conclusion, a new piece-wise constant strain model for multi-section soft robots has been presented which is based on the discretization of the continuous Cosserat model inheriting from it the fruitful geometrical and mechanical properties. The close relation between this model for soft robotics and the traditional model for rigid robotics is also highlighted. The PCS model has been extensively corroborated through simulation and experimental results of plane and out-of-plane multi-bending. Furthermore, the performances have been compared with the continuous Cosserat model showing comparable or even better results. It is worth to highlight, that a similar approach can be found in the context of recent finite element formulation for geometrically exact beam as in [26], making of this work a bridge between different engineering disciplines.