The softness distribution index: towards the creation of guidelines for the modeling of soft-bodied robots

Abstract

Modeling soft robots is not an easy task owing to their highly nonlinear mechanical behavior. So far, several researchers have tackled the problem using different approaches, each having advantages and drawbacks in terms of accuracy, ease of implementation, and computational burden. The soft robotics community is currently working to develop a unified framework for modeling. Our contribution in this direction consists of a novel dimensionless quantity that we call the softness distribution index (SDI). The SDI for a given soft body is computed based on the distribution of its structural properties. We show that the index can serve as a tool in the choice of a modeling technique among multiple approaches suggested in literature. At the moment, the investigation is limited to bodies performing planar bending. The aim of this work is twofold: (i) to highlight the importance of the distribution of the geometrical and material properties of a soft robotic link/body throughout its structure; and (ii) to demonstrate that a classification based on this distribution provides guidelines for the modeling.

Keywords

Modeling of soft robots softness stiffness soft-bodied robots soft robots

1. Introduction

Soft-bodied robots are often counterposed to rigid-linked robots. The former have gained huge popularity in recent years; the field is relatively young and since its appearance it has posed many challenges to researchers (Laschi et al., 2016; Lipson, 2014; Rus and Tolley, 2015). The latter belong to an older, deeply investigated field, in which researchers have solved major problems related to design and control; as a result, rigid-linked robots have found several industrial applications, in which both serial and parallel manipulators are used to perform various tasks, such as pick-and-place operations and handling of heavy metal sheets, at high speed and with high accuracy and precision.

Rigid-linked robots are usually provided with variable stiffness joints and actuators (Wolf et al., 2016), which improve the performances at the cost of increased complexity in the control. In these systems, named articulated soft robots, the joint stiffness is treated as a lumped parameter (LP), as explained by Albu-Schäffer and Bicchi (2016). In addition to the stiffness of the joints, the stiffness of the robotic links is taken into account in the study of elastodynamics: the eigenfrequencies of a manipulator depend on the structural properties of all its parts, and on how such parts are interconnected. Correct estimation and appropriate modeling of the stiffness of a robotic system are crucial to avoid undesired vibrations, which may jeopardize the performances in operative conditions. Therefore, the topic has been addressed by several authors in their works (to mention a very few of them, Pashkevich et al. (2009), Briot et al. (2009), Cammarata (2012), Zhang et al. (2015), Germain et al. (2015), and Rognant et al. (2010)). If the modeling of the stiffness plays such a key role for the so-called rigid-linked robots, it is easy to figure out that its importance is even greater for soft-bodied robots: in fact, they rely upon their low stiffness to perform a task, adapting their shape to the surrounding environment, often elegantly.

Unfortunately, the beauty of soft robots is also the beast that researchers have to fight. Different approaches have been proposed to tackle the problem of modeling soft-bodied robots. Among the most popular, we find the use of the finite element method (FEM). In fact, some authors use finite element commercial software to evaluate the mechanical response of their soft robots and actuators, to perform structural optimization or to develop their optimization method based on finite element simulations (e.g., Connolly et al., 2015; Elsayed et al., 2014; Moseley et al., 2016; Polygerinos et al., 2013; Runge et al., 2017; Suzumori et al., 2007). When using FEM, the stiffness appears in the so-called stiffness matrix, the size of which depends on the number and the kind of elements used. The computation of the inverse of this matrix is often a bottleneck: finite elements are indeed a powerful tool, but they come at a high computational cost. For this reason, other authors have implemented their own finite element code to make it suitable for real-time simulations; such huge work has required the efforts of a team and is part of a bigger framework (SOFA) presented in literature some years ago (Faure et al., 2012), followed by papers that expanded the work and prove its validity (such as Duriez, 2013; Ficuciello et al., 2018; Largilliere et al., 2015).

Other researchers based their approach on the Cosserat rod theory (Renda et al., 2017, 2014, 2018). In these works, the authors modeled the stiffness of soft bodies by means of a diagonal stiffness matrix; in some case, they discretize the structure of the body along its longitudinal axis and write a stiffness matrix for each segment. Recently, Grazioso et al. (2018) proposed a geometrically exact model that combines Cosserat rod theory and finite elements; in this work, as well as in previous approaches (such as that of Grazioso et al., 2016), the authors account for the stiffness by the same diagonal matrix. These works concern the modeling of continuum arms, that is, slender structures that can be discretized as a series of rods or beams.

Different approaches have also been found. For instance, some recur to the classical beam theory (Shapiro et al., 2015) or they start from the computation of the elastic strain energy (Connolly et al., 2017) to model soft pneumatic actuators, or use the elastica theory (Armanini, 2018; Armanini et al., 2017; Zhou et al., 2015) to describe the behavior of soft links to be used for soft robotics applications. Finally, some researchers do not develop any structural model for their soft robots. In these cases, the approach is purely experimental and the design and the development of the soft-bodied system is based on intuitions and qualitative considerations, and improved by trial and error.

All the approaches recalled above have both advantages and limitations in terms of accuracy, required computational burden, ease of implementation, and applicability. Researchers agree that the field still lacks a unified framework for modeling (Renda et al., 2017; Trimmer et al., 2015), and, more in general, that the development of a common approach for modeling, design, and control will be crucial for the success of the field (Bao et al., 2018). Despite the relevant progress in recent years, modeling soft robots still remains a challenging task. The difficulties lie in their highly nonlinear mechanical behavior. A first source of nonlinearity is geometrical and comes from the fact that they typically undergo large deflections, and in some cases, large strains. A second, no less important, source may come from the materials often used to build the body of the robot: polymers are characterized by nonlinear stress–strain curves, whose trends are captured by well-established models (e.g., Arruda and Boyce, 1993; Mooney, 1940; Ogden, 1972). Owing to the use of these materials (in addition to colloidal and granular matter), Wang and Iida (2015) used the term soft-matter robotics to highlight the role played by the material properties. In our opinion, there is a further aspect to consider: two soft bodies can behave very differently in terms of deformed shape, even though they are built with the same material, have the same size and undergo large displacements of the same order of magnitude, under the same external load. Very roughly speaking, some soft links can be considered more hyper-redundant than others.

In this article, we do not propose any new modeling technique. Instead, we explain why and how we have taken a step back and analyzed the problem of modeling soft-bodied robots under a different light. The focus of our work is not only on how soft a robotic link is, but also, and in particular, on the distribution of its structural properties, and on the role that such distribution plays when a modeling technique must be chosen among others. To refer to the distribution in a quantitative way, we define what we call the SDI for planar bending. We show that this index, which is a complex number, allows classes of soft bodies to be associated with one or more modeling techniques involving different numbers of parameters (and, therefore, having different computational cost).

To better share our view with the reader, we present the work following the stream of ideas and considerations that have led us to the results here reported. Therefore, Section 2 provides a discussion about stiffness and softness, from which we have found motivation for our work. Section 3 introduces the SDI. In Section 4, we briefly recall the modeling techniques that we use in this work, and that we relate to ranges of value of the index. Section 5 reports the study performed on a class of bodies, modeled by the techniques listed in Section 4; to complete the investigation, other bodies with particular distribution of the structural properties are considered in Section 6. In Section 7 we report few examples that show how the SDI can be used when performing optimization of a structure performing bending under a set of applied loads. Section 8 contains a discussion concerning the usefulness of the proposed index and highlights possible applications. The limitations of the current work are also clearly stated, to not mislead the reader. Conclusions follow.

2. Stiffness and softness

Stiffness is a property of structures and it depends on both the material and the geometry. It is a physical quantity defined in relation to an applied load and the consequent deflection along the considered direction; it is indeed a relative quantity. By the term soft we refer to a body, robot, or structure, having low stiffness. In this article, therefore, we often use the term softness as complementary to stiffness. By the term rigid, instead, we refer to a body having much greater stiffness than that of the surrounding environment (as is well known, the definition of rigid body is purely the result of abstraction).

2.1. Stiff or soft?

Let us now consider two systems: one made of two rigid links connected by an elastic rotational joint, counterposed to a single soft link. The difference between the two is clear: in the former, the compliance is localized at the joint; in the latter, the compliance is distributed and the relative position between the two end sections depends on the position and orientation of all the intermediate sections. Therefore, the compliance of the former can be modeled by one parameter, treating the joint as a torsional spring; for the soft link, such approach generally turns out to be inapplicable. These two examples represent the extrema of the path depicted in Figure 1: it is possible to consider infinite systems between these two, with gradually more and more distributed compliance. Indeed, one may ask whether there is precise, sharp boundary between what is soft and what is rigid, along this path. We all distinguish soft robotic links from rigid links based on our intuition, but can we draw a sharp line that separates these two classes?

A similar question was posed by Ananthasuresh (2013) about differences and similarities between compliant mechanisms (CMs) and rigid-body mechanisms (RBMs). In this work, Ananthasuresh highlighted that CMs have features in common with both RBMs and stiff structures: because they are continuum systems they can be seen as structures; at the same time, they do not differ from RBMs in terms of function as well as of modeling approach. In fact, the pseudo rigid-body model (PRBM) has been proposed (Howell and Midha, 1995; Howell et al., 1996), by which the continuum structure of the CM is transformed into its equivalent RBM with lumped springs. Later, Cao et al. (2015) proposed a unified synthesis approach without any prescription of the type of mechanisms obtained.

2.2. How to choose the modeling technique

CMs usually consist of thin leafsprings with uniform cross section or in relatively rigid rods connected by flexural pivots that provide mobility to the mechanism. As explained by Ananthasuresh (2013), in some cases finite elements are needed to model the mechanical response of the CMs to applied loads, especially if an accurate analysis of the mechanical stress must be performed. However, the PRBM well describes the deformation of most of the CMs and other CMs are studied by a spring–mass–lever (SML) model. An important remark is that models such as the PRBM and the SML are useful for synthesis purposes (e.g., Howell and Midha, 1994, 1996; Pucheta and Cardona, 2010). Therefore, the choice of a modeling technique is not merely a matter of convenience in terms of computational burden and accuracy; it also helps in synthesis, optimization, and to achieve a deeper understanding of how the system works. There is no reason to believe that these considerations do not hold for soft-bodied robots too.

In general, the bodies that constitute a soft robot can have irregular geometry and can be made of multiple materials. To provide a few examples, the actuator by Elsayed et al. (2014) and the soft tentacle by Martinez et al. (2013) are made of more than one polymeric rubber; the soft finger developed by Manti et al. (2015) has irregular geometry and is made of two materials; the robot that appears in the work by Largilliere et al. (2015) is a parallel machine whose soft limbs have variable cross section along their length. All these systems are capable of large displacements, but they differ in that they exhibit different deflected configurations under the same applied load. Some distributions can be effectively described by models based on few parameters, as happens for some CMs; others require more computationally expensive techniques. As the deformation of a body strongly depends on both the structural properties and their distribution, it seems useful to consider not only how soft a system is, but also how the softness is distributed throughout it, to select a modeling technique that describes the deformed configuration with the required accuracy and at the lowest computational cost. Let us consider again that a compliant link with a flexural pivot (such as that in Figure 1(b)) can be transformed into one with uniformly distributed compliance (such as that in Figure 1(f)) by a continuous transformation of its structure; if we cannot draw the aforementioned line that separates soft bodies from rigid bodies, how can we assess to what extent the methods developed for rigid-linked robots and CMs can be used for soft robots? It is worth investigating this aspect, to try to draw this line and provide quantitative, rather than qualitative, answers.

Therefore, we need a quantity that takes into account the variation of the structural properties along a given direction in the system. In the next section, we define such quantity and we discuss its meaning. As a considerable number of soft-bodied robots (such as soft continuum manipulators, soft fingers) perform their tasks by bending, the definition is based on this deformation mode.

3. SDI for planar bending

As already stated, a soft robot body can have irregular geometry and can be made of multiple materials. Figure 2 represents a beam in which both the cross section and the material properties vary along the longitudinal axis. Denoting by $x$ the coordinate along the longitudinal axis and by $L$ the length of the beam, the so-called flexural rigidity of the beam is the function

k_{f} (x) = E (x) J (x)

(1)

where $E (x)$ and $J (x)$ are the Young’s modulus and the moment of inertia, respectively, computed at the generic cross section $x$ in the interval $[0, L]$ . From now on, the beam will be denoted by $B_{i}$ , and it will be considered an element of the set $B$ , containing all the bodies performing planar bending.

Provided that $k_{f} (x)$ is a smooth function, its first derivative can be computed as follows:

α (x) = \frac{d k_{f} (x)}{dx} = E (x) \frac{dJ (x)}{dx} + J (x) \frac{dE (x)}{dx}

(2)

We denote by $x_{M}$ and $x_{m}$ the coordinates at which $α (x)$ has global maximum $α_{M}$ and global minimum $α_{m}$ , respectively. If $x_{M}$ and/or $x_{m}$ are not unique, they must be taken in order to minimize the norm $| x_{M} - x_{m} |$ . Based on Equations (1) and (2), we define the SDI for planar bending as the complex number

SDI = χ + i ψ

(3)

where it is

χ = {\begin{matrix} \begin{matrix} \frac{α_{M} - α_{m}}{\max {| α_{M} |, | α_{m} |}} \cdot \frac{L}{x_{M} - x_{m}} & if α (x) \neq const . \\ 0 & if α (x) = const . \end{matrix} \end{matrix}

(4)

and

ψ = \frac{k_{f_{M}} - k_{f_{m}}}{k_{f_{M}}} \cdot \frac{d_{M}}{L}

(5)

in which $d_{M}$ denotes the maximum dimension of the largest cross section (for example, the maximum diameter in the case of a circular cross section) and $k_{f_{M}}$ and $k_{f_{m}}$ denote the maximum and minimum value of $k_{f}$ , respectively.

The real part of the SDI accounts for the distribution of $k_{f}$ throughout the body: the first multiplying factor depends on the maximum and minimum increments of $k_{f}$ along $x$ , whereas the second depends on the coordinates at which these increments occur. Instead, the imaginary part accounts for $k_{f}$ and for the slenderness of the structure; by definition, it is always $k_{f_{M}} \geq k_{f_{m}}$ and, therefore, $ψ \geq 0$ .

The SDI can be computed for any body $B_{i}$ of the set $B$ ; it can be seen as a function $SDI : B \to C$ . This function is not injective (different bodies can be characterized by the same SDI) and not surjective in $C$ (in fact, there is no body having $ψ < 0$ ). An important property of the SDI is that it does not depend on the orientation of the $x$ -axis. This turns out to be fundamental: to be meaningful, the index must depend only on the structural characteristics of the body, and not on the provided mathematical description. The proof of the invariance under coordinate transformation is given in Appendix A.

From now on, the objective of the work will be to find a correlation between a body $B_{i}$ and a modeling technique $T_{ι}$ by means of the SDI, as outlined in Figure 3. In this article, we limit our investigation to four techniques that use different number of parameters.

4. Modeling techniques

As already stated, we consider four different approaches/techniques to study the static mechanical response of bodies:

finite element analysis (FEA);

nonlinear matrix structural analysis (NMSA);

elastica approach (EA);

LPs model (LPM).

In the following, we briefly recall these techniques and we explain why they have been selected among others for the purpose of this work. To avoid cumbersome mathematical expressions and excessively long discussion, the details are reported in the appendices at the end of the article, when needed.

4.1. Finite elements

In structural analysis, finite elements (FE) allow mechanical systems to be studied, however complex; that is, having irregular geometry, multiple materials, under any boundary conditions included distributed loads and occurrence of contact. They have been widely used in the last 60 years for a variety of engineering problems in which analytical models would be extremely difficult to handle, if not unworkable. The effectiveness and the consequent popularity of the method has led to the implementation of both commercial and open-source software, employed worldwide for analysis and optimization. To summarize, the method is based on the discretization of the investigated system in a number of elements having finite dimensions; the result of the discretization is called a mesh. For the generated mesh, a set of algebraic equations are written that correlates the vector of nodal displacements to the vector of nodal forces by means of the stiffness matrix.

A review on the topic goes far beyond the scope of this paper; we limit the discussion to the description of the settings adopted using a commercially available software (ABAQUS®). In all the performed simulations, the type and the number of elements has been selected in order to obtain reliable results. For each simulation, we have verified that no excessive distortion of elements had occurred and that a further increase of the fineness of the mesh did not lead to any significant change in the results obtained. Moreover, for the purpose of this work, all the simulations account for large deflections.

It is necessary to remark that in this work we do not implement a FEM model that needs to be validated, nor do we suggest the use of FEM as a general approach; our aim is to compare the results of FEA with those obtained by other techniques that involve far fewer parameters. From now on, therefore, we consider the deformation computed by FEA as the most accurate, but at the highest computational cost.

4.2. NMSA

Matrix structural analysis (MSA) can be considered the mother of the well-established FEM. For the interesting story of this technique and the way it finally turned into FEM, the reader can refer to Felippa (2001). Applications of the technique can be found in works addressing the elastodynamics of conventional rigid-linked manipulators (such as those cited in Section 1).

As stated previously, for the purposes of this work, we analyze soft bodies performing planar bending. Similarly to FEM, the body is discretized into a certain number of elements having finite dimensions; in this work, we use two-node beam elements, each node having three degrees of freedom (DOFs) on the plane (the translations along two orthogonal directions and the rotation along the axis orthogonal to the plane). As a result, each element is characterized by a $6 \times 6$ stiffness matrix. The size of the matrix modeling the stiffness of the entire structure depends on the total number of nodes of the mesh, as in FEM; however, the stiffness matrix in MSA is usually considerably smaller than that in FEM, owing to the low number of nodes in the mesh. Focusing on the modeling of the stiffness, the method shares similarities with both FEM and the Cosserat approach used for continuum manipulators: in fact, the stiffness matrices in MSA are formally identical to those written based on FEM, but the number of parameters involved is of the same order of magnitude of that appearing in Cosserat-based approaches (see, for instance, Renda et al., 2017), owing to the coarseness of the discretization.

A complete and detailed explanation of the method can be found in McGuire and Gallagher (1979). To model the dynamics of a system, a mass and a damping matrix are also defined; however, our investigation is limited to the static response of bodies. Therefore, when using MSA we will write the stiffness matrix only.

If the assumption of small displacements is valid, nodal displacements and nodal forces are related by the very well-known Hooke’s law, that is, by a linear relation. To account for large displacements, NMSA must be carried out. For this reason, in this work we implement a linear-incremental MSA: the external load is applied in incremental way in a predefined number of steps. Although each step is linear, the overall computation is nonlinear, because at each step the stiffness matrix is updated based on the displacements computed in the previous step. A detailed mathematical description is reported in Appendix B.

4.3. EA

Some authors have investigated the use of the EA to model the mechanical response of soft links or arms (Armanini, 2018), suggesting potential application in the field of soft robotics. Unlike FEM and MSA, which are based on algebraic equations written in matrix form, this approach involves differential equations. A brief discussion on the topic can be found, for instance, in Howell (2001).

The theory allows to account for large deflections of a beam subject to an applied load; the relation between the curvature and the bending moment is described by the equation

\frac{1}{ρ (s)} = \frac{d θ}{ds} = \frac{M (s)}{E (s) J (s)}

(6)

in which $s$ denotes the curvilinear coordinate along the axis of the beam. By $ρ$ and $θ$ we denote the radius of curvature and the rotation of the cross section located at $s$ , respectively, owing to the bending moment $M (s)$ . In Cartesian coordinates, $x$ is the axial coordinate and $y$ is the direction of the vertical deflection, thus the equation is written as

\frac{M (x)}{EJ} = \frac{\frac{d^{2} y}{d x^{2}}}{\sqrt{{(1 + {(\frac{dy}{dx})}^{2})}^{3}}}

(7)

For a beam undergoing small deflections, the derivative $\frac{dy}{dx}$ (that is, the rotation of the cross section) is small and can be neglected; therefore, the equation takes the form

\frac{M (x)}{EJ} = \frac{d^{2} y}{d x^{2}}

(8)

and it can be integrated twice over the length of the beam to compute the vertical deflection $y (x)$ under given boundary conditions. This equation is widely used for a variety of classical mechanical and civil engineering problems, in which the product $EJ$ is usually constant. It is worth to highlight that Equation (8) does not allow the displacement of the cross sections along $x$ to be computed; in fact, in the assumption of small vertical deflections, the horizontal deflections are neglected. This makes Equation (8) unsuitable for soft robotics applications, in which the systems typically undergo large deflections. Therefore, in this article, we only consider Equation (6). The Cartesian coordinates of the generic cross section are related to the curvilinear coordinate by the relations

\frac{dx}{ds} = \cos θ

(9)

\frac{dy}{ds} = \sin θ

(10)

The integration of Equation (6) may be a difficult task, depending on the loading and constraint conditions. In general, the equation is reworked, written in a form that contains elliptic integrals and solved by numerical techniques (see, for instance, Zhang and Chen, 2013).

4.4. LPs

LPs are used in a huge variety of domains other than robotics, such as the automotive industry, civil engineering, and biomechanics.

The advantage of using LPs is the low computational cost, especially compared with other techniques such as FEM. A lumped stiffness $k$ can be either constant or displacement-dependent; in any case, it consists of a single quantity that accounts for the distributed stiffness of one or more mechanical components. Lumped springs are successfully employed to model the behavior of CMs, as stated in Section 2.

However, in many cases LPs turn out to be insufficient to describe the deformed configuration of a system: for instance, the variable radius of curvature of a long, soft continuum arm that performs bending would be scarcely captured by few lumped springs.

In the following section, we model the planar bending of a body using few parameters, whose values are computed based on the structural properties.

5. A soft notched body made of uniform material

Keeping in mind the concept depicted in Figure 1, our investigation starts with a soft notched body made of one material having Young’s modulus $E$ . The body has length $L$ and rectangular cross section with constant width $b$ and variable height $h (x)$ , which ranges from a minimum $h_{m}$ to a maximum $h_{M}$ along the $x$ -axis, as represented in Figure 4:

\begin{matrix} h (x) = h_{M} + \frac{1}{2} (h_{M} - h_{m}) (\cos \frac{n π (x - \frac{L}{2} - \frac{L}{n})}{L} - 1) \cdot \\ H (x - \frac{L}{2} - \frac{L}{n}) H (\frac{L}{2} + \frac{L}{n} - x) \end{matrix}

(11)

in which $H (\cdot)$ denotes the Heaviside function. The parameter $n$ is used here to tune the length $L_{n}$ of the notch; it is, in fact, $L_{n} = 2 L / n$ . Increasing $n$ , a narrower notch is obtained; decreasing $n$ , the length of the notch tends to the length of the link, up to the limit condition $n = 2$ , in which it is $L_{n} = L$ . Figure 5 shows the geometries obtained for $n = 2, 4, 10$ . In this case, because the material properties do not change along the body, only the geometry plays a role in the distribution of the stiffness; the Young’s modulus reduces to a multiplying factor in the first derivative of $k_{f}$ :

α (x) = E \frac{dJ (x)}{dx}

(12)

Therefore, in this section we speak about the variation of $k_{f}$ or of the geometry along the $x$ -axis without distinction. Here and in the following, we set $E = 1$ MPa, a value in the range reported by Wang and Iida (2015) for soft-matter robotics.

By tuning the parameter $n$ we obtain more/less abrupt variations of the cross sections. Intuitively, we expect that a body characterized by a high value of $n$ behaves in a hinge-like way in response to an applied moment load; hence, the adoption of a LPM might provide accurate results at a very low computational cost. In contrast, a body with low $n$ has a more distributed softness. In this case, we expect that the use of a different technique is necessary to compute the relative displacement between the two end sections of the body with the same accuracy.

In this section, we compute the SDI for this class of bodies, we model the bodies by using the modeling techniques listed in Section 4 and, based on the results obtained, we derive the correlation between the SDI and the convenience of a modeling technique in terms of accuracy, computational burden and ease of implementation.

5.1. SDI for the notched body

The SDI as defined in Equations (3), (4), and (5) depends on $k_{f}$ , on its first derivative and on the aspect ratio of the body. We treat $n$ and the ratio $h_{m} / h_{M}$ as parameters which range in given intervals, to observe how the SDI varies depending on the geometry. The intervals and the values of the material and geometrical parameters adopted for this study are reported in Table 1.

Table 1.

Geometrical and material parameters of the notched body.

Parameter	Description	Value/interval
$n$	Tunes the length of the notch	$[2, 16]$
$L$	Length of the body	100 mm
$h_{M}$	Maximum height of the cross section	16 mm
$h_{m}$	Minimum height of the cross section	$[4, 14]$ mm
$b$	Width of the cross section	10 mm
$E$	Young’s modulus	1 MPa

To provide an example, Figure 6 shows the graph of the function $k_{f} (x)$ and the first derivative $α (x)$ for a notched body with $n = 10$ , $h_{m} = 8$ mm and other parameters as reported in Table 1. In this case, it is $x_{m} = 43.4$ mm, $x_{M} = 56.6$ mm, $α_{m} = - 534.2$ Nmm, and $α_{M} = 534.2$ Nmm; the maximum and minimum flexural rigidity are $k_{f_{M}} = 3.413 \times 10^{3}$ Nmm² and $k_{f_{m}} = 0.427 \times 10^{3}$ Nmm², respectively. Being $d_{M} = h_{M}$ , the computed SDI for this body is $SDI = 15.26 + i 0.14$ .

The values of $n$ and $h_{m} / h_{M}$ do not influence the trend of the functions $k_{f} (x)$ and $α (x)$ ; as far as it is $h_{m} < h_{M}$ , the trends are those shown in Figure 6. However, they greatly affect the values of the extrema of the functions and, hence, $χ$ and $ψ$ ; Figure 7 plots the functions $χ$ and $ψ$ in the intervals reported in Table 1.

We observe that $χ$ ranges between 2.86 and 29.63 in the considered domain and $ψ$ varies between 0.05 and 0.16. Combining the results shown in Figure 7, we gain a preliminary understanding of how the SDI maps the notched bodies to points belonging to the first quadrant of the complex plane; Figure 8 provides some examples of this map, showing nine bodies and their SDI on the complex plane. It can be seen that a body having nearly uniform cross section, body n.7, is characterized by lower values of $χ$ and $ψ$ compared with a body with a narrow and deep notch, body n.3.

5.2. Modeling

We study the deformation of this class of bodies owing to a constant external moment $\bar{M}$ by using the four modeling techniques listed in Section 4. The objective is to investigate the advantages and drawbacks of these techniques depending on the specific body, and to draw conclusions concerning the convenience of the considered approaches.

For the motivations reported in Section 4.1, we assume that the finite element technique provides the most accurate estimation of the deformation, among all the techniques; therefore, the results will be presented in terms of error with respect to those obtained by FEA. The aim of the work reported in this section is to assess whether it is possible to identify a particularly convenient modeling technique depending on the distribution of the softness in the notched body. We have analyzed a total of 21 bodies which differ in $n$ and $h_{m} / h_{M}$ : the former assumes the values in the set ${2, 4, 6, 8, 10, 12, 14}$ and the latter in ${0.25, 0.50, 0.75}$ . Before showing the results, we describe the implementation of the analyses. In all the cases, the system of coordinates is the same as introduced at the beginning of this section: $x$ denotes the abscissa along the longitudinal axis of the body and $y$ the direction of the vertical displacement. The body deflects in the plane $(x, y)$ . To obtain large displacement of the free end of the body, the applied moment $\bar{M}$ varies depending on the ratio $h_{m} / h_{M}$ . We have set the moment equal to 4.36 Nmm, 13.4 Nmm, and 39.3 Nmm for bodies having $h_{m} / h_{M}$ equal to 0.25, 0.5, and 0.75, respectively.

The reader should keep in mind that the investigation carried out here is limited to linear elastic material. Moreover, unlike FEM, the other methods do not account for shear effects or warping of the cross-section. Further discussion on these aspects is postponed to Section 8.

5.2.1. FEA

FEA is performed by means of the commercial software ABAQUS Standard®. In all the cases considered, the body is discretized in 5,712 hexahedral elements with reduced integration, for a total of 27,151 nodes. The kind of element selected has 20 nodes. We have selected hexahedral elements instead of tetrahedral, because in this case they allow for the creation of a structured mesh. By properly setting the element size, a subset of nodes can be located at the neutral plane of the body. This is necessary to extract results in terms of nodal displacements that can be compared with those obtained by the other modeling techniques adopted in this work, which are all relative to the neutral plane. All the simulations account for large deflections and have been performed using an iterative solver. The end moment load is gradually applied as a linear function of pseudo-time, in 100 steps. The material is considered linear elastic, with Young’s modulus $E = 1$ MPa and Poisson ratio $ν = 0.45$ . The applied moment varies depending on the body; it is taken such that a non-negligible horizontal displacement of the free end is obtained. As a general comment, we have noticed that for bodies with low $n$ and low $h_{m} / h_{M}$ it is easier for the simulation to converge, compared with the case in which a body with high $n$ and $h_{m} / h_{M}$ experiences the same displacement of the tip.

The results for each simulation are then exported and imported into the MATLAB® environment, in which the displacements of the nodes on the neutral plane are extracted and plotted, to be compared with the results obtained with the other modeling techniques.

5.2.2. NMSA

The structure is discretized into 6 beam elements having equal length $L_{e} = L / 6$ , numbering the nodes from 1 at $x = 0$ to 7 at $x = L$ . Since the body is made of one material, the Young’s modulus $E$ is the same for all elements. Adopting the same approach used by Renda et al. (2017), the area and moment of inertia of the cross section of the $i$ th element are calculated considering the average value of the function $h (x)$ in the interval $[(i - 1) L_{e}, i L_{e}]$ . Node 1 is fully constrained ( $u_{1} = 0$ ) and the end moment load is applied on node 7, as a constant increment $Δ M = \bar{M} / n_{s}$ . The number of steps is 300. The stiffness matrix is updated at each step, as explained in Appendix B.

5.2.3. EA

We study the deformation of a cantilever beam with variable cross section subject to end moment load $\bar{M}$ ; therefore, Equation (6) takes the form

\frac{d θ}{ds} = \frac{\bar{M}}{EJ (s)}

(13)

and the boundary conditions are as follows:

\begin{matrix} θ (0) = 0 \\ x (0) = 0 \\ y (0) = 0 \end{matrix}

(14)

We compute the deformed configuration of the beam by numerical integration of Equation (13) performed in MATLAB® using the built-in function integral. We discretize the interval $[0, L]$ into 200 linearly spaced points; we compute the rotation $θ (s_{i})$ for $i = 1, 2, \dots, 200$ as

θ (s_{i}) = \int_{0}^{s_{i}} \frac{\bar{M}}{EJ (\tilde{s})} d \tilde{s}

(15)

The Cartesian coordinates of the cross section $s_{i}$ are found by numerical integration of Equations (9) and (10).

5.2.4. LPs

The LPM uses the minimum number of parameters among the four considered techniques. In this work, we model the deformation of the notched body by means of the set of parameters defined as follows.

Let us denote by $x_{1} = L (\frac{1}{2} - \frac{1}{n})$ and $x_{2} = L (\frac{1}{2} + \frac{1}{n})$ the extrema of the notch (therefore, it is $| x_{2} - x_{1} | = L_{n}$ ). The cross section is uniform in the interval $[0, x_{1}]$ . The flexural stiffness of this part of the body is taken into account by the quantity

k_{U} = \frac{E J_{M}}{x_{1}}

(16)

where it is $J_{M} = \frac{{bh}_{M}^{3}}{12}$ . The quantity above relates the moment load and the rotation of the cross section at $x = x_{1}$ according to the classical beam theory. Similarly, the flexural stiffness of the notch is modeled by the parameter

k_{N} = \frac{E J_{m}}{l_{eq}}

(17)

The quantity $l_{eq}$ is the length given by the expression

l_{eq} = | x_{F} - x_{I} |

(18)

being $[x_{I}, x_{F}]$ the interval in which it is

h (x) \leq h_{m} + 0.267 (h_{M} - h_{m})

(19)

The multiplying factor 0.267 has been determined by imposing that a notch with $h_{m} / h_{M} = 0.25$ and $n = 16$ behaves like a flexural pivot, setting that the position of its free end section in the deformed configuration matches with that computed by FEA. It can be seen as a correction factor that accounts for the cosinusoidal profile of the notch (as discussed in more detail later, the choice of the numerical value of this factor is not trivial).

The stiffness coefficients defined previously are used to compute the coordinates of points $Q$ and $S$ in Figure 9, in which the body is represented as a broken line made of two segments. In fact, the angles

θ_{U} = \frac{\bar{M}}{k_{U}}

(20)

θ_{N} = \frac{\bar{M}}{k_{N}}

(21)

are the relative rotations between sections at points $P$ and $O$ , and between $R$ and $P$ , respectively; the relative rotation between $S$ and $R$ is equal to $θ_{U}$ as well, because the segments $\bar{OP}$ and $\bar{RS}$ have the same structural properties. Here, $P$ and $Q$ denote the extrema of the notch, whereas $Q$ is the midpoint.

As the vertical deflection of point $P$ due to the end moment $\bar{M}$ is

y_{P} = \frac{\bar{M} x_{1}^{2}}{2 E J_{U}}

(22)

the slope of the line passing through $O$ and $P$ can be calculated as $\arctan (y_{P} / x_{1}) = \bar{M} / 2 k_{U} = θ_{U} / 2$ . Analog considerations hold for segments $\bar{PQ}$ , $\bar{PR}$ , and $\bar{RS}$ . The three springs represented in Figure 9 have stiffness coefficients $k_{O} = k_{R} = 2 k_{U}$ and $k_{P} = 2 k_{N}$ . The coordinates of $Q$ and $S$ can be expressed as functions of $θ_{U}$ and $θ_{N}$ as follows:

\begin{matrix} x_{Q} = x_{1} \cos \frac{θ_{U}}{2} + \frac{L_{n}}{2} \cos (θ_{U} + \frac{θ_{N}}{4}) \\ y_{Q} = x_{1} \sin \frac{θ_{U}}{2} + \frac{L_{n}}{2} \sin (θ_{U} + \frac{θ_{N}}{4}) \end{matrix}

(23)

\begin{matrix} x_{S} = x_{R} + (L - x_{2}) \cos (θ_{U} + θ_{N} + \frac{θ_{U}}{2}) \\ y_{S} = y_{R} + (L - x_{2}) \sin (θ_{U} + θ_{N} + \frac{θ_{U}}{2}) \end{matrix}

(24)

being

\begin{matrix} x_{R} = x_{1} \cos \frac{θ_{U}}{2} + L_{n} \cos (θ_{U} + \frac{θ_{N}}{2}) \\ y_{R} = x_{1} \sin \frac{θ_{U}}{2} + L_{n} \sin (θ_{U} + \frac{θ_{N}}{2}) \end{matrix}

(25)

5.3. Results

The results obtained by using the four techniques are here reported and compared for the 21 bodies analyzed. Figure 10 shows the deformed configuration for the body with $n = 2$ and $h_{m} / h_{M} = 0.5$ as computed by FEA; from this configuration, only the nodes lying on the longitudinal axis of the body are extracted and plotted in Figure 11 together with the deformed axis computed by means of the other techniques. The procedure is repeated for the remaining 20 bodies.

It is evident from the plot that, for the considered body, all the techniques provide results that deviate from those given by FEA; however, the error in the position $p$ of the free end section largely varies from case to case. A discussion on the motivation for these discrepancies is given at the end of this section. At the moment, we only present the results for each case in terms of difference between the position $p_{FEA}^{(ij)}$ and the position $p_{ι}^{(ij)}$ computed by the technique $T_{ι}$ . The superscripts $i$ and $j$ here denote the values of $n$ and $h_{m} / h_{M}$ , respectively. The 2-norm of the displacement of the free end is given by

d_{FEA}^{(ij)} = | | p_{FEA}^{(ij)} - p_{0} | |

(26)

where $p_{0} = {L 0}^{T}$ is the position vector in the undeformed configuration. We compute the relative error for the technique $T_{ι}$ as

e_{ι}^{(ij)} = \frac{| | p_{ι}^{(ij)} - p_{FEA}^{(ij)} | |}{d_{FEA}^{(ij)}}

(27)

and $ι$ denotes EA, NMSA, or LPM.

The plots in Figure 12 show the errors relative to bodies with different $n$ and ratio $h_{m} / h_{M}$ equal to 0.25, 0.50, and 0.75, respectively. We observe that EA seems to provide the most accurate results in all the cases except for $h_{m} / h_{M} = 0.25$ . NMSA turns out to be more accurate than LPM for low values of $n$ and $h_{m} / h_{M} = 0.25, 0.50$ , and for any $n$ and $h_{m} / h_{M} = 0.75$ . Based on these results, it is possible to sketch lines to divide the complex plane in regions as shown in Figure 13. Each region includes a subset of bodies, identified by their SDI, and corresponds to the modeling technique that has provided the minimum error in the modeling. We observe that, as expected, the bodies with the narrower notch are better modeled by LPM than by NMSA; these bodies are characterized by high values of both $χ$ and $ψ$ .

We observe that the SDI, as defined in Section 3, effectively allows to identify regions of the complex plane that correspond to a modeling technique. In Figure 13, the full colored region indicates that the denoted technique has provided an error smaller than 10%; the other regions are relative to an error up to 20%. We also note that in some cases two techniques overlap. FEA can provide highly accurate results for any body, with the proper settings; hence, we can consider that it covers the entire space represented in the plot, and it is not displayed to make the figure more readable.

The LPM as implemented in this work has been developed to model hinge-like notched bodies: the lowest $e_{LPM}$ is obtained for notched bodies with $h_{m} / h_{M} = 0.25$ and high values of $n$ . The stiffness parameters adopted in the model are constant and they have been defined based on the classical linear beam theory. Therefore, as expected, they might not describe accurately the deformation of a body with gradually distributed compliance, which is better described by EA and NMSA; however, these techniques have their limitations as well. For instance, all the techniques as implemented here do not account for warping of the cross section, which is well captured by FEA, instead. As observed in the literature, EA is usually employed to model the deflection of slender, thin beams or leafsprings, having uniform cross section. For this kind of structure, warping is a negligible phenomenon; the technique produces accurate results in those cases. Consistently, bodies characterized by low values of $ψ$ (<0.15) are better described by EA compared with bodies with greater $ψ$ . NMSA is affected by the variability of the geometry: it is not easy to set the height of cross section of the elements in order to obtain accurate results for all the bodies. Therefore, the results are highly affected by the discretization, which should be performed carefully based on the specific case. For this reason, NMSA might have a limited applicability for the structural synthesis of soft mechanisms whose links have non-uniform stiffness.

It is indeed worth keeping in mind that the plot in Figure 13 is valid for the considered techniques in the way we have employed them, and that different implementations (e.g., an alternative set of LPs) would have provided different results. Another important aspect is related to the ratio between the amplitude of the displacement and the length of the body. The higher the ratio, the higher the error for any of the techniques that we have used in this article (later, we provide the reader with some insight concerning this aspect). However, the results that we have obtained are meaningful for two reasons: on the one hand, they show that the SDI can serve as a tool to relate a body to a technique; on the other hand, they provide at least a general guideline for the choice of a modeling technique. It must be clear that the SDI, at least at this stage of the work, does not provide directions about the appropriate modeling assumptions. Any technique, however accurate, leads to poor results if used under incorrect assumptions or without carefully checking the setup. This important matter is more widely addressed in Section 8.

5.4. A further investigation on notched bodies

So far, we have compared the results obtained by four techniques on bodies characterized by Equation (19). The definition of the SDI is based on the flexural stiffness, and therefore should have general validity. However, we repeat the study on bodies having a different profile. The aim is to assess to which extent the results in Figure 13 are valid in a case which differs from the one of cosinusoidal notch.

Therefore, here we consider bodies with rectangular cross section having height which varies according to the following expression:

\begin{matrix} h (x) = & h_{M} + f_{1} (x) \cdot H (\bar{x} - x) \cdot H (x - \bar{x} + a'_{1}) \\ + f_{2} (x) \cdot H (x - \bar{x}) \cdot H (- x + \bar{x} + a'_{2}) \end{matrix}

(28)

in which

f_{i} (x) = c \cdot h_{M} - 2 \bar{a} \sqrt{1 - \frac{{(x - \bar{x})}^{2}}{a_{i}^{2}}}

(29)

a_{i}' = a_{i} \sqrt{1 - {(\frac{c \cdot h_{M}}{2 \bar{a}})}^{2}}

(30)

for $i = 1, 2$ . The notch obtained in this way is elliptical, with smallest height at $x = \bar{x}$ . If it is chosen as $a_{1} \neq a_{2}$ , the notch is not symmetrical with respect to $\bar{x}$ . In fact, the two ellipses share only the semiaxis $\bar{a}$ :

\bar{a} = \frac{1}{2} [(1 + c) h_{M} - h_{m}]

(31)

Figure 14 clarifies the meaning of the coefficients appearing in the expressions above.

We report here the results concerning seven bodies that in our opinion are particularly meaningful. All the bodies considered are characterized by the numerical values listed in Table 2, while they differ for those in Table 3, where each body is denoted by a different capital letter.

Table 2.

Geometrical and material parameters common to all the notched bodies.

Parameter	Description	Value
$c$	Multiplying factor in Equation (29)	0.3
$L$	Length of the body	150 mm
$h_{M}$	Maximum height of the cross section	26 mm
$b$	Width of the cross section	15 mm
$E$	Young’s modulus	10 MPa
$\bar{x}$	Abscissa at which it is $h = h_{m}$	75 mm
$a_{2}$	Semiaxis of the ellipse at the right of $\bar{x}$	5.2 mm

Table 3.

Minimum height $h_{m}$ , semiaxis $a_{1}$ , applied moment $\bar{M}$ , and SDI for the seven bodies considered.

Body, $B$	$h_{m}$	$a_{1}$	$\bar{M}$	SDI
A.	6 mm	20.3 mm	160 Nmm	7.91+ $i$ 0.171
B.	6 mm	15.1 mm	144 Nmm	10.57+ $i$ 0.171
C.	6 mm	4.7 mm	210 Nmm	31.94+ $i$ 0.171
D.	13 mm	20.3 mm	800 Nmm	8.32+ $i$ 0.152
E.	13 mm	9.9 mm	800 Nmm	17.20+ $i$ 0.152
F.	18 mm	20.3 mm	880 Nmm	8.75+ $i$ 0.116
G.	18 mm	9.9 mm	1,040 Nmm	18.24+ $i$ 0.116

The four techniques are implemented as described for the previous set of bodies. The only difference consists of the numerical coefficient that appears in Equation (19): instead of 0.267, we use here 0.2. This choice is discussed at the end of this section.

Table 4 summarizes the results of this investigation. First, let us focus on EA: based on Figure 13, for bodies $B_{A}$ and $B_{B}$ we should expect an error between 10% and 20%, whereas for body $B_{D}$ we should expect an error close to 10%. Results confirm our expectations. This holds for all the other bodies as well: for $B_{F}$ and $B_{G}$ the error is less than 10%; for $B_{I}$ , as expected, we obtain an error slightly greater than 10%. It is worth noting, in particular, that the error is very large for $B_{C}$ : this result suggests the existence of a boundary of the hatched region for EA. The reader should not be surprised from this result. In fact, in Figure 12, we can see that for the body with $h_{m} / h_{M} = 0.25$ and $n = 14$ , EA is affected from an error slightly greater than 20%.

Table 4.

Error $e_{ι}$ for the seven bodies listed in Table 3, obtained by EA, NMSA, and LPM.

Body, $B$	EA	NMSA	LPM
A.	11.42%	72.30%	23.90%
B.	14.44%	74.52%	17.15%
C.	30.19%	73.06%	11.98%
D.	9.53%	26.37%	17.23%
E.	12.21%	27.02%	16.69%
F.	3.33%	7.03%	15.50%
G.	4.23%	6.81%	10.00%

Proceeding with the analysis, we see that also results obtained by NMSA are in accordance with Figure 13: there is no need for detailed comments about this set of results. In contrast, the case of LPM merits great attention. As we stated in Section 5.3, the boundaries of the regions of the complex plane that we can draw for a LPM strongly depend on the specific set of LPs that we define and use. It is true that LPM comes at relatively low computational cost; however, its drawback is that it is not always easy to assign the value of the parameters. In particular, in cases such as those investigated in this article, in which the flexural stiffness is not constant along the axis of the bodies, the identification of the equivalent length in Equation (18) may be tricky. A small difference in a coefficient, such as the numerical value in Equation (19), can relevantly increase or decrease the accuracy. To give a flavour of this fact, we point out that using the coefficient 0.25 instead of 0.20, allows the error $e_{LPM}$ for body $B_{C}$ to be reduced from 11.98% to 4.09%. Moreover, it is worth noting that for $B_{G}$ the error is 10%: in some cases, LPM may perform with such accuracy also for bodies with $ψ ~ 0.1$ . Consistently, going back to Figure 12, we note in fact that LPM provided slightly better results for bodies with $h_{m} / h_{M}$ equal to 0.75 than for those with 0.5. To conclude, we remark that the definition of the LPs turns out to be a critical aspect of the modeling. A detailed dissertation on the best criteria to assign a suitable set of parameters is beyond the scope of this work, and would require a long and deep investigation.

6. SDI for particular cases

The results found so far suggest that the SDI may allow bodies to be classified based on the distribution of the structural properties. However, further considerations must be made to prove its general validity: only the region of the complex plane shown in Figure (13) has been explored so far. In this section, we fill this gap by considering structures that turn out to be particularly meaningful for the purpose of this work. First, we discuss the case of bodies made of multiple materials; then, we introduce classes of bodies for which $χ = 0$ , $0 < χ < 1$ , and $χ < 0$ ; in addition, the particular case of a continuum arm with uniform stiffness is treated.

6.1. Body made of multiple materials

Let us consider a body made of two materials having different Young’s modulus: $E_{1} \neq E_{2}$ . Let us consider a cross section of the body in which both the materials are present. Instead of the actual cross section, it is possible to consider an equivalent cross section made of only one of the two materials (see Beer et al., 2002). If we want to obtain the equivalent section made of material with Young’s modulus $E_{1}$ , the transformation depends on the ratio

m = \frac{E_{2}}{E_{1}}

(32)

that is used as a multiplying factor to transform the area of the cross section in the direction parallel to the neutral axis: in other words, the width of the region with $E_{2}$ is contracted (if $m < 1$ ) or enlarged (if $m > 1$ ). Once the section has been transformed and the position of its neutral axis has been found, the equivalent moment of inertia $J_{e}$ can be computed and the flexural stiffness at that section can be expressed as $E_{1} J_{e}$ . Doing the same for all the cross sections, the function $k_{f_{e}} = E_{1} J_{e} (x)$ can be used to compute the SDI. Therefore, the case of multiple materials can be treated as that of a uniform material and variable cross section.

6.2. Bodies with $χ = 0$

A subset of bodies is characterized by $χ = 0$ . This condition occurs if it is $α (x) = α = const$ . An example is given by the class of bodies made of one material and with rectangular cross section with width $b$ and variable height

h (x) = \sqrt[3]{h_{1} x + h_{2}}

(33)

being $h_{1}$ and $h_{2}$ positive constant coefficients. In this case, it is

k_{f} (x) = \frac{Eb}{12} (h_{1} x + h_{2})

(34)

and, therefore,

α (x) = \frac{Eb h_{1}}{12} = α

(35)

Under this condition, the first factor in the definition of $χ$ is equal to zero; however, based on the provided definition, also the term $x_{M} - x_{m}$ is null. For this reason, $χ$ has been defined to circumvent this singularity, mapping this class of bodies to the imaginary axis of the complex plane. The modeling technique for this kind of body depends therefore only on the slenderness of the structure: bodies with high values of $ψ$ owing to high $d_{M} / L$ ratios are bulky; in most applications, these bodies can be treated as rigid. In contrast, if the body is a slender long arm, then it can be modeled by means of NMSA or EA. This is consistent with what is found in the literature: the stiffness of continuum arms and of slender soft links is modeled by the diagonal matrix used in Cosserat rod theory-based approaches, as mentioned in Section 1. A particular case is represented by the class of continuum arms with uniform cross section and single material. This case is discussed in the following.

6.3. Soft continuum arm with uniform cross section

Let us consider a continuum arm, $L_{a} = 150$ mm long, with circular cross section whose radius is $r_{a} = 2.5$ mm and made of linear elastic material with $E = 1$ MPa. The arm is clamped at one end section and the moment ${\bar{M}}_{a} = 0.8$ Nmm is applied at the free end section. As both the moment and the flexural stiffness are constant along the arm, the deformed configuration is circular. The deformation of the arm is computed by nonlinear FEA with a hexahedral mesh of 700 elements; similarly to the notched body, the body is considered symmetric with respect to the $XY$ plane. Using NMSA, the structure is discretized into 10 segments, all characterized by the same local stiffness matrix (82) with $E_{i} = E$ , $L_{i} = L_{a} / 10$ , $A_{i} = π r_{a}^{2}$ , and $J_{i} = π r_{a}^{4} / 4$ . The number of steps is 300. When using EA, the same settings as used for the notched body are adopted for consistency, although the differential equation is easy to integrate analitically in this case.

Figure 15 shows the arm, initially straight, in the deformed configuration as computed by FEA, in comparison with the deformation obtained by applying EA and NMSA. In this case, EA provides results that slightly differ from those computed by FEA: it is, in fact, $e_{EA} = 0.003$ . NMSA is affected by a greater error ( $e_{NMSA} = 0.025$ ), but describes the deformation of the arm with reasonable accuracy. LPM is not implemented in this case: it would be possible to compute the position of the free end section of the arm by defining an equivalent torsional spring as done by Howell and Midha (1994); however, such an approach would not model the distribution of the deformation along the longitudinal axis.

For this class of bodies, $α (x) = 0$ because $k_{f} (x) = const .$ Therefore, $χ = 0$ and $ψ = 0$ . Thus, the origin of the complex plane represents all the bodies with uniform stiffness, independently of their slenderness. This is the only case in which the choice of the modeling technique should be based on the stiffness only: the distribution of the structural properties is uniform.

6.4. Bodies with $0 < χ < 1$

Here, we consider the classes of bodies characterized by $0 < χ < 1$ . An example is given by conical bodies made of a single material. Being conical, the radius of their circular cross section varies linearly along the length. Denoting by $r_{0}$ and $r_{L}$ the radius at the first and the last section, respectively, the radius at the generic abscissa $x$ is

r (x) = r_{0} - \frac{r_{0} - r_{L}}{L} x

(36)

where $L$ is the length of the body and the origin of the coordinate axis is such that $r (0) = r_{0}$ . The flexural stiffness is the function

k_{f} (x) = E \frac{π}{4} r^{4} (x)

(37)

The first derivative of the stiffness is

α (x) = - E π \frac{r_{0} - r_{L}}{L} {(r_{0} - \frac{r_{0} - r_{L}}{L} x)}^{3}

(38)

We observe that $α (x)$ is a monotonically increasing function in $[0, L]$ . In fact, its first derivative

\frac{d α (x)}{dx} = 3 E π {(\frac{r_{0} - r_{L}}{L})}^{2} {(r_{0} - \frac{r_{0} - r_{L}}{L} x)}^{2}

(39)

is positive everywhere. Hence, it is $α (L) > α (0)$ . In addition, we note that:

for $r_{L} < r_{0}$ , $α (x) < 0$ $\forall x$ ; therefore, $α (L) = α_{M}$ , $α (0) = α_{m}$ , and $| α_{m} | > | α_{M} |$ , thus $\frac{α_{M} - α_{m}}{| α_{m} |} < 1$ ;

for $r_{L} > r_{0}$ , $α (x) > 0$ $\forall x$ ; therefore, $α (L) = α_{M}$ , $α (0) = α_{m}$ , and $| α_{M} | > | α_{m} |$ , thus $\frac{α_{M} - α_{m}}{| α_{M} |} < 1$ .

In both cases, $| x_{M} - x_{m} | = L$ . Hence, we conclude that it is always $0 < χ < 1$ . Figure 16 shows $χ$ and $ψ$ as functions of the ratio $r_{L} / r_{0}$ for conical bodies having $E = 10$ MPa, $r_{0} = 13$ mm, and $r_{L}$ ranging in the interval $[1, 26]$ mm. We can see that for $r_{L} / r_{0} = 1$ it is $SDI = 0 + i 0$ , consistently with what was reported in the previous section for a soft continuum arm with uniform cross section.

It is fairly intuitive that for this class of bodies, the most suitable choice in terms of modeling technique is represented either by EA or NMSA. According to Figure 13, by using EA we should compute the position of the free end section with an error up to 10%. To verify that this is true, we have used the techniques on 12 conical bodies. Material and geometrical parameters for this study are listed in Table 5, whereas the applied external moment at the free end section is reported in Table 6.

Table 5.

Geometrical and material parameters for conical bodies.

Parameter	Description	Value
$L$	Length of the body	200 mm
$r_{0}$	Radius at section $x = 0$	${7, 9, 11, 13}$ mm
$r_{L}$	Radius at section $x = L$	${2, 4, 6}$ mm
$E$	Young’s modulus	10 MPa

Table 6.

Moment $\bar{M}$ applied on conical bodies for varying $r_{L}$ and $r_{0}$ . The conical bodies are denoted by $C_{γ}$ with $γ$ ranging from 1 to 12.

$γ$	$r_{0}$	$r_{L}$	$\bar{M}$
	(mm)	(mm)	(Nmm)
1	7	2	12
2		4	50
3		6	200
4	9	2	16
5		4	120
6		6	280
7	11	2	16
8		4	150
9		6	300
10	13	2	20
11		4	150
12		6	340

Both techniques perform as expected. Figure 17 allows a comparison of the results given by EA, NMSA, and FEA for body $C_{6}$ . Table 7 summarizes the errors for the two techniques, as well as the SDI, for each body analyzed. We observe, in particular, that for $C_{10}$ NMSA performs with an error equal to 10.4%. This suggests that the boundary that we attempted to sketch for NMSA(0.1) in Figure 13 needs correction. It is worth pointing out that, for this set of analyses, we have used NMSA discretizing the structure into six elements, as done in Section 5.2.2. A finer discretization would lead to lower errors.

Table 7.

SDI, $e_{EA}$ , and $e_{NMSA}$ for conical bodies $C_{γ}$ .

$γ$	SDI	$e_{EA}$	$e_{NMSA}$
1	0.977+ $i$ 0.069	1.2%	5.9%
2	0.814+ $i$ 0.062	1.8%	2.0%
3	0.370+ $i$ 0.032	2.2%	2.5%
4	0.990+ $i$ 0.090	1.6%	8.3%
5	0.913+ $i$ 0.086	0.5%	3.5%
6	0.704+ $i$ 0.072	1.0%	2.1%
7	0.995+ $i$ 0.110	2.9%	8.2%
8	0.953+ $i$ 0.108	0.9%	4.3%
9	0.838+ $i$ 0.100	1.5%	1.9%
10	0.997+ $i$ 0.130	3.0%	10.4%
11	0.972+ $i$ 0.129	1.5%	4.1%
12	0.902+ $i$ 0.124	1.4%	1.7%

Conical bodies are not the only ones for which is it $0 < χ < 1$ . Another example is represented by the class of bodies with circular cross section whose radius varies exponentially:

r (x) = R_{0} e^{b \cdot x / L}

(40)

with $x \in [0, L]$ . Here $R_{0}$ is positive and the coefficient $b$ is either positive or negative. For these bodies, it is

k_{f} (x) = E \frac{π}{4} R_{0}^{4} e^{4 b \cdot x / L}

(41)

α (x) = \frac{b}{L} E π R_{0}^{4} e^{4 b \cdot x / L}

(42)

\frac{d α (x)}{dx} = 4 \frac{b^{2}}{L^{2}} E π R_{0}^{4} e^{4 b \cdot x / L}

(43)

As we observed for conical bodies, from Equation (43) we infer that $α (L) > α (0)$ , because $α (x)$ monotonically increases in $[0, L]$ . Moreover:

for $b < 0$ , $α (x) < 0$ everywhere; therefore, $α (L) = α_{M}$ , $α (0) = α_{m}$ , and $| α_{m} | > | α_{M} |$ (because $e^{- 4 | b |} < 1$ );

for $b > 0$ , $α (x) > 0$ everywhere; therefore, $α (L) = α_{M}$ , $α (0) = α_{m}$ , and $| α_{M} | > | α_{m} |$ (because $e^{4 b} > 1$ ).

We reach the same conclusions drawn for conical bodies. In addition, we can note that the case $b = 0$ is analogous to the case $r_{L} / r_{0} = 1$ discussed for the conical bodies: in both cases, the radius of the cross section is constant along the length. Figure 19 collects in the complex plane the results shown in Figure 16. In the same figure, we report the SDI computed for bodies with exponential profile with $R_{0} = 13$ mm and $b$ ranging in $[- 1, 1]$ , while keeping the same body length $L$ and Young’s modulus $E$ . For the sake of brevity, we do not report here the results of a set of analyses performed on this class of bodies, except the case $b = - 1$ for which the body has $SDI = 0.982 + i 0.127$ . For an applied moment $\bar{M} = 150$ Nmm, we obtain $e_{EA}$ =1.9% and $e_{NMSA}$ = 2.1%. Results are displayed in Figure 18.

6.5. Bodies with $χ < 0$

Finally, we provide an example of a class of bodies for which it is $χ < 0$ . Let us consider bodies made of a single material with variable rectangular cross section whose height is given by the function

\begin{matrix} h (x) & = h_{M} - (h_{M} - h_{m}) H (x - \frac{L}{3}) \\ + \frac{1}{2} (h_{M} - h_{m}) (\cos \frac{n' π (x - \frac{L}{3} - \frac{L}{n'})}{L} - 1) \cdot \\ H (\frac{L}{3} - x) H (x - \frac{L}{3} + \frac{L}{n'}) \end{matrix}

(44)

Again, $H (\cdot)$ is the Heaviside function and the parameter $n'$ is used to tune the distribution of the softness of the link. An example of bodies of this class is shown in Figure 20, in which $n' = 4$ , $h_{M} = 16$ mm, $h_{m} = 6$ mm, $b = 10$ mm, and $L = 100$ mm.

In this case, $χ = - 12.21$ and $ψ = 0.15$ . The deformation of bodies of this class can be effectively modeled by EA, same as the notched bodies characterized by $SDI = - χ + i ψ$ ; however, the partition of the complex plane as shown in Figure 13 is not symmetric with respect to the imaginary axis. The difference is that because the considered body has no notch, but only a gradual monotonic decreasing of the stiffness along the axis, there is no hinge-like region to model by means of a LP. Consistently with what was obtained for $χ > 0$ , instead, NMSA is suggested as a convenient technique for bodies having $ψ$ lower than ~0.11, such as the body represented in Figure 21, which differs from the previous in that $h_{m} = 12$ mm, resulting in $SDI = - 9.47 + i 0.09$ .

7. SDI for optimization

So far, we have observed that the SDI can lead to a consistent partition of the complex plane, allowing us to draw regions that we can link to one or more modeling techniques. This encouraging result pushes us to investigate another aspect: can the SDI be useful for design purposes? Can it play a role in the structural optimization of a soft link?

In this section, we present an optimization problem that has helped us to better assess the validity and the usefulness of the SDI. Thanks to its simplicity, the example reported here shows how the SDI can guide us in the optimization.

7.1. Statement of the problem

Let us design a link of length $L$ , clamped at one end and subject to a load applied at the other end as shown in Figure 22. By actuating the system with the moment $M_{a}$ we want not only to compensate for the deflection caused by the load $F_{v}$ , but also to displace the end section of the link to the point $Q$ denoted in figure. Then, the problem consists of finding a link whose end tip undergoes the prescribed displacement under the given loading and constraint conditions.

To work with a low number of design variables and simplify the discussion, we assume that the body is made of a single material that we treat as linear elastic. In addition, we prescribe that the body has circular cross section whose diameter can range between a minimum $d_{m}$ and a maximum $d_{M}$ (for instance, for reasons related to the manufacturing process).

7.2. Optimization

A very simple and intuitive model for this problem is represented in Figure 23: the body consists of two quasi-rigid links connected by a thin, circular beam that behaves as a flexural pivot and geometry of which has to be found. The compliance of the body can be considered concentrated between points $T$ and $S$ ; for this reason, in the following we refer to this model as a LPM. From now on, the diameter $d$ , the length $l$ of this beam, and the distance $D$ from the clamped end will be the design variables in the optimization problem. We collect them in the vector

v_{LP} = {D, l, d}

(45)

The two rigid links have given diameter $d_{M}$ . If $d_{M}$ is sufficiently greater than $d$ , the deformed configuration of the body can be easily obtained by computation of the deflection of the thin, deformable beam. The computation is performed with respect to the coordinate system shown in Figure 23.

Unlike the body in Section 5.2.4, here the beam is subject to a non-constant moment, owing to the force $F_{v}$ . We note that the moment depends both on the abscissa $x$ and on the position of the free end of the body. Therefore, here we follow an incremental procedure applying the load in a number $n_{s}$ of steps. If we assume that the section at point $T$ is clamped, the vertical displacement of point $S$ at the $s$ th step can be computed according to the classical beam theory as

y_{S_{s}} = y_{S_{s - 1}} + (m_{a} + f_{v} a \cos θ_{S_{s - 1}} + f_{v} l) \frac{l^{2}}{2 EJ} - \frac{f_{v} l^{3}}{6 EJ}

(46)

with $a = L - D - L_{p}$ , $J = π d^{4} / 64$ , $f_{v} = F_{v} / n_{s}$ , and $m_{a} = M_{a} / n_{s}$ . The rotation at point $S$ is computed as

θ_{S_{s}} = θ_{S_{s - 1}} + (m_{a} + f_{v} a \cos θ_{S_{s - 1}} + f_{v} l) \frac{l}{EJ} - \frac{f_{v} l^{2}}{2 EJ}

(47)

The displacement of $S$ along the $x$ -axis is neglected here. Therefore, $x_{S_{i}} = D + l$ for any $i$ . The coordinates of point $P$ at the $i$ th step are collected in the vector

p_{s} = {x_{S_{s}}, y_{S_{s}}}^{T} + {a \cos θ_{S_{s}}, a \sin θ_{S_{s}}}^{T}

(48)

The objective of the optimization is to minimize the distance between point $P$ in the deformed configuration and point $Q$ . Therefore, the objective function is written here as

ϕ (v) = | | p_{n_{s}} - q | |

(49)

Here, we perform optimization in MATLAB environment by the fmincon built-in function, which requires the assignment of an initial guess for the set of design variables. Based on our intuition, a reasonable initial guess is the set

v_{0_{LP}} = {\frac{2 L}{3}, \frac{L}{10}, d_{m}}

(50)

Moreover, the solution is searched accounting for the following set of constraints:

0 < D < L

(51)

0 < l < \frac{L}{2}

(52)

d_{m} \leq d \leq d_{M}

(53)

If we adopt the numerical values listed in Table 8, we find that the solution is $D = D^{(1)} = 109.2$ mm, $l = l^{(1)} = 52.2$ mm, and $d = d^{(1)} = 6$ mm. In the following, we denote the body having these dimensions as $B_{1}$ . The final value of the objective function turns out to be 7.3 mm: according to the model, the coordinates of point $P$ under the given loads would be $(176.3, 63.6)$ , against the prescribed $(170, 60)$ . Such error might seem acceptable, if we consider that soft robotics applications do not generally require high accuracy. Therefore, we might be tempted to accept the solution and build the body based on it. However, if we perform FEA using the commercial software on the body having the optimized geometry, we realize that the actual displacement of $P$ is quite far from the prescribed value: the coordinates of point $P$ in the deformed configuration are $(143.5, 72.1)$ according to FEA. Figure 24 shows this discrepancy, that is due to the approximations by which the model is affected (primarily, the displacement of point $S$ along the $x$ -axis is neglected). The error evaluated as defined by Equation (27) is equal to 37.0%.

Table 8.

Geometrical and material parameters in the optimization process.

Parallel	Description	Value
$L$	Length of the body	200 mm
$d_{M}$	Maximum diameter of the cross section	20 mm
$d_{m}$	Minimum diameter of the cross section	6 mm
$d_{m}'$	Minimum diameter of the cross section	6 mm
$F_{v}$	Applied force at end tip	−0.5 N
$M_{a}$	Applied moment at end tip	100 Nmm
$E$	Young’s modulus	50 MPa
$P$	Coordinates at initial configuration	$(200, 0)$ mm
$Q$	Prescribed coordinates	$(170, 60)$ mm
$Q'$	Prescribed coordinates	$(187, - 60)$ mm
$n_{s}$	Number of steps	40

As the problem seems to require a larger computational effort, let us try again using a different technique, setting the optimization problem by using NMSA. The objective function has still the form in Equation (49), but now it depends on a different set of design variables:

v_{NMSA} = {D, l_{1}, l_{2}, l_{3}, d_{1}, d_{2}, d_{3}}

(54)

whose meaning is clarified by Figure 25, which also shows the global enumeration of the nodes of the discretization. The structure is divided into five elements, serially connected; the nodes are enumerated from 1 to 6, being the first node clamped. As was done for the notched body, the load is applied incrementally in $n_{s}$ steps. The stiffness matrix of the structure depends on the design variables: in fact, the lengths of all the elements are unknown (the only condition is that the length of the fifth element must be equal to $L - D - \sum_{i} l_{i}$ ) as well as the areas and the moments of inertia of the second, third, and fourth elements.

The optimization is performed under the following constraints:

0 < D < \frac{2 L}{3}

(55)

0 < l_{i} < \frac{L}{3} for i = 1, 2, 3

(56)

d_{m} \leq d_{i} \leq d_{M} for i = 1, 2, 3

(57)

and from an initial guess based on the solution computed by LPM:

v_{0_{NMSA}} = {D^{(1)}, \frac{l^{(1)}}{3}, \frac{l^{(1)}}{3}, \frac{l^{(1)}}{3}, d_{ave}^{(1)}, d^{(1)}, d_{ave}^{(1)}}

(58)

where $d_{ave}^{(1)} = (d^{(1)} + d_{M}) / 2$ . The objective function takes the value $7 \times 10^{- 6}$ mm for $D = D^{(2)} = 104.8$ mm, $l_{1} = l_{1}^{(2)} = 14.2$ mm, $l_{2} = l_{1}^{(2)} = 49.9$ mm, $l_{3} = l_{1}^{(2)} = 18.3$ mm, $d_{1} = d_{1}^{(2)} = 12.5$ mm, $d_{2} = d_{3}^{(2)} = 6.0$ mm, and $d_{3} = d_{3}^{(2)} = 7.6$ mm. The result of this optimization will be hereafter denoted by $B_{2}$ . We note that, unlike in the previous case, the value of the objective function approaches zero. Moreover, as visible in Figure 26, the difference between the deformation computed by NMSA and by FEA is still non-negligible, but smaller. From FEA, we see that the coordinates of point $P$ in the deformed configuration are $(150.9, 62.8)$ . Therefore, by using NMSA, we obtain an error equal to 21.9%.

Based on the results reported previously, we might conclude that the LPM is not suitable at all, and that we should have recurred to NMSA from the beginning, possibly with a number of elements greater than five, to further reduce the error. However, this would be a mistake. We prove this with a simple example. Let us consider a similar problem, with a different prescribed point $Q'$ (reported in Table 8) and only the force $F_{v}$ applied. We search the solution by using the same LPM; under a new set of constraints

0 < D < L

(59)

0 < l < \frac{7 L}{100}

(60)

d_{m}' \leq d \leq d_{M}

(61)

and starting the search from the initial guess

v_{0} = {\frac{L}{3}, \frac{L}{20}, d'_{m}}

(62)

we find that the objective function takes the value $2.8 \times 10^{- 5}$ mm for $D = 55.0$ mm, $l = 13.5$ mm, and $d = 5.4$ mm. As done previously, let us perform FEA on a body having these dimensions. Figure 27 shows the comparison between the results obtained by FEA and by LPM. The deflection computed by LPM nearly overlaps with that obtained by FEA: the position of the free end is found to be $Q'$ and $(185.0, - 63.4)$ with the two methods, respectively, resulting in an error equal to 6.0%.

It is evident that, in the last case, the LPM represents a good candidate for the optimization process. Compared with NMSA, it is easier to implement and requires familiarity with the classical beam theory only. In addition, it makes the optimization process notably faster than that involving NMSA, owing to the lower number of design variables. In fact, the price that we pay for a more accurate model is related to the computational cost: to find the geometry of $B_{2}$ , we have performed a search in a seven-dimensional space, but for $B_{1}$ and $B_{3}$ we searched in a three-dimensional space. To provide the reader with at least a rough idea about the difference, we report that the geometries of $B_{1}$ and $B_{3}$ have been found in less than 0.05 seconds, whereas the search for $B_{2}$ has taken 4.8 seconds on the same machine. Although it may seem not of great importance here, we should consider that the examples presented above are extremely simple: the optimization is performed on a single body, prescribing the displacement of a single point, and keeping the number of design variables as low as possible; in addition, let us remark that we are accepting what probably is a local minimum as a solution, and we are not verifying whether it is the global minimum of the objective function. The matter assumes importance for more complex optimization problems, and becomes crucial when searching in high-dimensional spaces.

These results show that the suitability of the modeling technique should be evaluated for each specific case. Can the SDI help us to assess whether we have chosen an appropriate technique, without verifying the results by FEA or other computationally expensive methods?

7.3. The role of SDI in optimization

We now recapitulate all the results found in the previous section, and we compare the bodies $B_{1}$ , $B_{2}$ , and $B_{3}$ in terms of their SDI. Before starting, we note that for all the bodies under considerations, jump discontinuities of the cross-section occur along the longitudinal axis: from a rigorous mathematical viewpoint, this represents an issue in view of our purpose, because we want to compute the first derivative of the flexural stiffness (however, such discontinuity would not occur in a physical prototype, because perfectly sharp edges are a pure abstraction in mechanics). To obtain the SDI in these cases, we suggest to compute the function $α (x)$ recurring to finite differences. In the following, we discretize the interval $[0, L]$ with 1,000 evenly spaced points and we approximate

α (x_{i}) = \frac{Δ k_{f_{i}}}{Δ x} = \frac{k_{f_{i + 1}} - k_{f_{i}}}{x_{i + 1} - x_{i}}

(63)

Figure 28 shows the profile of the three bodies and the function $α (x)$ , while Table 9 lists the computed errors and the SDI for each body. Let us now recall the entire optimization process, to comment on it in view of these values. The process started by using a LPM, which returned as a result a body with $SDI = 7.6 + i 0.0992$ ; from the sketch in Figure 13, it seems that the use of a LPM is not suggested in this case, and that a more accurate technique should be employed. In fact, by implementing NMSA with a coarse discretization, we found a smaller error. In contrast, in the third optimization problem we found a body with $SDI = 31.4 + i 0.0995$ : for bodies having such high values of the real part of the SDI and $ψ ~ 0.1$ , LPMs may be a suitable option as shown in Section 5.4. The relatively small error (6.0%) seems a confirmation of this conclusion. Therefore, the SDI could have driven us in the appropriate choice of the technique. The role played by the SDI can be twofold, indeed: (i) helping in the modeling of an already existing body; and (ii) suggesting how to set up an optimization problem to find a body that fulfills prescribed requirements.

Table 9.

Summary: SDI and error computed for each body.

Body	SDI	Error
$B_{1}$	$7.6 + i 0.0992$	37.0% by LPM
$B_{2}$	$4.5 + i 0.0992$	21.9% by NMSA
$B_{3}$	$31.4 + i 0.0995$	6.0% by LPM

Concerning the latter, we provide in the following a further example of optimization, based on a different approach. The aim is to show that it is possible to choose a modeling technique a priori and use it for a constrained optimization whose constraints are written based on the SDI.

7.4. Optimization constraints based on SDI

Let us return again to the geometry used in Section 5. This time, the objective is to design a notched body that undergoes the vertical displacement $y_{P} = 35$ mm under an external moment $\bar{M}$ . Compared with the previous optimization problem, here we assign the technique to be used. As we have seen, each technique cannot be considered reliable in general: provided results are reasonably accurate depending on the specific distribution of the stiffness of the body.

In this optimization problem, we treat the parameters $n$ and $h_{m}$ as design variables, that is, we tune the length and the depth of the notch. The variables are collected in

w = {h_{m}, n} .

(64)

Numerical values of all other parameters are assigned as reported in Table 1 and we adopt the same system of coordinates used so far. We choose to implement the model using the EA, with the same settings described in Section 5.2.3; the applied moment is $\bar{M} = 8$ Nmm. Hereafter, we denote by $y_{f}$ the displacement along the vertical axis that the free end section of the body undergoes under the applied load. The objective here is to minimize the absolute value of the function

φ (w) = y_{f} - y_{P}

(65)

and we perform the optimization imposing the following constraints:

\frac{h_{M}}{10} < h < \frac{9 h_{M}}{10}

(66)

2 < n < 14

(67)

starting the search from

w_{0}^{(1)} = {\frac{3 h_{M}}{20}, 8}

(68)

The solver succeeds to find a local minimum while satisfying the assigned constraints: at the point $w^{(1)} = {4.84, 7.68}$ the value of the objective function is $1.7 \times 10^{- 7}$ mm.

However, we see from Figure 29a that there is a discrepancy between results obtained by EA and by FEA: the coordinates of the end point of the neutral axis of the body are $(84.96, 39.07)$ according to FEA, whereas we compute $(87.90, 35.00)$ by EA. The error is equal to 12.0%. Can we obtain a more accurate result? Owing to the nature of the prescribed requirement, it is fairly intuitive that the objective function is null at more than one point of the search space. Let us perform again the optimization process, starting the search from a different point:

w_{0}^{(2)} = {\frac{h_{M}}{4}, 11.5}

(69)

We find that the value of the objective function for $w^{(2)} = {4.06, 11.79}$ is $1.9 \times 10^{- 7}$ mm. As done before, we compare the results obtained by EA with those provided by FEA. In this case, we obtain the coordinates $(80.89, 42.55)$ by FEA, against $(87.78, 35.00)$ computed by EA. We conclude that $w_{0}^{(2)}$ is a more unfortunate initial guess than $w_{0}^{(1)}$ : by this second optimization, we obtain an error equal to 21.9%.

From what has been reported above, we draw two conclusions: (i) as expected, the problem admits more than one solution, i.e., there is more than one body fulfilling the prescribed requirement; and (ii) solutions may largely differ in terms of error $e_{EA}$ . Figure 30 discourages us from repeating the optimization with a different initial guess: the function $φ$ as computed by EA is characterized by a wavy trend and its intersection with the plane $z = 0$ consists of several curves (we show $φ$ , rather than its absolute value, for a mere convenience of display). Therefore, we try a different approach, involving the SDI. Let us add, indeed, an additional constraint to the set introduced above. We require that for the body the following inequality must hold:

ψ (w) < 0.15 \cdot H (χ (w) - 6)

(70)

and we perform again the optimization starting the search from $w_{0}^{(3)} = w_{0}^{(1)}$ . This time, we find that the objective function takes the value $1.2 \times 10^{- 7}$ mm at $w^{(3)} = {6.31, 4.00}$ while satisfying all the assigned constraints. The body is characterized by $χ = 5.85$ and $ψ = 0.15$ . The inequality constraint (70) has been written based on the sketched lines in Figure 13 to guarantee that the SDI of the body resulting from optimization lies inside the blue region. In other words, such constraint allows us to be confident about the choice of the modeling technique used for the optimization. According to Figure 13, in this case the error should be less than 10%. In fact, from FEA we find that the coordinates of the end point of the axis are $(86.70, 37.65)$ , against $(88.24, 35.00)$ , resulting in an error equal to 7.7%. Figure 29c shows the comparison between the results obtained by the two techniques.

It is also worth pointing out that the bodies in Figure 29a and (b) have SDI equal to $11.05 + i 0.1555$ and $16.58 + i 0.1573$ , respectively. Consistently, for the former we have obtained an error slightly greater than 10%; for the latter, an error around 20%. Concerning the latter, from Figure 12 we would predict an error around 18%. However, the external load applied here is greater than the one reported in Section 5. We remark, as previously, that the discrepancy between FEA and any technique tends to increase drastically as the external load increases, owing to the increasingly relevant effect of phenomena such as warping of the cross section.

8. Discussion

As stated in Section 2, the choice of a convenient modeling technique is crucial not only when analysis must be performed, but especially when it must be used for synthesis. For analysis, the choice depends of course on the specific phenomenon or condition that has to be modeled. For instance, finite elements are powerful when the mechanical stress and strain must be computed, or when contact between bodies (or even self-contact within the same body) occurs. However, if the objective of the analysis is the computation of the deformed configuration under an external load, or the preliminary design of a soft body, finite elements are not necessarily needed. We have shown that accurate results can be obtained, in some cases, by using alternative techniques, at a much lower computational cost and in a way at the reach of our intuition. In the case of synthesis, we believe that the choice of a technique has a non-negligible influence on our understanding of the working principles of soft robots and, in general, of all the systems. Not without reason, in fact, several soft robots and actuators have been built based on the researchers’ intuitions and qualitative considerations: in our opinion, this is a consequence of the lack of adequate modeling approaches.

The SDI proposed in this article allows to classify bodies based on the distribution of their structural properties. To summarize, it helps in the choice of a modeling technique based on such distribution; it can be considered as a tool to gain a hint on the number of parameters and on the technique needed to model the softness of a robotic link with a given accuracy. Let us consider a body such as that described in Section 5, with the same length but with two equal notches. The SDI would be equal to that computed with the same body with a single notch: this means that the SDI suggests the same technique for both the bodies. In fact, if the body with the single notch can be effectively modeled by LPM, then the body two notches also can: the difference is that an additional LP is needed in the model. If we consider, instead, a long continuum arm with uniform cross section except for the presence of a narrow notch, the real part of the SDI would assume a great value (depending on the length of the arm) and the imaginary part would be small. This is consistent with the plot in Figure 13: the arm should be modeled by EA or NMSA. Possibly, a hybrid technique including also a LP to model the notch could be adopted.

The SDI introduced in this article accounts for the distribution of the flexural stiffness. However, this is not the only important quantity to consider. Soft bodies for robotic applications may be, in general, intended not only to bend, but also to stretch, to undergo torsion, inflation, under combination of loads that act simultaneously (point forces, distributed loads, and moments). It will be necessary to define additional indices to fully classify soft-bodied robots and develop a comprehensive set of tools to be used for analysis and synthesis purposes. The path ahead is still long.

Another important aspect concerns the assumptions on which a given modeling technique is based. In this article, for instance, we have neglected the shear effects when implementing LPM, EA, and NMSA; however, in some cases it is needed to introduce shear correction factors that depend on the geometry of the body (as explained by Timoshenko beam theory). Another assumption, widely employed in soft robotics, is that segments of the body bend with constant curvature; in some works, researchers have shown that this assumption has led to acceptable results, but this does not mean that the assumption could always be adopted. In some works, it is assumed that the cross-section does not change with the deformation (see, e.g., Renda et al., 2017: in which the area and the moment of inertia of the cross section are constant). However, warping of the section, local deformations and high local strains may play a big role. In Section 7.3, we have seen that for body $B_{2}$ , NMSA with five elements performs with $e_{NMSA}$ = 21.9%. Based on the SDI of the body, we would have expected a better performance from NMSA. Unfortunately, NMSA as implemented in this work does not account for local peaks of strain, which occur at every abrupt change of the diameter of the cross section. A more sophisticated implementation, comprehensive of correction factors to account for this phenomenon, would have provided more accurate results.

To summarize, the assumptions on which a model can be based are numerous. A further aim of our future work will be to identify criteria that encourage or discourage the use of the commonly used assumptions, based on the distribution of the structural properties of the system. To reach this goal, the first derivative of $k_{f} (x)$ will be a valuable ally, together with its higher-order derivatives. Concerning the bending, this will allow a clearer definition of the partition of the complex plane (at the moment, merely sketched). We clarify that even a comprehensive set of SDIs will never serve as a substitute of the researcher’s skill in using the various techniques. The greatest caution should be exercised even when implementing the most accurate method: to give an example among several, finite elements with inappropriate aspect ratio can cause inaccurate results. The SDI represents our contribution to the creation of a framework to model soft-bodied robots: we believe, indeed, that such framework cannot be built without accounting for the distribution of the structural properties of the huge variety of soft bodies that can be used. A possible approach consists in discretizing the soft robotic structure in parts and compute the SDI for each part: in this way, one can operate a quantitative distinction between segments with highly distributed compliance from others whose modeling is possible with few parameters. We insist particularly on the fact that considerations made on the softness distribution of a robot, or body, should be quantitative and not only qualitative. Our intuition is powerful, indeed, but limited: it can push the development of the field up to an appreciable boundary, but mathematical tools are needed to cross the limit and achieve the design of efficient soft-bodied robots. We leave the reader with an open question, indeed: are we fully exploiting the potential offered by the structural properties of soft-bodied systems?

9. Conclusions

The main objective of this work has been to highlight the importance of the distribution of the structural properties in a body undergoing large deformation and that such distribution is related to the choice of a modeling technique. We have proposed a novel index, called SDI, that accounts for the structural properties of a body performing bending.

In order to carry out a quantitative investigation, we have modeled the bending of soft bodies recurring to different techniques and we have shown that it is possible to map bodies to the complex plane using the SDI. The complex plane has been divided into regions, each corresponding to one or more modeling techniques suitable for the body under investigation, showing that the provided definition of the SDI makes a classification possible.

The achievement of a unified framework for the modeling of soft robots will require further work from the soft robotics community. We believe that this framework should include all the modeling techniques employed so far by the researchers, possibly developing hybrid techniques that could combine the use of FE, MSA, Cosserat rod theory, elastica theory, LPs, and different assumptions (for instance, the piecewise constant curvature, frequently recalled in the literature), to be selected based on the SDI of the body. We do not claim that the proposed SDI represents a fully developed and tested tool to solve the problem completely; in contrast, we believe that additional efforts are required from the soft robotics community to achieve a complete classification of soft bodies based on the distribution of their properties, and the corresponding choice of a modeling approach. We believe that other indices should be defined to account for other typical deformation modes of soft robots (for instance, tension and compression, squeezing). In this article, we have provided evidence that it is possible to find a relation between bodies and modeling techniques by means of a mathematical tool.

Footnotes

Appendix A: Invariance of SDI under inversion of x -axis

Let us consider the function $k_{f}$ for a generic body at a given coordinate $\tilde{x}$ in a system of reference such as that in Figure 2. By inversion of the $x$ -axis and shifting the origin at the opposite end section, we operate the transformation $x \to x' = L - x$ and it is $k_{f} (x) = k'_{f} (x') = k'_{f} (L - x)$ .

By definition, the first derivative of $k_{f}$ at $\tilde{x}$ is

(71)

α (\tilde{x}) = lim_{Δ x \to 0} \frac{k_{f} (\tilde{x} + Δ x) - k_{f} (\tilde{x})}{Δ x}

Let $Δ x > 0$ and $\tilde{x}$ such that it holds

(72)

k_{f} (\tilde{x} - Δ x) < k_{f} (\tilde{x}) < k_{f} (\tilde{x} + Δ x)

that is, the function $k_{f}$ is increasing at $\tilde{x}$ (it is $α (\tilde{x}) > 0$ ). Using the aforementioned transformation, this condition can be rewritten as

(73)

k'_{f} (\tilde{x}' + Δ x') < k'_{f} (\tilde{x}') < k'_{f} (\tilde{x}' - Δ x')

Therefore, the difference quotient

(74)

\frac{{k'}_{f} (\tilde{x}' + Δ x') - {k'}_{f} (\tilde{x}')}{Δ x'}

is negative. Computing the limit for $Δ x' \to 0$ , we conclude that

(75)

α (\tilde{x}) = - α' (\tilde{x}')

If $x_{1}$ and $x_{2}$ are the coordinates at which $α (x)$ is maximum and minimum, respectively, it is $α (x_{1}) > α (x_{2})$ . Under the coordinate transformation and considering the condition (75), we have

(76)

α' (x_{1}') < α' (x_{2}')

that is, the function $k_{f}' (x')$ is maximum at $x'_{2}$ and minimum at $x'_{1}$ .

As it holds

(77)

\max {| α (x_{1}) |, | α (x_{2}) |} = \max {| α' (x'_{1}) |, | α' (x'_{2}) |}

according to the definition of $χ$ given $k' (x')$ , if the SDI is invariant it must be

(78)

\frac{α' (x_{2}') - α' (x_{1}')}{x_{2}' - x_{1}'} = \frac{α (x_{1}) - α (x_{2})}{x_{1} - x_{2}}

The condition above is satisfied. In fact, it is

(79)

α' (x_{2}') - α' (x_{1}') = - α (x_{2}) + α (x_{1})

and

(80)

x_{2}' - x_{1}' = L - x_{2} - L + x_{1} = x_{1} - x_{2}

It follows that $χ$ is invariant under the considered coordinate transformation. Moreover, because the maximum and minimum value of $k_{f}$ do not depend on the choice of the coordinate system, the imaginary part $ψ$ of the SDI is invariant. We conclude that the SDI does not depend on the orientation of the coordinate axis.

Appendix B: NMSA

Same as for the FEM, the order in which the DOFs are listed is arbitrary. In our work, following the common practise, we define the vector of nodal displacements of the $i$ th beam element in the local system of coordinates ( $ξ^{(i)}$ and $η^{(i)}$ ) (shown in Figure 31) as

(81)

u_{l_{i}} = {ξ_{1}^{(i)} η_{1}^{(i)} γ_{1}^{(i)} ξ_{2}^{(i)} η_{2}^{(i)} γ_{2}^{(i)}}^{T}

In this work, the nodal rotations $γ_{1}^{(i)}$ and $γ_{2}^{(i)}$ are considered positive if counter-clockwise.

Coherently, the local stiffness matrix $K_{l_{i}} \in R^{6 \times 6}$ is written as

(82)

K_{l_{i}} = \frac{E_{i}}{L_{i}} [\begin{matrix} A_{i} & 0 & 0 & - A_{i} & 0 & 0 \\ 0 & \frac{12 J_{i}}{L_{i}^{2}} & \frac{6 J_{i}}{L_{i}} & 0 & - \frac{12 J_{i}}{L_{i}^{2}} & \frac{6 J_{i}}{L_{i}} \\ 0 & \frac{6 J_{i}}{L_{i}} & 4 J_{i} & 0 & - \frac{6 J_{i}}{L_{i}} & 2 J_{i} \\ - A_{i} & 0 & 0 & A_{i} & 0 & 0 \\ 0 & - \frac{12 J_{i}}{L_{i}^{2}} & - \frac{6 J_{i}}{L_{i}} & 0 & \frac{12 J_{i}}{L_{i}^{2}} & - \frac{6 J_{i}}{L_{i}} \\ 0 & \frac{6 J_{i}}{L_{i}} & 2 J_{i} & 0 & - \frac{6 J_{i}}{L_{i}} & 4 J_{i} \end{matrix}]

in which $E_{i}$ is the Young’s modulus, $A_{i}$ the area of the cross section, $L_{i}$ is the length, and $J_{i}$ is the moment of inertia of the cross section of the $i$ th element.

The local nodal displacements are expressed in the global system of coordinates by means of the transformation

(83)

u_{i} = Q_{i} u_{l_{i}}

where $Q_{i}$ denotes the rotation matrix, that can be written as

(84)

Q_{i} = [\begin{matrix} {\tilde{Q}}_{i} & 0_{3 \times 3} \\ 0_{3 \times 3} & {\tilde{Q}}_{i} \end{matrix}]

in which it is

(85)

{\tilde{Q}}_{i} = [\begin{matrix} \cos (θ_{i}) & - \sin (θ_{i}) & 0 \\ \sin (θ_{i}) & \cos (θ_{i}) & 0 \\ 0 & 0 & 1 \end{matrix}]

being $θ_{i}$ the angle between the $ξ_{i}$ -axis and the global $x$ -axis, measured positive in the counter-clockwise direction.

The strain energy for the $i$ th element can be written as

(86)

V_{i} = \frac{1}{2} u_{l_{i}}^{T} K_{l_{i}} u_{l_{i}}

using the local quantities, or as

(87)

V_{i} = \frac{1}{2} u_{i}^{T} K_{i} u_{i}

as a function of the nodal displacements and the stiffness matrix expressed in global coordinates. As the strain energy is invariant under coordinate transformation, it must be

(88)

\frac{1}{2} u_{l_{i}}^{T} K_{l_{i}} u_{l_{i}} = \frac{1}{2} u_{i}^{T} K_{i} u_{i}

from which follows that

(89)

K_{i} = Q_{i} K_{l_{i}} Q_{i}^{T}

by substitution of Equation (83).

The stiffness matrix derived so far belongs to $R^{6 \times 6}$ . If $N$ is the total number of nodes of the structure, the size of the global stiffness matrix is $3 N \times 3 N$ ; therefore, each matrix $K_{i}$ must be expanded to match the global dimensions of the problem. Such expansion is performed by means of the Boolean localization $6 \times 3 N$ matrix

(90)

\begin{matrix} B_{i} = [B_{ι κ}^{(jk)}] = \\ {\begin{matrix} 1 & if (ι, κ) \in {(1, 3 j - 2), (2, 3 j - 1), (3, 3 j)} \\ 1 & if (ι, κ) \in {(4, 3 k - 2), (5, 3 k - 1), (6, 3 k)} \\ 0 & elsewhere \end{matrix} \end{matrix}

where $j$ and $k$ are the numbers of node 1 and 2, respectively, in the global enumeration of the discretization. The transformation has the form

(91)

K_{i}^{e} = B_{i}^{T} K_{i} B_{i}

and the full stiffness matrix in global coordinates is then computed as the sum of all the extended matrices:

(92)

K = \sum_{i} K_{i}^{e}

Similarly, all the nodal displacements in global coordinates are collected in a vector $u$ ; a vector $F$ contains the nodal forces. It must be stated now that the obtained stiffness matrix has determinant equal to zero (i.e., it is not invertible): this is due to the fact that all nodes are considered, so far, both free and constrained. In order to account for the constraint conditions, the order of the DOFs should be manipulated such that it is possible to separate free from constrained quantities and write the following equation:

(93)

[\begin{matrix} K_{ff} & K_{fc} \\ K_{cf} & K_{cc} \end{matrix}] {\begin{matrix} u_{f} \\ u_{c} \end{matrix}} = {\begin{matrix} F_{f} \\ F_{c} \end{matrix}}

where the subscripts f and c refer to free and constrained DOFs, respectively; the vector $F_{c}$ contains the components of the reaction forces and moments. In our work, it is $u_{c} = 0$ . Therefore, it is $K_{ff} u_{f} = F_{f}$ . For a given set of external nodal forces, the nodal displacements are computed as

(94)

u_{f} = K_{ff}^{- 1} F_{f}

which is the very well-known Hooke’s law. In order to account for large displacements, several approaches can be adopted (see McGuire and Gallagher, 1979); here, we apply the external load in a finite number of steps $n_{s}$ , as a constant increment $Δ F_{f} = F_{f} / n_{s}$ , computing the increment of nodal displacements at the $s$ th step as

(95)

Δ u_{f (s)} = K_{ff (s)}^{- 1} Δ F_{f}

where the stiffness matrix $K_{ff (s)}$ is updated based on the nodal displacements computed at step $(s - 1)$ th. In our implementation, the stiffness matrix is updated at each step. The low number of elements in the discretization helps to keep low the computational cost; on the other hand, a high number of steps is generally required.

In the initial configuration of the bodies considered in Sections 5.2.2 and 6.3, the orientation of the $i$ th local system of coordinates and the one of the global system are coincident. Therefore, at the first step, the rotation matrix in Equation (85) is the identity matrix. At all the other steps, the stiffness matrix is updated based on the computed displacements; at the end $s$ th step, the global coordinates of the $j$ th node are given by

(96)

\begin{matrix} x_{j}^{(s)} = x_{j}^{(s - 1)} + Δ x_{j}^{(s)} \\ y_{j}^{(s)} = y_{j}^{(s - 1)} + Δ y_{j}^{(s)} \end{matrix}

and the orientation of the $i$ th element is computed as

(97)

θ_{i}^{(s)} = \arctan \frac{y_{i + 1}^{(s)} - y_{i}^{(s)}}{x_{i + 1}^{(s)} - x_{i}^{(s)}}

being, in the global enumeration adopted, $j = i$ . This angle is used to update the stiffness matrix used in the subsequent step:

(98)

K_{i}^{(s + 1)} = Q_{i} (θ_{i}^{(s)}) K_{i}^{(s)} Q_{i} (θ_{i}^{(s)})^{T}

(for a matter of convenience, the subscript f has been dropped in the expression above).

Acknowledgements

We are sincerely grateful to the three anonymous reviewers for their helpful comments and suggestions.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was funded by GrowBot, the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement No. 824074.

ORCID iD

Giovanna A Naselli

References

Albu-Schäffer

Bicchi

(2016) Actuators for Soft Robotics. Cham: Springer International Publishing, pp. 499–530.

Ananthasuresh

G. K.

(2013). How Far are Compliant Mechanisms from Rigid-body Mechanisms and Stiff Structures? Heidelberg: Springer International Publishing, pp. 83–94.

Armanini

(2018) Instabilities and Dynamics of Elastic Rods in the Presence of Movable Constraints. PhD thesis, University of Trento.

Armanini

Dal Corso

Misseroni

Bigoni

(2017) From the elastica compass to the elastica catapult: an essay on the mechanics of soft robot arm. Proceedings of the Royal Society A 473(2198): 20160870.

Arruda

Boyce

(1993) A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids 41(2): 389–412.

Bao

Fang

Chen

, et al. (2018) Soft robotics: Academic insights and perspectives through bibliometric analysis. Soft Robotics 5(3): 229–241.

Beer

Johnston

DeWolf

(2002) Mechanics of Materials. New York: McGraw-Hill.

Briot

Pashkevich

Chablat

(2009) On the optimal design of parallel robots taking into account their deformations and natural frequencies. In: ASME 2009 International Design Engineering Technical Conferences and Computers and Information in Engineering ConferenceSan Diego, CA, 30 August-–2 September 2009 (Volume 7: 33rd Mechanisms and Robotics Conference, Parts A and B). New York: ASME, pp. 367–376.

Cammarata

(2012) On the stiffness analysis and elastodynamics of parallel kinematic machines. In: Kucuk

(ed.), Serial and Parallel Robot Manipulators. Rijeka: IntechOpen.

10.

Cao

Dolovich

Schwab

Herder

Zhang

(2015) Toward a unified design approach for both compliant mechanisms and rigid-body mechanisms: Module optimization. Journal of Mechanical Design 137(12): 122301.

11.

Connolly

Polygerinos

Walsh

Bertoldi

(2015) Mechanical programming of soft actuators by varying fiber angle. Soft Robotics 2(1): 26–32.

12.

Connolly

Walsh

Bertoldi

(2017) Automatic design of fiber-reinforced soft actuators for trajectory matching. Proceedings of the National Academy of Sciences 114(1): 51–56.

13.

Duriez

(2013) Control of elastic soft robots based on real-time finite element method. In 2013 IEEE International Conference on Robotics and Automation, pp. 3982–3987.

14.

Elsayed

Vincensi

Lekakou

, et al. (2014) Finite element analysis and design optimization of a pneumatically actuating silicone module for robotic surgery applications. Soft Robotics 1(4): 255–262.

15.

Faure

Duriez

Delingette

, et al. (2012) SOFA: A Multi-Model Framework for Interactive Physical Simulation. Berlin: Springer, pp. 283–321.

16.

Felippa

(2001) A historical outline of matrix structural analysis: A play in three acts. Computers & Structures 79(14): 1313–1324.

17.

Ficuciello

Migliozzi

Coevoet

Petit

Duriez

(2018) FEM-based deformation control for dexterous manipulation of 3D soft objects. In: IROS 2018 - IEEE/RSJ International Conference on Intelligent Robots and Systems, Madrid, Spain.

18.

Germain

Briot

Caro

Wenger

(2015) Natural frequency computation of parallel robots. Journal of Computational and Nonlinear Dynamics 10: 021004.

19.

Grazioso

Di Gironimo

Siciliano

(2018) A geometrically exact model for soft continuum robots: The finite element deformation space formulation. Soft Robotics, in press.

20.

Grazioso

Sonneville

Gironimo

Bauchau

Siciliano

(2016) A nonlinear finite element formalism for modelling flexible and soft manipulators. In: 2016 IEEE International Conference on Simulation, Modeling, and Programming for Autonomous Robots (SIMPAR), pp. 185–190.

21.

Howell

(2001) Compliant Mechanisms. New York: John Wiley & Sons, Inc.

22.

Howell

Midha

(1994) A method for the design of compliant mechanisms with small-length flexural pivots. Journal of Mechanical Design 116(1): 280–290.

23.

Howell

Midha

(1995) Parametric deflection approximations for end-loaded, large-deflection beams in compliant mechanisms. Journal of Mechanical Design 117(1): 156–165.

24.

Howell

Midha

(1996) A loop-closure theory for the analysis and synthesis of compliant mechanisms. Journal of Mechanical Design 118(1): 121–125.

25.

Howell

Midha

Norton

(1996) Evaluation of equivalent spring stiffness for use in a pseudo-rigid-body model of large-deflection compliant mechanisms. Journal of Mechanical Design 118(1): 126–131.

26.

Largilliere

Verona

Coevoet

Sanz-Lopez

Dequidt

Duriez

(2015) Real-time control of soft-robots using asynchronous finite element modeling. In ICRA 2015, Seattle, WA.

27.

Laschi

Mazzolai

Cianchetti

(2016) Soft robotics: Technologies and systems pushing the boundaries of robot abilities. Science Robotics 1(1): eaah3690.

28.

Lipson

(2014) Challenges and opportunities for design, simulation, and fabrication of soft robots. Soft Robotics 1(1): 21–27.

29.

Manti

Hassan

Passetti

D’Elia

Laschi

Cianchetti

(2015) A bioinspired soft robotic gripper for adaptable and effective grasping. Soft Robotics 2(3): 0009.

30.

Martinez

Branch

Fish

, et al. (2013) Robotic tentacles with three-dimensional mobility based on flexible elastomers. Advanced Materials 25(2): 205–212.

31.

McGuire

Gallagher

(1979) Matrix Structural Analysis. New York: John Wiley & Sons, Inc.

32.

Mooney

(1940) A theory of large deformation. Journal of Applied Physics 11: 582–592.

33.

Moseley

Florez

Sonar

Agarwal

Curtin

Paik

(2016) Modeling, design, and development of soft pneumatic actuators with finite element method. Advanced Engineering Materials 18(6): 978–988.

34.

Ogden

(1972) Large deformation isotropic elasticity - on the correlation of theory and experiment for incompressible rubberlike solids. Proceedings of the Royal Society London A: Mathematical, Physical and Engineering Sciences 326(1567): 565.

35.

Pashkevich

Chablat

Wenger

(2009) Stiffness analysis of overconstrained parallel manipulators. CoRR, abs/0904.0052.

36.

Polygerinos

Lyne

Wang

, et al. (2013) Towards a soft pneumatic glove for hand rehabilitation. In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 1512–1517.

37.

Pucheta

Cardona

(2010) Design of bistable compliant mechanisms using precision–position and rigid-body replacement methods. Mechanism and Machine Theory 45(2): 304–326.

38.

Renda

Boyer

Dias

JMM

Seneviratne

(2017) Discrete cosserat approach for multi-section soft robots dynamics. CoRR, abs/1702.03660.

39.

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.

40.

Renda

Giorgio Serchi

Boyer

Laschi

Dias

Seneviratne

(2018) A unified multi-soft-body dynamic model for underwater soft robots. The International Journal of Robotics Research 37(6): 648–666.

41.

Rognant

Courteille

Maurine

(2010) A systematic procedure for the elastodynamic modeling and identification of robot manipulators. IEEE Transactions on Robotics 26(6): 1085–1093.

42.

Runge

Peters

Raatz

(2017) Design optimization of soft pneumatic actuators using genetic algorithms. In: 2017 IEEE International Conference on Robotics and Biomimetics (ROBIO).

43.

Rus

Tolley

(2015) Design, fabrication and control of soft robots. Nature 521: 467.

44.

Shapiro

Gabor

Wolf

(2015) Modeling a hyperflexible planar bending actuator as an inextensible Euler–Bernoulli beam for use in flexible robots. Soft Robotics 2(2): 71–79.

45.

Suzumori

Endo

Kanda

Kato

Suzuki

(2007) A bending pneumatic rubber actuator realizing soft-bodied manta swimming robot. In Proceedings 2007 IEEE International Conference on Robotics and Automation, pp. 4975–4980.

46.

Trimmer

Vanderborght

Mengüç

Tolley

Schultz

(2015) Soft robotics as an emerging academic field. Soft Robotics 2(4): 131–134.

47.

Wang

Iida

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

48.

Wolf

Grioli

Eiberger

, et al. (2016) Variable stiffness actuators: Review on design and components. IEEE/ASME Transactions on Mechatronics 21(5): 2418–2430.

49.

Zhang

Chen

(2013) A comprehensive elliptic integral solution to the large deflection problems of thin beams in compliant mechanisms. Journal of Mechanisms and Robotics 5(2): 021006.

50.

Zhang

Zhao

Jin

(2015) Elastodynamic modeling and analysis for an exechon parallel kinematic machine. Journal of Manufacturing Science and Engineering 138(3): 031011.

51.

Zhou

Majidi

O’Reilly

(2015) Flexing into motion: A locomotion mechanism for soft robots. International Journal of Non-Linear Mechanics 74: 7–17.

The softness distribution index: towards the creation of guidelines for the modeling of soft-bodied robots

Abstract

Keywords

1. Introduction

2. Stiffness and softness

2.1. Stiff or soft?

2.2. How to choose the modeling technique

3. SDI for planar bending

4. Modeling techniques

4.1. Finite elements

4.2. NMSA

4.3. EA

4.4. LPs

5. A soft notched body made of uniform material

5.1. SDI for the notched body

5.2. Modeling

5.2.1. FEA

5.2.2. NMSA

5.2.3. EA

5.2.4. LPs

5.3. Results

5.4. A further investigation on notched bodies

6. SDI for particular cases

6.1. Body made of multiple materials

6.2. Bodies with χ = 0

6.3. Soft continuum arm with uniform cross section

6.4. Bodies with 0 < χ < 1

6.5. Bodies with χ < 0

7. SDI for optimization

7.1. Statement of the problem

7.2. Optimization

7.3. The role of SDI in optimization

7.4. Optimization constraints based on SDI

8. Discussion

9. Conclusions

Footnotes

Appendix A: Invariance of SDI under inversion of x -axis

Appendix B: NMSA

Acknowledgements

Funding

ORCID iD

References

6.2. Bodies with $χ = 0$

6.4. Bodies with $0 < χ < 1$

6.5. Bodies with $χ < 0$