A mechanical model for a system of coupled elastic filaments for terrestrial robot locomotion

Abstract

Multi-legged robots are a class of devices offering potential solutions in several operational scenarios that involve the deployment in unstructured and uncertain environments, in which trade-offs between specialization and robustness need to be considered. Several aquatic and terrestrial organisms evolutionarily converged to robust locomotion solutions. An important example is the emergence of collective beating dynamics in coupled arrays of flexible protrusions, which is at the core of the locomotion mechanics of small swimmers, and of terrestrial walkers with flexible, elongated bodies. Here, we formulate the dynamics of a system of flexible elastic filaments coupled through a solid medium, intended to be a simplified model for the locomotion mechanism for a legged terrestrial robot. The legs’ coupling is modeled via linear elastic lumped elements, and metachronal wave patterns are enforced via ad hoc leg actuation. Simulation results show that this model predicts the persistence of wave patterns that result in locomotion across a flat terrain.

1. Introduction

The design of locomotion systems for robot devices is dictated by the environment in which they operate. Solutions for robots deployed in well-structured environments, performing highly specialized tasks, are optimized for the related narrow operational conditions. On the other hand, there are several applications, including non-destructive environmental inspection [1–3], commercial implementations [4,5], military and defence applications [6–8], and operation in hazardous and generally not accessible environments [9–11], that are characterized by unstructured and/or unknown environments, for which robust locomotion solutions with low specialization with respect to the workspace are necessary. Further examples include a specific design to tackle rough terrains [12] and sensor-based methods used to efficiently navigate through unstructured environments [13].

Robust locomotion systems have been developed by several biological organisms through evolutionary adaptation processes. In particular, several small-scale organisms evolved hair-like beating protrusions called cilia or flagella for mobility in fluid environments [14,15]. The core structure of these protrusions is the axoneme, which consists of a cylindrical arrangement of elastic filaments connected by motor proteins acting as actuators [16]. In fluidic environments, micro-organisms actuate axonemal protrusions synchronously in different patterns, resulting in locomotion due to coupling with the surrounding fluids. Complex beating patterns such as synchronization, phase locking, and metachronal waves can emerge from synchronous oscillations of the protrusions [17], resulting in efficient swimmer-type locomotion. Metachronal coordination in an array of filaments is a synchronous pattern in which all protrusions oscillate with a constant phase shift among contiguous ones, resulting in traveling waves through the array. This has important implications for multi-legged terrestrial locomotion [18,19], since from a design standpoint, arrays of flexible protrusions used as legs can be a robust solution to navigate rough terrains as they can adapt and morph to suit a variety of structures. Furthermore, metachronal coordination is a suitable gait for multi-legged robot locomotion [20,21].

In the operation of autonomous robots, minimization of energy consumption is a crucial constraint. Therefore, it is important to understand what type of actuation is minimally necessary for inducing coordination patterns that are suitable for locomotion. Previous work [22] has shown that emergent synchronous patterns are aided by hydrodynamic and basal coupling among the filaments. Therefore, instead of purely relying on internal forces, microorganisms can actuate the locomotion protrusions through hydrodynamic and basal coupling. This suggests that pure basal coupling can achieve beating patterns in ciliary and flagellar arrays in the absence of fluid interaction forces, with important implications for terrestrial locomotion, and ultimately for the energy management of autonomous multi-legged walkers.

Here, we extend the work in [23] by presenting the mechanical model of a terrestrial planar robot consisting of an elongated rigid body, representing the payload, and of a locomotion system consisting of an array of flexible filaments coupled at the base through elastic elements hosted in the body. The stiff filaments connect the body to the terrain, therefore acting as legs. A strong formulation of the system is obtained by modeling the filaments as bending-stiff and inextensible, and a weak form is projected along a finite-dimensional subspace of the original solution space for numerical simulation. The filaments are actuated with time-varying moments designed to sustain metachronal coordination. The resulting motion mimics those observed in biological systems approximated by the model. The objective is to present a mechanical model for a class of terrestrial robots with robust locomotion mechanics and to illustrate in simulation the possibility of inducing metachronal coordinated patterns without internal axonemal actuation on each leg. This is the foundation for current and future work, focused on predicting the conditions for the emergence of coordinated patterns exploitable for locomotion, and on the understanding of the minimal number of actuation inputs that guarantee metachronal coordination and its sustainability through basal coupling.

2. Mechanical model

In this section, we present the mechanical model of a slender-legged robot in plane motion. In this simplified setting, the body is rigid and the legs are modeled as bending-stiff inextensible filaments.

The mechanical system, schematized in Figure 1, includes N flexible legs. A global coordinate frame with unit basis vectors $\widehat{\mathit{i}}$ and $\widehat{\mathit{j}}$ is adopted with origin placed in such a way that the $\widehat{\mathit{i}}$ basis vector is along the axis of symmetry of the body/payload. The $\widehat{\mathit{j}}$ basis vector aligns with the line joining the extremities of the legs attached to the body in the undeformed configuration.

Figure 1.

Schematics of the mechanical model.

Each leg, modeled as a bending-stiff inextensible elastic filament, connects to contiguous legs and to the body through linear springs as schematized. Referring to Figure 1, in which they are represented in their undeformed configurations, there is a total of $2N+1$ springs, partitioned into the set $\{S_{{\alpha }_1},\dots ,S_{{\alpha }_{N+1}}\}$ across contiguous pairs of legs, and into the set $\{S_{{\beta }_1},\dots ,S_{{\beta }_N}\}$ connecting one extremity of the legs to the body. The elements of the two sets have, respectively, stiffness α and β, and undeformed length $L_{\alpha }$ and $L_{\beta }$ .

In designing walking robots, one of the main challenges is modeling the contact mechanics with the terrain. This is especially true for the class of devices considered here, for which the bio-inspired locomotion mechanics is aimed at robustness with respect to the environment/terrain. However, in order to focus on the induction of metachronal coordination as locomotion mechanism, we restrict the analysis to a flat terrain parallel to the basis vector $\widehat{\mathit{j}}$ , and described by $x=X_t$ , where $X_t$ is the perpendicular distance along the axis $\widehat{\mathit{i}}$ , and x is the coordinate along the same axis.

2.1. Body mechanical model and Lagrangian coordinates

By referring to Figure 1, the robot’s body is modeled as a rigid U-shaped object of length ${\tau }_L$ , breadth ${\tau }_B$ , thickness ${\tau }_T$ , and clearance ${\tau }_C$ . We will also refer to this object as the payload, since in this simplified model it lumps the corresponding operational load of the robot. Given that the payload is a rigid object in a two-dimensional space, we use the generalized coordinates ${\mathit{q}}_b(t)=(X_b(t),Y_b(t),\theta (t))$ , where $X_b(t)$ and $Y_b(t)$ are the coordinates of the center of mass, and $\theta (t)$ is the orientation. The total payload mass is m.

The geometric description of the body is completed by defining the length ${\tau }_L=2{\tau }_T+(N+1)L_{\alpha }$ , the breadth ${\tau }_B={\tau }_T+L_{\beta }+{\tau }_C$ , and the area ${\mathrm{\Pi }}_A={\tau }_L{\tau }_T+2({\tau }_B-{\tau }_T){\tau }_T$ on the plane of motion. Assuming that at the initial time the supporting and connecting springs between the legs and the body are unstretched, the coordinates of the center of mass in the undeformed configuration are

$$\begin{array}{ccc}\multicolumn{1}{c}{X_b(0)=}& \multicolumn{1}{c}{\frac{{\tau }_T}{2{\mathrm{\Pi }}_A}\left(2({\tau }_C-L_{\beta })({\tau }_B-{\tau }_T)-{\tau }_L(2L_{\beta }+{\tau }_T)\right)}& \multicolumn{1}{c}{}\end{array}$$

(1a)

$$\begin{array}{ccc}\multicolumn{1}{c}{Y_b(0)=}& \multicolumn{1}{c}{0}& \multicolumn{1}{c}{}\end{array}$$

(1b)

In the adopted global frame, let ${\mathit{r}}_{B_i}(t)$ be the position vector of point $B_i$ , $i=1,\dots ,N+2$ at time t. In terms of the introduced Lagrangian coordinates, we have

$$\begin{array}{cc}\begin{array}{cc}\multicolumn{1}{c}{{\mathit{r}}_{B_i}(t)}& \multicolumn{1}{c}{={\mathit{r}}_C(t)+\mathit{R}(\theta (t)){\mathit{r}}_{B_{i}C}(0)}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=\left[\begin{array}{c}\multicolumn{1}{c}{X_b(t)}\\ \multicolumn{1}{c}{Y_b(t)}\end{array}\right]+\left[\begin{array}{cc}\multicolumn{1}{c}{\cos \theta (t)}& \multicolumn{1}{c}{-\sin \theta (t)}\\ \multicolumn{1}{c}{\sin \theta (t)}& \multicolumn{1}{c}{\cos \theta (t)}\end{array}\right]{\mathit{r}}_{B_{i}C}(0)}\end{array}& \end{array}$$

(2)

where ${\mathit{r}}_C(t)$ is the position of the center of mass in the same global frame, $\mathit{R}$ is the rotation matrix, and ${\mathit{r}}_{B_{i}C}(0)$ is the relative position of point $B_i$ with respect to the center of mass in the reference configuration, which is time invariant due to the rigid body constraint. The representation of ${\mathit{r}}_{B_{i}C}$ in terms of the body’s generalized coordinates is

$${\mathit{r}}_{B_{1}C}(0)=\left(-X_b(0),-\frac{N+1}{2}{L}_{\alpha }\right)$$

(3a)

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{r}}_{B_{i}C}(0)=\left(-X_b(0)-L_{\beta },-\left(\frac{N+1}{2}-(i-1)\right)L_{\alpha }\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{\mathrm{for}2\le i\le N+1}& \multicolumn{1}{c}{}\end{array}$$

(3b)

$${\mathit{r}}_{B_{N+2}C}(0)=\left(-X_b(0),\frac{N+1}{2}{L}_{\alpha }\right)$$

(3c)

2.2. Legs as stiff elastic filaments

The structure and bending behaviour of a cilium or flagellum can be effectively approximated via an elastic rod [24]. Mechanical modeling of the plane kinematics of this class of structures is the object of relatively extensive literature, such as [14,25], particularly in the context of beating patterns induced by coupling with highly viscous fluids. Considering that elastic rod models can describe large deformations and displacements undergone by thin filament-like structures with bending stiffness, we use them to model the system’s legs, which couple the body with the terrain and exert actions that result in motion.

As depicted in Figure 2, the planar section of this structure comprises two parallel elastic filaments separated by a constant distance a. These filaments are rigidly coupled at their base and are constrained from sliding relative to each other. The coordinate parameter s denotes the arc length along the neutral axis. All points along the neutral axis, with respect to the global coordinate frame, are defined using the position vector $\mathit{r}(s)=(X(s),Y(s))$ . Consequently, the filaments can be described using the position vectors ${\mathit{r}}_{\alpha }(s)=\mathit{r}(s)-(a∕2)\mathit{n}(s)$ and ${\mathit{r}}_{\beta }(s)=\mathit{r}(s)+(a∕2)\mathit{n}(s)$ , where $\mathit{n}(s)$ is a unit normal vector at point s on the neutral axis.

Figure 2.

Schematic of the simplified two-dimensional axonemal structure.

Given this geometry, the kinematics can be described via the neutral axis. By defining the related unit tangent vector $\mathit{t}(s)$ , the local geometry can be differentially characterized by the Frenet–Serret relations [26]

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{r}}^{\prime }=\mathit{t},{\mathit{t}}^{\prime }=\kappa \mathit{n},{\mathit{n}}^{\prime }=-\kappa \mathit{t}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(4)

Here, the prime symbol denotes derivatives with respect to the arc length s, and κ is the curvature. Since the tangent and normal are unit vectors, and by using the relations in (4), we have $\kappa =\parallel {\mathit{r}}^{\mathrm{\prime \prime }}\parallel =\sqrt{{X^{″}}^2+{Y^{″}}^2}$ . Additionally, $\psi (s)$ is the angle between the tangent vector and basis vector $\widehat{\mathit{i}}$ .

The filament can bend under the action of external forces and moments at the boundaries/extremities, and as a consequence of distributed actions that model the macroscopic effect of the interplay between molecular motors and elastic connections between the filaments constituting an axoneme [27]. This system of distributed forces, schematized in Figure 2, is equivalent to an internal couple of magnitude $a\mathit{f}$ that is dual to the relative sliding displacement $\mathrm{\Delta }(s)$ towards the distal end of the elastic filaments, given by to $\mathrm{\Delta }(s)=a{\int }_0^s\kappa (\sigma )\mathrm{d\sigma }$ [14].

Let ${\mathit{r}}_i(s_i,t)=(X_i(s_i,t),Y_i(s_i,t))$ be the generalized coordinates describing the neutral axis of the $i\mathrm{th}$ leg, where $s_i\in [0,\ell ]$ and ℓ is the length. The kinetic energy of the $i\mathrm{th}$ leg is

$${\mathcal{K}}_i=\frac{1}{2}\underset{0}{\overset{\ell }{\int }}\rho {\stackrel{·}{\mathit{r}}}_i\cdot {\stackrel{·}{\mathit{r}}}_{i}d{s}_i$$

(5)

where ρ is the mass density per unit line, and the superimposed dot denotes time derivative. The strain energy, governed by the curvature and bending of the system through internal forces [14,25], combined with the gravitational potential energy, yields the total potential energy of the $i\mathrm{th}$ leg

$${\mathcal{E}}_i=\underset{0}{\overset{\ell }{\int }}\left(\frac{K_M}{2}{\kappa }_i^2+f_i{\mathrm{\Delta }}_i+\mathrm{\rho g}(X_t-{\mathit{r}}_i\cdot \widehat{\mathit{i}})\right)d{s}_i$$

(6)

where $K_M$ is the bending stiffness of the filaments and $f_i=\parallel {\mathit{f}}_i\parallel$ is the magnitude of internal couple ${\mathit{f}}_i$ acting on the $i\mathrm{th}$ leg. Partial integration of (6) leads to [14]

$${\mathcal{E}}_i=\underset{0}{\overset{\ell }{\int }}\left(\frac{K_M}{2}{\kappa }_i^2-a{F}_i{\kappa }_i+\mathrm{\rho g}(X_t-{\mathit{r}}_i\cdot \widehat{\mathit{i}})\right)d{s}_i$$

(7)

where $F_i(s)=-{\int }_s^{\ell }{f}_i(\sigma )\mathrm{d\sigma }$ . The variation of the potential energy along the shape ${\mathit{r}}_i(s,t)$ gives

$$\begin{array}{cccc}\multicolumn{1}{c}{\delta {\mathcal{E}}_i=}& \multicolumn{1}{c}{{\int }_0^{\ell }\left(K_M{\mathit{r}}_i^{\mathrm{IV}}-{\left(\frac{a{F}_i}{\parallel {\mathit{r}}_i^{″}\parallel }{\mathit{r}}_i^{″}\right)}^{″}-\mathrm{\rho g}\widehat{\mathit{i}}\right)\cdot \delta {\mathit{r}}_{i}d{s}_i}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{+\left(K_M-\frac{a{F}_i}{\parallel {\mathit{r}}_i^{″}\parallel }\right){\mathit{r}}_i^{″}\cdot \delta {\mathit{r}}_i^{\prime }{|}_0^{\ell }-\left(K_M{\mathit{r}}_i^{\mathrm{III}}-{\left(\frac{a{F}_i}{\parallel {\mathit{r}}_i^{″}\parallel }{\mathit{r}}_i^{″}\right)}^{\prime }\right)\cdot \delta {\mathit{r}}_i{|}_0^{\ell }}& \multicolumn{1}{c}{}\end{array}$$

(8)

The details of the derivations are given in Appendix 1.

The inextensibility constraint is imposed on the filaments modeling the legs, under the hypothesis that the kinematics is purely described by the curvature [14]. Locally, the constraint is formalized as ${\mathit{r}}_i^{\prime }\cdot {\mathit{r}}_i^{\prime }=1$ [14,25]. In the work by Camalet and Jülicher [14], the constraint is enforced through the use of a Lagrange multiplier in the action functional. Here, we consider instead the penalty-based method by Bayo et al. [28]. Unlike the Lagrange multiplier approach, penalty methods do not introduce additional variables, thereby reducing the computational complexity at the cost of not-exact enforcement. The inextensibility of the $i\mathrm{th}$ leg in integral form is written as:

$$\underset{0}{\overset{\ell }{\int }}\left({\mathit{r}}_i^{\prime }\cdot {\mathit{r}}_i^{\prime }-1\right)d{s}_i=0.$$

(9)

The penalty method will be integrated in the weak form to enforce this constraint.

The Rayleigh dissipation function

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{R}}_i=\frac{\zeta }{2}\underset{0}{\overset{\ell }{\int }}{\stackrel{·}{\mathit{r}}}_i\cdot {\stackrel{·}{\mathit{r}}}_{i}d{s}_i}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(10)

is used to model viscous damping in the legs, where ζ is a constant dissipation parameter [29]. The damping force acting on the $i\mathrm{th}$ leg is the negative velocity gradient of the dissipation function [29], yielding

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{F}}_{D_i}=-\zeta \underset{0}{\overset{\ell }{\int }}{\stackrel{·}{\mathit{r}}}_{i}d{s}_i.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(11)

Consequently, the virtual work done by the damping forces on the $i\mathrm{th}$ leg is

$$\delta {\overline{W}}_{D_i}=-\zeta \underset{0}{\overset{\ell }{\int }}{\stackrel{·}{\mathit{r}}}_i\cdot \delta {\mathit{r}}_{i}d{s}_i.$$

(12)

As illustrated in Figure 3, each filament representing a leg is also subjected to spring reaction forces and a time-varying actuation moment at its base. Additionally, normal and frictional forces acting at the tip model contact interactions with the terrain. Let $\mathrm{\Delta }{x}_{{\alpha }_i}=L_{\alpha }-{\gamma }_{{\alpha }_i}$ be the change in length of spring $S_{{\alpha }_i}$ , where

$${\gamma }_{{\alpha }_i}=\{\begin{array}{cc}\multicolumn{1}{c}{\parallel {\mathit{r}}_1(0,t)-{\mathit{r}}_{B_1}(t)\parallel }& \multicolumn{1}{c}{\mathrm{for}i=1}\\ \multicolumn{1}{c}{\parallel {\mathit{r}}_i(0,t)-{\mathit{r}}_{i-1}(0,t)\parallel }& \multicolumn{1}{c}{\mathrm{for}2\le i\le N}\\ \multicolumn{1}{c}{\parallel {\mathit{r}}_N(0,t)-{\mathit{r}}_{B_{N+2}}(t)\parallel }& \multicolumn{1}{c}{\mathrm{for}i=N+1.}\end{array}$$

(13)

Figure 3.

Depiction of external forces and moments acting on the $i\mathrm{th}$ leg.

Similarly, let $\mathrm{\Delta }{x}_{{\beta }_i}=L_{\beta }-{\gamma }_{{\beta }_i}$ be the change in length of the spring $S_{{\beta }_i}$ with

$${\gamma }_{{\beta }_i}=\{\begin{array}{cc}\multicolumn{1}{c}{\parallel {\mathit{r}}_1(0,t)-{\mathit{r}}_{B_2}(t)\parallel }& \multicolumn{1}{c}{\mathrm{for}i=1}\\ \multicolumn{1}{c}{\parallel {\mathit{r}}_i(0,t)-{\mathit{r}}_{B_{i+1}}(t)\parallel }& \multicolumn{1}{c}{\mathrm{for}2\le i\le N-1}\\ \multicolumn{1}{c}{\parallel {\mathit{r}}_N(0,t)-{\mathit{r}}_{B_{N+1}}(t)\parallel }& \multicolumn{1}{c}{\mathrm{for}i=N.}\end{array}$$

(14)

The spring reaction force ${\mathit{F}}_{S_i}$ acting at the base/attachment point of the $i\mathrm{th}$ filament is

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{F}}_{S_1}=}& \multicolumn{1}{c}{\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_1}}{{\gamma }_{{\alpha }_1}}({\mathit{r}}_1(0,t)-{\mathit{r}}_{B_1}(t))+\frac{\beta \mathrm{\Delta }{x}_{{\beta }_1}}{{\gamma }_{{\beta }_1}}({\mathit{r}}_1(0,t)-{\mathit{r}}_{B_2}(t))+\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_2}}{{\gamma }_{{\alpha }_2}}({\mathit{r}}_1(0,t)-{\mathit{r}}_2(0,t))}& \multicolumn{1}{c}{}\end{array}$$

(15a)

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{F}}_{S_i}=}& \multicolumn{1}{c}{\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_i}}{{\gamma }_{{\alpha }_i}}({\mathit{r}}_i(0,t)-{\mathit{r}}_{i-1}(0,t))+\frac{\beta \mathrm{\Delta }{x}_{{\beta }_i}}{{\gamma }_{{\beta }_i}}({\mathit{r}}_i(0,t)-{\mathit{r}}_{B_{i+1}}(t))+\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_{i+1}}}{{\gamma }_{{\alpha }_{i+1}}}({\mathit{r}}_i(0,t)-{\mathit{r}}_{i+1}(0,t)}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathrm{for}2\le i\le N-1}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(15b)

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{F}}_{S_N}=}& \multicolumn{1}{c}{\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_N}}{{\gamma }_{{\alpha }_N}}({\mathit{r}}_N(0,t)-{\mathit{r}}_{N-1}(0,t))+\frac{\beta \mathrm{\Delta }{x}_{{\beta }_N}}{{\gamma }_{{\beta }_N}}({\mathit{r}}_N(0,t)-{\mathit{r}}_{B_{N+1}}(t))++\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_{N+1}}}{{\gamma }_{{\alpha }_{N+1}}}({\mathit{r}}_N(0,t)-{\mathit{r}}_{B_{N+2}}(t)).}& \multicolumn{1}{c}{}\end{array}$$

(15c)

Actuating moments are applied at the bases of all filaments to induce clockwise rotations at a desired frequency, thus enabling the robot to move forward through contact constraints with the terrain. A constant phase difference is set between contiguous legs, throughout its movement across the terrain. A frequency-locked motion of the legs with constant phase difference between contiguous ones corresponds to a metachronal gait [17,21]. Here, we induce it by actuating each leg with a suitable actuating moment to illustrate the ability of the proposed model to reproduce this type of motion. Ongoing and future developments include the study of the emergence of metachronal coordination based on the system’s parameter and operational condition, and the study of minimal actuation requirements to sustain locomotion-inducing gait patterns, leveraging on the coupling between legs. Let the desired frequency of rotation be ω. The moment applied at the base of the $i\mathrm{th}$ leg is

$$M_{0_i}(t)=K_p\left(-\omega -{\stackrel{·}{\psi }}_i(0,t)\right)$$

(16)

where $K_p$ is a proportional gain and ${\stackrel{·}{\psi }}_i(0,t)$ is the rate of change of the local tangent angle at the base. Here, we also refer to $\psi (0,t)$ as the basal angle. The negative sign for ω is used to induce clockwise legs’ rotation.

When the tip of the $i\mathrm{th}$ leg engages with the terrain, a normal force $F_{{\ell }_{x_i}}(t)$ and a frictional force $F_{{\ell }_{y_i}}(t)$ are applied at the point of contact to model the action of the terrain on the system. These forces are implemented as penalty functions when the terrain contact conditions are met [30]

$$F_{{\ell }_{x_i}}(t)=\{\begin{array}{cc}\multicolumn{1}{c}{ϵ\left(X_t-X_i(\ell ,t)\right)}& \multicolumn{1}{c}{\mathrm{if}{X}_i(\ell ,t)\ge X_t}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{\mathrm{otherwise}}\end{array}$$

(17a)

$$F_{{\ell }_{y_i}}(t)=\{\begin{array}{cc}\multicolumn{1}{c}{\nu \left(-{Ẏ}_i(\ell ,t)\right)}& \multicolumn{1}{c}{\mathrm{if}{X}_i(\ell ,t)\ge X_t}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{\mathrm{otherwise}.}\end{array}$$

(17b)

Here, for simplicity, the penalty functions are based on a linear Hooke model, in contrast to a Kelvin–Voigt model or a constraint-based approach in the form of Lagrange multipliers [30]. Combining the contact forces into a vector ${\mathit{f}}_{{\ell }_i}=(F_{{\ell }_{x_i}},F_{{\ell }_{y_i}})$ , the virtual work done by all the external forces and moments on the $i\mathrm{th}$ leg is

$$\delta {\overline{W}}_{{\mathrm{Ext}}_i}=M_{0_i}\delta {\psi }_i(0,t)+{\mathit{F}}_{S_i}\cdot \delta {\mathit{r}}_i(0,t)+{\mathit{F}}_{{\ell }_i}\cdot \delta {\mathit{r}}_i(\ell ,t).$$

(18)

From the expression of the variation of the curvature $\delta {\kappa }_i={\mathit{n}}_i\cdot \delta {\mathit{r}}_i^{″}$ [14,31], by using (4) and ${\kappa }_i(s)={\psi }_i^{\prime }(s)$ we obtain

$$\begin{array}{ccc}\multicolumn{1}{c}{\delta {\psi }_i={\mathit{n}}_i\cdot \delta {\mathit{r}}_i^{\prime }=\frac{{{\mathit{r}}_i}^{\mathrm{\prime \prime }}}{\parallel {\mathit{r}}_i^{\mathrm{\prime \prime }}\parallel }\cdot \delta {\mathit{r}}_i^{\prime }.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(19)

Also from ${\mathit{t}}_i=(X_i^{\prime },Y_i^{\prime })$ , ${\mathit{n}}_i=(-Y_i^{\prime },X_i^{\prime })$ [26], we have

$$\begin{array}{ccc}\multicolumn{1}{c}{\delta {\psi }_i=-Y_i^{\prime }\delta X_i^{\prime }+X_i^{\prime }\delta Y_i^{\prime }.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(20)

Hence, (18) is rewritten as:

$$\begin{array}{ccc}\multicolumn{1}{c}{\delta {\overline{W}}_{{\mathrm{Ext}}_i}}& \multicolumn{1}{c}{=M_{0_i}\frac{{\mathit{r}}_i^{″}(0,t)}{\parallel {\mathit{r}}_i^{″}(0,t)\parallel }\cdot \delta {\mathit{r}}_i^{\prime }(0,t)+{\mathit{f}}_{S_i}\cdot \delta {\mathit{r}}_i(0,t)+{\mathit{f}}_{{\ell }_i}\cdot \delta {\mathit{r}}_i(\ell ,t).}& \multicolumn{1}{c}{}\end{array}$$

(21)

2.3. Governing equations for the legs

The equations of motion for the legs are derived using the extended Hamilton’s principle, which involves the first variation of the Lagrangian ${\mathcal{L}}_i={\mathcal{K}}_i-{\mathcal{E}}_i$ and the virtual work done by the damping forces, external forces, and moments. By considering two arbitrary time instants $t_1$ and $t_2$ at which the generalized coordinates are prescribed, the variational principle is stated as follows:

$$\underset{t_1}{\overset{t_2}{\int }}\left(\delta {\mathcal{K}}_i-\delta {\mathcal{E}}_i+\delta {\overline{W}}_{D_i}+\delta {\overline{W}}_{{\mathrm{ext}}_i}\right)\mathrm{dt}=0.$$

(22)

The detailed calculations for the variation of the kinetic energy are presented in Appendix 1. Substituting into (22) the expressions for the variations of the kinetic energy and the potential energy, along with the virtual work done by the external forces and moments, we obtain the equations of motion of the $i\mathrm{th}$ leg as

$$\rho {\stackrel{\mathrm{··}}{\mathit{r}}}_i+\zeta {\stackrel{·}{\mathit{r}}}_i+K_M{\mathit{r}}_i^{\mathrm{IV}}-{\left(\frac{a{F}_i}{\parallel {\mathit{r}}_i^{″}\parallel }{\mathit{r}}_i^{″}\right)}^{\mathrm{\prime \prime }}-\mathrm{\rho g}\widehat{\mathit{i}}=\mathbf{0}$$

(23)

together with the inextensibility constraint

$$\underset{0}{\overset{\ell }{\int }}\left({\mathit{r}}_i^{\prime }\cdot {\mathit{r}}_i^{\prime }-1\right)d{s}_i=0.$$

(24)

The boundary conditions are

$$s=\ell :\{\begin{array}{c}\multicolumn{1}{c}{{\mathit{F}}_{{\ell }_i}=-K_M{\mathit{r}}_i^{\mathrm{III}}(\ell ,t)+{\left(\frac{a{F}_i}{\parallel {\mathit{r}}_i^{″}\parallel }{\mathit{r}}_i^{″}\right)}^{\prime }{|}_{s=\ell }}\\ \multicolumn{1}{c}{0=K_M\parallel {\mathit{r}}_i^{″}(\ell ,t)\parallel }\end{array}$$

(25a)

$$s=0:\{\begin{array}{c}\multicolumn{1}{c}{{\mathit{F}}_{S_i}=K_M{\mathit{r}}_i^{\mathrm{III}}(0,t)-{\left(\frac{a{F}_i}{\parallel {\mathit{r}}_i^{″}\parallel }{\mathit{r}}_i^{″}\right)}^{\prime }{|}_{s=0}}\\ \multicolumn{1}{c}{M_{0_i}=-K_M\parallel {\mathit{r}}_i^{″}(0,t)\parallel +a{F}_i(0,t).}\end{array}$$

(25b)

Here, ${\mathit{F}}_{{\ell }_i}$ , ${\mathit{F}}_{S_i}$ , and $M_{0_i}$ are, respectively, the contact forces, spring reaction forces, and the external moment applied on the $i\mathrm{th}$ leg. All boundary conditions are dynamic, or von Neumann, accounting for the actions of the restoring springs, the actuating moment at the base, and for the contact-with-terrain condition at the free-end, with the assumption that the contact moment is zero, that is, the terrain does no work on the rotation of the extremity in contact. Notice that the spring forces couple the legs to the body. The initial conditions will be specified in the simulation section.

2.4. Weak form of the legs’ equations and finite element formulation

The numerical solution of the boundary value problems governing the legs’ dynamics is obtained by first deriving a weak formulation of the equations and then applying the finite element method to find an approximate solution. By multiplying (23) with the virtual displacement $\delta {\mathit{r}}_i$ and integrating it over the spatial domain $[0,\ell ]$ , after integrating by parts and by substituting the boundary terms, we obtain

$$\begin{array}{cc}\begin{array}{cc}\multicolumn{1}{c}{\underset{0}{\overset{\ell }{\int }}}& \multicolumn{1}{c}{(\left(\rho {\stackrel{\mathrm{··}}{\mathit{r}}}_i+\zeta {\stackrel{·}{\mathit{r}}}_i-\mathrm{\rho g}\widehat{\mathit{i}}\right)\cdot \delta {\mathit{r}}_i+\left(K_M-\frac{a{F}_i}{\parallel {\mathit{r}}_i^{″}\parallel }\right){\mathit{r}}_i^{″}\cdot \delta {\mathit{r}}_i^{″})d{s}_i}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{-{\mathit{F}}_{{\ell }_i}\cdot \delta {\mathit{r}}_i(\ell ,t)-{\mathit{F}}_{S_i}\cdot \delta {\mathit{r}}_i(0,t)-M_{0_i}\frac{{\mathit{r}}_i^{″}(0,t)}{\parallel {\mathit{r}}_i^{″}(0,t)\parallel }\cdot \delta {\mathit{r}}_i^{\prime }(0,t)=0.}\end{array}& \end{array}$$

(26)

Details of the related calculations are presented in Appendix 1. By using

$$\begin{array}{ccc}\multicolumn{1}{c}{\frac{{\mathit{r}}_i^{″}}{\parallel {\mathit{r}}_i^{″}\parallel }\cdot \delta {\mathit{r}}_i^{\prime }=\delta {\psi }_i=-Y_i^{\prime }\delta X_i^{\prime }+X_i^{\prime }\delta Y_i^{\prime }}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(27)

and ${\mathit{r}}_i=(X_i,Y_i)$ , we get

$$\begin{array}{c}\multicolumn{1}{c}{{\int }_0^{\ell }(\rho \left({Ẍ}_i\delta X_i+\mathrm{Ÿ\delta }{Y}_i\right)+\zeta \left({Ẋ}_i\delta X_i+{Ẏ}_i\delta Y_i\right)-\mathrm{\rho g\delta }{X}_i+(K_M-\frac{a{F}_i}{\sqrt{X_i^{″2}+Y_i^{″2}}})\left(X_i^{″}\delta X_i^{″}+Y_i^{″}\delta Y_i^{″}\right))d{s}_i}\\ \multicolumn{1}{c}{-F_{{\ell }_{x_i}}(t)\delta X_i(\ell ,t)-F_{{\ell }_{y_i}}(t)\delta Y_i(\ell ,t)-F_{S_{x_i}}(t)\delta X_i(0,t)-F_{S_{y_i}}(t)\delta Y_i(0,t)}\\ \multicolumn{1}{c}{-M_{0_i}(t)\left(-Y_i^{\prime }(0,t)\delta X_i^{\prime }(0,t)+X_i^{\prime }(0,t)\delta Y_i^{\prime }(0,t)\right)=0}\end{array}$$

(28)

where $F_{S_{x_i}}$ and $F_{S_{y_i}}$ are the components of the spring reaction force in the global frame. Therefore, there are in total $2N$ trial solution fields $(X_i,Y_i)$ to which we associate the $2N$ test function fields $(\delta X_i,\delta Y_i)$ . In the context of mechanical systems, the test functions are virtual displacements [32].

The spatial domain $[0,\ell ]$ is discretized by n nodes with coordinates $s_{i_j}\in [0,\ell ]$ , for $i=1,\dots ,N$ and $j=1,\dots ,n$ . Trial functions and test functions are approximated by the expansions

$$\begin{array}{ccc}\multicolumn{1}{c}{X_i(s_i,t)=\sum _{k=1}^n\left({\varphi }_{{\alpha }_k}(s_i)a_{i_k}(t)+{\varphi }_{{\beta }_k}(s_i)b_{i_k}(t)\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(29a)

$$\begin{array}{ccc}\multicolumn{1}{c}{\delta X_i(s_i,t)=\sum _{k=1}^n\left({\varphi }_{{\alpha }_k}(s_i)\delta a_{i_k}(t)+{\varphi }_{{\beta }_k}(s_i)\delta b_{i_k}(t)\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(29b)

$$\begin{array}{ccc}\multicolumn{1}{c}{Y_i(s_i,t)=\sum _{k=1}^n\left({\varphi }_{{\alpha }_k}(s_i)c_{i_k}(t)+{\varphi }_{{\beta }_k}(s_i)d_{i_k}(t)\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(29c)

$$\begin{array}{ccc}\multicolumn{1}{c}{\delta Y_i(s_i,t)=\sum _{k=1}^n\left({\varphi }_{{\alpha }_k}(s_i)\delta c_{i_k}(t)+{\varphi }_{{\beta }_k}(s_i)\delta d_{i_k}(t)\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(29d)

where ${\varphi }_{{\alpha }_1},\dots ,{\varphi }_{{\alpha }_n}$ and ${\varphi }_{{\beta }_1},\dots ,{\varphi }_{{\beta }_n}$ are Lagrange polynomials in $H^2$ , satisfying the conditions ${\int }_0^{\ell }{({\varphi }_{{\alpha }_j}^{″})}^2\mathrm{ds}<\infty$ and ${\int }_0^{\ell }{({\varphi }_{{\beta }_j}^{″})}^2\mathrm{ds}<\infty$ for all $j=1,\dots ,n$ , so that integrals in (28) involving second derivatives $X_i^{″}$ and $Y_i^{″}$ are well defined [32]. Furthermore, the coefficient sets $a_{i_1},\dots ,a_{i_n}$ , $b_{i_1},\dots ,b_{i_n}$ , $c_{i_1},\dots ,c_{i_n}$ , and $d_{i_1},\dots ,d_{i_n}$ are the unknown time-varying amplitudes. The trial solutions evaluated at the node locations give $X_i(s_{i_k},t)=a_{i_k}(t)$ , $X_i^{\prime }(s_{i_k},t)=b_{i_k}(t)$ , $Y_i(s_{i_k},t)=c_{i_k}(t)$ , and $Y_i^{\prime }(s_{i_k},t)=d_{i_k}(t)$ due to the following interpolating properties of the Lagrange polynomials [32]

$${\varphi }_{{\alpha }_j}(s_{i_k})=\{\begin{array}{cc}\multicolumn{1}{c}{1}& \multicolumn{1}{c}{k=j}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{k\ne j}\end{array},{\varphi }_{{\alpha }_j}^{\prime }(s_{i_k})=0;j,k=1,\dots ,n$$

(30a)

$${\varphi }_{{\beta }_j}^{\prime }(s_{i_k})=\{\begin{array}{cc}\multicolumn{1}{c}{1}& \multicolumn{1}{c}{k=j}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{k\ne j}\end{array},{\varphi }_{{\beta }_j}(s_{i_k})=0;4\mathrm{ptj},k=1,\dots ,n.$$

(30b)

Note that the set of Lagrange polynomials used for both trial solutions $X_i$ and $Y_i$ is the same since they have to satisfy the same continuity requirements. The Lagrange polynomials are collected into ${\mathit{\varphi }}_{\alpha }=[{\varphi }_{{\alpha }_1},\dots ,{\varphi }_{{\alpha }_n}]$ , ${\mathit{\varphi }}_{\beta }=[{\varphi }_{{\beta }_1},\dots ,{\varphi }_{{\beta }_n}]$ , and $\mathit{\varphi }=[{\mathit{\varphi }}_{\alpha },{\mathit{\varphi }}_{\beta }]$ . Similarly, the unknown time-varying coefficients are collected into ${\mathit{a}}_i=[a_{i_1},\dots ,a_{i_n}]$ , ${\mathit{b}}_i=[b_{i_1},\dots ,b_{i_n}]$ , ${\mathit{c}}_i=[c_{i_1},\dots ,c_{i_n}]$ , ${\mathit{d}}_i=[d_{i_1},\dots ,d_{i_n}]$ , ${\mathit{q}}_{X_i}=[{\mathit{a}}_i,{\mathit{b}}_i]$ , ${\mathit{q}}_{Y_i}=[{\mathit{c}}_i,{\mathit{d}}_i]$ , and ${\mathit{q}}_i=[{\mathit{q}}_{X_i},{\mathit{q}}_{Y_i}]$ . Therefore, the trial functions and the test functions can be compactly compactly rewritten as:

$$X_i(s_i,t)=\mathit{\varphi }{\mathit{q}}_{X_i}^T$$

(31a)

$$\delta X_i(s_i,t)=\delta {\mathit{q}}_{X_i}{\mathit{\varphi }}^T$$

(31b)

$$Y_i(s_i,t)=\mathit{\varphi }{\mathit{q}}_{Y_i}^T$$

(31c)

$$\delta Y_i(s_i,t)=\delta {\mathit{q}}_{Y_i}{\mathit{\varphi }}^T.$$

(31d)

Substituting the projections (31) in the weak form (28), we obtain the finite element discrete form which holds for all virtual displacements $\delta {\mathit{q}}_i$

$$\rho {\mathit{M}}_i{\stackrel{\mathrm{··}}{\mathit{q}}}_i+\zeta {\mathit{M}}_i{\stackrel{·}{\mathit{q}}}_i+{\mathit{N}}_i{\mathit{q}}_i-{\mathit{U}}_i=0$$

(32)

where ${\mathit{q}}_i$ is a list of length $4n$ collecting the unknown coefficients, ${\mathit{M}}_i$ is the mass matrix of dimensions $4n\times 4n$ , ${\mathit{N}}_i$ is a $4n\times 4n$ nonlinear matrix that accounts for nonlinearities in the dynamics of the leg, and ${\mathit{U}}_i$ is the input vector of length $4n$ that depends on the external forces and moments acting on the leg. Further details about the FE matrices can be found in Appendix 1.

2.5. Inextensibility constraint

The constraint is enforced using the penalty method formulated by Bayo et al. [28], where a fictitious mass–spring–damper system based on the constraints is added to the system’s Lagrangian. We first introduce the projected constraint along the finite element basis functions by substituting (31) into (9). For the $i\mathrm{th}$ filament, we obtain the following:

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathcal{V}}_i=}& \multicolumn{1}{c}{\underset{0}{\overset{\ell }{\int }}\left({\mathit{r}}_i^{\prime }.{\mathit{r}}_i^{\prime }-1\right)d{s}_i=\underset{0}{\overset{\ell }{\int }}(X_i^{\prime 2}+Y_i^{\prime 2}-1)d{s}_i}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{={\mathit{q}}_i\left[\begin{array}{cc}\multicolumn{1}{c}{{\int }_0^{\ell }\left({\mathit{\varphi }}^{\prime T}{\mathit{\varphi }}^{\prime }\right)d{s}_i}& \multicolumn{1}{c}{\mathbf{0}}\\ \multicolumn{1}{c}{\mathbf{0}}& \multicolumn{1}{c}{{\int }_0^{\ell }\left({\mathit{\varphi }}^{\prime T}{\mathit{\varphi }}^{\prime }\right)d{s}_i}\end{array}\right]{\mathit{q}}_i^T-\ell .}\end{array}$$

(33)

Since the constraint ${\mathcal{V}}_i({\mathit{q}}_i)$ is holonomic in the finite element nodal coordinates ${\mathit{q}}_i$ , the penalty system adds the following term to the discretized system’s dynamics [28]:

$${\mathit{L}}_i=\lambda \left(\frac{d^2{\mathcal{V}}_i}{d{t}^2}+2\mathrm{\eta \mu }\frac{d{\mathcal{V}}_i}{\mathrm{dt}}+{\eta }^2{\mathcal{V}}_i\right){\frac{\partial {\mathcal{V}}_i}{\partial {\mathit{q}}_i}}^T$$

(34)

where λ, η, and μ are penalty parameters, with λ being a large constant, and η and μ representing the natural frequency and damping ratio of a fictitious penalty oscillator that enforces the constraint. After adding the penalty term to the legs’ dynamics, the finite element discrete form becomes

$$\rho {\mathit{M}}_i{\stackrel{\mathrm{··}}{\mathit{q}}}_i+\zeta {\mathit{M}}_i{\stackrel{·}{\mathit{q}}}_i+{\mathit{N}}_i{\mathit{q}}_i+{\mathit{L}}_i-{\mathit{U}}_i=0.$$

(35)

In total, there are $4n\times N$ equations describing the dynamics of N legs.

2.6. Body dynamics

Since the body of the robot is a planar rigid object with generalized coordinates ${\mathit{q}}_b(t)=(X_b(t),Y_b(t),\theta (t))$ , the kinetic energy is

$${\mathcal{K}}_b=\frac{1}{2}m\left({Ẋ}_b^2+{Ẏ}_b^2\right)+\frac{1}{2}{\mathcal{I}}_b{\stackrel{·}{\theta }}^2$$

(36)

where

$${\mathcal{I}}_b=\frac{1}{12}m\left(\left({\tau }_L^2+{\tau }_B^2\right)-\left({({\tau }_L-2{\tau }_T)}^2+{({\tau }_B-{\tau }_T)}^2\right)\right)$$

(37)

is the moment of inertia of the payload about its center of mass, with ${\tau }_L$ , ${\tau }_B$ , and ${\tau }_T$ being the length, breadth, and thickness as in Figure 1. The gravitational potential energy is given by

$${\mathcal{E}}_b=\mathrm{mg}\left(X_t-X_b\right).$$

(38)

A Rayleigh dissipation function

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{R}}_b=\frac{\zeta }{2}{\stackrel{·}{\mathit{q}}}_b\cdot {\stackrel{·}{\mathit{q}}}_b=\frac{\zeta }{2}({Ẋ}_b^2+{Ẏ}_b^2+{\ell }^2{\stackrel{·}{\theta }}^2)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(39)

is used also for the payload to model a stabilizing external linear viscous damping. Here, ζ is assumed to be the same dissipation constitutive parameter used for the legs. Hence, the generalized damping force vector acting on the body is given by

$${\mathit{Q}}_{D_b}=-\frac{\partial {\mathcal{R}}_b}{\partial {\stackrel{·}{\mathit{q}}}_b}.$$

(40)

The spring reaction force ${\mathit{F}}_{S_{\mathrm{bod}{y}_i}}$ acting at point $B_i$ (see Figure 1) on the payload is given by

$${\mathit{F}}_{S_{\mathrm{bod}{y}_i}}=\{\begin{array}{cc}\multicolumn{1}{c}{\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_1}}{{\gamma }_{{\alpha }_1}}({\mathit{r}}_{B_1}(t)-{\mathit{r}}_1(0,t))}& \multicolumn{1}{c}{\mathrm{for}i=1}\\ \multicolumn{1}{c}{\frac{\beta \mathrm{\Delta }{x}_{{\beta }_{i-1}}}{{\gamma }_{{\beta }_{i-1}}}({\mathit{r}}_{B_i}(t)-{\mathit{r}}_{i-1}(0,t))}& \multicolumn{1}{c}{\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathrm{for}2\le i}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{i\le N+1}\\ \multicolumn{1}{c}{}\end{array}}\\ \multicolumn{1}{c}{\frac{\alpha \mathrm{\Delta }{x}_{{\alpha }_{N+1}}}{{\gamma }_{{\alpha }_{N+1}}}({\mathit{r}}_{B_{N+2}}(t)-{\mathit{r}}_N(0,t))}& \multicolumn{1}{c}{\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathrm{for}}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{i=N+2.}\end{array}}\end{array}$$

(41)

The generalized spring force is

$${\mathit{Q}}_{S_b}=\sum _{i=1}^{N+2}{\mathit{F}}_{{S_b}_i}.\frac{\partial {\mathit{r}}_{B_i}}{\partial {\mathit{q}}_b}.$$

(42)

We consider the possibility of contact between the body and the terrain at points $P_1$ and $P_2$ , as shown in Figure 1, since at least one of these two points is the closest to the flat terrain. By using the body’s generalized coordinates, the positions in the global frame are:

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{r}}_{P_1}(t)=}& \multicolumn{1}{c}{\left[\begin{array}{c}\multicolumn{1}{c}{X_b(t)}\\ \multicolumn{1}{c}{Y_b(t)}\end{array}\right]}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{+\left[\begin{array}{cc}\multicolumn{1}{c}{\cos \theta (t)}& \multicolumn{1}{c}{-\sin \theta (t)}\\ \multicolumn{1}{c}{\sin \theta (t)}& \multicolumn{1}{c}{\cos \theta (t)}\end{array}\right]\left[\begin{array}{c}\multicolumn{1}{c}{-X_b(0)+{\tau }_C}\\ \multicolumn{1}{c}{\frac{N+1}{2}{L}_{\alpha }+{\tau }_B}\end{array}\right]}\end{array}$$

(43a)

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathit{r}}_{P_2}(t)=}& \multicolumn{1}{c}{\left[\begin{array}{c}\multicolumn{1}{c}{X_b(t)}\\ \multicolumn{1}{c}{Y_b(t)}\end{array}\right]}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{+\left[\begin{array}{cc}\multicolumn{1}{c}{\cos \theta (t)}& \multicolumn{1}{c}{-\sin \theta (t)}\\ \multicolumn{1}{c}{\sin \theta (t)}& \multicolumn{1}{c}{\cos \theta (t)}\end{array}\right]\left[\begin{array}{c}\multicolumn{1}{c}{-X_b(0)+{\tau }_C}\\ \multicolumn{1}{c}{-\frac{N+1}{2}{L}_{\alpha }-{\tau }_B}\end{array}\right].}\end{array}$$

(43b)

At the contact point, the action of the terrain on the body is then modeled by the following penalty functions [30]:

$$F_{N_{P_{1_x}}}(t)=\{\begin{array}{cc}\multicolumn{1}{c}{ϵ(X_t-r_{P_{1_x}}(t))}& \multicolumn{1}{c}{\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathrm{if}{r}_{P_{1_x}}(t)\ge X_t}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathrm{and}\theta (t)\le 0}\end{array}}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{\mathrm{otherwise}}\end{array}$$

(44a)

$$F_{N_{P_{2_x}}}(t)=\{\begin{array}{cc}\multicolumn{1}{c}{ϵ(X_t-r_{P_{2_x}}(t))}& \multicolumn{1}{c}{\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathrm{if}{r}_{P_{2_x}}(t)\ge X_t}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathrm{and}\theta (t)>0}\end{array}}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{\mathrm{otherwise}}\end{array}$$

(44b)

where $r_{P_{1_x}}$ and $r_{P_{2_x}}$ are the $\widehat{\mathit{i}}$ components of the vectors ${\mathit{r}}_{P_1}$ and ${\mathit{r}}_{P_2}$ , respectively. The penalty functions $F_{N_{P_{1_x}}}$ and $F_{N_{P_{2_x}}}$ are the $\widehat{\mathit{i}}$ components of the force vectors ${\mathit{F}}_{N_{P_1}}=(F_{N_{P_{1_x}}},0)$ and ${\mathit{F}}_{N_{P_2}}=(F_{N_{P_{2_x}}},0)$ that act on points $P_1$ and $P_2$ , respectively. The generalized normal force vector is

$${\mathit{Q}}_{N_b}={\mathit{F}}_{N_{P_1}}\cdot \frac{\partial {\mathit{r}}_{P_1}}{\partial {\mathit{q}}_b}+{\mathit{F}}_{N_{P_2}}\cdot \frac{\partial {\mathit{r}}_{P_2}}{\partial {\mathit{q}}_b}.$$

(45)

The vector of generalized forces acting on the robot’s body is

$${\mathit{Q}}_b={\mathit{Q}}_{D_b}+{\mathit{Q}}_{S_b}+{\mathit{Q}}_{N_b}.$$

(46)

Hence, the Euler–Lagrange equations describing the dynamics of the payload are

$$\frac{d}{\mathrm{dt}}\left(\frac{\partial {\mathcal{L}}_b}{\partial {\stackrel{·}{\mathit{q}}}_b}\right)-\frac{\partial {\mathcal{L}}_b}{\partial {\mathit{q}}_b}={\mathit{Q}}_b$$

(47)

which results in three coupled ordinary differential equations (ODEs) describing the evolution of the coordinates $X_b(t)$ , $Y_b(t)$ , and $\theta (t)$ , coupled with the legs’ dynamics via the spring reaction forces. In total, the discretized robot’s dynamics is described by $4\mathrm{nN}+3$ ODEs.

3. Simulation results

The set of simulations presented in this section aims to illustrate the main features of the terrestrial locomotion system presented above. Specifically, results show that metachronal coordination induced by external actuation can be used for locomotion, demonstrating that the model is suitable for autonomous walking robot design.

The full set of material and geometric parameters used in the simulation is listed in Table 1.

Table 1.

List of simulation parameters.

Parameter	Value
Length ℓ of each leg	$1\mathrm{m}$
Distance a between filaments	$0.05\mathrm{m}$
Linear mass density ρ of each leg	$1\mathrm{kg}{m}^{-1}$
Bending rigidity $K_M$ of each filament	$50\mathrm{kg}{m}^3{s}^{-2}$
Mass m of the robot’s body	$0.01\mathrm{kg}$
Damping coefficient ζ	$5\mathrm{N}s{m}^{-1}$
Spring constants α and β	$500\mathrm{N}{m}^{-1}$
Initial length $L_{\alpha }$	$1\mathrm{m}$
Initial length $L_{\beta }$	$0.5\mathrm{m}$
Thickness of the robot’s body ${\tau }_T$	$1\mathrm{m}$
Clearance of the robot’s body ${\tau }_C$	$0.2\mathrm{m}$
Value of $X_t$	$1\mathrm{m}$
Frequency of rotation of the legs ω	$0.5{\mathrm{s}}^{-1}$
Proportional gain $K_p$	500
Penalty constant ϵ	1000
Penalty constant ν	200
Penalty constant λ	$10,000$
Natural frequency η	175
Damping ratio μ	50

In the initial configuration, contiguous legs are phase shifted by $90\circ$ with legs 1, 5, and 9 in contact with the terrain. The distributed actuation forces of the elastic filaments are set to zero to isolate the effect of actuating moments applied at the legs’ bases. A relatively high leg-bending rigidity $K_M$ is assigned to facilitate the locomotion by transferring the terrain reactions to the robot’s body. Each leg is discretized into five uniformly spaced finite element nodes in the range $[0,\ell ]$ . The finite element basis functions and their derivatives used for the legs’ spatial discretization are plotted in Figure 4.

Figure 4.

Plot of finite element basis functions: (a) plot of ${\mathit{\varphi }}_{\alpha }$ , (b) plot of ${\mathit{\varphi }}_{\beta }$ , (c) plot of ${\mathit{\varphi }}_{\alpha }^{\prime }$ and (d) plot of ${\mathit{\varphi }}_{\beta }^{\prime }$ .

The system of equations is generated and solved using a code written in Mathematica 13 with simulation time $t_f=50\mathrm{s}$ . The in-built Mathematica solver, NDSolve, is used to find the solutions to the system of ODEs (35) using local convergence control and adaptive time steps.

Figure 5 shows snapshots of the robot configuration at different instants along with the position of the robot and the envelopes of the tips of the legs. The simulation illustrates the robot’s capability to successfully traverse flat terrain through the use of sustained metachronal coordination between the legs. The legs bend when in contact with the terrain, due to the flexibility of the adopted filament model. The envelope, obtained by joining the leg tips with straight lines, resembles a forward propagating wave due to the imposed coordinated motion of the legs. For this set of simulations, the springs’ stiffnesses were tuned to maintain a smooth gait cycle for the robot. The springs’ stiffness can be determined by trade-off considerations between coupling efficiency between legs, total actuation effort, and oscillatory motion of the body.

Figure 5.

Snapshots of the animation of the robot at different time instants: (a) $t=0s$ , (b) $t=10s$ , (c) $t=20s$ , (d) $t=30s$ , (e) $t=40s$ and (f) $t=50s$ .

The plot of the time histories $\psi (0,t)$ of the orientations of the legs’ bases is depicted in Figure 6. The graph shows that the implemented actuating moments propel the robot forward while maintaining the desired $90\circ$ phase difference between adjacent legs.

Figure 6.

Evolution of the basal angles.

The periodic nature of the motion is further illustrated by plotting the coordinates of the tips of the filaments in Figure 7. The emergent periodic patterns observed in the time evolution plots related to the rotation and displacement of the legs can be attributed to the metachronal coordination between filaments imposed at the legs’ bases.

Figure 7.

Time histories of the tip displacements of the legs: (a) X-coordinate of the tips. (b) Y-coordinate of the tips.

The inextensibility condition imposed via the penalty function method discussed above is plotted in Figure 8 for all the filaments. The tuned parameters for the penalty function result in errors of order of magnitude of $10^{-7}$ throughout the simulation; this is negligible compared to the filament length of $1\mathrm{m}$ . This confirms the effectiveness of the penalty function and the tuned parameters in enforcing the inextensibility constraint.

Figure 8.

Evolution of the inextensibility error of the filaments.

4. Conclusion

We presented a mechanical model for the locomotion system of a class of multi-legged terrestrial robots. The locomotion system consists of an array of flexible filaments, coupled at their bases through connecting elastic elements. This structure is inspired by axonemal filaments forming cilia and flagella, which among other functions are associated with motility of small biological organisms. Different coordinated motions and beating patterns can emerge. For terrestrial locomotion, metachronal coordination is of particular relevance, as it is the dominant locomotion mechanism used by a variety of terrestrial organisms such as millipedes and centipedes. We demonstrated in simulation that the evolution of the system results in locomotion of a payload by inducing metachronal coordination on the legs through ad hoc actuation consisting of moments applied at the legs’ bases.

This work is a first step towards the development of an autonomous robotic system with robust locomotion in unstructured environments. Several open questions are the subject of ongoing and future research. These include the prediction of the parameter configurations corresponding to the emergence of coordination patterns for locomotion. In addition, the influence of the legs’ coupling constitutive model needs to be investigated in the context of minimization of actuation energy, to determine minimal actuation conditions compatible with sustained coordinated beating patterns in the legs.