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

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 $ℝ^{6}$ in se(3).

Converts $ℝ^{3}$ in so(3).

$∊ ℝ$ Time.

X $∊ ℝ$ Reference arc-length parametrization.

ϕ ∊S¹ Angle of revolution.

$(e_{1}, e_{2}, e_{3})$ $∊ ℝ^{3} \times ℝ^{3} \times ℝ^{3}$ Ambient reference frame.

$(e_{r}, e_{ϕ}, e_{3})$ $(X, ϕ, t) ∊ ℝ^{3} \times ℝ^{3} \times ℝ^{3}$ Frame attached to the point (X,ϕ).

$(a, b, - e_{ϕ})$ $(X, ϕ, t) ∊ ℝ^{3} \times ℝ^{3} \times ℝ^{3}$ Director orthogonal frame.

g(X,ϕ,t)∊SE(3) Configuration matrix.

ξ(X, t)∊se(3) Local deformation twist vector.

η(X, t)∊se(3) Local velocity twist vector.

ℱ] $(X, t) ∊ s e (3) *$ Wrench vector.

R(X,ϕ,t) ∊ SO(3) Orientation matrix.

$r$ $(X, ϕ, t) ∊ ℝ^{3}$ Position vector.

$g$ $(X, t) ∊ ℝ^{3}$ Linear strain.

$k$ $(X, t) ∊ ℝ^{3}$ Curvature vector.

$v$ $(X, t) ∊ ℝ^{3}$ Linear velocity.

w $(X, t) ∊ ℝ^{3}$ Angular velocity vector.

c $(X, t) ∊ ℝ^{3}$ Added mass load.

r $(X, t) ∊ ℝ^{3}$ Radius.

z $(X, t) ∊ ℝ$ Altitude.

θ(X,t) ∊ S¹ Fibre angle.

λ $(X, t) ∊ ℝ^{+}$ Tangential strain.

β $(X, t) ∊ ℝ$ Thickness strain.

μ $(X, t) ∊ ℝ$ Curvature function.

v_a $(X, t) ∊ ℝ$ Tangential velocity.

v_b $(X, t) ∊ ℝ$ Perpendicular velocity.

w $(X, t) ∊ ℝ$ Angular velocity function.

N_X $(X, t) ∊ ℝ$ Internal tangential force along a.

H $(X, t) ∊ ℝ$ Internal thickness force.

M_X $(X, t) ∊ ℝ$ Internal torque along $- e_{ϕ}$ .

$N_{ϕ}$ $(X, t) ∊ ℝ$ Internal tangential force along $- e_{ϕ}$ .

$M_{ϕ}$ $(X, t) ∊ ℝ$ Internal torque along a.

f_a $(X, t) ∊ ℝ$ External tangential force.

f_b $(X, t) ∊ ℝ$ External perpendicular force.

l $(X, t) ∊ ℝ$ External torque.

u $(t) ∊ ℝ$ Mantle position.

v $(t) ∊ ℝ$ 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 fishlike 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 resemblant 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 essentially relies on an apriori kinematics of the shell deformation, 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 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.

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 capturing the mechanics of the deformation of the soft tissues during propulsion. Ideally, these models would be analytical models allowing us to idealize the principles of this mode of locomotion. Here we propose a numerical coupling 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 terrrestrial 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 measures 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 ρ.

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:

𝒞 : = {(X^{1}, X^{2}) ∊ B \mapsto (r, b) ∊ ℝ^{3} \times S^{2}}

(1)

where $B \subset ℝ^{2}$ defines the material domain of the mid-surface, r represents the field of position of the points of B and b stands for the field of unit vectors attached to the shell fibres, i.e., the directors living in the two dimensional unit sphere S². In other words, every point of a geometrically exact Reissner shell can translate in space (r) and perform both pitch and roll (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 (e₁,e₂,e₃), where e₃ supports the symmetry axis. Denoting by ϕ the angle of revolution and by X the material coordinate along the meridian curves, the orthogonal basis (e_r,e_ϕ,e₃) fixed to the material point of abscissa X and meridian ϕ is defined by:

g_{1} ∊ S E (3) :

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

\exp ({\tilde{e}}_{3} ϕ) = (\begin{matrix} \cos (ϕ) & - \sin (ϕ) & 0 \\ \sin (ϕ) & \cos (ϕ) & 0 \\ 0 & 0 & 1 \end{matrix})

is a rotation ϕ around the axis e₃ and $r : (X, t) ∊ (0, L) \times$ $[0, \infty) \mapsto r (X, t) =$ ${(\cos (ϕ) r (X, t), \sin (ϕ) r (X, t), z (X, t))}^{T}$ is a curve that define the profile of the shell and for which, r(.) and z(.) 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 (e_r, e₃,−e_ϕ) adding the transformation

g_{2} ∊ S E (3) : g_{2} = (\begin{matrix} \exp ({\tilde{e}}_{r} π / 2) & 0 \\ 0 & 1 \end{matrix})

Now let us call θ(X, t) the angle between e₃ and the shell fibre located at any X along the ϕ-meridian, then the so-called director orthogonal frame (a,b,−e_ϕ) is defined at each instant t, by:

\begin{array}{l} g_{3} ∊ S E (3) : (X) ∊ [0, L] \mapsto \\ g_{3} (X) = (\begin{matrix} \exp (- {\tilde{e}}_{ϕ} θ) & 0 \\ 0 & 1 \end{matrix}) \end{array}

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

g ∊ S E (3) : g (X, ϕ) = g_{1} g_{2} g_{3} = (\begin{matrix} R & r \\ 0 & 1 \end{matrix})

where $R = \exp ({\tilde{e}}_{3} ϕ) \exp ({\tilde{e}}_{r} π / 2) \exp (- {\tilde{e}}_{ϕ} θ)$ has been defined.

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

\bar{𝒞} : = {(X, ϕ) ∊ \bar{B} \mapsto g ∊ S E (3)}

(2)

where, referring to the more general context of (1), X ¹ = X and X² have been replaced by ϕ $(\bar{B} \subset ℝ \times S^{1})$ with $X^{2} = r^{o} (X) ϕ$ , 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, ϕ, t)$ is represented by two vector field: $\hat{ξ_{1}} (X, t) = g^{- 1} \partial g / \partial X = g^{- 1} g^{'}$ and $\hat{ξ_{2}} (X, t) = g^{- 1} \partial g / \partial (r^{o} ϕ) = g^{- 1} g^{≀}$ . The hat represents the isomorphism between the twist vector space $ℝ^{6}$ and the Lie algebra se(3). Below their components are specified:

\begin{matrix} \hat{ξ_{1}} = (\begin{matrix} {\tilde{k}}_{1} & g_{1} \\ 0 & 0 \end{matrix}) ∊ s e (3) \\ \hat{ξ_{2}} = (\begin{matrix} {\tilde{k}}_{2} & g_{2} \\ 0 & 0 \end{matrix}) ∊ s e (3) \\ ξ_{1} = {(k_{1}^{T}, g_{1}^{T})}^{T} = {(0,0, μ, λ, β,0)}^{T} ∊ ℝ^{6} \\ ξ_{2} = {(k_{2}^{T}, g_{2}^{T})}^{T} = {(\frac{\sin (θ)}{r^{o}}, \frac{\cos (θ)}{r^{o}},0,0,0, \frac{- r}{r^{o}})}^{T} ∊ ℝ^{6} \end{matrix}

where the curvature μ and the longitudinal and transversal strains (λ, β) have been defined as:

\begin{matrix} μ = θ^{'} \\ λ = \cos (θ) r^{'} + \sin (θ) z^{'} \\ β = \cos (θ) z^{'} - \sin (θ) r^{'} \end{matrix}

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

\begin{matrix} \hat{η} = (\begin{matrix} \tilde{w} & v \\ 0 & 0 \end{matrix}) ∊ s e (3) \\ η = {(w^{T}, v^{T})}^{T} = {(0,0, w, v_{a}, v_{b},0)}^{T} ∊ ℝ^{6} \end{matrix}

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

\begin{matrix} w = \dot{θ} \\ v_{a} = \cos (θ) \dot{r} + \sin (θ) \dot{z} \\ v_{b} = \cos (θ) \dot{z} - \sin (θ) \dot{r} \end{matrix}

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

\dot{g} = g \hat{η}

(3)

that outlining the components becomes:

\begin{array}{l} \dot{θ} = w \\ \dot{r} = \cos (θ) v_{a} - \sin (θ) v_{b} \\ \dot{z} = \sin (θ) v_{a} + \cos (θ) v_{b} \end{array}

(4)

2.2 Compatibility Equations

It has been shown above that $g^{'} = g \hat{ξ_{1}}$ . By taking the derivative of this equation with respect to time and recalling that $\dot{g} = g \hat{η}$ , we obtain the following compatibility equation between the velocity and the deformation variables: $\hat{{\dot{ξ}}_{1}} = \hat{η^{'}} + \hat{ξ_{1}} \hat{η} - \hat{η} \hat{ξ_{1}}$ . In terms of twist vectors this can be written as:

Figure 2.

Axisymmetric shell kinematics (left). Profile view, beam-like parametrization (right).

\begin{array}{l} {\dot{ξ}}_{1} = η^{'} + a d_{ξ_{1}} (η) \end{array}

(5)

where

a d_{ξ_{1}} = (\begin{matrix} {\tilde{k}}_{1} & 0 \\ {\tilde{g}}_{1} & {\tilde{k}}_{1} \end{matrix})

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

Outlining the components we obtain:

\begin{array}{l} \dot{μ} = w^{'} \\ \dot{λ} = v_{a^{'}} + β w - μ v_{b} \\ \dot{β} = v_{b^{'}} - λ w + μ v_{a} \end{array}

(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 R³ is induced on the surface in any of its points. In our case, this is defined in each point (X,ϕ) 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₁, g₂)(X,ϕ). For our revolution surface it can be simply expressed as:

h = (\begin{matrix} h_{11} (X) & 0 \\ 0 & h_{22} (X) \end{matrix}) = (\begin{matrix} λ^{2} + β^{2} & 0 \\ 0 & \frac{r^{2}}{r^{o 2}} \end{matrix})

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_{α γ} = g_{α} \cdot {\tilde{k}}_{γ} 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{matrix} k_{11} (X) & 0 \\ 0 & k_{22} (X) \end{matrix}) = (\begin{matrix} - μ λ & 0 \\ 0 & - \frac{r \sin (θ)}{r^{o 2}} \end{matrix})

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

e = \frac{1}{2} (h - h^{o}) = \frac{1}{2} (\begin{matrix} λ^{2} + β^{2} - 1 & 0 \\ 0 & \frac{r^{2}}{r^{o 2}} - 1 \end{matrix})

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{matrix} g_{1} \cdot b - g_{1}^{o} \cdot b^{o} \\ g_{2} \cdot b - g_{2}^{o} \cdot b^{o} \end{matrix}) = (\begin{matrix} β \\ 0 \end{matrix})

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

d = k - k^{o} = (\begin{matrix} μ^{o} - μ λ & 0 \\ 0 & \frac{\sin (θ^{o})}{r^{o}} - \frac{r \sin (θ)}{r^{o 2}} \end{matrix})

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., β^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{array}{l} \frac{\partial (j n^{α})}{j \partial X^{α}} + \bar{n} = 2 ρ h \ddot{r} \\ \frac{\partial (j m^{α})}{j \partial X^{α}} + \frac{\partial r}{\partial X^{α}} \times n^{α} + \bar{m} = 2 ρ h b \times \ddot{b} \end{array}

(7)

where $j = \sqrt{\det (h)} = r / r^{o} \sqrt{λ^{2} + β^{2}}$ , while the vectors n^α and m^α are, respectively, the resultant of internal stress forces and couples applied by the left part (x < X^α) into the right part (x ≥ X ^α) of the shell, across the section X ^α, normalized with the surface Jacobian j. $\bar{n}$ and $\bar{m}$ are the external force and couple per unit of mid-surface area. For the repeated α 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 (r, 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:

ℳ \dot{η} = 1 / j (j ℱ_{i}^{1}) - a d_{ξ_{α}}^{*} (ℱ_{i}^{α}) + ℱ_{e} + a d_{η}^{*} (ℳ η)

(8)

Here $ℱ_{i}^{α}$ is the wrench of internal forces, $ℱ_{e}$ is the external wrench of distributed applied forces, M is the screw inertia matrix and

a d_{ξ_{α}}^{*} = (\begin{matrix} {\tilde{k}}_{α}^{T} & {\tilde{g}}_{α}^{T} \\ 0 & {\tilde{k}}_{α}^{T} \end{matrix}) a d_{η}^{*} = (\begin{matrix} {\tilde{w}}^{T} & {\tilde{v}}^{T} \\ 0 & {\tilde{w}}^{T} \end{matrix})

are the co-adjoint maps.

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

ℱ_{i}^{1} (X, t) = {(0,0, M_{X}, N_{X}, H,0)}^{T} ∊ R^{6}

(9)

ℱ_{i}^{2} (X, t) = {(M_{ϕ},0,0,0,0, - N_{ϕ})}^{T} ∊ R^{6}

(10)

ℱ_{e} (X, t) = {(0,0, l, f_{a}, f_{b},0)}^{T} ∊ R^{6}

(11)

and the screw inertia matrix is equal to: $ℝ^{6} \otimes ℝ^{6} э$ $ℳ = d i a g (ρ J, ρ J_{b}, ρ J,2 ρ h,2 ρ h,2 ρ h)$ .

In the equations above J and Jb are the second moment of the cross sectional line equal to J =h ²/3, J_b ~ 0.

By outlining the components, we obtain:

\begin{array}{l} ρ J \dot{w} = 1 / j {(j M_{X})}^{'} + λ H - β N_{X} - \frac{\cos (θ)}{r^{o}} M_{ϕ} + l \\ 2 ρ h {\dot{v}}_{a} = 1 / j {(j N_{X})}^{'} - μ H - \frac{\cos (θ)}{r^{o}} N_{ϕ} + f_{a} + 2 ρ h w v_{b} \\ 2 ρ h {\dot{v}}_{b} = 1 / j {(j H)}^{'} + μ N_{X} + \frac{\sin (θ)}{r^{o}} N_{ϕ} + f_{b} - 2 ρ h w v_{a} \end{array}

(12)

Figure 3.

Internal (dark) and external (red) loads exerted in (X,ϕ)

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{matrix} {\tilde{n}}^{β α} = \frac{2 E h}{1 - ν^{2}} H^{β α γ δ} e_{γ δ}, {\tilde{m}}^{β α} = \frac{2 E h J}{1 - ν^{2}} H^{β α γ δ} d_{γ δ} \\ {\tilde{q}}^{α} = 2 G h h^{o α β} s_{β}, \end{matrix}

(13)

where E is the Young modulus, G the shear modulus and ν is the Poisson modulus; h^oαβ are the elements of $h^{o - 1} \equiv I^{- 1} = I$ and H^βαγδ define a four order (Hook-like) tensor given by:

H^{β α γ δ} = [ν (h^{o β α} h^{o γ δ}) + 1 / 2 (1 - ν) (h^{o β γ} h^{o α δ} + h^{o β δ} h^{o α γ})] .

and, in our case, resulting in:

H^{β α γ δ} = {\begin{matrix} 1 & β α γ δ = 1111,2222 \\ ν & β α γ δ = 1122,2211 \\ 0 & o t h e r w i s e \end{matrix}}

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{matrix} {\tilde{n}}^{β α} = \frac{2 E h}{1 - ν^{2}} H^{β α γ δ} e_{γ δ} + \frac{6 υ h}{1 - ν^{2}} H^{β α γ δ} {\dot{e}}_{γ δ} \\ {\tilde{m}}^{β α} = \frac{2 E h J}{1 - ν^{2}} H^{β α γ δ} d_{γ δ} + \frac{6 υ h J}{1 - ν^{2}} H^{β α γ δ} {\dot{d}}_{γ δ} \\ {\tilde{q}}^{α} = 2 G h h^{o α β} s_{β} + 2 υ h h^{o α β} {\dot{s}}_{β} \end{matrix}

(14)

where the Young modulus and the shear modulus have been replaced by the shear viscosity constants 3ν and ν, 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}}^{α β}$ ${\tilde{m}}^{α β}$ and ${\tilde{q}}^{α}$ respectively, stand for the components of the effective resultant traction, couple and shearing stresses¹, which are related to the stresses of (7) through the relations:

n^{α} = ({\tilde{n}}^{β α} + λ_{μ}^{β} {\tilde{m}}^{α μ}) \frac{\partial r}{\partial X^{β}} + ({\tilde{q}}^{α} + λ_{μ}^{3} {\tilde{m}}^{α μ}) b

(15)

m^{α} = b \times {\tilde{m}}^{α} = b \times ({\tilde{m}}^{β α} \frac{\partial r}{\partial X^{β}} + {\tilde{m}}^{3 α} b)

(16)

where the functions λ_μ^β and λ_μ ³ are defined by the following equation:

b_{, X^{μ}} = λ_{μ}^{β} \frac{\partial r}{\partial X^{β}} + λ_{μ}^{3} b

(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:

\begin{matrix} M_{X} = - λ {\tilde{m}}^{11}, M_{ϕ} = - \frac{r}{r^{o}} {\tilde{m}}^{22}, \\ N_{X} = λ {\tilde{n}}^{11} - μ {\tilde{m}}^{11}, N_{ϕ} = \frac{r}{r^{o}} {\tilde{n}}^{22} - \frac{\sin (θ)}{r^{o}} {\tilde{m}}^{22}, \\ H = {\tilde{q}}^{1} + β {\tilde{n}}^{11} . \end{matrix}

(18)

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

\begin{array}{l} N_{X} = \frac{2 E h}{1 - ν^{2}} [λ (e_{11} + ν e_{22}) - J μ (d_{11} + ν d_{22})] \\ + \frac{6 υ h}{1 - ν^{2}} [λ ({\dot{e}}_{11} + ν {\dot{e}}_{22}) - J μ ({\dot{d}}_{11} + ν {\dot{d}}_{22})] \\ N_{ϕ} = \frac{2 E h}{1 - ν^{2}} [\frac{r}{r^{o}} (e_{22} + ν e_{11}) - J \frac{\sin (θ)}{r^{o}} (d_{22} + ν d_{11})] \\ + \frac{6 υ h}{1 - ν^{2}} [\frac{r}{r^{o}} ({\dot{e}}_{22} + ν {\dot{e}}_{11}) - J \frac{\sin (θ)}{r^{o}} ({\dot{d}}_{22} + ν {\dot{d}}_{11})] \\ H = 2 h β [G + \frac{E}{1 - ν^{2}} (e_{11} + ν e_{22})] \\ + 2 h \dot{β} [υ + \frac{6 υ}{1 - ν^{2}} ({\dot{e}}_{11} + ν {\dot{e}}_{22})] \\ M_{X} = - \frac{2 E h J}{1 - ν^{2}} λ (d_{11} + ν d_{22}) - \frac{6 υ h J}{1 - ν^{2}} λ ({\dot{d}}_{11} + ν {\dot{d}}_{22}) \\ M_{ϕ} = - \frac{2 E h J}{1 - ν^{2}} \frac{r}{r^{o}} (d_{22} + ν d_{11}) - \frac{6 υ h J}{1 - ν^{2}} \frac{r}{r^{o}} ({\dot{d}}_{22} + ν {\dot{d}}_{11}) \end{array}

(19)

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.

ℱ_{e} (X, t) = {(0,0,0, d_{a},0,0)}^{T} + {(0,0,0, c^{T})}^{T} .

(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) = - ρ_{w} C_{a} (v_{a} + \sin (θ) v) | v_{a} + \sin (θ) v |

(21)

where $ρ_{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 e₃ (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.

\begin{array}{l} c (X, t) = - \frac{\partial (ρ F (v + R^{T} \dot{u}))}{\partial t} = \\ - ρ F (\dot{v} + R^{T} \ddot{u}) + w \times ρ F (v + R^{T} \dot{u}) \end{array}

(22)

where $F (X) ∊ ℝ^{3} \otimes ℝ^{3}$ is a tensor that incorporates the geometric and hydrodynamics factors. It is equal to F =diag(0,2h B_b,0), where B_b is the hydrodynamic correction coefficient for the added mass model. The vector $\dot{u}$ is the swimming velocity of the mantle with respect to the fixed Galilean frame (e₁, e₂, e₃) 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,

(ρ_{w} V + ρ V_{m}) \dot{v} = - \frac{1}{2} ρ_{w} C_{d} A_{r e f} v | v | - B \dot{v} + q ρ_{w} \dot{V}

(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{α_{0}}{2 - α_{0}} m

with

α_{0} = \frac{2 (1 - e^{2})}{e^{3}} (\frac{1}{2} l n \frac{1 + e}{1 - e} - e) e = 1 - {(b / a)}^{2}

and

m = ρ_{w} V + ρ 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 $ρ_{w} \dot{V}$ occurring within the robot due to the collapse of the elastic chamber.

Figure 4.

Depiction of the elements accounted for in the dynamics model of the pulsed-jet propulsion

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 ρ_{w} \dot{V} = - \frac{ρ_{w} C_{f} \dot{V} | \dot{V} |}{A_{n}}

(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:

\begin{matrix} A_{n} (t) = π r {(L, t)}^{2} A_{r e f} (t) = π m a x {(r (\cdot, t))}^{2} \\ V_{m} = - \int_{0}^{z (L)} 4 π r^{o} h d z = - 4 π h \int_{0}^{L} r^{o} z^{o}' d X \end{matrix}

while the mantle inner volume and his time derivative are:

\begin{matrix} V (t) = - \int_{0}^{z (L)} π r^{2} d z = - π \int_{0}^{L} r^{2} z^{'} d X \\ \dot{V} (t) = - π \int_{0}^{L} (2 r \dot{r} z^{'} + r^{2} \dot{z}') d X \end{matrix}

where,

\begin{matrix} z^{'} = \sin (θ) λ + \cos (θ) β \dot{r} = \cos (θ) v_{a} - \sin (θ) v_{b} \\ \dot{z}' = \sin (θ) (\dot{λ} - w β) + \cos (θ) (\dot{β} + w λ) \end{matrix}

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₃-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₁(ϕ,t) with the following

\begin{matrix} \bar{g_{1}} ∊ S E (3) : \\ (ϕ, X) \mapsto \bar{g_{1}} (ϕ, X) = (\begin{matrix} \exp ({\tilde{e}}_{3} ϕ) & (r + u) \\ 0 & 1 \end{matrix}) \end{matrix}

where 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}} g_{2} g_{3}$ and the corresponding velocity twist vector field $\bar{η} (X, t)$ is

{\bar{g}}^{- 1} \bar{g} = \hat{\bar{η}} = (\begin{matrix} \tilde{w} & (v + R^{T} \dot{u}) \\ 0 & 0 \end{matrix}) ∊ s e (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{η} (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 show in figure 5.

Figure 5.

Scheme of the coupling between the shell model and the jet propulsion model for every time step

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{array}{l} \dot{θ} = w \\ \dot{r} = \cos (θ) v_{a} - \sin (θ) v_{b} \\ \dot{z} = \sin (θ) v_{a} + \cos (θ) v_{b} + v \\ \dot{μ} = w^{'} \\ \dot{λ} = v_{a^{'}} + β w - μ v_{b} \\ \dot{β} = v_{b^{'}} - λ w + μ v_{a} \\ \dot{w} = \frac{{(j M_{X})}^{'} / j + λ H - β N_{X} - \cos (θ) M_{ϕ} / r^{o} + l}{ρ J} \\ {\dot{v}}_{a} = \frac{{(j N_{X})}^{'} / j - μ H - \cos (θ) N_{ϕ} / r^{o} + f_{a}}{2 ρ h} + w v_{b} \\ {\dot{v}}_{b} = \frac{{(j H)}^{'} / j + μ N_{X} + \sin (θ) N_{ϕ} / r^{o} + f_{b} - 2 ρ h w v_{a}}{2 h (ρ + ρ B_{b})} \\ \dot{u} = v \\ \dot{v} = \frac{q ρ_{w} \dot{V} - \frac{1}{2} ρ_{w} C_{d} A_{r e f} v | v |}{ρ_{w} V + ρ V_{m} + B} \end{array}

(25)

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:

η (0) = η_{-} = 0

(26)

where $η_{-}$ 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, +}, H (L) = f_{b, +}, M_{X} (L) = - l_{+},

where $f_{a, +}$ $f_{b, +}$ and l ₊ denote the eventual external forces and torque along a, b and −e_ϕ 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 (θ) - f_{b} \sin (θ)$ 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 distribution of the actuation is assumed with the form: P(X)=(P−Po/L)^* X + Po (figure 6).

Figure 6.

A sketch plot of the two linear distribution of the external actuation used in the simulations

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 1.

Families of Variable Parameters for the External Actuation

Parameter	Value 1	test #	Value 2	test #	Value 3	test #
E and	110 Pa	1-12	1100 Pa	13-24	11000 Pa	25-36
P	0.3 N		3 N		30 N
T and	1 s	1-4	1.5 s	5-8	1.5 s	9-12
		13-16		17-20		21-24
T_c	0.5 s	25-28	0.5 s	29-32	1 s	33-36
P_o	P	odd	P/2	even
	test #	test #
Ellipsoid	1-2, 5-6,	9-10, 13-14
	17-18,	25-26,
	21-22	29-30
	33-34
Cone	3-4, 7-8,	11-12,
		15-16
	19-20,	27-28,
	23-23	31-32
	35-36

Table 2.

Cone and Ellipsoid Parameters

Parameter	Value
L	80 cm
h	2 mm
R_c	5 cm
Rn	15 cm
a	50 cm
b	20 cm

Table 3.

Load Parameters

Parameter	Value
μ	50 Pa • s
v	0
ρ	1.08 kg/dm³
_ρw	1.02 kg/dm³
C_a	0.01
B_b	1
C_d	1
C_f	1

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 $(ξ_{1} (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, counterbalanced 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 = 15, 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 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 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.

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.

Figure 8.

Snapshots of the ellipsoidal mantle configuration with E = 110e2 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).

Figure 9.

Velocity profiles of the best performers for the three stiffnesses under the external actuation

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 β^o, see 2.3) with a desired time-varying strain field (h^d, k^d and β^d) obtained by the following constrain inputs: $β^{d} (X, t) = 0$ $λ^{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 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.

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 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 e best performer, i.e., a conical mantle with E = 110e2Pa under the T_c = 0.5, T = 1s, are presented. In figure 15 his of 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 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 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.

Table 4.

Families of Variable Parameters for the Internal Actuation

Parameter	Value 1	test #	Value 2	test #	Value 3	test #
E	110 Pa	1-6	1100 Pa	7-12	11000 Pa	13-18
T and	1 s	1-2, 7-8	1.5 s	3-4, 9-10	1.5 s	5-6, 11-12
T_c	0.5 s	13-14	0.5 s	15-16	1 s	17-18
Ellipsoid	odd
Cone	even

Figure 12.

Desired radius functions for the conical (up) and ellipsoidal (down) mantle

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.

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).

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. Furthermore, 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.

Figure 15.

Velocity profiles of the best performers for the three stiffnesses under the internal actuation

In figure 17 the results of the two geometries (conical and ellipsoidal) are isolated. The conical mantle performances 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 = 0.5s w.r.t. the other routines.

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_c⁼¹ is clearly selected.

Figure 17.

Bar plot of the simulations for an ellipsoidal (top) and conical (down) geometry for the internal actuation

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_c = 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 performances 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 fas st jea tt ling 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 rocketlike 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.

Footnotes

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].

Appendix

In this appendix, the derivation of the constitutive equation (19) is outlined. With the equations (9) and (10), we have defined: $n^{1} = N_{X} a + H b$ $n^{2} = N_{ϕ} e_{ϕ}$ $m^{1} = - M_{X} e_{ϕ}$ and $m^{2} = M_{ϕ} a$ . Thanks to equations (15) and (16), we find the following vector equations:

(27)

\begin{array}{l} N_{X} a + H b = ({\tilde{n}}^{11} + λ_{1}^{1} {\tilde{m}}^{11} + λ_{2}^{1} {\tilde{m}}^{12}) r' + \\ ({\tilde{n}}^{21} + λ_{1}^{2} {\tilde{m}}^{11} + λ_{2}^{2} {\tilde{m}}^{12}) r^{≀} + ({\tilde{q}}^{1} + λ_{1}^{3} {\tilde{m}}^{11} + λ_{2}^{3} {\tilde{m}}^{12}) b \\ N_{ϕ} e_{ϕ} = ({\tilde{n}}^{12} + λ_{1}^{1} {\tilde{m}}^{21} + λ_{2}^{1} {\tilde{m}}^{22}) r' + \\ ({\tilde{n}}^{22} + λ_{1}^{2} {\tilde{m}}^{21} + λ_{2}^{2} {\tilde{m}}^{22}) r^{≀} + ({\tilde{q}}^{2} + λ_{1}^{3} {\tilde{m}}^{21} + λ_{2}^{3} {\tilde{m}}^{22}) b \\ - M_{X} e_{ϕ} = b \times ({\tilde{m}}^{11} r' + {\tilde{m}}^{21} r^{≀} + {\tilde{m}}^{31} b) \\ M_{ϕ} a = b \times ({\tilde{m}}^{12} r' + {\tilde{m}}^{22} r^{≀} + {\tilde{m}}^{32} b) \end{array}

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 $b' = {\tilde{k}}_{1} b = - μ a$ , and $b^{'} = {\tilde{k}}_{2} b = - (\sin (θ) / r^{0}) e_{ϕ}$ from equation (17) we have:

(28)

\begin{array}{l} - μ a = λ_{1}^{1} r' + λ_{1}^{2} r^{≀} + λ_{1}^{3} b \\ - (\sin (θ) / r^{0}) e_{ϕ} = λ_{2}^{1} r' + λ_{2}^{2} r^{≀} + λ_{2}^{3} b \end{array}

from which we can immediately find: λ₁²=λ₂¹=λ₂³ = 0. Finally, replacing r’ with λa + βb and $r^{≀}$ with (r / r⁰)e_ϕ equation (27) becomes:

(29)

\begin{array}{l} N_{X} a + H b = ({\tilde{n}}^{11} + λ_{1}^{1} {\tilde{m}}^{11}) (λ a + β b) + ({\tilde{q}}^{1} + λ_{1}^{3} {\tilde{m}}^{11}) b \\ N_{ϕ} e_{ϕ} = ({\tilde{n}}^{22} + λ_{2}^{2} {\tilde{m}}^{22}) (r / r^{0}) e_{ϕ} \\ - M_{X} e_{ϕ} = λ {\tilde{m}}^{11} e_{ϕ} \\ M_{ϕ} a = - (r / r^{0}) {\tilde{m}}^{22} a \end{array}

and equation (28) becomes:

(30)

\begin{array}{l} - μ a = λ_{1}^{1} (λ a + β b) + λ_{1}^{3} b \\ - (\sin (θ) / r^{0}) e_{ϕ} = λ_{2}^{2} (r / r^{0}) e_{ϕ} \end{array}

At this point, from the last equivalence, we find: λ ₁¹=-μ/λ, λ₂² = - sin(θ) / r and λ₁³ = (μβ) /λ. Substituting these results into (29), we have:

(31)

\begin{array}{l} N_{X} a + H b = (λ {\tilde{n}}^{11} - μ {\tilde{m}}^{11}) a + ({\tilde{q}}^{1} + β {\tilde{n}}^{11}) b \\ N_{ϕ} e_{ϕ} = (\frac{r}{r^{0}} {\tilde{n}}^{22} - \frac{\sin (θ)}{r^{0}} {\tilde{m}}^{22}) e_{ϕ} \\ - M_{X} e_{ϕ} = λ {\tilde{m}}^{11} e_{ϕ} \\ M_{ϕ} a = - (r / r^{0}) {\tilde{m}}^{22} a \end{array}

which leads to equation (18).

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

(32)

\begin{array}{l} {\tilde{n}}^{11} = \frac{2 E h}{1 - ν^{2}} (e_{11} + ν e_{22}) + \frac{6 υ h}{1 - ν^{2}} ({\dot{e}}_{11} + ν {\dot{e}}_{22}) \\ {\tilde{n}}^{22} = \frac{2 E h}{1 - ν^{2}} (e_{22} + ν e_{11}) + \frac{6 υ h}{1 - ν^{2}} ({\dot{e}}_{22} + ν {\dot{e}}_{11}) \\ {\tilde{m}}^{11} = \frac{2 E h J}{1 - ν^{2}} (d_{11} + ν d_{22}) + \frac{6 υ h J}{1 - ν^{2}} ({\dot{d}}_{11} + ν {\dot{d}}_{22}) \\ {\tilde{m}}^{22} = \frac{2 E h J}{1 - ν^{2}} (d_{22} + ν d_{11}) + \frac{6 υ h J}{1 - ν^{2}} ({\dot{d}}_{22} + ν {\dot{d}}_{11}) \\ {\tilde{q}}^{1} = 2 G h β + 2 υ h \dot{β} \end{array}

Which, substituted into equation (18), returns the constitutive equation (19).

References

Trueman

E. R.

Packard

(1968). Motor performances of some cephalopods. J. of Exp. Bio. 49: 495–507.

Bartol

I. K.

Patterson

M. R.

Mann

(2001) Swimming mechanics and behavior of the shallow-water brief squid Lolliguncula brevis. J. of Exp. Bio. 204: 3655–3682.

Sfakiotakis

Kazakidi

Pateromichelakis

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.

Gosline

J. M.

DeMont

M. E.

(1985) Jet-propelled swimming in squids. Sci. Am. 252: 96–103.

Krieg

Mohseni

(2012) New perspectives on collagen fibers in the squid mantle. J. Morphol. 273: 586–95.

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.

Krueger

P. S.

Gharib

(2005) Thrust augmentation and vortex ring evolution in a fully pulsed jet AIAA Journal. 43: 792–801.

Krieg

Mohseni

(2008) Thrust characterization of a bio-inspired vortex ring generator for locomotion of underwater robots. IEEE J. of Ocean Eng. 33: 123–132.

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

Laschi

(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

Laschi

(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

Laschi

(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

Mohseni

(2013) Modelling circulation, impulse, and kinetic energy of starting jets with non-zero radial velocity. J. of Fluid Mech. 719: 488–526.

17.

Giorelli

Giorgio

Serchi F.

Laschi

(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.

Antman

S. S.

(2005) Nonlinear problems of elasticity, 2nd edn. New York: Springer.

20.

Anderson

E. J.

Quinn

DeMont

(2001) Hydrodynamics of locomotion in the squid Loligo pealei. J. of Fluid Mech. 436: 249–266.

21.

Cosserat

(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

(1945) The effect of transverse shear deformation on the bending of elastic plates. ASME Journal of Applied Mechanics. 12: 68–77.

24.

Boyer

Porez

Leroyer

(2010) Poincaré-Cosserat rquations for the Lighthill three-dimensional large amplitude elongated body theory: application to robotics. Journal of Nonlinear Science. 20: 47–79.

25.

Boyer

Porez

Leroyer

Visonneau

(2008) Fast dynamics of an eel-like robot-comparisons With Navier-Stokes simulations. Robotics, IEEE Transactions on. 24: 1274–1288.

26.

Renda

Giorelli

Calisti

Cianchetti

Laschi

(2014) Dynamic model of a multibending soft robot arm driven by cables. Robotics, IEEE Transactions on. 30: 1109–1122.

27.

Linn

Lang

Tuganov

(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

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

Giorgio

Serchi F.

Boyer

Laschi

(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.