# Structural Dynamics of a Pulsed-Jet Propulsion System for Underwater Soft Robots

Federico Renda, Francesco Giorgio Serchi, Frédéric Boyer, Cecilia Laschi

► To cite this version:

Federico Renda, Francesco Giorgio Serchi, Frédéric Boyer, Cecilia Laschi. Structural Dynamics of a Pulsed-Jet Propulsion System for Underwater Soft Robots. International Journal of Advanced Robotic Systems, 2015, 12, pp.68. 10.5772/60143. hal-01201680

HAL Id: hal-01201680

https://imt-atlantique.hal.science/hal-01201680v1

Submitted on 8 Oct 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Structural Dynamics of a Pulsed-Jet Propulsion System for Underwater Soft Robots

Regular Paper

Federico Renda<fnref n="1" />,<fnref n="2" />*, Francesco Giorgio Serchi<fnref n="2" />, Frederic Boyer<fnref n="3" /> and Cecilia Laschi<fnref n="2" />

Received 30 October 2014; Accepted 03 February 2015

DOI: 10.5772/60143

© 2015 Author(s). Licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## Abstract

This paper entails the study of the pulsed-jet propulsion inspired by cephalopods in the frame of underwater bioinspired robotics. This propulsion routine involves a sequence of consecutive cycles of inflation and collapse of an elastic bladder, which, in the robotics artefact developed by the authors, is enabled by a cable-driven actuation of a deformable shell composed of rubber-like materials. In the present work an all-comprehensive formulation is derived by resorting to a coupled approach that comprises of a model of the structural dynamics of the cephalopod-like elastic bladder and a model of the pulsed-jet thrust production. The bladder, or mantle, is modelled by means of geometrically exact, axisymmetric, nonlinear shell theory, which yields an accurate estimation of the forces involved in driving the deformation of the structure in water. By coupling these results with those from a standard thrust model, the behaviour of the vehicle propelling itself in water is derived. The constitutive laws of the shell are also exploited as control laws with the scope of replicating the muscle activation routine observed in cephalopods. The model is employed to test various shapes, material properties and actuation routines of the mantle. The results are

compared in terms of speed performance in order to identify suitable design guidelines. Altogether, the model is tested in more than 50 configurations, eventually providing useful insight for the development of more advanced vehicles and bringing evidence of its reliability in studying the dynamics of both man-made cephalopod-inspired robots and live specimens.

**Keywords** Dynamics, Continuum Robots, Soft Robots, Biologically-Inspired Robots

## Nomenclature

* 0 Variable in the reference configuration.

* Derivative with respect to time $t$.

* ' Derivative with respect to $X$.

* Converts $\mathbb{R}^6$ in $se(3)$.

* Converts $\mathbb{R}^3$ in $so(3)$.

* $t \in \mathbb{R}$ Time.

<footnote n="1">1 Khalifa University Robotics Institute, Khalifa University of Science Technology and Research (KUSTAR), Abu Dhabi, United Arab Emirates</footnote>

<footnote n="2">2 The BioRobotics Institute, Scuola Superiore Sant'Anna, Pisa, Italy</footnote>

<footnote n="3">3 Research Institute of Communication and Cybernetics of Nantes, Ecole des Mines de Nantes, France</footnote>

<footnote>*Corresponding author(s) E-mail: federico.renda@kustar.ac.ae</footnote>
• $X \in \mathbb{R}$ Reference arc-length parametrization.

• $\phi \in S^1$ Angle of revolution.

• $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3) \in \mathbb{R}^3 \times \mathbb{R}^3 \times \mathbb{R}^3$ Ambient reference frame.

• $(\mathbf{e}_r, \mathbf{e}_\phi, \mathbf{e}_3)(X, \phi, t) \in \mathbb{R}^3 \times \mathbb{R}^3 \times \mathbb{R}^3 \text{ Frame attached to the point } (X, \phi)$.

• $(\mathbf{a}, \mathbf{b}, -\mathbf{e}_\phi)(X, \phi, t) \in \mathbb{R}^3 \times \mathbb{R}^3 \times \mathbb{R}^3 \text{ Director orthogonal frame.}$

• $g(X, \phi, t) \in SE(3)$ Configuration matrix.

• $\xi(X,t) \in se(3)$ Local deformation twist vector.

• $\eta(X,t) \in se(3)$ Local velocity twist vector.

• $\mathcal{F}(X,t) \in se(3)*$ Wrench vector.

• $R(X, \phi, t) \in SO(3)$ Orientation matrix.

• $\mathbf{r}(X, \phi, t) \in \mathbb{R}^3$ Position vector.

• $\mathbf{g}(X,t) \in \mathbb{R}^3$ Linear strain.

• $\mathbf{k}(X,t) \in \mathbb{R}^3$ Curvature vector.

• $\mathbf{v}(X,t) \in \mathbb{R}^3$ Linear velocity.

• $\mathbf{w}(X,t) \in \mathbb{R}^3$ Angular velocity vector.

• $c(X,t) \in \mathbb{R}^3$ Added mass load.

• $\mathbf{r}(X,t) \in \mathbb{R}^3$ Radius.

• $z(X,t) \in \mathbb{R}$ Altitude.

• $\theta(X,t) \in S^1$ Fibre angle.

• $\lambda(X,t) \in \mathbb{R}^+$ Tangential strain.

• $\beta(X,t) \in \mathbb{R}$ Thickness strain.

• $\mu(X,t) \in \mathbb{R}$ Curvature function.

• $v_a(X,t) \in \mathbb{R}$ Tangential velocity.

• $v_b(X,t) \in \mathbb{R}$ Perpendicular velocity.

• $w(X,t) \in \mathbb{R}$ Angular velocity function.

• $N_X(X,t) \in \mathbb{R}$ Internal tangential force along a.

• $H(X,t) \in \mathbb{R}$ Internal thickness force.

• $M_X(X,t) \in \mathbb{R}$ Internal torque along $-\mathbf{e}_\phi$.

• $N_\phi(X,t) \in \mathbb{R}$ Internal tangential force along $-\mathbf{e}_\phi$.

• $M_\phi(X,t) \in \mathbb{R}$ Internal torque along a.

• $f_a(X,t) \in \mathbb{R}$ External tangential force.

• $f_b(X,t) \in \mathbb{R}$ External perpendicular force.

• $l(X,t) \in \mathbb{R}$ External torque.

• $u(t) \in \mathbb{R}$ Mantle position.

• $v(t) \in \mathbb{R}$ Mantle velocity.

# 1. Introduction

Cephalopods, such as squids and octopuses, often propel themselves in water by resorting to a sequence of cyclic contractions and expansions of a soft cavity of their body, commonly referred to as the mantle [1]. While squid may rely on fin-assisted swimming [2] and octopuses are observed to use arm sculling [3], here we will be dealing exclusively with the pulsed-jet mode of propulsion. During each pulsation cycle, the mantle inflates, ingesting ambient water, and then abruptly contracts, thus expelling a slug of water which, by reaction, generates the forward thrust [4]. The contraction of the mantle is ensured by a network of circular muscles symmetrically arranged all around the body [5]. Like other aquatic animals, cephalopod swimming locomotion occurs by discontinuous bursts of acceleration associated with the expulsion of water from the mantle cavity [6]. However, cephalopod locomotion is in many ways different from finned or caudal flapping fish-like propulsion for a number of reasons. First, the very nature of the slug of fluid expelled across the nozzle is known to give rise to a vortex ring which, in turn, is shown to provide a critical contribution to the production of thrust [7]. The benefit provided by pulsed-jet propulsion has gained the attention of several research groups, eventually leading to the design and development of a number of underwater vehicles that rely on this locomotion strategy (e.g., [8, 9, 10]). In addition, the shape change associated with the shrinkage of the mantle during pulsation was found to participate in further increasing the thrust generated via the re-capture of the energy from the ambient flow [11]. Because of these peculiar features, lately cephalopods have represented an important source of inspiration for the development of a new kind of underwater thrusters composed of soft materials. A first example of a pulsed-jet propelled soft robot inspired by the octopus is presented in [12], where an underwater thruster actuated via cable transmission is described (fig. 1). The prototype presented in [12] and later revised in [13] and [14] is capable of propelling itself in water by performing a thrust production routine analogous to that of cephalopods. These vehicles are composed of an elastic, hollow shell, somewhat similar to a bladder which undergoes periodic phases of collapse and inflation. During these stages the vehicle respectively expels fluid and successively refills the bladder, thus performing a sequence of pulsed-jets in a closely resembling fashion to what cephalopods do when swimming [15]. By propelling itself via the actual collapse of the collapsible bladder, the cephalopod-inspired underwater vehicles not only benefit from the advantages provided by vortex-ring aided thrust production, [16], but they also capitalize on the positive feedback that added-mass recovery has on thrust during shape change, [11]. The
realization of this new kind of vehicle was enabled by
addressing new design solutions in order to account for the
use of soft materials as well as introducing new models and
control strategies (e.g., [17] and [13]). The existing models
of the actuation mechanism of the soft-bodied thruster of
[12] is based on highly simplified assumptions and essen-
tially relies on an apriori kinematics of the shell deforma-
tion, see [13]. However, the need has arisen to evaluate the
mechanics inherent to the deformation of the elastic shell
throughout the pulsation sequence. As an example, it is of
major concern to quantify the stresses generated within the
shell thickness in order to predict the speed of passive
inflation of the mantle after the fluid expulsion phase. This
not only yields invaluable information for the improved
mechanical design and choice of the constituent materials,
but also aids in determining the optimal pulsating routines
and hence defining optimized control laws for the vehicle.

<figure><img src="image_2.png" /><figcaption>Figure 1. The first cephalopod-inspired, jet-propelled soft robot [12]: (a) side, (b) frontal and (c) underneath view of the prototype. In (a), the numbers respectively refer to: 1. the outflow nozzle, 2. the ingestion valve (operating like the pallial valves of cephalopods), 3. the artificial, silicone mantle of the vehicle and 4. the power supply wire.</figcaption></figure>

In order to study this mode of swimming mechanisms in
detail and translate these on to the robotic artefacts, we
need to state the locomotion models of these jet-propelled
soft animals as well as derive quantitative tools for captur-
ing the mechanics of the deformation of the soft tissues
during propulsion. Ideally, these models would be analyt-
ical models allowing us to idealize the principles of this
mode of locomotion. Here we propose a numerical cou-
pling between a structural model of the mantle, based on
the framework of the geometrically exact theory of shells
in finite transformations [18] [19], and a model for thrust
production via jet propulsion based on the analogy of this

mode of thrust production with that of rockets, as was done
earlier by [20].

In the geometrically exact approach the shell is taken as a
Cosserat medium, i.e., a continuous assembly of rigid
micro-solids whose rigid overall motions are considered
from the beginning without any approximations [21]. This
aspect is crucial to tackle locomotion problems where we
not only need to model the internal strains occurring in the
body but also the net rigid overall motions in space. While
in Cosserat beams the micro-solids are the cross sections of
the beam to each of which a rigid reference frame (three
unit vectors) is attached, in the case of shells, the micro-
solids are rigid fibres transversally attached to the mid-
surface of the shell to which one can only attach a single
vector referred to as the 'director' [22]. Based on this
Cosserat model, several internal kinematics can be adopted
depending on whether the shell is thin or not. In the first
case, the directors remain perpendicular to the mid-surface
while in the second case they can rotate freely with respect
to the mid-surface with two additional degrees of freedom
which, in turn, induce two further strain fields named
'transverse shearing'. The first kinematics correspond to the
so-called Kirchhoff model of shells while the second
correspond to the Reissner model [23]. Geometrically exact
beam theories have been recently applied to continuous
(hyper-redundant) and soft robotics in the context of
underwater and terrestrial locomotion of fish [24] and
snakes [25] and for manipulation of octopus like arms [26].
In this new article, we use the geometrically exact shell
model to address the issue of cephalopod jet propelling.
The model of the mantle proposed here is based on this
second theory which leads to partial differential equations
of minimum order. Furthermore, taking inspiration from
actual cephalopods, the shell will be taken as axisymmetric
[20]. As it is customary in nonlinear structural dynamics,
the model is derived by first defining the shell kinematics.
From these kinematics, we will build a set of strain meas-
ures and will derive the dynamic balance equations in their
Cauchy form, i.e., in terms of internal stresses. Finally, this
picture will be completed in terms of the constitutive laws
in section 2.5. All these developments will be achieved in
the case of an axisymmetric shell. Afterwards, the jet
propulsion model will be described and the coupling
between the two models will be formalized. This represents
a first step in modelling the complex fluid-structure
interaction problem, which will be further investigated in
a future work. The coupled model is then employed for
simulating a range of scenarios where a variety of mantle
shapes and activation sequences are taken into account and
the results from these simulations are illustrated.

## 2. Mantle Model

The mantle is a piece of tissue forming a cavity which opens
into ambient water through an orifice (see Figure 2). In the
remainder of this article, this will be modelled as an
axisymmetric shell of thickness $2h$ composed of a soft hyperelastic isotropic material of density $\rho$.

Adopting the Reissner model (that takes into account shear deformation through the thickness of the shell), the configuration space of the shell can be first defined by:

$$C := \{(X^1, X^2) \in B \mapsto (\mathbf{r}, \mathbf{b}) \in \mathbb{R}^3 \times \mathbb{S}^2\} \quad (1)$$

where $B \subset \mathbb{R}^2$ defines the material domain of the mid-surface, $\mathbf{r}$ represents the field of position of the points of $B$ and $\mathbf{b}$ stands for the field of unit vectors attached to the shell fibres, i.e., the directors living in the two dimensional unit sphere $\mathbb{S}^2$. In other words, every point of a geometrically exact Reissner shell can translate in space ($\mathbf{r}$) and perform both pitch and roll ($\mathbf{b}$). As a result of the degree of symmetry of the body accounted for, the configuration space can be reduced further. In the next section, details of this reduction is argued step-by-step.

## 2.1 Kinematics of Axisymmetric Shells

Mathematically, an axisymmetric surface or 'surface of revolution', is obtained by rotating a planar curve or 'profile' around a fixed axis named the symmetry axis. This rotation changes the 'profile curve' into any of the meridian curves that constitute the shell. The ambient Euclidean space is endowed with a fixed base of orthogonal unit vectors ($\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$), where $\mathbf{e}_3$ supports the symmetry axis. Denoting by $\phi$ the angle of revolution and by $X$ the material coordinate along the meridian curves, the orthogonal basis ($\mathbf{e}_r, \mathbf{e}_\phi, \mathbf{e}_3$) fixed to the material point of abscissa $X$ and meridian $\phi$ is defined by:

$$g_1 \in SE(3):$$

$$(\phi, X) \in [0, 2\pi[ \times [0, L] \mapsto g_1(\phi, X) = \begin{pmatrix} \exp(\tilde{\mathbf{e}}_3\phi) & \mathbf{r} \\ 0 & 1 \end{pmatrix}$$

where $\exp$ is the exponential in $SO(3)$ and the tilde is the usual isomorphism between a vector of $\mathbb{R}^3$ and the corresponding skew-symmetric matrix. Then

$$\exp(\tilde{\mathbf{e}}_3\phi) = \begin{pmatrix} \cos(\phi) & -\sin(\phi) & 0 \\ \sin(\phi) & \cos(\phi) & 0 \\ 0 & 0 & 1 \end{pmatrix}$$

is a rotation $\phi$ around the axis $\mathbf{e}_3$ and $\mathbf{r}: (X,t) \in (0, L) \times [0, \infty) \to \mathbf{r}(X,t) = (\cos(\phi)r(X,t), \sin(\phi)r(X,t), z(X,t))^T$ is a curve that define the profile of the shell and for which, $r(\cdot)$ and $z(\cdot)$ are two smooth functions that define the radius and the

altitude of the point $X$ on the profile (see Figure 2). For the sake of convenience, we introduce another reference frame ($\mathbf{e}_r, \mathbf{e}_3, -\mathbf{e}_\phi$) adding the transformation

$$g_2 \in SE(3): g_2 = \begin{pmatrix} \exp(\tilde{\mathbf{e}}_r, \pi/2) & 0 \\ 0 & 1 \end{pmatrix}$$

Now let us call $\theta(X,t)$ the angle between $\mathbf{e}_3$ and the shell fibre located at any $X$ along the $\phi$-meridian, then the so-called *director orthogonal frame*($\mathbf{a}, \mathbf{b}, -\mathbf{e}_\phi$) is defined at each instant $t$, by:

$$g_3 \in SE(3): (X) \in [0, L] \mapsto$$

$$g_3(X) = \begin{pmatrix} \exp(-\tilde{\mathbf{e}}_\phi \theta) & 0 \\ 0 & 1 \end{pmatrix}$$

Finally, putting them all together, the shell configuration matrix is

$$g \in SE(3): g(X, \phi) = g_1 g_2 g_3 = \begin{pmatrix} R & \mathbf{r} \\ 0 & 1 \end{pmatrix}$$

where $R = \exp(\tilde{\mathbf{e}}_3\phi)\exp(\tilde{\mathbf{e}}_r\pi/2)\exp(-\tilde{\mathbf{e}}_\phi\theta)$ has been defined.

As a result, we can now introduce the following definition of the configuration space of an axisymmetric shell:

$$\bar{C} := \{(X, \phi) \in \bar{B} \mapsto g \in SE(3)\} \quad (2)$$

where, referring to the more general context of (1), $X^1 = X$ and $X^2$ have been replaced by $\phi(\bar{B} \subset \mathbb{R} \times S^1)$ with $X^2 = r^o(X)\phi$, and $r^o(X)$ the value of $r(X)$ in the reference configuration of the shell (before any deformation).

The tangent plane on the surface $g(X, \phi, t)$ is represented by two vector field: $\hat{\xi}_1(X, t) = g^{-1}\partial g / \partial X = g^{-1}g'$ and $\hat{\xi}_2(X, t) = g^{-1}\partial g / \partial(r^o\phi) = g^{-1}g^l$. The hat represents the isomorphism between the twist vector space $\mathbb{R}^6$ and the Lie algebra $se(3)$. Below their components are specified:

$$\hat{\xi}_{1} = \begin{pmatrix} \tilde{\mathbf{k}}_{1} & \mathbf{g}_{1} \\ 0 & 0 \end{pmatrix} \in se(3)$$

$$\hat{\xi}_{2} = \begin{pmatrix} \tilde{\mathbf{k}}_{2} & \mathbf{g}_{2} \\ 0 & 0 \end{pmatrix} \in se(3)$$

$$\xi_1 = (\mathbf{k}_1^T, \mathbf{g}_1^T)^T = (0, 0, \mu, \lambda, \beta, 0)^T \in \mathbb{R}^6$$
$$ \xi_2 = (\mathbf{k}_2^T, \mathbf{g}_2^T)^T = \left( \frac{\sin(\theta)}{r^o}, \frac{\cos(\theta)}{r^o}, 0, 0, 0, \frac{-r}{r^o} \right)^T \in \mathbb{R}^6 $$

where the curvature $\mu$ and the longitudinal and transversal strains ($\lambda, \beta$) have been defined as:

$$ \mu = \theta' \\ \lambda = \cos(\theta)r' + \sin(\theta)z' \\ \beta = \cos(\theta)z' - \sin(\theta)r' $$

The time evolution of the configuration curve $g$ is represented by the twist vector field $\eta(X,t) \in \mathbb{R}^6$ defined by $\hat{\eta} = g^{-1}\partial_g / \partial t = g^{-1}\dot{g}$. Let us specify the component of

$$ \hat{\eta} = \begin{pmatrix} \tilde{\mathbf{w}} & \mathbf{v} \\ 0 & 0 \end{pmatrix} \in se(3) \\ \eta = (\mathbf{w}^T, \mathbf{v}^T)^T = (0, 0, w, v_a, v_b, 0)^T \in \mathbb{R}^6 $$

where the angular velocity $w$ and the longitudinal and transversal velocity ($v_a, v_b$) have been defined as:

$$ w = \dot{\theta} \\ v_a = \cos(\theta)\dot{r} + \sin(\theta)\dot{z} \\ v_b = \cos(\theta)\dot{z} - \sin(\theta)\dot{r} $$

Finally, from $\hat{\eta} = g^{-1}\dot{g}$, the kinematic equation is

$$ \dot{g} = g\hat{\eta} \qquad (3) $$

that outlining the components becomes:

$$ \begin{aligned} \dot{\theta} &= w \\ \dot{r} &= \cos(\theta)v_a - \sin(\theta)v_b \\ \dot{z} &= \sin(\theta)v_a + \cos(\theta)v_b \end{aligned} \qquad (4) $$

## 2.2 Compatibility Equations

It has been shown above that $g' = g \hat{\xi}_1$. By taking the derivative of this equation with respect to time and recalling that $\dot{g} = g \hat{\eta}$, we obtain the following compatibility equation

<figure><img src="image_3.png" /><figcaption>Figure 2. Axisymmetric shell kinematics (left). Profile view, beam-like parametrization (right).</figcaption></figure>

between the velocity and the deformation variables: $\hat{\xi}_1 = \eta' + \hat{\xi}_1 \hat{\eta} - \hat{\eta} \hat{\xi}_1$. In terms of twist vectors this can be written as:

$$ \dot{\xi}_1 = \eta' + ad_{\xi_1}(\eta) \qquad (5) $$

where

$$ ad_{\xi_1} = \begin{pmatrix} \tilde{\mathbf{k}}_1 & 0 \\ \tilde{\mathbf{g}}_1 & \tilde{\mathbf{k}}_1 \end{pmatrix} $$

is the adjoint map that represents the action of the Lie algebra on itself.

Outlining the components we obtain:

$$ \begin{aligned} \dot{\mu} &= w' \\ \dot{\lambda} &= v_a' + \beta w - \mu v_b \\ \dot{\beta} &= v_b' - \lambda w + \mu v_a \end{aligned} \qquad (6) $$

## 2.3 Strain Measures

First, let us introduce two quadratic forms named: $h$ and $k$, that represent the first and the second fundamental forms of a Reissner shell, respectively.

The first fundamental form of a surface is a quadratic form that determines how the Euclidean metric of $\mathbb{R}^3$ is induced on the surface in any of its points. In our case, this is defined in each point $(X, \phi)$ of the surface by a tensor the components of which are deduced from the scalar products of all the vectors of the field of basis (said natural basis) $(g_1, g_2)(X, \phi)$. For our revolution surface it can be simply expressed as:
$$h = \begin{pmatrix} h_{11}(X) & 0 \\ 0 & h_{22}(X) \end{pmatrix} = \begin{pmatrix} \lambda^2 + \beta^2 & 0 \\ 0 & \frac{r^2}{r^{o2}} \end{pmatrix}$$

which, due to the axisymmetry, depends only on X.

While the first fundamental form defines the scalar product of any tangent vectors to the mid surface, the second fundamental form defines the curvature of the surface in any tangent direction (defined by an unit tangent vector). For a Reissner shell the components of k(X) are such that $k_{\alpha\gamma} = g_{\alpha} \cdot \tilde{k}_{\gamma} \cdot b$.

Contrary to the traditional first and second fundamental form of a surface, this 'special' form takes into account the effect of the shear between two material elements, where the two overlap if the director b points in the direction normal to the mid-surface of the shell (i.e., no shear strain).

For an axisymmetric Reissner shell, we have:

$$k = \begin{pmatrix} k_{11}(X) & 0 \\ 0 & k_{22}(X) \end{pmatrix} = \begin{pmatrix} -\mu\lambda & 0 \\ 0 & -\frac{r \sin(\theta)}{r^{o2}} \end{pmatrix}$$

In accordance with [18], the following strain tensor field has been adopted:

$$e = \frac{1}{2}(h - h^o) = \frac{1}{2} \begin{pmatrix} \lambda^2 + \beta^2 - 1 & 0 \\ 0 & \frac{r^2}{r^{o2}} - 1 \end{pmatrix}$$

which describes the membrane strain state in the mid-surface. For what concerns the shear strain state, we use the following strain vector:

$$s = \begin{pmatrix} \mathbf{g}_1 \cdot \mathbf{b} - \mathbf{g}_1^o \cdot \mathbf{b}^o \\ \mathbf{g}_2 \cdot \mathbf{b} - \mathbf{g}_2^o \cdot \mathbf{b}^o \end{pmatrix} = \begin{pmatrix} \beta \\ 0 \end{pmatrix}$$

while the flexural strain state is parameterized by using the following tensor field:

$$d = k - k^o = \begin{pmatrix} \mu^o - \mu\lambda & 0 \\ 0 & \frac{\sin(\theta^o)}{r^o} - \frac{r \sin(\theta)}{r^{o2}} \end{pmatrix}$$

In all the above definitions, the upper index o represents a field when it is evaluated in the reference relaxed configuration. Furthermore, by taking X as a material coordinate, the Euclidean curvilinear abscissa along the corresponding meridian when it is in the reference configuration, we have $h_{11}^o=1$, while it is natural to consider that there is no

transverse shearing in the reference resting configuration, i.e., $\beta^o=0$.

## 2.4 Dynamics

The dynamic model of the shell is given by the balance of kinetic momenta, i.e., by Newton's laws or a variational principle. In any case, this model takes the form of a set of partial differential equations (p.d.e.'s) which govern the time evolution of the system (the shell, in this case) on its configuration space. With the definition (1) of the configuration space of a shell (not necessarily axisymmetric), these p.d.e.'s have been derived in [18] as follows:

$$\begin{gathered} \frac{\partial(j\mathbf{n}^\alpha)}{j\partial X^\alpha} + \bar{\mathbf{n}} = 2\rho h \dot{\mathbf{r}} \\ \frac{\partial(j\mathbf{m}^\alpha)}{j\partial X^\alpha} + \frac{\partial \mathbf{r}}{\partial X^\alpha} \times \mathbf{n}^\alpha + \bar{\mathbf{m}} = 2\rho h \mathbf{b} \times \ddot{\mathbf{b}} \end{gathered} \quad (7)$$

where $j=\sqrt{\det(h)}=r/r^o\sqrt{\lambda^2 + \beta^2}$, while the vectors $\mathbf{n}^\alpha$ and $\mathbf{m}^\alpha$ are, respectively, the resultant of internal stress forces and couples applied by the left part $(x<X^a)$ into the right part $(x \ge X^a)$ of the shell, across the section $X^a$, normalized with the surface Jacobian $j$. $\bar{\mathbf{n}}$ and $\bar{\mathbf{m}}$ are the external force and couple per unit of mid-surface area. For the repeated $\alpha$ the Einstein convention has to be used as in the rest of the paper. As expected, these equations give the time evolution of all the pairs $(\mathbf{r}, \mathbf{b})$ over the shell as a function of the external load and the internal stress.

With respect to the local reference frame, equation (7) can be written, in a geometric notation, as:

$$\mathcal{M}\dot{\eta} = 1/j(j\mathcal{F}_i^1) - ad_{\xi\alpha}^*(\mathcal{F}_i^\alpha) + \mathcal{F}_e + ad_\eta^*(\mathcal{M}\eta) \quad (8)$$

Here $\mathcal{F}_i^\alpha$ is the wrench of internal forces, $\mathcal{F}_e$ is the external wrench of distributed applied forces, $\mathcal{M}$ is the screw inertia matrix and

$$ad_{\xi\alpha}^* = \begin{pmatrix} \tilde{\mathbf{k}}_{\alpha}^T & \tilde{\mathbf{g}}_{\alpha}^T \\ 0 & \tilde{\mathbf{k}}_{\alpha}^T \end{pmatrix} \qquad ad_{\eta}^* = \begin{pmatrix} \tilde{\mathbf{w}}^T & \tilde{\mathbf{v}}^T \\ 0 & \tilde{\mathbf{w}}^T \end{pmatrix}$$

are the co-adjoint maps.

Due to the axisymmetry, the internal and external wrench fields take the particular form [19]:

$$\mathcal{F}_i^1(X,t) = (0,0,M_X,N_X,H,0)^T \in \mathbb{R}^6 \quad (9)$$

$$\mathcal{F}_i^2(X,t) = (M_\phi, 0, 0, 0, 0, -N_\phi)^T \in \mathbb{R}^6 \quad (10)$$
$$ \mathcal{F}_c(X,t) = (0,0,l, f_a, f_b, 0)^T \in \mathbb{R}^6 \quad (11) $$

and the screw inertia matrix is equal to: $\mathbb{R}^6 \otimes \mathbb{R}^6 \ni \mathcal{M} = \text{diag}(\rho J, \rho J_b, \rho J, 2\rho h, 2\rho h, 2\rho h)$.

In the equations above $J$ and $J_b$ are the second moment of the cross sectional line equal to $J=h^2/3$, $J_b \sim 0$.

By outlining the components, we obtain:

$$ \begin{aligned} \rho J \dot{w} &= 1/j(jM_X)' + \lambda H - \beta N_X - \frac{\cos(\theta)}{r^\alpha} M_\phi + 1 \\ 2\rho h \dot{v}_a &= 1/j(jN_X)' - \mu H - \frac{\cos(\theta)}{r^\alpha} N_\phi + f_a + 2\rho h w v_b \quad (12) \\ 2\rho h \dot{v}_b &= 1/j(jH)' + \mu N_X + \frac{\sin(\theta)}{r^\alpha} N_\phi + f_b - 2\rho h w v_a \end{aligned} $$

<figure><img src="image_4.png" /><figcaption>Figure 3. Internal (dark) and external (red) loads exerted in $(X, \phi)$</figcaption></figure>

## 2.5 Constitutive Equations

According to [18], for a hyper-elastic isotropic material the general constitutive equations of a shell are defined in terms of the strain measures as follows:

$$ \begin{aligned} \tilde{n}^{\beta\alpha} &= \frac{2Eh}{1-v^2} H^{\beta\alpha\gamma\delta} e_{\gamma\delta}, \quad \tilde{m}^{\beta\alpha} = \frac{2EhJ}{1-v^2} H^{\beta\alpha\gamma\delta} d_{\gamma\delta} \\ \tilde{q}^{\alpha} &= 2Gh h^{\alpha\beta} s_{\beta}, \end{aligned} \quad (13) $$

where $E$ is the Young modulus, $G$ the shear modulus and $\nu$ is the Poisson modulus; $h^{\alpha\beta}$ are the elements of $h^{\alpha-1} = I^{-1} = I$ and $H^{\beta\alpha\gamma\delta}$ define a four order (Hook-like) tensor given by:

$$ H^{\beta\alpha\gamma\delta} = [\nu(h^{\alpha\beta\alpha}h^{\gamma\delta}) + 1/2(1-\nu)(h^{\alpha\beta\gamma}h^{\gamma\alpha\delta} + h^{\alpha\beta\delta}h^{\gamma\alpha\gamma})]. $$

and, in our case, resulting in:

$$ H^{\beta\alpha\gamma\delta} = \begin{cases} 1 & \beta\alpha\gamma\delta = 1111, 2222 \\ \nu & \beta\alpha\gamma\delta = 1122, 2211 \\ 0 & \text{otherwise} \end{cases} $$

We extend this formulation to a viscoelastic constitutive model based on the Kelvin-Voigt model, which simply adds, to the elastic term, a viscous contribution linearly proportional to the rate of strain:

$$ \begin{aligned} \tilde{n}^{\beta\alpha} &= \frac{2Eh}{1-v^2} H^{\beta\alpha\gamma\delta} e_{\gamma\delta} + \frac{6\nu h}{1-v^2} H^{\beta\alpha\gamma\delta} \dot{e}_{\gamma\delta} \\ \tilde{m}^{\beta\alpha} &= \frac{2EhJ}{1-v^2} H^{\beta\alpha\gamma\delta} d_{\gamma\delta} + \frac{6\nu hJ}{1-v^2} H^{\beta\alpha\gamma\delta} \dot{d}_{\gamma\delta} \\ \tilde{q}^{\alpha} &= 2Gh h^{\alpha\beta} s_{\beta} + 2\nu h h^{\alpha\beta} \dot{s}_{\beta} \end{aligned} \quad (14) $$

where the Young modulus and the shear modulus have been replaced by the shear viscosity constants $3\nu$ and $\nu$, respectively, following the same procedure developed in [27] for a beam-like structure (novel analysis for the shell case is currently in progress).

The $\tilde{n}^{\beta\alpha}$, $\tilde{m}^{\beta\alpha}$ and $\tilde{q}^{\alpha}$ respectively, stand for the components of the effective resultant traction, couple and shearing stresses<fnref n="1" />, which are related to the stresses of (7) through the relations:

$$ \mathbf{n}^{\alpha} = (\tilde{n}^{\beta\alpha} + \lambda_{\mu}^{\beta} \tilde{m}^{\alpha\mu}) \frac{\partial \mathbf{r}}{\partial X^{\beta}} + (\tilde{q}^{\alpha} + \lambda_{\mu}^{3} \tilde{m}^{\alpha\mu}) \mathbf{b} \quad (15) $$

$$ \mathbf{m}^{\alpha} = \mathbf{b} \times \tilde{\mathbf{m}}^{\alpha} = \mathbf{b} \times \left( \tilde{\mathbf{m}}^{\beta\alpha} \frac{\partial \mathbf{r}}{\partial X^{\beta}} + \tilde{\mathbf{m}}^{3\alpha} \mathbf{b} \right) \quad (16) $$

where the functions $\lambda_{\mu}^{\beta}$ and $\lambda_{\mu}^{3}$ are defined by the following equation:

$$ \mathbf{b}_{,X^{\mu}} = \lambda_{\mu}^{\beta} \frac{\partial \mathbf{r}}{\partial X^{\beta}} + \lambda_{\mu}^{3} \mathbf{b} \quad (17) $$

By comparing (9) and (10) with (15) and (16), after some algebra (Appendix) the following relations between the effective resultant stress-couple and our axialsymmetric resultant stress-couple are derived:

<footnote n="1">1 The word 'effective' here is meant to underline that these functions can be directly related to the constitutive equations of the three dimensional theory [18].</footnote>
$$M_X = -\lambda \tilde{m}^{11}, M_\phi = -\frac{r}{r^0} \tilde{m}^{22},$$

$$N_X = \lambda \tilde{n}^{11} - \mu \tilde{m}^{11}, N_\phi = \frac{r}{r^0} \tilde{n}^{22} - \frac{\sin(\theta)}{r^0} \tilde{m}^{22}, \quad (18)$$

$$H = \tilde{q} + \beta \tilde{n}^{11}.$$

With these relations and the constitutive equations (14), the constitutive equations for our internal stresses (i.e., $N_X(X)$, $N_\phi(X)$, $H(X)$, $M_X(X)$ and $M_\phi(X)$) are given by:

$$N_X = \frac{2Eh}{1-v^2} \left[ \lambda(e_{11} + v e_{22}) - J\mu(d_{11} + v d_{22}) \right] + \frac{6\nu h}{1-v^2} \left[ \lambda(\dot{e}_{11} + v \dot{e}_{22}) - J\mu(\dot{d}_{11} + v \dot{d}_{22}) \right]$$

$$N_\phi = \frac{2Eh}{1-v^2} \left[ \frac{r}{r^0}(e_{22} + v e_{11}) - J \frac{\sin(\theta)}{r^0}(d_{22} + v d_{11}) \right] + \frac{6\nu h}{1-v^2} \left[ \frac{r}{r^0}(\dot{e}_{22} + v \dot{e}_{11}) - J \frac{\sin(\theta)}{r^0}(\dot{d}_{22} + v \dot{d}_{11}) \right] \quad (19)$$

$$H = 2h\beta \left[ G + \frac{E}{1-v^2}(e_{11} + v e_{22}) \right] + 2h\dot{\beta} \left[ \nu + \frac{6\nu}{1-v^2}(\dot{e}_{11} + v \dot{e}_{22}) \right]$$

$$M_X = -\frac{2EhJ}{1-v^2} \lambda(d_{11} + v d_{22}) - \frac{6\nu hJ}{1-v^2} \lambda(\dot{d}_{11} + v \dot{d}_{22})$$

$$M_\phi = -\frac{2EhJ}{1-v^2} \frac{r}{r^0}(d_{22} + v d_{11}) - \frac{6\nu hJ}{1-v^2} \frac{r}{r^0}(\dot{d}_{22} + v \dot{d}_{11})$$

In order to avoid the polar singularity ($r^o=0$) in the constitutive law, the mantle is assumed to be connected to a small rigid spherical cap axisymmetric with respect to the mantle axis which crosses the cap in the cap's pole (see Figure 2). The inertia of the cap is assumed to be negligible.

## 2.6 External Loads

The external loads taken into account are the drag and the added mass exerted by the fluid. Below, the expression of all the external loads are shown.

$$\mathcal{F}_e(X,t) = (0,0,0,d_a,0,0)^T + (0,0,0,\mathbf{c}^T)^T. \quad (20)$$

### 2.6.1 Drag Load

The drag load is proportional to the square of the velocity and directed in the opposite direction. In particular, the drag force was taken to act tangentially to the shell, as expressed by eq. (20). The magnitude of the drag load is also determined by the geometry of the director X and by hydrodynamics phenomena expressed by empirical coefficients.

The equation (21) shows the resultant expression used in this model.

$$d_a(X,t) = -\rho_w C_a (v_a + \sin(\theta)v) |v_a + \sin(\theta)v| \quad (21)$$

where $\rho_w$ is the water density, $C_a$ is an empirical hydrodynamic coefficients which incorporates the geometric and hydrodynamics factors in the viscosity model, and $v$ is the swimming velocity of the mantle along $\mathbf{e}_3$ (see section 3).

### 2.6.2 Added Mass Load

The added mass, representing the load locally exerted by pressure, acts in the direction normal to the surface of the shell and is proportional to the acceleration. As in the case of the drag load, the magnitude is also determined by the geometry of the director X and by hydrodynamics phenomena expressed, in part, by correction coefficients. The problem accounted for here, ideally requires a distinct formulation for the added mass as far as the external and internal domains are concerned. While an expression for the added mass terms is straightforward for the external flow case, treatment of the condition inside the shell is significantly more complicated because of the closed geometry and the time-varying shape of the domain accounted for. As a first approximation this problem is disregarded and the following simplified formulation is adopted.

The equation (22) shows the resultant added mass vector used in this model.

$$\mathbf{c}(X,t) = -\frac{\partial(\rho\mathbf{F}(\mathbf{v} + R^T\dot{\mathbf{u}}))}{\partial t} = -\rho\mathbf{F}(\dot{\mathbf{v}} + R^T\dot{\mathbf{u}}) + \mathbf{w} \times \rho\mathbf{F}(\mathbf{v} + R^T\dot{\mathbf{u}}) \quad (22)$$

where $F(X) \in \mathbb{R}^3 \otimes \mathbb{R}^3$ is a tensor that incorporates the geometric and hydrodynamics factors. It is equal to $F = \text{diag}(0,2h B_b, 0)$, where $B_b$ is the hydrodynamic correction coefficient for the added mass model. The vector $\dot{\mathbf{u}}$ is the swimming velocity of the mantle with respect to the fixed Galilean frame $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$ equal to $[0,0,v]^T$.

## 3. Jet Propulsion Model

The propulsion modelling is based on the standard one-dimensional momentum equation for a neutrally buoyant, rigid body translating in water,

$$(\rho_w V + \rho V_m) \dot{v} = -\frac{1}{2} \rho_w C_d A_{ref} v |v| - B \dot{v} + q \rho_w \dot{V} \quad (23)$$

where $V(t)$ is the mantle inner volume and $V_m$ is the volume of the elastic material composing the external shell of the mantle. $v=\dot{u}$ is the swimming velocity (where $u$ is the mantle position). The first right-hand side (RHS) term represents the drag, $C_d$ being the drag coefficient and $A_{ref}(t)$ a reference area of the mantle; the second RHS term is the added mass, with $B$ being the axial added mass coefficient of an ellipsoid of revolution ([28]) expressed by:
$$B = \frac{\alpha_0}{2 - \alpha_0} m$$

with

$$\alpha_0 = \frac{2(1 - e^2)}{e^3} \left( \frac{1}{2} \ln \frac{1 + e}{1 - e} - e \right) \qquad e = 1 - (b/a)^2$$

and

$$m = \rho_w V + \rho V_m$$

where a and b respectively represent the major and minor semi-axis of the shell.

The third term in eq. (23) is the thrust, given by the speed of the outflow q, with respect to the mantle, across the nozzle-exit area $A_n(t)$ and the variation of mass $\rho_w \dot{V}$ occurring within the robot due to the collapse of the elastic chamber.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Depiction of the elements accounted for in the dynamics model of the pulsed-jet propulsion</figcaption></figure>

The outflow speed is given by:

$$q = - \frac{C_f \dot{V}}{A_n}$$

where $C_f$ defines a flow loss coefficient at the nozzle entrance, which is taken to vary between 0.6 and 1 [29].

The thrust term is thus rewritten,

$$q \rho_w \dot{V} = - \frac{\rho_w C_f \dot{V} |\dot{V}|}{A_n} \qquad (24)$$

In order to study the thrust generated by the deformation of the axisymmetric shell, the propulsion model presented above was adapted to the geometry of the elastic mantle-like shell. Thus, in terms of the configuration variables $r(X,t)$ and $z(X,t)$, the nozzle-exit area, the reference area and the volume of the mantle are:

$$A_n(t) = \pi r(L,t)^2 \qquad A_{ref}(t) = \pi \max(r(\cdot,t))^2$$

$$V_m = - \int_{0}^{z^{(L)}} 4\pi r^0 h dz = -4\pi h \int_{0}^{L} r^0 z^{\alpha'} dX$$

while the mantle inner volume and his time derivative are:

$$V(t) = - \int_{0}^{z^{(L)}} \pi r^2 dz = -\pi \int_{0}^{L} r^2 z' dX$$

$$\dot{V}(t) = -\pi \int_{0}^{L} (2r\dot{z}' + r^2 \dot{z}') dX$$

where,

$$z' = \sin(\theta)\lambda + \cos(\theta)\beta \qquad \dot{r} = \cos(\theta)v_a - \sin(\theta)v_b$$

$$\dot{z}' = \sin(\theta)(\dot{\lambda} - w\beta) + \cos(\theta)(\dot{\beta} + w\lambda)$$

Given the orientation chosen in the kinematics (fig. 2) the orientation of the mantle surface comes out to be negative, which leads to the minus sign in the integral above.

## 4. Pulsed-Jet Dynamic Model

The previous formulation can be used to address the study of cephalopod locomotion as well as the dynamics of a cephalopod-like underwater robot such as those presented in [12, 13] and [14]. We consider that the mantle is jet propelled along the e_3-axis while the axisymmetric external loads cannot generate a net displacement in another dimensions of the Lie group SE(3). To that end, the kinematics developed in section 2.1 have to be slightly modified by replacing the transformation $g_1(\phi,t)$ with the following

$$\bar{g}_1 \in SE(3):$$

$$(\phi, X) \mapsto \bar{g}_1(\phi, X) = \begin{pmatrix} \exp(\bar{e}_3\phi) & (\mathbf{r} + \mathbf{u}) \\ 0 & 1 \end{pmatrix}$$

where $\mathbf{u}(t)$ is the position vector of the mantle equal to $(0,0,u(t))^T$, which depends only on $t$. Consequently, we have $\bar{g} = \bar{g}_1 \bar{g}_2 \bar{g}_3$ and the corresponding velocity twist vector field $\bar{\eta}(X,t)$ is

$$\bar{g}^{-1} \bar{g} = \bar{\eta} = \begin{pmatrix} \tilde{\mathbf{w}} & (\mathbf{v} + R^T \dot{\mathbf{u}}) \\ 0 & 0 \end{pmatrix} \in se(3)$$

that define the new kinematic equations (first three lines of (25)). Due to the different nature of the models to be coupled (the shell model is formulated relative to the body
while the jet propulsion model is earth-fixed) this velocity field just introduced has been used only for the kinematics equation. This implies that $\bar{\eta}(X,t)$ do not enter either the compatibility equation (5) nor the dynamic equation (8), which refers to the internal dynamic while the global mantle dynamic is governed by the equations (23).

The coupling between the two is manifested in the equations of the model presented above. To summarize, at every time step the shell model uses the velocity variable ($v$) of the jet propulsion model to calculate the drag and added mass loads on the mantle, while the jet propulsion model uses the surface configuration variable of the shell ($r$ and $z$) to calculate the geometry of the mantle and consequently the drag, added mass and thrust loads. This concept is shown in figure 5.

<figure><img src="image_6.png" /><figcaption>Figure 5. Scheme of the coupling between the shell model and the jet propulsion model for every time step</figcaption></figure>

Now we can state the system of second-order partial differential equations by gathering the (modified) kinematics equations, the compatibility equations (6), the mantle dynamic equations (12) and the jet propulsion models (23), in the state-space form $\dot{x} = f(x, x', x'', t)$ as follows:

$$
\begin{align}
\dot{\theta} &= w \nonumber \\
\dot{r} &= \cos(\theta)v_a - \sin(\theta)v_b \nonumber \\
\dot{z} &= \sin(\theta)v_a + \cos(\theta)v_b + v \nonumber \\
\dot{\mu} &= w' \nonumber \\
\dot{\lambda} &= v_a' + \beta w - \mu v_b \nonumber \\
\dot{\beta} &= v_b' - \lambda w + \mu v_a \nonumber \\
\dot{w} &= \frac{(jM_X')' / j + \lambda H - \beta N_X - \cos(\theta) M_\phi / r^o + l}{\rho J} \tag{25} \\
\dot{v}_a &= \frac{(jN_X')' / j - \mu H - \cos(\theta) N_\phi / r^o + f_a}{2\rho h} + w v_b \nonumber \\
\dot{v}_b &= \frac{(jH')' / j + \mu N_X + \sin(\theta) N_\phi / r^o + f_b - 2\rho h w v_a}{2h(\rho + \rho B_b)} \nonumber \\
\dot{u} &= v \nonumber \\
\dot{v} &= \frac{q \rho_w \dot{V} - \frac{1}{2} \rho_w C_d A_{ref} v |v|}{\rho_w V + \rho V_m + B} \nonumber
\end{align}
$$

where the second-order partial differentiation is due to the viscous component of the constitutive equation.

In these forward dynamics, the state vector $x$ is infinite dimensional since all its components (along with those of $f$) are some functions of the profile abscissa $X$. As a result, the above state equation has to be first space-discretised on a grid of nodes along $[0, L]$ before being time integrated using explicit or implicit time integrators starting from the initial state $x(0)$. In this grid, all the space derivatives appearing in the $f$ vector can be approximated by finite difference schemes, while the external and internal stress forces and couples are respectively given by (21), (22) and the constitutive law (19). The system (25) has been solved in Matlab® through an *ad hoc* explicit second-order finite-difference scheme. On an AMD Phenom(TM) II X4 965 processor, at 784 Mhz, 3.25 GB of RAM it took almost 27sec for one second of simulation and a space resolution of 1.25% of the total length.

## 4.1 Boundary Conditions

At the boundary overlooking the rigid cap, we have the following static condition:

$$ \eta(0) = \eta_{-} = 0 \qquad (26) $$

where $\eta_{-}$ is fixed by the boundaries of the cap in which the mantle is clamped. At the other end of the mantle, we have the following natural boundary conditions:

$$ N_X(L) = f_{a,+}, \quad H(L) = f_{b,+}, \quad M_X(L) = -l_{+}, $$

where $f_{a,+}$, $f_{b,+}$ and $l_{+}$ denote the eventual external forces and torque along $a$, $b$ and $-e_\phi$, respectively, applied onto the sharp boundaries of the mantle orifice.

## 5. Simulation Results

The model developed in the previous sections has been used to analyse the properties of a jet-propelled soft robot in order to provide a support in the process of mechanical design and control-algorithm formulation of a cephalopod-inspired robot. Two actuation modes are implemented: a more artificial (robot-like) external actuation and a more biological (muscle-like) internal actuation.

### 5.1 External Actuation

The rhythmic external actuation can be modelled by taking the radial (along $e_r$) force $f_r(X,t) = f_a \cos(\theta) - f_b \sin(\theta)$ as a $T$-periodic function with two phases. In the first phase $[0, T_c]$, i.e., the contraction phase, $f_r(X,t) = P(X)$. In the second phase, i.e., the relaxation phase $(T_c, T]$, $f_r(\cdot, t) = 0$, the mantle passively recovers its resting shape thanks to the internal restoring stresses. From now on, a linear spatial
<figure><img src="image_7.png" /><figcaption>**Figure 6.** A sketch plot of the two linear distribution of the external actuation used in the simulations</figcaption></figure>

distribution of the actuation is assumed with the form: $P(X)=(P-P_o/L)*X + P_o$ (figure 6).

For the external actuation mode, an extensive simulation analysis has been conducted by varying the shape of the mantle, the material stiffness and the spatiotemporal properties of the external actuation. The two geometries accounted for are the conical and the ellipsoidal shape. The other families of parameters with their respective values used in the simulations are listed in Table 1. Every test is a mathematical combination of all these possibilities. The test number corresponding to a specific value is reported in Table 1. In Table 2 and 3 the fixed geometrical and load parameters are reported.

<table><caption>**Table 1.** Families of Variable Parameters for the External Actuation</caption><thead><tr><th>Parameter</th><th>Value 1</th><th>test #</th><th>Value 2</th><th>test #</th><th>Value 3</th><th>test #</th></tr></thead><tbody><tr><td>E and P</td><td>110 Pa<br>0.3 N</td><td>1-12</td><td>1100 Pa<br>3 N</td><td>13-24</td><td>11000 Pa<br>30 N</td><td>25-36</td></tr><tr><td>T and T_c</td><td>1 s<br>0.5 s</td><td>1-4<br>25-28</td><td>1.5 s<br>0.5 s</td><td>5-8<br>29-32</td><td>1.5 s<br>1 s</td><td>9-12<br>21-24<br>33-36</td></tr><tr><td>P_o</td><td>P</td><td>odd</td><td>P/2</td><td>even</td><td></td><td></td></tr><tr><td>Ellipsoid</td><td>1-2, 5-6, 9-10, 13-14</td><td>test #</td><td>test #</td><td></td><td></td><td></td></tr><tr><td></td><td>17-18,<br>21-22</td><td>25-26,<br>29-30</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>33-34</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Cone</td><td>3-4, 7-8,<br>19-20,<br>23-23</td><td>11-12,<br>15-16<br>27-28,<br>31-32</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>35-36</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 2.** Cone and Ellipsoid Parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>L</td><td>80 cm</td></tr><tr><td>h</td><td>2 mm</td></tr><tr><td>R_c</td><td>5 cm</td></tr><tr><td>R_w</td><td>15 cm</td></tr><tr><td>a</td><td>50 cm</td></tr><tr><td>b</td><td>20 cm</td></tr></tbody></table>

<table><caption>**Table 3.** Load Parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>μ</td><td>50 Pa · s</td></tr><tr><td>ν</td><td>0</td></tr><tr><td>ρ</td><td>1.08 kg/dm³</td></tr><tr><td>ρ_w</td><td>1.02 kg/dm³</td></tr><tr><td>C_a</td><td>0.01</td></tr><tr><td>B_b</td><td>1</td></tr><tr><td>C_d</td><td>1</td></tr><tr><td>C_f</td><td>1</td></tr></tbody></table>

Ultimately, 36 simulations were performed. The comparison among the various tests is established solely in terms of maximum mantle velocity (v) attained after a specified interval. In this respect, energetic analysis based on estimation of the local deformation (ε₁(X)) is left for future development. In figure 7 a bar plot of the maximum velocities is shown. The three highest peaks correspond to the ellipsoidal mantles with the highest stiffness. The results suggest that, of the geometries taken into consideration, an ellipsoid of revolution gives a highest rate of volume contraction (that corresponds to a higher thrust) augmented by the high reactivity of the stiffness, counter-balanced by an higher actuation load (see Table 1). In figure 8 few snapshots of the best performing test, i.e., an ellipsoidal mantle with E = 110e2Pa under the T_c = 0.5, T = 1.5s, are presented. In figure 9 his velocity profile is shown, together with the best performer of the two other stiffnesses discussed below. The maximum velocity reached by this specimen is 2.2m/s.

In figure 10 the results for the cases with the three different stiffness value are portrayed separately. This underlines that, in the case of the lowest stiffness, the best time routine for the external actuation is T = 1.5, T_c = 1 while for the medium stiffness value it is T = 1.5, T_c = 0.5 and for the highest one there are a sort of *ex aequo* between T = 1, T_c = 0.5 and T = 1.5, T_c = 0.5.

In figure 11 the results of the two geometries (conical and ellipsoidal) are compared. Since two consecutive bars differ
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Bar plot of all the 36 simulations for the external actuation. The three highest peaks correspond to the ellipsoidal mantles with the highest stiffness.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Snapshots of the ellipsoidal mantle configuration with $E = 110e^2$ Pa under the $T_c = 0.5$, $T = 1.5$ s actuation routine, the best performer for the external actuation. From top left: $t=0$ (a), $t=1$ (b), $t=3$ (c), $t=5$ (d), $t=7$ (e) and $t=9$ (f).</figcaption></figure>

only for the spatial distribution of the actuation (figure 6), in figures 11 it is demonstrated the extent by which this parameter influences the performances of the ellipsoidal geometry and how little, instead, it affects the performances of the conical geometry. This is due to the fact that, in the conical case, the linear distribution of the activation

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Velocity profiles of the best performers for the three stiffnesses under the external actuation</figcaption></figure>

results in a shrinkage of the nozzle area which, in turn, significantly enhances the thrust, since the low exit speed of the outflow is the most significant limitation of the conical geometry.

## 5.2 Internal Actuation

Let us reiterate that the squid is modelled by an axisymmetric shell internally actuated by a network of circular muscles organized in rings around the shell axis [5]. In order to address the study of cephalopod locomotion with an internal actuation we have replaced the reference strains of the shell ($h^o, k^o$ and $\beta^o$, see 2.3) with a desired time-varying strain field ($h^d, k^d$ and $\beta^d$) obtained by the following constrain inputs: $\beta^d(X,t)=0$, $\lambda^d(X,t)=1$ and $r^d(X,t)=f(X,t)$ where $r^d(X,t)$ is the desired value of the radius of the axisymmetric strip of squid which passes through the point of abscissa X along the shell profile. Since the muscles can only contract, the rhythmic muscular activity can be modelled by taking $f(X,t)$ as a T-periodic function with two phases. In the first phase, i.e., the contraction phase, it decreases over a short part of $[0,T]$; in the second phase, i.e., the relaxation phase, $f(X,t)=r^o(X)$, and the mantle recover its resting shape passively thanks to the internal restoring stresses. In these conditions, the constitutive law (19) is partly used as a kind of linear control law where the elastic coefficients stand for proportional control gains. This approach ideally enables us to expand the analysis beyond isotropic materials and consider more general constitutive laws of the type (19) where the constant passive elastic coefficients are replaced by active time varying ones which model the nonlinear muscular activity. In spite of all these possibilities, we restrict our investigations to the case of (19) with constant coefficients.

As was performed earlier for the external actuation, the conical and the ellipsoidal geometries are accounted for, while the other families of parameters with their different values used in the simulations are listed in Table 4. The two desired radius functions for the conical and ellipsoidal mantle are illustrated in figure 12.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Bar plot of the simulations with $E = 110$ Pa (top), $E = 1100$ Pa (middle) and $E = 11000$ Pa (down) for the external actuation. Under the red square the best time routine for each case, for the lowest stiffness we have $T = 1.5$, $T_c = 1$ while for the middle one $T = 1.5$, $T_c = 0.5$ and for the highest one there are a sort of *ex aequo* between $T = 1$, $T_c = 0.5$ and $T = 1.5$, $T_c = 0.5$.</figcaption></figure>

In this case, only 18 simulations were performed, since the spatial distribution of the actuation is no longer taken into consideration. In figure 13 a bar plot of the maximum velocities reached is shown, proving that the best speed

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Bar plot of the simulations for an ellipsoidal (top) and conical (down) geometry for the external actuation. The influence of the two spatial distributions of the actuation (figure 6) is much more evident for the ellipsoidal case than for the conical one.</figcaption></figure>

<table><caption>**Table 4.** Families of Variable Parameters for the Internal Actuation</caption><thead><tr><th>Parameter</th><th>Value 1</th><th>test #</th><th>Value 2</th><th>test #</th><th>Value 3</th><th>test #</th></tr></thead><tbody><tr><td>E</td><td>110 Pa</td><td>1-6</td><td>1100 Pa</td><td>7-12</td><td>11000 Pa</td><td>13-18</td></tr><tr><td>T and</td><td>1 s</td><td>1-2, 7-8</td><td>1.5 s</td><td>3-4, 9-10</td><td>1.5 s</td><td>5-6, 11-12</td></tr><tr><td>T<sub>c</sub></td><td>0.5 s</td><td>13-14</td><td>0.5 s</td><td>15-16</td><td>1 s</td><td>17-18</td></tr><tr><td>Ellipsoid</td><td>odd</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Cone</td><td>even</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

performances are those observed with the conical mantle with high stiffness. This is probably due to the fact that the desired radius function is linear (constant) for most of the mantle profile and thus more attuned with the linear conical profile. In figure 14 few snapshots of the best performer, i.e., a conical mantle with $E = 110e2$ Pa under the $T_c = 0.5$, $T = 1$ s, are presented. In figure 15 his velocity profile is shown, together with the best performer of the two other stiffness discussed below. The velocity profile of the highest stiffness decade toward zero because, after six
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Desired radius functions for the conical (up) and ellipsoidal (down) mantle</figcaption></figure>

iterations of the actuation routine, the mantle kinetic energy is left to dissipate in water. The maximum velocity reached by this specimen is 4.4m/s.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Bar plot of all the 18 simulations for the internal actuation. The highest peaks correspond to the conical mantles with the highest stiffness.</figcaption></figure>

In figure 16 the results for the three different stiffness are depicted separately. This brings evidence that there is not an actual best performing routine, except in the case of low stiffness, where the $T=1.5$, $T_c=1$ stands out clearly. Fur-

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Snapshots of the conical mantle configuration with $E=110e2 Pa$ under the $T_c=0.5$, $T=1$ s actuation routine, the best performer for the internal actuation. From top left: $t=0.5$ (a), $t=1$ (b), $t=1.75$ (c) and $t=2$ (d).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Velocity profiles of the best performers for the three stiffnesses under the internal actuation</figcaption></figure>

thermore, comparing figure 16 with figure 10, we observe that the internal routine is capable of achieving better performances, with the only exception of the case with $E=1100Pa$ where the two actuations have similar results.

In figure 17 the results of the two geometries (conical and ellipsoidal) are isolated. The conical mantle performances
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Bar plot of the simulations with *E* = 110 Pa (top), *E* = 1100 Pa (middle) and *E* = 11000 Pa (down) for the internal actuation. There is not a real favourite routine except in the case of low stiffness, where the *T* = 1.5, *T*<sub>c</sub> = 1 is clearly selected.</figcaption></figure>

are almost independent of the actuation routine, while the
ellipsoidal mantle, in test number 9, benefits from a
particular resonance, hence doubling the efficacy for the
case of T = 1.5, T_c = 0.5s w.r.t. the other routines.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Bar plot of the simulations for an ellipsoidal (top) and conical (down) geometry for the internal actuation</figcaption></figure>

**6. Discussion**

While a rigorous validation of the model goes beyond the
scope of the present work and has been thoroughly
undertaken by the authors in [30], a broad comparison
between the current results and existing experiments
performed with swimming cephalopods as well as with the
cephalopod-inspired vehicle of figure 1 can be established.
The results from the model of figure 7, which more closely
resembles an actual cephalopod, demonstrate that, given a
Young modulus of *E* =11000 and a pulsation routine with
*T*<sub>c</sub>=0.5s and *T* =1.5s, a swimming speed in the range of 2.75
and 0.625 bdl/s (body lengths per second) is achieved. This
estimate is consistent with measurements on cephalopods
performed by [1] and [31]. In the case of *Loligo vulgaris*, the
common squid, a large variability in swimming perform-
ances is found: specimen are known to range in size
between 20 and 70 cm (mantle length) and to perform
sustained swimming at about 1.4 m/s [1] and escape
manoeuvre speeds as high as 2.0 m/s [31] or even 4.2 m/s
[1]. This suggests that fast jetting squid can easily be
observed to swim at speed beyond 5 bdl/s; however, for
steadily swimming specimens, a value of about 1.0 bdl/s
was found [32] coincident with the peak of swimming efficiency. As for the *Octopus vulgaris*, i.e., the common octopus, specimen are found to grow to a mantle length of about 25 cm and swim at an average speed of 0.18 m/s, which is 0.7 bdl/s.

Finally, a last comparison shall be established with the cephalopod-inspired vehicle of [14], where an artificial mantle of rubber-like materials was complemented with the actuators capable of replicating a pulsed-jet routine analogous to that of actual cephalopods; see figure 1. This vehicle weighs 335 g, has a mantle maximum extension in the axial direction of 16 cm and a mantle capacity of 35 ml, making it fairly comparable to an average octopus. During experimental testing, the vehicle was found to swim at an average speed of 0.135 m/s (0.84 bdl/s) with a burst speed at 0.2 m/s (1.25 bdl/s) at a jet frequency of 1.6 Hz.

The major limitation and potential source of error of the model presented lies in the simplification associated with the treatment of the fluid-structure interaction which, in the present case, was dealt with by resorting to the coupling scheme illustrated in figure 5. While the solution adopted here allows for fast predictions of the swimming performance of cephalopods or cephalopod-like vehicles, a more rigorous approach would require to take into account the continuous, nonlinear pressure distribution arising within and outside the shell during the pulsation routine. The external part of the fluid problem has already been discussed by Anderson and DeMont [20] with recourse to a slender body potential flow model, but acknowledgement of the internal and external flow along with the elastic body structural deformations requires significant revision of the existing models and represents ongoing work. A final remark concerns the assumption of axisymmetry which constitutes an additional limitation over the range of possible geometries that could be taken into consideration in the study of the cephalopod-inspired vehicles.

## 7. Conclusion

In this paper, a general model of an elastic chamber resembling the mantle of a cephalopod is formulated. The model is based on a geometrically exact formulation of a generic axisymmetric elastic shell coupled with a rocket-like jet propulsion thrust model. The coupling of these models enables the capture of the deformation that the mantle undergoes during sequences of collapse and inflation with the purpose of replicating the pulsating propulsion routine of living cephalopods or cephalopod-inspired underwater robots. A broad range of scenarios can be simulated where various shapes, degrees of stiffness and activation routines are combined and tested. This in turn provides essential information with relevance both to the study of the biomechanics of swimming cephalopods and on the design and control of cephalopod-inspired aquatic machines.

An example of how the model developed herein can be used both for biological and robotics studies is reported in

section 5.1 and 5.2 where the results from a sequence of 54 simulations are analysed in terms of the maximum speed achieved by the mantle in order to establish the best performing combination of the parameters examined. These results are divided according to their more robotic-like or more muscle-like actuation suggesting that, with the former, an ellipsoid of revolution composed of an elastic material of Young modulus $E=110e2Pa$ and a spatially linear activation with a routine $T_c=0.5, T=1.5s$ performs best. In the latter case, on the contrary, a conical mantle with $E=110e2Pa$ and $T_c=0.5, T=1s$ swims faster.

These results demonstrate the aptness of the tool developed in studying both the structural and fluid dynamics aspect of pulsed-jet propulsion and prove that, once properly tuned, the model can be exploited to investigate the fine-scale dynamics of this as-yet poorly explored mode of aquatic locomotion.

## 8. Acknowledgements

This article is a revised and expanded version of a paper entitled 'Structural Dynamics and Propulsion Modelling of a Pulsed-Jet Underwater Soft Robot' presented at Bio-inspired Robotics, Frascati, 14-15 May 2014.

This work was jointly supported by the EU Commission in the frame of the CFD-OctoProp Project FP7 European Reintegration Grant and the Foundation Grant project PoseiDRONE of the Cassa di Risparmio di Livorno.

## 9. References

[1] Trueman E. R., Packard A. (1968). Motor performances of some cephalopods. J. of Exp. Bio. 49: 495-507.

[2] Bartol I. K., Patterson M. R., Mann R. (2001) Swimming mechanics and behavior of the shallow-water brief squid *Loliguncula brevis*. J. of Exp. Bio. 204: 3655-3682.

[3] Sfakiotakis M., Kazakidi A., Pateromichelakis N., Tsakiris D. P. (2013) Octopus-inspired eight-arm robotic swimming by sculling movements. Robotics and Automation (ICRA), 2013 IEEE International Conference on. 6-10 May 2013. pp. 5155,5161.

[4] Gosline J. M., DeMont M. E. (1985) Jet-propelled swimming in squids. Sci. Am. 252: 96-103.

[5] Krieg M., Mohseni K. (2012) New perspectives on collagen fibers in the squid mantle. J. Morphol. 273: 586-95.

[6] Anderson E. J., Edwin DeMont M. (2000) The mechanics of locomotion in squid *Loligo pealei*: locomotory function and unsteady hydrodynamics of the jet and intramantle pressure. J. of Exp. Bio. 203: 2851-2863.
[7] Krueger P. S., Gharib M. (2005) Thrust augmentation and vortex ring evolution in a fully pulsed jet AIAA Journal. 43: 792-801.

[8] Krieg M., Mohseni K. (2008) Thrust characterization of a bio-inspired vortex ring generator for locomotion of underwater robots. IEEE J. of Ocean Eng. 33: 123-132.

[9] Ruiz L. A., Whittlesey R. W., Dabiri J. O. (2011) Vortex-enhanced propulsion J. of Fluid Mech. 668: 5-32.

[10] Moslemi A. A., Krueger P. S. (2010) Propulsive efficiency of a biomorphic pulsed-jet underwater vehicle. Bioin. and Biomi. 5: 1-14.

[11] Weymouth G. D., Triantafyllou M. S. (2013) Ultrafast escape of a deformable jet-propelled body. J. of Fluid Mech. 721: 367-385.

[12] Giorgio Serchi F., Arienti A., Laschi C. (2013) Biomimetic vortex propulsion: toward the new paradigm of soft unmanned underwater vehicles. Mechatronics, IEEE/ASME Transactions on. 18: 484-493.

[13] Giorgio Serchi F., Arienti A., Laschi C. (2013) An elastic pulsed-jet thruster for soft unmanned underwater vehicles. IEEE International Conference on Robotics and Automation (ICRA). 2013 May 6-10; Karlsruhe.

[14] Giorgio Serchi F., Arienti A., Laschi C. (2013) A soft unmanned underwater vehicle with augmented thrust capability. MTS/IEEE Proceedings Oceans. 2013 Sep 23-27; San Diego.

[15] Gosline J. M., Steeves J. D., Hrman A. D., Demont M. E. (1983) Patterns of circular and radial mantle muscle activity in respiration and jetting of the squid *Loligo opalescens*. J. of Exp. Bio. 104: 97-109.

[16] Krieg M., Mohseni K. (2013) Modelling circulation, impulse, and kinetic energy of starting jets with non-zero radial velocity. J. of Fluid Mech. 719: 488-526.

[17] Giorelli M., Giorgio Serchi F., Laschi C. (2013) Forward speed control of a pulsed-jet soft-bodied underwater vehicle. Proceedings of the MTS/IEEE OCEANS Conference. 2013 Sep 21-27; San Diego.

[18] Simo J. C., Fox D. D. (1989) On stress resultant geometrically exact shell model. Part I: formulation and optimal parametrization. Journal Computer Methods in Applied Mechanics and Engineering. 72: 267-304.

[19] S. S. Antman (2005) Nonlinear problems of elasticity, 2nd edn. New York:Springer.

[20] Anderson E. J., Quinn W., DeMont M. (2001) Hydrodynamics of locomotion in the squid *Loligo pealei*. J. of Fluid Mech. 436: 249-266.

[21] Cosserat E., Cosserat F. (1909) Théorie des corps déformables. Paris:Hermann.

[22] Green A. E., Naghdi P. M., Wainwright W. L. (1965) A general theory of a Cosserat surface. Archive for Rational Mechanics and Analysis. 20: 287-308.

[23] Reissner E. (1945) The effect of transverse shear deformation on the bending of elastic plates. ASME Journal of Applied Mechanics. 12: 68-77.

[24] Boyer F., Porez M., Leroyer A. (2010) Poincaré-Cosserat equations for the Lighthill three-dimensional large amplitude elongated body theory: application to robotics. Journal of Nonlinear Science. 20: 47-79.

[25] Boyer F., Porez M., Leroyer A., Visonneau M. (2008) Fast dynamics of an eel-like robot-comparisons With Navier-Stokes simulations. Robotics, IEEE Transactions on. 24: 1274-1288.

[26] Renda F., Giorelli M., Calisti M., Cianchetti M., Laschi C. (2014) Dynamic model of a multibending soft robot arm driven by cables. Robotics, IEEE Transactions on. 30: 1109-1122.

[27] Linn J., Lang H., Tuganov A. (2013) Geometrically exact Cosserat rods with Kelvin-Voigt type viscous damping. Mech. Sci. 4: 79-96.

[28] Fossen T. I. (1999) Guidance and control of ocean vehicles. Wiley.

[29] Johnson W., Soden P. D., Trueman E. R. (1972) A study in jet propulsion: an analysis of the motion of the squid, Loligo Vulgaris. J. of Exp. Bio. 56: 155-165.

[30] Renda F., Giorgio Serchi F., Boyer F., Laschi C. (2015) Propulsion and elastodynamic model of an underwater shell-like soft robot. Robotics and Automation (ICRA), 2013 IEEE International Conference on. 2015 May 26-30; Seattle: in Preparation.

[31] O'Dor R. K., Webber D. M. (1986) The constraints on cephalopods: why squid aren't fish. Canadian Journal of Zoology. 64: 1591-1605.

[32] Anderson E. J., Grosenbaugh M. A. (2005) Jet flow in steadily swimming adult squid. J. of Exp. Bio. 208: 1125-1146.
Appendix

In this appendix, the derivation of the constitutive equation (19) is outlined. With the equations (9) and (10), we have defined: **n**<sup>1</sup>=N<sub>X</sub>**a** + H**b**, **n**<sup>2</sup>=N<sub>φ</sub>**e**<sub>φ</sub>, **m**<sup>1</sup>=-M<sub>X</sub>**e**<sub>φ</sub> and **m**<sup>2</sup>=M<sub>φ</sub>**a**. Thanks to equations (15) and (16), we find the following vector equations:

$$
\begin{align}
N_X \mathbf{a} + H \mathbf{b} &= (\tilde{n}^{11} + \lambda_1^1 \tilde{m}^{11} + \lambda_2^1 \tilde{m}^{12}) \mathbf{r}' + \nonumber \\
&\quad (\tilde{n}^{21} + \lambda_1^2 \tilde{m}^{11} + \lambda_2^2 \tilde{m}^{12}) \mathbf{r}' + (\tilde{q}^1 + \lambda_1^3 \tilde{m}^{11} + \lambda_2^3 \tilde{m}^{12}) \mathbf{b} \nonumber \\
N_\phi \mathbf{e}_\phi &= (\tilde{n}^{12} + \lambda_1^1 \tilde{m}^{21} + \lambda_2^1 \tilde{m}^{22}) \mathbf{r}' + \nonumber \\
&\quad (\tilde{n}^{22} + \lambda_1^2 \tilde{m}^{21} + \lambda_2^2 \tilde{m}^{22}) \mathbf{r}' + (\tilde{q}^2 + \lambda_1^3 \tilde{m}^{21} + \lambda_2^3 \tilde{m}^{22}) \mathbf{b} \tag{27} \\
-M_X \mathbf{e}_\phi &= \mathbf{b} \times (\tilde{m}^{11} \mathbf{r}' + \tilde{m}^{21} \mathbf{r}' + \tilde{m}^{31} \mathbf{b}) \nonumber \\
M_\phi \mathbf{a} &= \mathbf{b} \times (\tilde{m}^{12} \mathbf{r}' + \tilde{m}^{22} \mathbf{r}' + \tilde{m}^{32} \mathbf{b}) \nonumber
\end{align}
$$

Equation (27) can be simplified by noticing, from equation (14), that: $\tilde{n}^{12} = \tilde{n}^{21} = \tilde{m}^{12} = \tilde{m}^{21} = \tilde{q}^2 = 0$. Furthermore, since $\mathbf{b}' = \tilde{\mathbf{k}}_1 \mathbf{b} = -\mu \mathbf{a}$ and $\mathbf{b}' = \tilde{\mathbf{k}}_2 \mathbf{b} = -(\sin(\theta)/r^0) \mathbf{e}_\phi'$, from equation (17) we have:

$$
\begin{equation}
\begin{aligned}
-\mu \mathbf{a} &= \lambda_1^1 \mathbf{r}' + \lambda_1^2 \mathbf{r}' + \lambda_1^3 \mathbf{b} \\
&\quad - (\sin(\theta)/r^0) \mathbf{e}_\phi = \lambda_2^1 \mathbf{r}' + \lambda_2^2 \mathbf{r}' + \lambda_2^3 \mathbf{b}
\end{aligned}
\tag{28}
\end{equation}
$$

from which we can immediately find: $\lambda_1^2 = \lambda_1^1 = \lambda_2^3 = 0$. Finally, replacing $\mathbf{r}'$ with $\lambda \mathbf{a} + \beta \mathbf{b}$ and $\mathbf{r}'$ with $(r/r^0)\mathbf{e}_\phi'$, equation (27) becomes:

$$
\begin{align}
N_X \mathbf{a} + H \mathbf{b} &= (\tilde{n}^{11} + \lambda_1^1 \tilde{m}^{11}) (\lambda \mathbf{a} + \beta \mathbf{b}) + (\tilde{q}^1 + \lambda_1^3 \tilde{m}^{11}) \mathbf{b} \nonumber \\
N_\phi \mathbf{e}_\phi &= (\tilde{n}^{22} + \lambda_2^2 \tilde{m}^{22}) (r/r^0) \mathbf{e}_\phi \nonumber \\
-M_X \mathbf{e}_\phi &= \lambda \tilde{m}^{11} \mathbf{e}_\phi \nonumber \\
M_\phi \mathbf{a} &= -(r/r^0) \tilde{m}^{22} \mathbf{a} \tag{29}
\end{align}
$$

and equation (28) becomes:

$$
\begin{equation}
\begin{aligned}
-\mu \mathbf{a} &= \lambda_1^1 (\lambda \mathbf{a} + \beta \mathbf{b}) + \lambda_1^3 \mathbf{b} \\
-(\sin(\theta)/r^0) \mathbf{e}_\phi &= \lambda_2^2 (r/r^0) \mathbf{e}_\phi
\end{aligned}
\tag{30}
\end{equation}
$$

At this point, from the last equivalence, we find: $\lambda_1^1 = -\mu / \lambda$, $\lambda_2^2 = -\sin(\theta)/r$ and $\lambda_1^3 = (\mu\beta)/\lambda$. Substituting these results into (29), we have:

$$
\begin{align}
N_X \mathbf{a} + H \mathbf{b} &= (\lambda \tilde{n}^{11} - \mu \tilde{m}^{11}) \mathbf{a} + (\tilde{q}^1 + \beta \tilde{n}^{11}) \mathbf{b} \nonumber \\
N_\phi \mathbf{e}_\phi &= \left( \frac{r}{r^0} \tilde{n}^{22} - \frac{\sin(\theta)}{r^0} \tilde{m}^{22} \right) \mathbf{e}_\phi \tag{31} \\
-M_X \mathbf{e}_\phi &= \lambda \tilde{m}^{11} \mathbf{e}_\phi \nonumber \\
M_\phi \mathbf{a} &= -(r/r^0) \tilde{m}^{22} \mathbf{a} \nonumber
\end{align}
$$

which leads to equation (18).

Now we are in the position of deducing equation (19) by
calculating the remaining terms of (14).

$$
\begin{align}
\tilde{n}^{11} &= \frac{2Eh}{1-\nu^2} (e_{11} + \nu e_{22}) + \frac{6\nu h}{1-\nu^2} (\dot{e}_{11} + \nu \dot{e}_{22}) \nonumber \\
\tilde{n}^{22} &= \frac{2Eh}{1-\nu^2} (e_{22} + \nu e_{11}) + \frac{6\nu h}{1-\nu^2} (\dot{e}_{22} + \nu \dot{e}_{11}) \nonumber \\
\tilde{m}^{11} &= \frac{2EhJ}{1-\nu^2} (d_{11} + \nu d_{22}) + \frac{6\nu hJ}{1-\nu^2} (\dot{d}_{11} + \nu \dot{d}_{22}) \tag{32} \\
\tilde{m}^{22} &= \frac{2EhJ}{1-\nu^2} (d_{22} + \nu d_{11}) + \frac{6\nu hJ}{1-\nu^2} (\dot{d}_{22} + \nu \dot{d}_{11}) \nonumber \\
\tilde{q}^1 &= 2Gh\beta + 2\nu h\dot{\beta} \nonumber
\end{align}
$$

Which, substituted into equation (18), returns the constit-
utive equation (19).