A unified multi-soft-body dynamic model for underwater soft robots

Abstract

A unified formulation that accounts for the dynamics of a general class of aquatic multi-body, soft-structured robots is presented. The formulation is based on a Cosserat formalism where the description of the ensemble of geometrical entities, such as shells and beams, gives rise to a multi-soft-body system capable of simulating both manipulation and locomotion. Conceived as an advanced tool for a priori hardware development, n-degree-of-freedom dynamics analysis and control design of underwater, soft, multi-body, vehicles, the model is validated against aquatic locomotion experiments of an octopus-inspired soft unmanned underwater robot. Upon validation, the general applicability of the model is demonstrated by predicting the self-propulsion dynamics of a diverse range of new viable combinations of multi-soft-body aquatic system.

Keywords

Soft robotics bio-inspired robotics marine robotics dynamics

1. Introduction

Fostered by the growing needs of the marine and maritime industry to perform increasingly daunting tasks in always more forbidding environments, a renovated effort is being made nowadays to endow underwater robots with enhanced manoeuvring capabilities. On the one hand this has entailed the revision of traditional systems (Elvander and Hawkes, 2012; Vaganay et al., 2009; Vasilescu et al., 2010), or by improving navigation and positioning systems, i.e. by combining data collected jointly from doppler velocity logger (DVL), GPS, pressure-depth sensors, synthetic aperture sonars, and multi-beam echo sounders (Hover et al., 2012). On the other hand, alternative design criteria have been taken in consideration by capitalizing on the study of water-dwelling organisms.

In recent times underwater robotics has largely benefited from the growing fascination for bioinspired aquatic locomotion and, motivated by the abundance of outstanding feats that aquatic animals display, has started to pave the way for the development of new vehicles capable of feats yet to be seen in commercially available unmanned underwater vehicles (UUVs). Hovering, short radius turning, fast start/slowdown, and low-speed manoeuvring are just few examples that highlight how the design of underwater robots can profit massively from the investigation of the swimming strategies, hydrodynamics, and physiology of aquatic animals.

Several examples exist of aquatic organisms which have been taken as the source of inspiration for designing a trustworthy robotic replica. These include flagellates (Abbott et al., 2009), turtles (Licht et al., 2004), eels (Yu et al., 2012), and, of course, fish. The finned and caudal flapping of fish (e.g. Conte et al., 2010; Saimek and Li, 2001), has gathered the most recognition in the scientific community, in part because of the sound understanding of the underlying physics involved in their locomotion (Colgate and Lynch, 2004).

The locomotion of aquatic organisms and their robotics counterparts commonly involves periodic oscillatory deformations of one or more body parts that, in turn, give rise to the unsteady hydrodynamics responsible for generating thrust. The actuation mechanisms that enable these deformations in the bioinspired water-dwelling robots has, in most cases, entailed the replacement of continuously deforming bodies by reducing the number of degrees of freedom (DOFs) with a finite sequence of rigid links and joints. Alternatively, the compliant nature and the multi-DOF bending capability of the biological counterpart has been accounted for by resorting to continuously deforming soft structures and actuators (Marchese et al., 2014). The effort in designing structurally compliant underwater robots is fostered, to a large extent, by the acknowledgement that safe physical human–robot interaction and manoeuvring in highly perturbed, unstructured scenarios can effectively be handled via the recourse to soft-bodied components rather than by making the control finer (Mortl et al., 2012; Wang and Iida, 2015; Woodman et al., 2012). Indeed, commercial underwater robots are safe traveling in open stretches of sea, but suffer from non-negligible limitations when navigating close to the seabed or in close proximity to submerged structures where unintended impacts must be prevented consistently. The exploitation of soft-bodied vehicles that benefit of the assets from bioinspired propulsion and manipulation systems can provide a viable solution to complex tasks that existing remotely operated underwater vehicles (ROVs) and autonomous underwater vehicles (AUVs) are unfit for such as operations in current-perturbed domains, performing maintenance over the hull of ships and harbors and working in synergy with divers.

This has encouraged the authors to design and develop an innovative class of soft-bodied, bioinspired, underwater robots (Figure 1). These consist of octopus-resembling machines endowed with a number of continuous manipulators (element 4 in Figure 2), and a central unit devoted to thrust generation (element 1 in Figure 2), which essentially defines an underwater multi-limbed soft vehicle (Giorgio-Serchi et al., 2017). The robot is composed by as much as 90% in volume by elastomeric materials and actuation is provided by electric motors and cable transmission thus enabling the robot to profit from a high overall degree of structural compliance. In analogy with its biological source of inspiration, this kind of robot is capable of performing basic manipulation, legged locomotion, and waterborne propulsion. On one_ hand, this kind of design offers a number of assets due to its structure and mode of actuation that have been covered at length by Giorgio-Serchi et al. (2016, 2015), as far as propulsion is concerned, and by Renda et al. (2014) for manipulation. On the other hand, however, the morphology of the robot also requires an ad hoc formulation in order to treat the dynamics of this flexible multi-body system. This represent the focus of the present work. The mechanical system of Figure 1, with its articulated configuration and its combination of flexible and rigid components represents the state-of-the-art paradigm for developing and validating an advanced mathematical framework capable of dealing with such complexity.

Fig. 1.

(a) The soft, multi-body aquatic robot developed by the authors. (b) During testing at sea.

Fig. 2.

A schematic of the soft, multi-body aquatic vehicle developed by the authors. Numbers refer to: (1) pulsed-jet thruster; (2) the nozzle; (3) the cables that drive the shell collapse; (4) the continuum manipulators; (5) the actuators of the manipulators; (6) the actuator of the shell; and (7) the cable that drives manipulator actuation.

With the advancement in bioinspired robots, increasingly sophisticated mathematical models have been developed with the scope of accounting for the growing complexity of such systems (Krieg et al., 2015). As far as aquatic robots are concerned, the swimming routine commonly entails caudal or finned flapping and whole-body undulatory oscillations. An extensive literature exists that accounts for the dynamics entailed with these kinds of swimming routines, as well as with the associated accurate flow features. Whereas a rigorous treatment of these kinds of systems requires the solution of a full hydroelastic problem, physically sound and less numerically intensive approaches have been found to be well suited for the purpose of control and design optimization. These have entailed the reduction of the problem to a coupling between the body, regarded either as a series of rigid links or as continuously deforming beams, and the fluid as a quasi-inviscid one, where reactive inertial terms are computed via potential flow theory and resistive viscous terms are derived from empirically determined coefficients.

This approach has been applied in recent years to the swimming of fish via an extended Lighthill model that expands the large-amplitude elongated body theory (LAEBT) of (Lighthill, 1970) to the case of self-propelled three-dimensional swimming (Boyer et al., 2010). This latter formulation has been employed to encompass the case of non-quiescent ambient flows, such as in the Von Karman vortex street from which a dead fish is shown to extract energy in order to passively propel itself upstream (Candelier et al., 2013). This model accounts for most of fish morphologies (anguilliform, carangiform) by treating the fish body as a nonlinear Cosserat beam in finite transformations (deformations and rotations), i.e. an infinite set of rigid cross-sections (thus modeling, in a way, the fish vertebrae) regarded as continuously stacked along the vertebral axis of the animal. It can be shown (Candelier et al., 2011), by exploiting the tapered shape of the body, that the fluid forces exerted on a beam cross-section only depend on the fluid velocities and acceleration of the water slice that prolongs the beam cross-section in the fish surrounding. Based on this remark, it becomes possible to extend the Newton–Euler-based approach of rigid discrete multi-body systems dynamics to the case of continuous systems where the cross-sections label stand for the body index of the discrete case. In this context, the inverse (computed torque) algorithm of Luh et al. (1980) has been extended to the locomotion of continuous elongated systems in Boyer et al. (2006). Remarkably, the resulting dynamics approach exploits the topology of these continuous systems to design fast dynamics algorithms where the usual recursions of the Newton–Euler algorithms are replaced by ordinary differential equations (o.d.e.s) that are solved forward and/or backward along the beams axisin a global time loop. Furthermore, compared with other approaches, such as those based on the floating frame (Canavin and Likins, 1977), they can naturally tackle the finite deformations observed in soft animals, an advantage that is crucial in the context of this article.

Whereas modern modeling paradigms essentially pertain to the multi-body or continuous approach, the need arises to reconcile these supposedly divergent perspectives into a more general view capable of encompassing both counterparts. In the present work, of which a preliminary version has been presented in Renda et al. (2018), the dichotomy between the multi-rigid-body and the single-continuous-body paradigms is relaxed by expressing the whole-body dynamics of an octopus-like robot via a multi-soft-body formulation. To do so, we expand on the state-of-the-art geometrical models of archetypal elements, such as beams and shells, to construct a unified framework where various appendages are allowed to participate to the dynamics of a single entity. As a mean of validating this construct, the model is employed to replicate the robot depicted in Figures 1 and 2 and compare the simulated with the experimentally observed dynamics during aquatic self-propulsion. The versatility of the mathematical framework introduced is demonstrated by extending it to account for a set of diverse geometrical configurations and actuation routines.

The first section of this paper entails an extended general description of the mathematical frames exploited throughout. In Section 3 this is employed to illustrate the modeling formalism adopted for beams and shell-shaped soft bodies, which are later considered in a unified system (Section 5). The model thus formulated is then validated for the case of a four-limbed, self-propelling, soft robot closely resembling one of the vehicles developed previously by the authors (Section 6). Eventually, in Section 7 we demonstrated how this unified model for multi-soft-body vehicles can be exploited with the purpose of exploring innovative design paradigms by predicting the locomotion performances of two distinct soft underwater vehicles.

2. Model description

The basic structure of the model is made by a rigid body, called rigid root-body, and several soft appendages attached to it through one extreme or boundary of the soft body (beams and axisymmetric shells in this work). The rigid root-body is not kinematically attached to any hard frame, but instead it is free to move in the 3D space, whereas the soft appendages are not connected to each other preventing the realization of closed-loop mechanisms. This kind of multi-body system is said to have a star structure (Selig, 2007).

Let us call $(o, e_{1}, e_{2}, e_{3})$ the inertial frame of the Euclidean space and $(o, E_{1}, E_{2}, E_{3})$ the reference orthogonal frame attached to the rigid root-body whose inertial motion defines the net (rigid overall) motion of the entire system. The configuration space of the root-body is $g_{r} \in SE (3)$ , which maps the inertial onto the reference mobile frame.

Each of the soft appendages are modeled as Cosserat medium, which can be intuitively considered as a continuous staking of a rigid small solid named “microstructure” along one (beam) or two (shell) material dimensions. As a result, the configuration space of such a medium can be intrinsically defined as the set of maps (Figure 3):

C = {g : (X^{1}, X^{2}, \dots, X^{p}) \mapsto g (X^{1}, X^{2}, \dots, X^{p}) \in SE (3)}

(1)

where g is a field of rigid transformations mapping the inertial frame onto the frames attached to each of the microstructures that constitute the medium. This definition holds for beams (p=1) and shells (p=2).

Fig. 3.

Schematic description of the general configuration of a multi-soft-body system.

2.1. Multi-soft-body configuration space

In our definition of the configuration space of the multi-soft-body system, the net motion is separated from the deformation of the soft subsystems and the configuration space of the system is defined as $SE (3) \times C \times C \times \dots \times C$ where SE(3) denotes the configuration space of the reference rigid root-body and each copy $C$ of (1) denotes a field of transformations mapping a dummy frame connecting the reference rigid body and a soft subsystem onto its microstructures. In this context, $g_{rs}$ , $g_{rb} \in SE (3)$ denote the constant transformations that map the reference frame onto the dummy frames attached to shells and beams $((O_{s}, E_{1 s}, E_{2 s}, E_{3 s}))$ and $((O_{b}, E_{1 b}, E_{2 b}, E_{3 b}))$ (Figure 4(b)), respectively, whereas $g_{s} (X_{s}, ϕ)$ and $g_{b} (X_{b}) \in SE (3)$ denote the set of maps of (1) that map the dummy frame onto the microstructure frames of the two-dimensional axisymmetric shell, parameterized by the curvilinear abscissa $X_{s}$ and the revolution angle around the symmetry axis $ϕ$ , and to the material element of the one-dimensional beam, parameterized by the curvilinear abscissa $X_{b}$ , respectively. With these definitions at hand, the connections between the rigid and soft subsystem are defined as locked joints through the equations $g_{r} g_{rb} g_{b} (0) = g (0)$ and $g_{r} g_{rs} g_{s} (0, ϕ) = g (0, ϕ)$ for beams and shells, respectively.

Fig. 4.

Schematic description of the multi-soft-body configuration and the map’s hierarchy between the frames used in the model.

It is worth noting that the configuration spaces of the constitutive bodies (soft or rigid) sharing as common structure the Lie group SE(3), their dynamic models can be encompassed in a common framework that is presented in the subsequent developments.

3. Cosserat model for soft robotics

In this section, a brief description, based on the authors previous works (Renda et al., 2014, 2015c), of the kinematics and dynamics of soft robot arms (SRAs) and soft shell mantles (SSMs) for underwater soft robotics is given.

3.1. Kinematics

In the Cosserat theory, according to Equation (1), the configuration of a micro-solid of a soft body with respect to the inertial frame at a certain time is characterized by a position vector u and a material orientation matrix R, parameterized by the material abscissas, that are $ϕ \in [0, 2 π [$ (the angle of revolution of the axisymmetric surface) and $X_{s} \in [0, L_{s}]$ (the abscissa along the meridian) for the SSM; and $X_{b} \in [0, L_{b}]$ (the abscissa along the robot arm) for the SRA, the subscripts _s and _b denote the shell and beam, respectively. Thus, the configuration space is defined as a curve $g (\cdot) : X_{b} \mapsto g (X_{b}) \in SE (3)$ or a surface $g (\cdot) : (X_{s}, ϕ) \mapsto g (X_{s}, ϕ) \in SE (3)$ , with

g = (\begin{matrix} R & u \\ 0 & 1 \end{matrix})

As described in the previous section, the map g is the composition of three transformations, $g_{r}$ , $g_{rs}$ , and $g_{s}$ for shells and correspondingly $g_{r}$ , $g_{rb}$ , and $g_{b}$ for beams, giving $g = g_{r} g_{rs} g_{s}$ or $g = g_{r} g_{rb} g_{b}$ (Figure 5), with

\begin{matrix} g_{r} = (\begin{matrix} R_{r} & u_{r} \\ 0 & 1 \end{matrix}), g_{rs} = (\begin{matrix} R_{rs} & u_{rs} \\ 0 & 1 \end{matrix}), \\ g_{rb} = (\begin{matrix} R_{rb} & u_{rb} \\ 0 & 1 \end{matrix}), g_{s} = (\begin{matrix} R_{s} & u_{s} \\ 0 & 1 \end{matrix}), g_{b} = (\begin{matrix} R_{b} & u_{b} \\ 0 & 1 \end{matrix}) \end{matrix}

Fig. 5.

Sketch of the kinematics which show the geometrical meaning of the elements g, $ξ$ and $η$ . The reference frames on the figure are those used in the model.

Furthermore, exploiting the axisymmetry of the shell, the transformation $g_{s}$ can be specified as the following combination:

\begin{matrix} g_{s} = g_{1} g_{2} g_{3} = (\begin{matrix} \exp (ẽ_{3} ϕ) & u_{s} \\ 0 & 1 \end{matrix}) (\begin{matrix} \exp (ẽ_{r} π ∕ 2) & 0 \\ 0 & 1 \end{matrix}) \\ (\begin{matrix} \exp (- ẽ_{ϕ} θ) & 0 \\ 0 & 1 \end{matrix}) \end{matrix}

where $\exp$ is the exponential in SO(3), the tilde is the isomorphism between a vector of $ℝ^{3}$ and the corresponding skew-symmetric matrix ( $\in s o (3)$ ) and $θ (X_{s})$ is the angle between $E_{3 s}$ and the shell microstructure located at any $X_{s}$ along the $ϕ$ -meridian (Figure 5(a)). In this case $u_{s} (\cdot)$ takes the form: $u_{s} (X_{s}) = {(\cos (ϕ) r (X_{s}), \sin (ϕ) r (X_{s}), z (X_{s}))}^{T}$ for which $r (.)$ and $z (.)$ are two smooth functions that define the radius and the altitude of the point $X_{s}$ on the profile (g₂ is a constant transformation chosen by the authors for the sake of convenience).

Based on these kinematics, the strain state of the beam is defined by the vector field along the curve $g (\cdot) : X_{b} \mapsto g (X_{b})$ given by $\hat{ξ} (\cdot) : X_{b} \mapsto \hat{ξ} (X_{b}) = g^{- 1} ∂g ∕ \partial X_{b} = g_{b}^{- 1} \partial g_{b} ∕ \partial X_{b} = g_{b}^{- 1} g_{b}^{'} = {\hat{ξ}}_{b} (X_{b}) \in s e (3)$ and the two vector fields on the surface $g (\cdot) : (X_{s}, ϕ) \mapsto g (X_{s}, ϕ)$ , defined by the two vector fields: ${\hat{ξ}}_{1} (\cdot) : X_{s} \mapsto {\hat{ξ}}_{1} (X_{s}) = g^{- 1} ∂g ∕ \partial X_{s} = g_{s}^{- 1} \partial g_{s} ∕ \partial X_{s} = g_{s}^{- 1} g_{s}^{'} = {\hat{ξ}}_{s 1} (X_{s})$ and ${\hat{ξ}}_{2} (\cdot) : X_{s} \mapsto {\hat{ξ}}_{2} (X_{s}) = g^{- 1} ∂g ∕ ∂ϕ = g_{s}^{- 1} \partial g_{s} ∕ ∂ϕ = {\hat{ξ}}_{s 2} (X_{s})$ (where the ${hat}^{^}$ represents the isomorphism between the twist vector space $ℝ^{6}$ and the Lie algebra $s e (3)$ ). After simple algebra based on the above transformation kinematics, the components of these fields are specified in the microstructure frames as

\begin{matrix} ξ_{b} = (\begin{matrix} k_{b} \\ q_{b} \end{matrix}), ξ_{s 1} = (\begin{matrix} k_{s 1} \\ q_{s 1} \end{matrix}) \\ = (\begin{matrix} 0 \\ 0 \\ θ^{'} \\ \cos (θ) r^{'} + \sin (θ) z^{'} \\ \cos (θ) z^{'} - \sin (θ) r^{'} \\ 0 \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ μ \\ λ \\ β \\ 0 \end{matrix}), \\ ξ_{s 2} = (\begin{matrix} k_{s 2} \\ q_{s 2} \end{matrix}) = (\begin{matrix} \sin (θ) \\ \cos (θ) \\ 0 \\ 0 \\ 0 \\ - r \end{matrix}) \in ℝ^{6} \end{matrix}

where $q_{s} (X_{s})$ , $q_{b} (X_{b})$ represents the linear strains, and $k_{s} (X_{s})$ , $k_{b} (X_{b})$ the angular strain. It worth noting that the two fields ${\hat{ξ}}_{s 1}$ and ${\hat{ξ}}_{s 2}$ depend only on $X_{s}$ due to the shell’s axisymmetry. The time evolution of the configuration curve $g (\cdot) : X_{b} \mapsto g (X_{b})$ andsurface $g (\cdot) : (X_{s}, ϕ) \mapsto g (X_{s}, ϕ)$ is represented by the twist vector fields $\hat{η} (\cdot) : X_{b} \mapsto \hat{η} (X_{b})$ and $\hat{η} (\cdot) : X_{s} \mapsto \hat{η} (X_{s}) \in s e (3)$ defined by $\hat{η} (X_{b}) = g {(X_{b})}^{- 1} ∂g (X_{b}) ∕ ∂t = g^{- 1} ġ$ and $\hat{η} (X_{s}) = g {(X_{s}, ϕ)}^{- 1} ∂g (X_{s}, ϕ) ∕ ∂t = g^{- 1} ġ$ , respectively. It can be shown that these two vector fields can be written as $\hat{η} (X_{b}) = {\hat{η}}_{b} +^{b} {\hat{η}}_{r}$ and $\hat{η} (X_{s}) = {\hat{η}}_{s} +^{s} {\hat{η}}_{r}$ , where we have defined ${\hat{η}}_{b} = g_{b}^{- 1} ġ_{b}$ , ${\hat{η}}_{s} = g_{s}^{- 1} ġ_{s}$ and ${\hat{η}}_{r} = g_{r}^{- 1} ġ_{r}$ . Note that this decomposition is consistent with our definition of the multi-soft-body configuration space that separates the fields of inertial transformations on each of the subsystems into a rigid component $g_{r}$ and a shape component ( $g_{b} and g_{s}$ ). These twist can be detailed in terms of their components in the microstructure frames as

\begin{matrix} η_{b} = (\begin{matrix} w_{b} \\ v_{b} \end{matrix}), η_{s} = (\begin{matrix} w_{s} \\ v_{s} \end{matrix}) \\ = (\begin{matrix} 0 \\ 0 \\ \dot{θ} \\ \cos (θ) ṙ + \sin (θ) \dot{z} \\ \cos (θ) \dot{z} - \sin (θ) ṙ \\ 0 \end{matrix}) \\ = (\begin{matrix} 0 \\ 0 \\ Ω \\ V_{x} \\ V_{y} \\ 0 \end{matrix}), η_{r} = (\begin{matrix} w_{r} \\ v_{r} \end{matrix}) \in ℝ^{6} \end{matrix}

where $v_{b} (X_{b})$ , $v_{s} (X_{s})$ , $v_{r}$ and $w_{b} (X_{b})$ , $w_{s} (X_{s})$ , $w_{r}$ are the linear and angular velocity at a given instant, respectively. The back-superscripts ^s and ^b represent an element of the Lie algebra $s e (3)$ of $g_{r}$ expressed in the microstructure frame of the shell and of the beam, respectively, while the absence of any back-superscripts indicate that the Lie algebra element is expressed in the frame to which it is related. To exploit the axisymmetry of the SSM, we limit $g_{r}$ of $g (X_{s}, ϕ) = g_{r} g_{rs} g_{s}$ to its translational component about the axis of symmetry $E_{3 s}$ , leading to $η_{r} = {(0, 0, 0, (0, 0, V_{E_{3 s}}) R_{rs}^{T})}^{T}$ (i.e. the velocity of the multi-soft-body system in the $E_{3 s}$ direction) in that case. Remembering this distinction, the twist versions ${{}^{s}η}_{r}$ and ${{}^{b}η}_{r}$ can be calculated as ${Ad}_{g_{s}^{- 1} g_{rs}^{- 1}} η_{r}$ and ${Ad}_{g_{b}^{- 1} g_{rb}^{- 1}} η_{r}$ , where $Ad$ is the Adjoint map, defined as

{Ad}_{g} = (\begin{matrix} R & 0 \\ \tilde{u} R & R \end{matrix})

3.1.1. Strain measures

There are different ways to measure the strain of a continuous media, we choose the most commonly used in the specialized literature for the beam (Simo, 1985) and shell (Simo and Fox, 1989), respectively.

For the SRA, the strains are defined as the difference between the X -rate of g in the deformed configuration $ξ$ , and on the reference configuration $ξ^{⋆}$ (denoted with a ^⋆). In particular, the components of $k - k^{⋆}$ measure the torsion and the bending state in the two directions. Similarly, the components of $q - q^{⋆}$ represent the longitudinal strain (extension, compression) and the two shear strains.

For the SSM, in accordance with Simo and Fox (1989) as described in Renda et al. (2015c), the strain tensor field which describes the membrane strain state in the mid-surface is $e (\cdot) : X_{s} \mapsto e (X_{s}) = 1 ∕ 2 (h - h^{⋆}) \in ℝ^{2} \otimes ℝ^{2}$ where h is the first fundamental form of the Reissner shell equal to $h = diag (λ^{2} + β^{2}, r^{2})$ . Thus, we have $e = (1 ∕ 2) diag (λ^{2} + β^{2} - 1, r^{2} - r^{⋆ 2})$ , in which we have defined $h_{11}^{⋆} = 1$ . For what concerns the shear strain state, we have $s (X_{s}) = β - β^{⋆}$ . Finally, the flexural strain state is parametrized by the tensor field $d (\cdot) : X_{s} \mapsto d (X_{s}) = k - k^{⋆} \in ℝ^{2} \otimes ℝ^{2}$ , where k is the second fundamental form equal to $k = diag (- μλ, - r \sin (θ))$ . Thus, we have $d = diag (μ^{⋆} - μλ, r^{⋆} \sin (θ^{⋆}) - r \sin (θ))$ . Furthermore, it is natural to consider that there is no transverse shearing in the reference resting configuration, i.e. $β^{⋆} = 0$ .

3.2. Compatibility equations

We have seen above that $g {(X_{b})}^{'} = g {\hat{ξ}}_{b}$ and $g {(X_{s}, ϕ)}^{'} = g {\hat{ξ}}_{s 1}$ . By taking the derivative of these equations with respect to time and recalling that $ġ (X_{b}) = g ({\hat{η}}_{b} +^{b} {\hat{η}}_{r})$ and $ġ (X_{s}, ϕ) = g ({\hat{η}}_{s} +^{s} {\hat{η}}_{r})$ , we obtain ${\dot{ξ}}_{b} = {(η_{b} +^{b} η_{r})}^{'} + {ad}_{ξ_{b}} (η_{b} +^{b} η_{r})$ and ${\dot{ξ}}_{s 1} = {(η_{s} +^{s} η_{r})}^{'} + {ad}_{ξ_{s 1}} (η_{s} +^{s} η_{r})$ , where $ad$ is the adjoint map defined as

{ad}_{ξ} = (\begin{matrix} \tilde{k} & 0 \\ \tilde{q} & \tilde{k} \end{matrix})

Then, we can simplify these equations by noting that ${({}^{b}η r)}^{'} = {({Ad}_{g_{b}^{- 1} g_{rb}^{- 1}} η_{r})}^{'} = {({Ad}_{g_{b}^{- 1} g_{rb}^{- 1}})}^{'} η_{r} = - {Ad}_{g_{b}^{- 1} g_{rb}^{- 1}} ({ad}_{{{}^{r}ξ}_{b}} η_{r}) = - {ad}_{ξ_{b}} {{}^{b}η}_{r}$ (the same holds for ${{}^{s}η}_{r}$ ), leading to the following compatibility equations between velocity and deformation variables:

\begin{matrix} {\dot{ξ}}_{b} = η_{b}^{'} + {ad}_{ξ_{b}} η_{b} \end{matrix}

(2)

\begin{matrix} {\dot{ξ}}_{s 1} = η_{s}^{'} + {ad}_{ξ_{s 1}} η_{s} \end{matrix}

(3)

which remarkably depend only on the “shape” component of the multi-soft-body configuration space.

3.3. Dynamics

The partial differential equations (p.d.e.s) describing the evolution of a Cosserat rod and shell (not necessarily axisymmetric) have been derived by Reissner (1990) and exploited in Simo (1985) and (Simo and Fox, 1989), respectively, for nonlinear finite element analysis. More recently, they have been used in the context of continuous and soft robotics in Boyer et al. (2006) Candelier et al. (2013), Renda et al. (2014), and Renda et al. (2015c). Boyer and Primault (2005); Boyer and Renda (2016) showed that these beam (respectively, shell) p.d.e.s, together with their boundary conditions, can be derived directly from an extension of a variational calculus on Lie groups historically introduced by Poincaré (1901):

\begin{matrix} M_{b} [{\dot{η}}_{b} + ({}^{b}{\dot{η}}_{r}] = {(F_{bi} - F_{ba})}^{'} - {ad}_{ξ_{b}}^{T} (F_{bi} - F_{ba}) + {\bar{F}}_{be} \\ + {ad}_{(η_{b} +^{b} η_{r})}^{T} M_{b} (η_{b} +^{b} η_{r}) \\ η_{b} (0) = 0, (F_{bi} - F_{ba}) (L_{b}) = F_{be} (L_{b}) \end{matrix}

(4)

\begin{matrix} M_{s} [{\dot{η}}_{s} + ({}^{s}{\dot{η}}_{r}] = \frac{{[r^{⋆} (F_{si}^{1} - F_{sa}^{1})]}^{'}}{r^{⋆}} - {ad}_{ξ_{sα}}^{T} (F_{si}^{α} - F_{sa}^{α}) \\ + {\bar{F}}_{se} + {ad}_{(η_{s} +^{s} η_{r})}^{T} M_{s} (η_{s} +^{s} η_{r}) \\ η_{s} (0) = 0, (F_{si}^{1} - F_{sa}^{1}) (L_{s}) = F_{se} (L_{s}) \end{matrix}

(5)

where $F_{bi} (X_{b})$ and $F_{si}^{α} (X_{s})$ are the wrenches of internal elastic forces in the surface directions given by $q_{b}$ and $q_{sα}$ ( $α$ running over ${1, 2}$ ), $F_{ba} (X_{b}, t) and F_{sa}^{α} (X_{s}, t)$ are the internal actuation loads $, {\bar{F}}_{be} (X_{b}) and {\bar{F}}_{se} (X_{s})$ are the external wrench of distributed applied forces $, M_{b} (X_{b})$ and $M_{s} (X_{s})$ are the screw inertia matrix. For the repeated $α$ , the Einstein convention has to be used as in the rest of the paper.

Let us specify the angular and linear components of the internal and external wrenches (for the axisymmetric shell refer to Renda et al. (2015c) and Antman (2006)):

\begin{matrix} F_{bi} = (\begin{matrix} M_{bi} \\ N_{bi} \end{matrix}), F_{si}^{1} = (\begin{matrix} M_{si}^{1} \\ N_{si}^{1} \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ M_{X} \\ N_{X} \\ H \\ 0 \end{matrix}), \\ F_{si}^{2} = (\begin{matrix} M_{si}^{2} \\ N_{si}^{2} \end{matrix}) = (\begin{matrix} M_{ϕ_{x}} \\ M_{ϕ_{y}} \\ 0 \\ 0 \\ 0 \\ - N_{ϕ} \end{matrix}) \in ℝ^{6} \end{matrix}

\begin{matrix} F_{ba} = (\begin{matrix} M_{ba} \\ N_{ba} \end{matrix}), F_{sa}^{1} = (\begin{matrix} M_{sa}^{1} \\ N_{sa}^{1} \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ L_{za} \\ F_{xa} \\ F_{ya} \\ 0 \end{matrix}), \\ F_{sa}^{2} = (\begin{matrix} M_{sa}^{2} \\ N_{sa}^{2} \end{matrix}) = (\begin{matrix} L_{xa} \\ L_{ya} \\ 0 \\ 0 \\ 0 \\ - F_{za} \end{matrix}) \in ℝ^{6} \end{matrix}

\begin{matrix} {\bar{F}}_{be} = (\begin{matrix} m_{b} \\ n_{b} \end{matrix}), {\bar{F}}_{se} = (\begin{matrix} m_{s} \\ n_{s} \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ l \\ f_{x} \\ f_{y} \\ 0 \end{matrix}) \in ℝ^{6} \end{matrix}

where $N_{si}^{α} (X_{s})$ , $N_{bi} (X_{b})$ and $M_{si}^{α} (X_{s})$ , $M_{bi} (X_{b})$ are the internal elastic force and torque vectors, respectively, $N_{sa}^{α} (X_{s}, t)$ , $N_{ba} (X_{b}, t)$ and $M_{sa}^{α} (X_{s}, t)$ , $M_{ba} (X_{b}, t)$ are the internal actuation force and torque, whereas $n_{s} (X_{s})$ , $n_{b} (X_{b})$ and $m_{s} (X_{s})$ , $m_{b} (X_{b})$ are the external force and torque for unit of $X_{s}$ , $X_{b}$ , $r^{⋆} ϕ$ . The screw inertia matrices are equal to $M_{b} = diag (I_{b}, J_{b}, J_{b}, A, A, A) ρ_{b}$ and $M_{s} = diag (J_{s}, I_{s}, J_{s}, 2 h_{s}, 2 h_{s}, 2 h_{s}) ρ_{s} \in ℝ^{6} \otimes ℝ^{6}$ , where $ρ_{b}$ and $ρ_{s}$ are the body densities, $A (X_{b})$ is the section area equal to $A = π h_{b}^{2}$ ( $h_{b} (X_{b})$ being the cross-section radius) $, h_{s}$ is the half of the shell thickness and $J_{s} (X_{s})$ $, J_{b} (X_{b})$ $, I_{s} (X_{s})$ $, I_{b} (X_{b})$ are the second moment of inertia of the micro-solid equal to $J_{b} = π h_{b}^{4} ∕ 4$ $, J_{s} = h_{s}^{2} ∕ 3$ $, I_{b} = π h_{b}^{4} ∕ 2$ , and $I_{s} ~ 0$ .

As for the compatibility equations, we have $({}^{b}{\dot{η}}_{r}) = {Ad}_{g_{b}^{- 1} g_{rb}^{- 1}} {\dot{η}}_{r} + \dot{({Ad}_{g_{b}^{- 1} g_{rb}^{- 1}})} η_{r} =^{b} {\dot{η}}_{r} - {ad}_{η_{b}} {{}^{b}η}_{r}$ (the same holds for ${{}^{s}η}_{r}$ ), hence the left-hand side of Equations (4) and (5) become $M_{b} ({\dot{η}}_{b} +^{b} {\dot{η}}_{r} - {ad}_{η_{b}} {{}^{b}η}_{r})$ and $M_{s} ({\dot{η}}_{s} +^{s} {\dot{η}}_{r} - {ad}_{η s} {{}^{s}η}_{r})$ .

The imposed internal actuation wrenches ( $F_{sa}^{α} (X_{s}, t)$ $, F_{ba} (X_{b}, t)$ ) represents the input of the model and depends directly on time. It can be thought of as the action of the muscle fiber of the body in the case of a living organism or the result of embedded cable-driven actuation as in Renda et al. (2014).

3.4. Constitutive equations

A linear viscoelastic constitutive equation, based on the Kelvin Voigt model, is chosen. Following Linn et al. (2013) for the SRA and Simo and Fox (1989) for the SSM we respectively obtain

F_{bi} = Σ (ξ_{b} - ξ_{b}^{⋆}) + ϒ ({\dot{ξ}}_{b})

(6)

\begin{matrix} N_{X} = \frac{2 E h_{s}}{1 - ν^{2}} [λ (e_{11} + \frac{ν}{r^{⋆ 2}} e_{22}) - J_{s} μ (d_{11} + \frac{ν}{r^{⋆ 2}} d_{22})] \\ + \frac{6 υ h_{s}}{1 - ν^{2}} [λ (ė_{11} + \frac{ν}{r^{⋆ 2}} ė_{22}) - J_{s} μ (ḋ_{11} + \frac{ν}{r^{⋆ 2}} ḋ_{22})] \\ N_{ϕ} = \frac{2 E h_{s}}{1 - ν^{2}} [r (\frac{e_{22}}{r^{⋆ 4}} + \frac{ν}{r^{⋆ 2}} e_{11}) - J_{s} \sin (θ) (\frac{d_{22}}{r^{⋆ 4}} + \frac{ν}{r^{⋆ 2}} d_{11})] \\ + \frac{6 υ h_{s}}{1 - ν^{2}} [r (\frac{ė_{22}}{r^{⋆ 4}} + \frac{ν}{r^{⋆ 2}} ė_{11}) - J_{s} \sin (θ) (\frac{ḋ_{22}}{r^{⋆ 4}} + \frac{ν}{r^{⋆ 2}} ḋ_{11})] \\ H = 2 h_{s} β [G + \frac{E}{1 - ν^{2}} (e_{11} + \frac{ν}{r^{⋆ 2}} e_{22})] \\ + 2 h_{s} \dot{β} [υ + \frac{3 υ}{1 - ν^{2}} (ė_{11} + \frac{ν}{r^{⋆ 2}} ė_{22})] \\ M_{X} = - \frac{2 E h_{s} J_{s}}{1 - ν^{2}} λ (d_{11} + \frac{ν}{r^{⋆ 2}} d_{22}) - \frac{6 υ h_{s} J_{s}}{1 - ν^{2}} λ (ḋ_{11} + \frac{ν}{r^{⋆ 2}} ḋ_{22}) \\ M_{ϕ_{x}} = - \frac{2 E h_{s} J_{s}}{1 - ν^{2}} r (\frac{d_{22}}{r^{⋆ 4}} + \frac{ν}{r^{⋆ 2}} d_{11}) - \frac{6 υ h_{s} J_{s}}{1 - ν^{2}} r (\frac{ḋ_{22}}{r^{⋆ 4}} + \frac{ν}{r^{⋆ 2}} ḋ_{11}) \\ M_{ϕ_{y}} = 0 \end{matrix}

(7)

where $Σ (X)$ and $ϒ (X) \in ℝ^{6} \otimes ℝ^{6}$ are the screw stif- fness matrix and the screw viscosity matrix, equal to $Σ = diag (G I_{b}, E J_{b}, E J_{b}, EA, GA, GA)$ , $ϒ = diag (I_{b}, 3 J_{b}, 3 J_{b}, 3 A, A, A) υ$ , E being the Young modulus (different between the shell and the beam), G is the shear modulus (equal to $G = E ∕ 2 (1 + ν)$ for an isotropic material with Poisson ratio $ν$ ), and $υ$ is the shear viscosity modulus.

3.5. External loads

The external loads taken into account are those exerted by the fluid (i.e. drag, added mass, buoyancy, and thrust) in addition to the gravity load. Mathematically, we have

n_{b} = s_{b} + b_{b} + d_{b} + a_{b}

(8)

n_{s} = d_{s} + a_{s} + t_{s}

(9)

where $s_{b} (X_{b})$ is the gravity, $b_{b} (X_{b})$ is the buoyancy, $t_{s}$ is the thrust load, $d_{s} (X_{s})$ , $d_{b} (X_{b})$ are the drag, and $a_{s} (X_{s})$ , $a_{b} (X_{b})$ are the added mass, whereas $m_{b} = m_{s} = {(0, 0, 0)}^{T}$ .

An exhaustive derivation and interpretation of the fluid force model for the SSM has been presented in Renda et al. (2015b), based on the usual model of net external forces exerted on a rigid rocket, uniformly distributed overthe mantle. In this formulation, certain terms which participate in the definition of the total propulsive thrust are neglected. These concern, in particular, the internal pressure contribution associated with the dynamics of the cavity collapse, referred to by Krieg and Mohseni (2015) as total jetting force (see also Anderson and DeMont (2000)), as well as the positive feedback from the added mass variation of the collapsible shell (Giorgio-Serchi and Weymouth, 2016a,b). Here only the final equation is reported, whereas the SSM gravity and buoyancy loads are accounted for together with the rigid root-body so that the shell symmetry is not broken. For the SRA, the fluid force models have been originally derived in Boyer et al. (2006) and then introduced in a soft robotics context in Renda et al. (2014).

Gravity and buoyancy are simply the product between the mass per unit of $X_{b}$ of the robot arm of the water and the gravity acceleration, respectively,

s_{b} + b_{b} = (ρ_{b} - ρ_{w}) A R^{T} G

(10)

where G is the gravity acceleration vector, equal to $G = {(0, 0, - 9.81)}^{T}$ .

The drag load vector is proportional to the square of the velocity vector and is directed in the opposite direction. The amplitude of the drag load is also determined by the geometry of section $X_{s}$ , $X_{b}$ and by hydrodynamics phenomena expressed by empirical coefficients. The equations are

d_{b} = - ρ_{w} D_{b} (v_{b} +^{b} v_{r}) ∥ v_{b} +^{b} v_{r} ∥

(11)

d_{s} = - R_{s}^{T} {(0, 0, \frac{ρ_{w} C_{d} A_{ref} V_{E_{3 s}} | V_{E_{3 s}} |}{2 A_{m}})}^{T}

(12)

where $D_{b} (X_{b}) \in ℝ^{3} \otimes ℝ^{3}$ is equal to $D_{b} = diag (\frac{1}{2} π C_{bx}, C_{by}, C_{bz}) h_{b}$ for circular cross-sections of radius $h_{b}$ ; $C_{bx}$ , $C_{by}$ , $C_{bz}$ being the empirical hydrodynamic coefficients; $C_{d}$ is the net drag coefficient of the SSM, $A_{ref}$ is the reference area equal to $π {(max (r (X)))}^{2}$ and $A_{m}$ is the total surface of the SSM equal to $A_{m} = \int_{0}^{L_{s}} 2 π r^{⋆} d X_{s}$ .

The added mass load vector is proportional to the acceleration vector and is directed in the opposite direction. The amplitude is also determined by the geometry of section $X_{s}$ , $X_{b}$ and by hydrodynamics phenomena expressed in part by correction coefficients. The equations are

\begin{matrix} a_{b} = - \frac{d [ρ_{w} F_{b} (v_{b} +^{b} v_{r})]}{dt} \\ = - ρ_{w} F_{b} [{\dot{v}}_{b} + ({}^{b}{\dot{v}}_{r})] - (w_{b} +^{b} w_{r}) \times ρ_{w} F_{b} (v_{b} +^{b} v_{r}) \end{matrix}

(13)

\begin{matrix} a_{s} = - B_{s} ρ_{s} 2 h_{s} [{\dot{v}}_{s} + ({}^{s}{\dot{v}}_{r})] \\ = - B_{s} ρ_{s} 2 h_{s} (\begin{matrix} {\dot{V}}_{x} + \sin (θ) {\dot{V}}_{E_{3 s}} + Ω \cos (θ) V_{E_{3 s}} \\ {\dot{V}}_{y} + \cos (θ) {\dot{V}}_{E_{3 s}} - Ω \sin (θ) V_{E_{3 s}} \\ 0 \end{matrix}) \end{matrix}

(14)

where $B_{s}$ is the net added mass coefficient and $F_{b} (X_{b}) \in ℝ^{3} \otimes ℝ^{3}$ is a tensor that incorporates the geometric and hydrodynamics factors, equal to $F_{b} = diag (0, A B_{b}, A B_{b})$ , $B_{b}$ being the hydrodynamic correction coefficient.

The thrust load is

t_{s} = - R_{s}^{T} {(0, 0, \frac{ρ_{w} C_{f} \dot{U} | \dot{U} |}{A_{n} A_{m}})}^{T}

(15)

where $A_{n}$ is the nozzle area equal to $A_{n} = A_{no}$ for the outflow and equal to $A_{n} = A_{ni}$ for the inflow, $C_{f}$ defines a flow loss coefficient at the nozzle entrance, which is taken to vary between 0.6 and 1 (Johnson et al., 1972), and U is the mantle inner volume equal to $U = \int_{0}^{L_{s}} π r^{2} \sqrt{λ^{2} + β^{2}} d X_{s}$ .

4. Rigid root-body model

We here seek the dynamic model of the system net motions controlled by the shape deformations of the soft subsystems (Boyer and Porez, 2015). To do that, the kinematics and dynamics of the rigid root-body that connects the soft bodies in a star system are presented, together with the reaction wrenches due the soft appendages. This leaves to the following equation:

M_{r} {\dot{η}}_{r} = F_{ria} + F_{re} + {ad}_{η_{r}}^{T} M_{r} η_{r}

(16)

where $M_{r}$ is the screw inertia matrix and $F_{ria}$ , $F_{re}$ are the internal elastic/actuation (due to the soft appendages) and external wrench loads acting on the root body, respectively. This dynamic equation has to be supplemented with the following kinematic model:

\dot{g_{r}} = g_{r} {\hat{η}}_{r}

(17)

which plays the role of a reconstruction equation. For our soft unmanned underwater vehicle (SUUV), the rigid root-body is composed by four rectangular parallelepiped bars in a pyramidal configuration (Figure 6), hence, the inertia matrix takes the form

M_{r} = \sum_{β = 1}^{4} {Ad}_{g_{rrβ}^{- 1}}^{T} M_{rr} {Ad}_{g_{rrβ}^{- 1}}

where $M_{rr} = diag (m_{rr} (L_{rr}^{2} + h_{rr}^{2}) ∕ 12, m_{rr} (L_{rr}^{2} + h_{rr}^{2}) ∕ 12$ , $m_{rr} (2 h_{rr}^{2}) ∕ 12, m_{rr}, m_{rr}, m_{rr})$ is the principal inertia matrix of a single parallelepiped ( $m_{rr}$ $, L_{rr}$ , and $h_{rr}$ being its mass, length, and thickness, respectively) and $g_{rrβ}$ ( $β$ running from 1 to 4) is the constant rigid transformation between the frame $(O, E_{1}, E_{2}, E_{3})$ and the principal inertia frame of the bar $β$ (centered in the middle of the bar and aligned with the long axis).

Fig. 6.

illustrative scheme of the soft unmanned underwater vehicle (SUUV) kinematics, where $(e_{1}, e_{2}, e_{3})$ is the Euclidean fixed frame, $(E_{1}, E_{2}, E_{3})$ represent the root-body net-motion and is the soft bodies reference frame, and $(x, y, z)$ is the microstructure frame.

4.1. Root-body external loads

The external loads comes from the interaction of the rigid root-body with the environment. In our case, they are equal to

F_{re} = F_{r_{s}} + F_{r_{b}} + (\begin{matrix} 0 \\ d_{r} + a_{r} \end{matrix})

(18)

where $F_{r_{s}}$ is the gravity, $F_{r_{b}}$ is the buoyancy, $d_{r}$ is the drag load, and $a_{r}$ is the added mass.

Similarly to Section 3.5, for the external loads we have

F_{r_{s}} + F_{r_{b}} = (1 - ρ_{w} ∕ ρ_{r}) M_{r} {Ad}_{g_{r}}^{- 1} {(0, 0, 0, G^{T})}^{T}

(19)

d_{r} = - \frac{1}{2} ρ_{w} C_{r} v_{r} ∥ v_{r} ∥

(20)

a_{r} = - \frac{d (m_{r} B_{r} v_{r})}{dt} = - m_{r} B_{r} {\dot{v}}_{r} - w_{r} \times m_{r} B_{r} v_{r}

(21)

where $ρ_{r}$ and $m_{r}$ are the root-body density and total mass, respectively, whereas $C_{r}$ ( $[m^{2}]$ ) and $B_{r}$ are the empirical hydrodynamic coefficients.

4.2. Root-body reaction loads

For what concern the reaction loads due to the attached soft appendages ( $F_{ria}$ ), it has to be equal and opposite to the load that support the motion of the soft body through the junction (Featherstone, 2014). Thus, as for any kinematic tree, we have

\begin{matrix} F_{ria} = \sum_{β = 1}^{4} \int_{0}^{L_{b}} {Ad}_{g_{bβ}^{- 1} g_{rbβ}^{- 1}}^{T} \\ ({\bar{F}}_{beβ} - M_{b} [{\dot{η}}_{bβ} + ({}^{b \dot{β}}{\dot{η}}_{r})] + {ad}_{(η_{bβ} +^{bβ} η_{r})}^{T} M_{b} (η_{bβ} +^{bβ} η_{r})) \\ d X_{b} + {Ad}_{g_{bβ}^{- 1} g_{rbβ}^{- 1} (L_{b})}^{T} F_{beβ} (L_{b}) \\ + \int_{0}^{L_{s}} \int_{0}^{2 π} {Ad}_{g_{s}^{- 1} g_{rs}^{- 1}}^{T} \\ ({\bar{F}}_{se} - M_{s} [{\dot{η}}_{s} + ({}^{s}{\dot{η}}_{r})] + {ad}_{(η_{s} +^{s} η_{r})}^{T} M_{s} (η_{s} +^{s} η_{r})) \\ r^{⋆} dϕ d X_{s} + \int_{0}^{2 π} {Ad}_{g_{s}^{- 1} g_{rs}^{- 1} (L_{s}, ϕ)}^{T} F_{se} (L_{s}, ϕ) r^{⋆} (L_{s}) dϕ \end{matrix}

where $F_{be} (L_{b})$ and $F_{se} (L_{s}, ϕ)$ are the beam and shell boundary external loads.

Making use of the dynamic equations (4) and (5) (and their boundary condition), the right-hand side of the equation above can be derived by integrating the internal and actuation loads of the soft bodies leading to

\begin{matrix} F_{ria} = \sum_{β = 1}^{4} {Ad}_{g_{bβ}^{- 1} g_{rbβ}^{- 1} (0)}^{T} (F_{biβ} (0) - F_{baβ} (0)) \\ + \int_{0}^{2 π} {Ad}_{g_{s}^{- 1} g_{rs}^{- 1} (0, ϕ)}^{T} (F_{si}^{1} (0, ϕ) - F_{sa}^{1} (0, ϕ)) r^{⋆} dϕ \\ + \int_{0}^{L_{s}} \int_{0}^{2 π} (1 - ρ_{w} ∕ ρ_{s}) {Ad}_{g_{s}^{- 1} g_{rs}^{- 1}}^{T} M_{s} \\ {Ad}_{g}^{- 1} {(0, 0, 0, G^{T})}^{T} r^{⋆} dϕ d X_{s} \end{matrix}

(22)

where the neglected gravity and buoyancy loads of the SSM have been recovered.

5. Multi-soft-body dynamic model

In Section 3, the kinematics, compatibility, and dynamic equations for two type of soft bodies, i.e. the SRA and the SSM, have been given, whereas in Section 4, the rigid root-body kinematics and dynamics governed by the reaction wrenches of the soft appendage have been developed. Finally, here we are able to present the multi-soft-body system model and outline the solution algorithm that leads to the complete motion of the underwater soft vehicle.

5.1. Star system dynamic model

The final system of equations is composed by the o.d.e.s of the root-body and the second-order p.d.e.s of the soft bodies. The system of o.d.e.s for the root-body is composed by the kinematic equation (17) and the dynamic equation (16), endowed with the reaction load (22) and external loads (19), (20), and (21). The system of p.d.e.s for the SSM and the SRA is composed by the kinematics equation presented in Section 3.1, the compatibility equations (3) and (2), and the dynamic equations (5) and (4), in turn complemented with the internal elastic stresses (7), (6) and the external loads (10), (11), (12), (13), (14), and (15). Finally, in the state form $ẋ = f (x, x^{'}, x^{″}, t)$ , the star system dynamic model is

\begin{matrix} \begin{matrix} ġ_{r} = g_{r} {\hat{η}}_{r} \\ ġ_{b} = g_{b} {\hat{η}}_{b} \\ ġ_{s} = g_{s} {\hat{η}}_{s} \\ {\dot{ξ}}_{b} = η_{b}^{'} + {ad}_{ξ_{b}} η_{b} \\ {\dot{ξ}}_{s 1} = η_{s}^{'} + {ad}_{ξ_{s 1}} η_{s} \\ {\dot{η}}_{r} = M_{r}^{- 1} (F_{ria} + F_{re} + {ad}_{η_{r}}^{T} M_{r} η_{r}) \\ {\dot{η}}_{b} = M_{b}^{- 1} [{(F_{bi} - F_{ba})}^{'} - {ad}_{ξ_{b}}^{T} (F_{bi} - F_{ba}) \\ + {\bar{F}}_{be} + {ad}_{(η_{b} +^{b} η_{r})}^{T} M_{b} (η_{b} + {{}^{b}η}_{r})] - ({}^{b}{\dot{η}}_{r}) \\ {\dot{η}}_{s} = M_{s}^{- 1} [\frac{{[r^{⋆} (F_{si}^{1} - F_{sa}^{1})]}^{'}}{r^{⋆}} - {ad}_{ξ_{sα}}^{T} (F_{si}^{α} - F_{sa}^{α}) \\ + {\bar{F}}_{se} + {ad}_{(η_{s} +^{s} η_{r})}^{T} M_{s} (η_{s} + {{}^{s}η}_{r})] - {}^{s}{\dot{η}}_{r} \end{matrix} \end{matrix}

(23)

In Figure7 a diagram of the time integration loop is shown. The input of the model, directly function of time, are the actuation loads $F_{ba} (X_{b}, t)$ and $F_{sa} (X_{s}, t)$ of the SRAs and SSM, respectively. They enter in the dynamic equation of the soft appendages ((4), (5)) as well as in the expressions of the reaction wrenches acting on the rigid root-body (22). At each time step, given the current status (composed by $g_{r}$ , $g_{s}$ , $g_{b}$ , $ξ_{s}$ , $ξ_{b}$ , $η_{r}$ , $η_{s}$ , and $η_{b}$ ) and the actuation input, the kinematics equations, compatibility equations, and dynamics equations are computed, which returns the time derivative of the status vector. It is worth noting here that the dynamic of the soft bodies (5) and (4) take advantage of the calculation of the root-body dynamic (16) that returns ${\dot{η}}_{r}$ to calculate the relative acceleration ${\dot{η}}_{s}$ and ${\dot{η}}_{b}$ , as indicate in Figure 7. Furthermore, the system (23) is infinite dimensional since all its components are some functions of the profile abscissas $X_{s}$ and $X_{b}$ . As a result, it has to be first space-discretized on a grid of nodes before being time integrated using explicit or implicit time integrators starting from the initial state. In this grid, all the space derivatives appearing in the p.d.e.s system can be approximated by finite difference schemes, with the following boundary conditions given in (5) and (4).

Fig. 7.

Diagram of the time integration loop algorithm. On the left, the current status of the state variables is plug into the multi-soft-body model where $g_{r}$ transforms the inertial frame into the reference frame of the rigid body, $g_{s}$ transforms the dummy frame of the shell into the microstructure frames of the shell, $g_{b}$ does the same job for the beam, $ξ_{s}$ and $ξ_{b}$ represent the strain fields for the shell and the beam, respectively, $η_{r}$ is the velocity screw vector of the rigid body, and, finally, $η_{s}$ , $η_{b}$ are the velocity fields of the shell and beam, respectively. On the right, the time evolution of the state variables is the output of the multi-soft-body model.

The algorithm has been implemented in MATLAB®. The numerical scheme used is a decentralized (for the SRA) and centralized (for the SSM) space differentiation finite difference method, based on a fourth-order Runge–Kutta time integration with variable time step (by means of the MATLAB® ode45 function) (Renda et al., 2014, 2015b). A spatial distribution of one material point for every 5 mm (for SRA) and 1 mm (for the SSM) was adopted. Our implementation of the MATLAB® code is available on GitHub (github.com/federicorenda/Unified-Multi-soft-body-Dynamics) under the permissive BSD 3-clause license.

6. Experimental results

Having defined the modeling framework valid for an arbitrary multi-soft-body system, it is necessary to assess the degree of accuracy of such formulation with respect to experimental data. With this purpose in mind, the data collected from experimental trials of the vehicle depicted in Figures 1 and 2 are employed. Upon assessment of the degree of accuracy of the formulation presented, the model can be employed for innovative design exploration, as later demonstrated in Section 7, design optimization and control purposes.

Tests on the vehicle were performed in a controlled environment to provide the basis for model verification. The tests entailed the robot moving along a straight track inside a working space with the shape of a rectangular box delimited by eight markers, see Figure 8. By making use of two cameras and three additional markers (see element [2] in Figure 8) fitted on the central part of the robot, a three-dimensional reconstruction of the body position and orientations is derived via direct linear transformation.

Fig. 8.

Two screenshots from the digital camera recordings taken during robot testing at (a) the start and (b) end point of the experiment. The robot is portrayed with an arm twirled around a screwdriver [1]; the variable buoyancy module [3] and the LED markers for 3D tracking [2] are visible.

The experiments were performed in a 1150 mm long, 590 mm wide, and 500 mm deep tank filled with fresh water. The tests consist in recording the displacement of the body as it propels itself from one end of the tank to the other. Recordings are performed with a digital camera at 25 fps and later processed with an image tracking software and a Savitzky–Golay low-pass filter eventually yielding the displacement and velocity in the surge direction.

The robot is allowed to travel along a straight line inside the tank by letting the motor revolve at a quasi-constant angular velocity. The overall body of the vehicle is slightly negatively buoyant, thus requiring the use of an inflatable buoyancy module (see element [3] in Figure 8) fitted to the dorsal part of the robot to achieve a condition of consistent neutral buoyancy. Tests were repeated at motor angular velocity ranging from 5 to 15 rad/s, i.e. from 0.8 to 2.4 pulsations per second (pps). The vehicle is allowed to translate along the surge direction only, whereas the motor is supplied with a constant voltage. The result of the pulsed-jet mode of propulsion generates a quasi-sinusoidal velocity signal.

In addition, from the recordings of electric current, the tension of the cable bundle at the crank is derived; this is taken as the input force in the elastodynamics model. The electric current supplied to the motor throughout the pulsation cycle was measured in the case of a forcedly stationary prototype. This grants that the dynamic effect of the external flow during vehicle displacement could not affect the load acted upon the shell during actuation. A short section of the recordings of the cable tension pattern during actuation of the prototype is shown in Figure 9.

Fig. 9.

Cables tension of a forcedly not moving prototype with a motor frequency of 2 pps (pulsations per second).

6.1. Underwater locomotion comparison

The separated model for the SRA and the SSM have been experimentally validated in Renda et al. (2014) and Renda et al. (2015b), respectively (see also Renda et al. (2015a) for the SSM in steady-state condition). Taking advantage of this fact, the same mechanical and dynamical parameters obtained there are applied here to the new geometries. The geometry of the SSM has been represented with an half sphere with a radius $T_{s} = 31$ mm truncated with a rigid cup on the section of angle $π ∕ 4$ and merged with a cylinder of length 94 mm, for a total profile length of $L_{s} = 118$ mm. The shell half thickness is homogeneously equal to $h_{s} = 1$ mm. The SRA has a conical shape with a radius linearly decreasing from $h_{b_{\max}} = 10$ to $h_{b_{\min}} = 2$ mm and a length of $L_{b} = 245$ mm. These descriptions, together with the natural assumption of no shear in the reference configuration define the initial condition for the kinematics maps $g_{b}$ and $g_{s}$ . The mechanical, dynamical, and geometrical parameters are summarized for the SRA and the SSM in Table 1. For what concerns the rigid root-body, the geometrical and inertia parameters have been measured directly whereas the dynamical parameters have been calibrated manually to minimize the experimental error, which includes the effect of the buoyancy module. Each of the four bars has a length equal to $L_{rr} = 141$ mm, and a width equal to $h_{rr} = 5$ mm. The total mass of the root-body is 204 g, hence we assume the mass of a single bar to be $m_{rr} = 51$ g. The initial condition for the map $g_{r}$ is simply the identity I₄ (i.e. the Euclidean fixed frame $(o, e_{1}, e_{2}, e_{3})$ and the reference frame $(o, E_{1}, E_{2}, E_{3})$ coincide at t=0). The dynamical and geometrical parameters of the root-body are summarized in Table 2. The constant maps $g_{rs}$ and $g_{rbβ}$ are fixed and determined by the design of the prototype. In our case we have (reference unit [mm]):

\begin{matrix} g_{rb 1} = (\begin{matrix} 0.70 & - 0.71 & 0.06 & 76 \\ 0.70 & 0.71 & - 0.06 & 76 \\ 0.09 & 0 & 0.99 & - 90 \\ 0 & 0 & 0 & 1 \end{matrix}), \\ g_{rs} = (\begin{matrix} 0 & 0 & 1 & 90 \\ 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & - 50 \\ 0 & 0 & 0 & 1 \end{matrix}), \\ g_{rr 1} = (\begin{matrix} 0.45 & - 0.71 & - 0.54 & 38 \\ 0.45 & 0.71 & - 0.54 & 38 \\ 0.77 & 0 & 0.64 & - 45 \\ 0 & 0 & 0 & 1 \end{matrix}) \end{matrix}

whereas $g_{rb 2}$ , $g_{rb 3}$ , $g_{rb 4}$ and $g_{rr 2}$ , $g_{rr 3}$ , $g_{rr 4}$ are obtained by repetitively rotating $g_{rb 1}$ and $g_{rr 1}$ , respectively, by $π ∕ 2$ with respect to the axis E₃.

Table 1.

Parameters of the SRA and the SSM.

Beam parameter	Value	Shell parameter	Value
Young modulus E	110 kPa	Young modulus E	40 kPa
Poisson ratio $ν$	0	Poisson ratio $ν$	0
Shear viscosity modulus $υ$	300 Pa $\cdot$ s	Shear viscosity modulus $υ$	500 Pa $\cdot$ s
Mass density $ρ_{b}$	$1080 kg ∕ m^{3}$	Mass density $ρ_{s}$	7415 kg/m³
Hydrodynamic coefficient $C_{bx}$	0.01	Drag coefficient $C_{d}$	1.7
Hydrodynamic coefficient $C_{by}$	2.5	Added mass coefficient $B_{s}$	1.1
Hydrodynamic coefficient $C_{bz}$	2.5	Flow loss coefficient $C_{f}$	1
Hydrodynamic coefficient $B_{b}$	1.5	Length $L_{s}$	118 mm
Length $L_{b}$	245 mm	Thickness $h_{s}$	1 mm
Maximum radius $h_{b_{\max}}$	10 mm	Cylinder radius $T_{s}$	31 mm
Minimum radius $h_{b_{\min}}$	2 mm	Outflow area $A_{no}$	491 mm²
		Inflow area $A_{ni}$	$1400 m m^{2}$
		Left limit actuation $X_{s_{1}}$	25 mm
		Right limit actuation $X_{s_{2}}$	110 mm

Table 2.

Parameters of the rigid root-body for the Poseidron, Quadropus, and Monopus.

Parameter Poseidron	Value	Parameter Quadropus and Monopus	Value
Parallelepiped mass $m_{rr}$	51 g	Body mass $m_{rr}$	204 g
Drag coefficient $C_{r}$	0.02 m²	Drag coefficient $C_{r}$	$0.02 m^{2}$
Added mass coefficient $B_{r}$	0	Added mass coefficient $B_{r}$	0
Parallelepiped length $L_{rr}$	141 mm	Body length $L_{rr}$	112 mm
Parallelepiped thickness $h_{rr}$	5 mm	Body radius $h_{rr}$	20 mm

To derive a force input for the model, readings from the motor encoder are employed. Given a torque constant of 6.6 mNm/A, a motor maximum efficiency of 79% and a gearhead transmission efficiency of 73%, the estimate of time-varying experimental torque output, and hence model force input, is computed from current data as depicted in Figure 9. The harmonic oscillations depicted in Figure 9 are associated with the stages of inflation and deflation of the elastic shell which result in the periodic pull and release exerted upon the cables, as discussed at length in Giorgio-Serchi et al. (2016). The input of the model is the rhythmic actuation of the SSM provided by the cables. This can be modeled by taking the radial force $f_{r} (X_{s}, t)$ as a T -periodic function with two phases: a contraction phase and a relaxation phase, with the period T being the inverse of the motor frequency of the experiments. The actuation pressure $f_{r} (X_{s}, t)$ has been taken equal to zero for all the $X_{s}$ except for a central interval of application $[X_{s_{1}}, X_{s_{2}}]$ (equal to 85 mm in our case). At these points, $f_{r} (X_{s}, t)$ has been modeled as the ratio of the sum of the cable tension ( $F (t)$ , Figure 9) and the surface of the area of application, giving $f_{r} (X_{s}, t) = - F (t) ∕ \int_{X_{s_{1}}}^{X_{s_{2}}} 2 π r^{⋆} d X_{s}$ . The SRAs have not been actuated in these experiments, i.e. $F_{ba} (X_{b}, t) = 0$ . It can be shown that to have a radial pressure on the SSM with distribution $f_{r} (X_{s}, t)$ the internal actuation wrench takes the form

F_{sa} (X_{s}, t) = (\begin{matrix} R_{s}^{T} (\begin{matrix} - \sin (ϕ) \\ \cos (ϕ) \\ 0 \end{matrix}) \frac{f_{r}}{r^{⋆}} (\int_{X_{s}}^{X_{s_{2}}} z r^{⋆} d X_{s} - z \int_{X_{s}}^{X_{s_{2}}} r^{⋆} d X_{s}) \\ R_{s}^{T} (\begin{matrix} \cos (ϕ) \\ \sin (ϕ) \\ 0 \end{matrix}) \frac{f_{r}}{r^{⋆}} \int_{X_{s}}^{X_{s_{2}}} r^{⋆} d X_{s} \end{matrix})

(24)

where the lower bound of the integral therein becomes equal to $X_{s_{1}}$ for all $X_{s} \leq X_{s_{1}}$ and $X_{s_{2}}$ for all $X_{s} \geq X_{s_{2}}$ . Furthermore, to meet the experimental constraint of traveling on a straight line, the pitch and lift motions have been controlled.

The results of the comparison for three motor frequencies (1.89, 1.51, and 1.26 pps) are shown in Figure 10. The distance between the two values has been evaluated with respect to the mean swimming velocity ${\bar{V}}_{E_{3 s}}$ , thus the error is defined as: $e = | {\bar{V}}_{E_{3 s}} - {\bar{V}}^{*} | ∕ {\bar{V}}^{*}$ where ${\bar{V}}^{*}$ is the average experimental swimming speed. The error for the three experiments are 15.5 %, 18.4 %, and 22.9 %, respectively, from the fastest (1.89 pps) to the slowest (1.26 pps) case.

Fig. 10.

Real and simulated swimming velocity of the soft robot for tests performed at 1.89, 1.51, and 1.26 pps.

At present, the main source of error pertains to the hydrodynamic loads prediction that is largely affected by the difficulty to estimate in closed-form solution the contribution from the time-varying shape variations inside and around the SSM body. Inertial and viscous effect during expulsion and suction of fluid from the cavity do represent prominent terms in the dynamics of the shell. However, given the accuracy of the validation, it is reasonable to expect that, for the range of actuation frequencies investigated, neglecting these terms may represent an acceptable assumption.

7. Exploration of alternative designs

In this section, we explore the capabilities of the model to predict the dynamics of new conceptual prototype which are based on different arrangements of the baseline reference structures. In particular, the behavior of an underwater soft robot, referred to as Quadropus due to its four-limbed structure, with four SRAs stacked at the back of the SSM is analyzed, and the navigation capabilities of a vehicle, called Monopus, with one single SRA used as a steering mechanism are shown. These examples demonstrate the capability of the model to deal with simultaneous actuation of different modules during six-dimensional underwater swimming scenarios. Whereas these analyses are not meant to be conclusive, we aim at demonstrating the flexibility of the present model to treat a broad range of geometrical configurations and actuation routines and how these can be used to infer critical design parameters relevant to robot design such as power consumption and control optimization. Ultimately, this kind of analysis may be of value to study aquatic living organisms to derive an accurate biomechanical characterization of their swimming strategies.

In Figure 11, the new configurations of the soft bodies and the rigid root-body for the Quadropus and the Monopus is illustrated. From a modeling perspective, the description of these new morphological configuration is reflected in the value of the constant transformations $g_{rr}$ , $g_{rs}$ , $g_{rbβ}$ (Quadropus), and $g_{rb}$ (Monopus), which are the only variations that need to be implemented to model the new systems. The values used in the simulations are (reference unit [mm]):

\begin{matrix} g_{rb} = (\begin{matrix} 0 & 0 & 1 & - 110 \\ 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}), g_{rb 1} = (\begin{matrix} 0 & 0 & 1 & - 112 \\ 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & - 10 \\ 0 & 0 & 0 & 1 \end{matrix}), \end{matrix}

\begin{matrix} g_{rs} = (\begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}), g_{rr} = (\begin{matrix} 0 & 0 & 1 & - 84 \\ 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}) \end{matrix}

whereas $g_{rb 2}$ , $g_{rb 3}$ , and $g_{rb 4}$ (Quadropus) are obtained by repetitively rotating $g_{rb 1}$ by $π ∕ 2$ with respect to the axis E₁ (Figure 11).

Fig. 11.

illustrative scheme of the kinematics of the new structures explored. The new structures are obtained simply by modifying the value of the constant transformations $g_{rs}$ and $g_{rb}$ and the number of SRAs.

For both of these designs, the parameters that define the SRA and SSM are kept as in the previous section and summarized in Table 1. The rigid root-body, on the other hand, is now a solid cone 112 mm long, with a base radius of 20 mm and the same mass and hydrodynamic coefficients of the prototype root-body (Table 2). Hence, the only variation introduced in terms of parameters lies in the different rotational inertia of the cone geometry with respect to the four-bar system.

7.1. Quadropus dynamics

The Quadropus represents an octopus-like body that is able to exploit the combined effects of pulsed-jetting as well as the sculling of its tentacles. Pulsed-jet activation is implemented as a force acting along the circumference of the shell and given by $F_{sa}^{2} (X_{s}, t) = {(0, 0, 0, 0, 0, - F_{za} (X_{s}, t))}^{T}$ , where $F_{za} (X_{s}, t)$ is a time-dependent triangular-shaped function with a frequency of 1.51 pps and a negative amplitude (i.e. contracting force) from 0 to $- 15 ∕ r^{⋆} (X_{s})$ . To qualitatively estimate the swimming performance of the Quadropus, the locomotion velocity profile is compared for different actuation of the SRAs. As a matter of fact, the Quadropus design can benefit of impulsive acceleration thanks to the coordinated actuation of the SRAs at the beginning of the sequence. This swimming mode is used by the octopus during fast escape maneuvers and has been thoroughly studied and implemented in a octopus inspired prototype in Sfakiotakis et al. (2014). To reproduce such behavior in our model, the following actuation wrench for each SRAs is used:

F_{ba} (X_{b}, t) = (\begin{matrix} 0 \\ 0 \\ (m_{baz} - \frac{m_{baz}}{L_{b} t ∕ Δt} X_{b}) H (L_{b} t ∕ Δt - X_{b}) \\ 0 \\ 0 \\ 0 \end{matrix})

(25)

for $0 < t \leq Δt$ , whereas $F_{ba} (X_{b}, t) = {(0, 0, 0, 0, 0, 0)}^{T}$ otherwise. Here, $H (\cdot)$ is the Heaviside step function, $m_{baz}$ is the maximum torque load experienced by the $X_{b} = 0$ section, and $Δt$ is the time of recruitment of all the sections. Hence, the function $F (X_{b}, t) = (m_{baz} - \frac{m_{baz}}{L_{b} t ∕ Δt} X_{b})$ assigns a linear decreasing torque load to each section $X_{b}$ whereas the function $H (L_{b} t ∕ Δt - X_{b})$ recruits a portion of the sections (from the base to the tip) proportional to the portion of time elapsed before the total recruitment time $Δt$ .

The employment of different amount of torque in the actuation of the SRMs for $m_{baz}$ respectively equal to 0, 5, and 7 ([mN-m]) yields different swimming speeds shown in Figure 12, where the high acceleration of the Quadropus with active arms can be appreciated. A sequence of snapshots of the fast escape maneuver are presented in Figure 13.

Fig. 12.

Comparison between the swimming performance of the Quadropus under different SRAs actuations.

Fig. 13.

Few snapshots of the fast escaping maneuver of the Quadropus with active arms ( $m_{baz} = 5$ [mN-m]).

Decomposition of the force contributions to the whole-body dynamics from the actuated components can be observed in Figure 14. The force exerted by the arms at the three actuation values can be appreciated during the initial phase of arm expansion ( $1.0 < t < 2.0$ ) followed by the thrust peak coincident with the arm sudden closure. From this point onwards, the arms remain unactuated and the effect of their passive elasticity is projected on the root-body dynamics as a periodic retarding effect. This is due to alternate phases of stretching and compression of the back-bone under the co-axial contribution of the thrust force generated by the pulsed jetting. This analysis highlights how force characterization in such complex morphological structures can be performed in a segregated manner.

Fig. 14.

Comparison of force exchange between the arms and root-body in the three actuation scenarios. The load exerted by the pulsating shell is also reported.

We estimate efficiency of this system by using the definition of Maertens et al. (2015) for fish. This is the quasi-propulsive efficiency $η_{QP} = \frac{\bar{P_{out}}}{\bar{P_{in}}}$ and it is defined by the ratio between the useful work $P_{out}$ and the expended energy $P_{in}$ , i.e. the power required to overcome the resisting viscous forces at the terminal speed recorded

\begin{matrix} P_{out} = - \int_{0}^{L_{s}} \int_{0}^{2 π} (0, 0, V_{E_{3 s}}) R_{s} d_{s} r^{⋆} dϕ d X_{s} \\ - \int_{0}^{L_{b}} {(v_{b} +^{b} v_{r})}^{T} d_{b} d X_{b} - v_{r}^{T} d_{r} \end{matrix}

(26)

and the power required to actuate the pulsed jetting with

P_{in} = \int_{0}^{L_{s}} \int_{0}^{2 π} F_{sa}^{2 T} {\dot{ξ}}_{s 2} r^{⋆} dϕ d X_{s}

(27)

Because the analysis must be performed at steady state, the initial contribution from the sculling arms does not participate in the estimation of the quasi-propulsive efficiency. The comparison between actuation power and steady-swimming power is portrayed in Figure 15. Based on these, an estimate efficiency of 31% is inferred, which falls within the range of values observed for fish in Maertens et al. (2015).

Fig. 15.

Comparison between the input and output power during swimming at steady-state conditions.

7.2. Monopus dynamics

A further example of the capabilities of the presented model is provided by the case of a sperm-like body actuated by the synergic actuation of a pulsed-jetting shell and by the flapping of a rear-pointing soft manipulator. The SRA, in this case, is studied as a mean of steering actuator for the system. This example lends itself to the analysis of the turning moments generated by the actuation of differential parts within the system. In this case, the SSM of the Monopus is actuated as in the previous example with a frequency of 1.51 pps, whereas the SRA is bent in different directions to produce a 6D underwater turning motion. The actuation wrench for the SRA takes the form

\begin{matrix} F_{ba} (X_{b}, t) = (\begin{matrix} 0 \\ 0 \\ F (X_{b}, t) H (L_{b} t ∕ Δ t_{1} - X_{b}) \\ 0 \\ 0 \\ 0 \end{matrix}), \\ F_{ba} (X_{b}, t) = (\begin{matrix} 0 \\ F (X_{b}, t) H (\frac{L_{b}}{Δ t_{3}} (t - Δ t_{1} - Δ t_{2}) - X_{b}) \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}) \end{matrix}

(28)

respectively for $0 < t \leq Δ t_{1}$ and $Δ t_{1} + Δ t_{2} < t \leq Δ t_{1} + Δ t_{2} + Δ t_{3}$ , whereas $F_{ba} (X_{b}, t) = {(0, 0, 0, 0, 0, 0)}^{T}$ otherwise. This implies that the SRA is first gradually bent and released, from the base to the tip, around the local z -axis and than around the local y -axis. This results in a smooth and natural rotation, as reported in the sequence of Figure 16.

Fig. 16.

Few snapshots of the 6D underwater locomotion of the Monopus ( $m_{baz} = 5$ [mN-m]).

The dynamics of the Monopus is decomposed with the purpose of inferring baseline metrics that are commonly employed in manoeuvring tests of aquatic vehicles. With reference to the screw parameters (Murray et al., 1994):

\begin{matrix} r = \begin{matrix} \frac{∥ {\tilde{w}}_{r} v_{r} ∥}{(L_{s} + L_{b}) ∥ w_{r} ∥^{2}} & for & w \neq 0 \end{matrix} \\ m = \begin{matrix} ∥ w_{r} ∥ & for & w_{r} \neq 0 \end{matrix} \end{matrix}

(29)

These enable the derivation of the non-dimensional time-dependent radius of curvature defined as the distance from the axis of rotation scaled by the total length of the Monopus ( $r [-]$ ), as shown in Figure 17. Similarly, the magnitude of angular velocity ( $m [rad / s]$ ) can be derived to yield the values depicted in Figure 17. The turning radius varies according to the speed of the body, thus explaining the initial values of about $1$ and the latter decrease to $~ 0.3$ , which falls within the range commonly encountered for sea-dwelling organisms (Domenici et al., 2004). As the body approaches a straight swimming direction, the radius of curvature tends to infinity.

Fig. 17.

Turning radius and speed of angular velocity of the Monopus for two values of the SRA actuation.

8. Conclusion

In this work, a geometrically exact model for underwater soft robots is presented. The model is capable of representing a group of soft bodies connected together via a rigid root-body; this is done by taking into account the geometrical nonlinearity of the soft bodies (treated here as a Cosserat medium) along with the elastic responses, the mechanical actuation and the inertial loads exchanged by the interconnected bodies. This model is of general applicability since the dynamic equations of the soft bodies are derived from the unique knowledge of the Lagrangian density by means of a continuous extension of Poincaré’s equations. To the best of the authors’ knowledge, this is the first example in the literature where such a multi-soft-body geometrically exact model has been presented and used for the a priori evaluation of possible robotics design.

The work described in this manuscript constitutes the first milestone in our road map to the modeling and control of soft robots inspired from flexible aquatic organisms such as cephalopods. The difficulty to exhaustively address the coupled fluid–solid interactions acting upon the shell currently represents the major limitation of this formulation.

In particular, the model of the internal and external pressure due to the soft robot motion is simplified using the well-known solution of a one-dimensional momentum equation for a neutrally buoyant, rigid body translating in water (as explained in Renda et al. (2015b)) that drastically neglects potentially significant dynamical terms associated with body-shape modification. Indeed, the formal definition of the total propulsive force is simplified. Although it does not prevent the model from capturing the overall dynamics of the vehicle in surge motion, further refinements are sought for in order to attain a more sound representation of the complex physics involved in the system. Evidence is emerging that a substantial contribution to the total thrust of jetting bodies lies in the linked internal–external fluid dynamics of the collapsing shell. To acknowledge these terms in detail, a fully coupled fluid–structure interaction solver is mandatory. However, for the purpose of controller design and quick dynamics investigation, a fast and computationally inexpensive alternative must be looked into. This can be found from a coupled model that encompasses a quasi-analytic solutions of unsteady potential flows of radially varying slender, axisymmetric bodies (see Karamcheti (1966) and Anderson et al. (2001)) and an integral description of the internal pressure of the deforming cavity body, as in (Krieg and Mohseni, 2015). In this formulation, a prescribed kinematics of the shape-changing body enables the estimation of the added-mass variation as well as supplying the condition for accurate estimation of the forces generated by the expulsion of fluid across the nozzle. Despite the limited degree of reproducibility of the experiments, the model has been experimentally compared with a real multi-soft-body prototype with satisfactory results and then used to explore the design space of underwater soft robots characterized by different morphologies and actuation strategies. We show the varied range of dynamical analysis that can be performed on the newly designed conceptual prototypes by deriving quasi-propulsive efficiency of a four-arm pulsed-jetting octopus-like body and the time-dependent radius of curvature of a sperm-like vehicle. These demonstrate the capability of this mathematical formulation to represent an unlimited range of possible designs as well as to perform a priori evaluation of their maneuvering capabilities and swimming performances.

The authors believe that the versatility, accuracy and conceptual simplicity of the model presented here make this approach one of the most suitable in the frame of mobile soft robotics.

Footnotes

Appendix

Funding

This work was partially supported by Foundation Grant project PoseiDRONE of the Cassa di Risparmio di Livorno. Dr. F. Giorgio-Serchi is supported by the Natural Environmental Research Council (grant number NE/P003966/1) and by the SMMI-HEIF grant.

ORCID iD

Federico Renda

Jorge Dias

References

Abbott

Cosentino Lagomarsino

Zhang

Dong

Nelson

(2009) How should microrobots swim? The International Journal of Robotics Research 28(11–12): 1434–1447.

Anderson

DeMont

(2000) The mechanics of locomotion in the squid Loligo pealei: locomotory function and unsteady hydrodynamics of the jet and intramantle pressure. Journal of Experimental Biology 203(18): 2851–2863.

Anderson

Quinn

DeMont

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

Antman

(2006) Nonlinear Problems of Elasticity (Applied Mathematical Sciences, vol. 107). New York: Springer.

Boyer

Porez

(2015) Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects. Bioinspiration and Biomimetics 10(2): 025007.

Boyer

Porez

Khalil

(2006) Macro-continuous computed torque algorithm for a three-dimensional eel-like robot. IEEE Transactions on Robotics 22(4): 763–775.

Boyer

Porez

Leroyer

(2010) Poincaré cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: Application to robotics. Journal of Nonlinear Science 20(1): 47–79.

Boyer

Primault

(2005) The Poincaré–Chetayev equations and flexible multibody systems. Journal of Applied Mathematics and Mechanics 69(6): 925–942.

Boyer

Renda

(2016) Poincaré’s equations for Cosserat media: Application to shells. Journal of Nonlinear Science. DOI: 10.1007/s00332-016-9324-7.

10.

Canavin

Likins

(1977) Floating reference frames for flexible spacecraft. Journal of Spacecraft and Rockets 14(12): 724–732.

11.

Candelier

Boyer

Leroyer

(2011) Three-dimensional extension of Lighthill’s large-amplitude elongated-body theory of fish locomotion. Journal of Fluid Mechanics 674: 196–226.

12.

Candelier

Porez

Boyer

(2013) Note on the swimming of an elongated body in a non-uniform flow. Journal of Fluid Mechanics 716: 616–637.

13.

Colgate

Lynch

(2004) Mechanics and control of swimming: a review. IEEE Journal of Oceanic Engineering 29(3): 660–673.

14.

Conte

Modarres-Sadeghi

Watts

Hover

Triantafyllou

(2010) A fast-starting mechanical fish that accelerates at 40 m s². Bioinspiration and Biomimetics 5(3): 035004.

15.

Domenici

Standen

Levine

(2004) Escape manoeuvres in the spiny dogfish (squalus acanthias). Journal of Experimental Biology 207(13): 2339–2349.

16.

Elvander

Hawkes

(2012) ROVs and AUVs in support of marine renewable technologies. In: Oceans, 2012, pp. 1–6.

17.

Featherstone

(2014) Rigid Body Dynamics Algorithms. Springer US.

18.

Giorgio-Serchi

Arienti

Corucci

Giorelli

Laschi

(2017) Hybrid parameter identification of a multi-modal underwater soft robot. Bioinspiration and Biomimetics. DOI: 10.1088/1748-3190/aa5ccc.

19.

Giorgio-Serchi

Arienti

Laschi

(2016) Underwater soft-bodied pulsed-jet thrusters: actuator modelling and performance profiling. The International Journal of Robotics Research 35: 1308–1329.

20.

Giorgio-Serchi

Renda

Calisti

Laschi

(2015) Thrust depletion at high pulsation frequencies in underactuated, soft-bodied, pulsed-jet vehicles. In: OCEANS 2015, Genova, pp. 1–6.

21.

Giorgio-Serchi

Weymouth

(2016a) Drag cancellation by added-mass pumping. Journal of Fluid Mechanics 798: R3.

22.

Giorgio-Serchi

Weymouth

(2016b) Underwater soft robotics, the benefit of body-shape variations in aquatic propulsion. Biosystems and Biorobotics 17: 37–46.

23.

Hover

Eustice

Kim

Englot

Johannsson

Kaess

Leonard

(2012) Advanced perception, navigation and planning for autonomous in-water ship hull inspection. The International Journal of Robotics Research 31(12): 1445–1464.

24.

Johnson

Soden

Trueman

(1972) A study in jet propulsion: An analysis of the motion of the squid, Loligo vulgaris. Journal of Experimental Biology 56(1): 155–165.

25.

Karamcheti

(1966) Principles of ideal-fluid aerodynamics. Krieger.

26.

Krieg

Mohseni

(2015) Pressure and work analysis of unsteady, deformable, axisymmetric, jet producing cavity bodies. Journal of Fluid Mechanics 769: 337–368.

27.

Krieg

Sledge

Mohseni

(2015) Design considerations for an underwater soft-robot inspired from marine invertebrates. Bioinspiration & Biomimetics 10(6): 065004.

28.

Licht

Polidoro

Flores

Hover

Triantafyllou

(2004) Design and projected performance of a flapping foil auv. IEEE Journal of Oceanic Engineering 29(3): 786–794.

29.

Lighthill

(1970) Aquatic animal propulsion of high hydromechanical efficiency. Journal of Fluid Mechanics 44: 265–301.

30.

Linn

Lang

Tuganov

(2013) Geometrically exact Cosserat rods with Kelvin–Voigt type viscous damping. Mechanical Sciences 4(1): 79–96.

31.

Luh

JYS

Walker

Paul

RPC

(1980) On-line computational scheme for mechanical manipulator. Transactions of the ASME, Journal of Dynamical Systems 102: 69–76.

32.

Maertens

Triantafyllou

Yue

DKP

(2015) Efficiency of fish propulsion. Bioinspiration and Biomimetics 10: 046013.

33.

Marchese

Onal Cagdas

Rus

(2014) Autonomous soft robotic fish capable of escape maneuvers using fluidic elastomer actuators. Soft Robotics 1(1): 75–673.

34.

Mortl

Lawitzky

Kucukyilmaz

Sezgin

Basdogan

Hirche

(2012) The role of roles: Physical cooperation between humans and robots. The International Journal of Robotics Research 31(13): 1656–1674.

35.

Murray

Sastry

(1994) A Mathematical Introduction to Robotic Manipulation. London: Taylor & Francis.

36.

Poincaré

(1901) Sur une forme nouvelle des equations de la mecanique. Compte Rendu de l’Academie des Sciences de Paris 132: 369–371.

37.

Reissner

(1990) Special issue: Frontiers in computational mechanics on a one-dimensional theory of finite bending and stretching of elastic plates. Computers and Structures 35(4): 417–423.

38.

Renda

Giorelli

Calisti

Cianchetti

Laschi

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

39.

Renda

Giorgio-Serchi

Boyer

Laschi

(2015a) Locomotion and elastodynamics model of an underwater shell-like soft robot. In: 2015 IEEE International Conference on Robotics and Automation (ICRA), pp. 1158–1165.

40.

Renda

Giorgio-Serchi

Boyer

Laschi

(2015b) Modelling cephalopod-inspired pulsed-jet locomotion for underwater soft robots. Bioinspiration and Biomimetics 10(5): 055005.

41.

Renda

Giorgio-Serchi

Boyer

Laschi

(2015c) Structural dynamics of a pulsed-jet propulsion system for underwater soft robots. International Journal of Advanced Robotic Systems 12(6): 60143.

42.

Renda

Giorgio-Serchi

Boyer

Laschi

Dias

Seneviratne

(2018) A multi-soft-body dynamic model for underwater soft robots. In: Robotics Research. Cham: Springer, pp. 143–160.

43.

Saimek

(2001) Motion planning and control of a swimming machine. In: Proceedings of the 2001 American Control Conference, vol. 1, pp. 125–130.

44.

Selig

(2007) Geometric Fundamentals of Robotics (Monographs in Computer Science). New York: Springer.

45.

Sfakiotakis

Kazakidi

Chatzidaki

Evdaimon

Tsakiris

(2014) Multi-arm robotic swimming with octopus-inspired compliant web. In: 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2014), pp. 302–308.

46.

Simo

(1985) A finite strain beam formulation. the three-dimensional dynamic problem. Part I. Computer Methods in Applied Mechanics and Engineering 49(1): 55–70.

47.

Simo

Fox

(1989) On stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization. Computer Methods in Applied Mechanical Engineering 72(3): 267–304.

48.

Vaganay

Gurfinkel

Elkins

Jankins

Shurn

(2009) Hovering autonomous underwater vehicle-system design improvements and performance evaluation results. In: International Symposium on Unmanned Untethered Submarine Technology.

49.

Vasilescu

Detweiler

Doniec

, et al. (2010) AMOUR V: A hovering energy efficient underwater robot capable of dynamic payloads. The International Journal of Robotics Research 29(5): 547–570.

50.

Wang

Iida

(2015) Deformation in soft-matter robotics: A categorization and quantitative characterization. IEEE Robotics Automation Magazine 22(3): 125–139.

51.

Woodman

Winfield

Harper

Fraser

(2012) Building safer robots: Safety driven control. The International Journal of Robotics Research 31(13): 1603–1626.

52.

Ding

Yang

Tan

Wang

Zhang

(2012) On a bio-inspired amphibious robot capable of multimodal motion. IEEE/ASME Transactions on Mechatronics 17(5): 847–856.