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.

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?

Fig. 1.

From rigid to soft robotic system, passing through intermediate systems with gradually distributed compliance, and viceversa. The dashed arrow represents the giant leap from a rigid system provided with a rotational joint (a) to a continuum, compliant system (b). Full arrows denote the transformation from (b) to continuum systems such as that in (d) passing through an infinite number of systems (c) in which the structural properties are more and more uniformly distributed; continuing with this transformation (e), the final result is a continuum, soft link (f) unequivocally. The same reasoning applies in the reverse way: from the distributed compliance of system (f), systems with more and more concentrated compliance are obtained up to the hinge-like mechanism (b), and a final giant leap consists of jumping to the non-continuum rigid-linked robot (a).

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\left(x\right)=E\left(x\right)J\left(x\right)$$

(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 ${\mathcal{B}}_i$ , and it will be considered an element of the set $\mathcal{B}$ , containing all the bodies performing planar bending.

Fig. 2.

Beam with variable cross section and material along the longitudinal axis $X$ .

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

$$\alpha \left(x\right)=\frac{d{k}_f(x)}{\mathrm{dx}}=E\left(x\right)\frac{\mathrm{dJ}\left(x\right)}{\mathrm{dx}}+J\left(x\right)\frac{\mathrm{dE}\left(x\right)}{\mathrm{dx}}$$

(2)

We denote by $x_M$ and $x_m$ the coordinates at which $\alpha (x)$ has global maximum ${\alpha }_M$ and global minimum ${\alpha }_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

$$\mathrm{SDI}=\chi +i\psi$$

(3)

where it is

$$\chi =\{\begin{array}{c}\begin{array}{cc}\frac{{\alpha }_M-{\alpha }_m}{\max \{|{\alpha }_M|,|{\alpha }_m|\}}·\frac{L}{x_M-x_m}\hfill & \mathrm{if}\alpha (x)\ne \mathrm{const}.\hfill \\ \multicolumn{1}{c}{0}& \mathrm{if}\alpha (x)=\mathrm{const}.\hfill \end{array}\hfill \end{array}$$

(4)

and

$$\psi =\frac{k_{f_M}-k_{f_m}}{k_{f_M}}·\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}\ge k_{f_m}$ and, therefore, $\psi \ge 0$ .

The SDI can be computed for any body ${\mathcal{B}}_i$ of the set $\mathcal{B}$ ; it can be seen as a function $\mathrm{SDI}:\mathcal{B}\to \mathbb{C}$ . This function is not injective (different bodies can be characterized by the same SDI) and not surjective in $\mathbb{C}$ (in fact, there is no body having $\psi <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 ${\mathcal{B}}_i$ and a modeling technique ${\mathcal{T}}_{\iota }$ 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.

Fig. 3.

The SDI associates the body ${\mathcal{B}}_i$ with the pair $(\chi ,\psi )$ denoting a point on the complex plane $\mathbb{C}$ ; the further step is to associate such point with a modeling technique, ${\mathcal{T}}_{\iota }$ , which allows the deformation of the body to be modeled with the required accuracy and with the fewest 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}{\rho (s)}=\frac{d\theta }{\mathrm{ds}}=\frac{M\left(s\right)}{E\left(s\right)J\left(s\right)}$$

(6)

in which $s$ denotes the curvilinear coordinate along the axis of the beam. By $\rho$ and $\theta$ 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\left(x\right)}{\mathrm{EJ}}=\frac{\frac{d^{2}y}{d{x}^2}}{\sqrt{{(1+{\left(\frac{\mathrm{dy}}{\mathrm{dx}}\right)}^2)}^3}}$$

(7)

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

$$\frac{M\left(x\right)}{\mathrm{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\left(x\right)$ under given boundary conditions. This equation is widely used for a variety of classical mechanical and civil engineering problems, in which the product $\mathrm{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{\mathrm{dx}}{\mathrm{ds}}=\cos \theta$$

(9)

$$\frac{\mathrm{dy}}{\mathrm{ds}}=\sin \theta$$

(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{array}{c}h(x)=h_M+\frac{1}{2}(h_M-h_m)(\cos \frac{n\pi (x-\frac{L}{2}-\frac{L}{n})}{L}-1)·\hfill \\ H(x-\frac{L}{2}-\frac{L}{n})H(\frac{L}{2}+\frac{L}{n}-x)\hfill \end{array}$$

(11)

in which $H(·)$ 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=2L/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$ :

$$\alpha \left(x\right)=E\frac{\mathrm{dJ}\left(x\right)}{\mathrm{dx}}$$

(12)

Fig. 4.

Notched soft body made of one material and with variable cross section along the $x$ -axis. Dimensions are normalized with respect to $L$ . The distance from the neutral plane is denoted by $h/(2L)$ . Here $L\prime =L_n/L$ is the normalized length of the notch.

Fig. 5.

Bodies with variable cross section characterized by different $n$ . The smaller $n$ , the more distributed compliance along the axis of the body.

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 $\alpha (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, ${\alpha }_m=-534.2$ Nmm, and ${\alpha }_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 $\mathrm{SDI}=15.26+i0.14$ .

Fig. 6.

Notched body (top) with $n=10$ and $h_m=8$ mm, its function $k_f(x)$ (middle), and the derivative $\alpha (x)$ (bottom), whose maximum and minimum are denoted by red dots.

The values of $n$ and $h_m/h_M$ do not influence the trend of the functions $k_f(x)$ and $\alpha (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, $\chi$ and $\psi$ ; Figure 7 plots the functions $\chi$ and $\psi$ in the intervals reported in Table 1.

Fig. 7.

Real part of the SDI versus $n$ and $h_m/h_M$ (top) and imaginary part versus $h_m/h_M$ (bottom).

We observe that $\chi$ ranges between 2.86 and 29.63 in the considered domain and $\psi$ 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 $\chi$ and $\psi$ compared with a body with a narrow and deep notch, body n.3.

Fig. 8.

The SDI is computed for nine bodies whose notches are determined by different $n$ and $h_m/h_M$ (left) and denoted in the complex plane by a dot (right).

5.2. Modeling

We study the deformation of this class of bodies owing to a constant external moment $\overline{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 $\overline{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 $\nu =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 ( ${\mathbf{u}}_1=0$ ) and the end moment load is applied on node 7, as a constant increment $\mathrm{\Delta }M=\overline{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 $\overline{M}$ ; therefore, Equation (6) takes the form

$$\frac{d\theta }{\mathrm{ds}}=\frac{\overline{M}}{\mathrm{EJ}\left(s\right)}$$

(13)

and the boundary conditions are as follows:

$$\begin{array}{c}\theta \left(0\right)=0\hfill \\ \multicolumn{1}{c}{x\left(0\right)=0}\\ \multicolumn{1}{c}{y\left(0\right)=0}\end{array}$$

(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 $\theta \left(s_i\right)$ for $i=1,2,\dots ,200$ as

$$\theta \left(s_i\right)={\int }_0^{s_i}\frac{\overline{M}}{\mathrm{EJ}\left(\tilde{s}\right)}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{{\mathrm{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_{\mathrm{eq}}}$$

(17)

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

$$l_{\mathrm{eq}}=|x_F-x_I|$$

(18)

being $[x_I,x_F]$ the interval in which it is

$$h\left(x\right)\le 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

$${\theta }_U=\frac{\overline{M}}{k_U}$$

(20)

$${\theta }_N=\frac{\overline{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 ${\theta }_U$ as well, because the segments $\overline{\mathrm{OP}}$ and $\overline{\mathrm{RS}}$ have the same structural properties. Here, $P$ and $Q$ denote the extrema of the notch, whereas $Q$ is the midpoint.

Fig. 9.

The beam is represented in two segments, denoted by the dashed broken line; the coordinates of point $S$ are calculated accounting for the relative rotation between segments, which depend on the lumped torsional springs at points $O$ , $P$ , and $R$ .

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

$$y_P=\frac{\overline{M}{x}_1^2}{2E{J}_U}$$

(22)

the slope of the line passing through $O$ and $P$ can be calculated as $\arctan (y_P/x_1)=\overline{M}/2k_U={\theta }_U/2$ . Analog considerations hold for segments $\overline{\mathrm{PQ}}$ , $\overline{\mathrm{PR}}$ , and $\overline{\mathrm{RS}}$ . The three springs represented in Figure 9 have stiffness coefficients $k_O=k_R=2k_U$ and $k_P=2k_N$ . The coordinates of $Q$ and $S$ can be expressed as functions of ${\theta }_U$ and ${\theta }_N$ as follows:

$$\begin{array}{c}{x}_Q=x_1\cos \frac{{\theta }_U}{2}+\frac{L_n}{2}\cos ({\theta }_U+\frac{{\theta }_N}{4})\hfill \\ \multicolumn{1}{c}{y_Q=x_1\sin \frac{{\theta }_U}{2}+\frac{L_n}{2}\sin ({\theta }_U+\frac{{\theta }_N}{4})}\end{array}$$

(23)

$$\begin{array}{c}{x}_S=x_R+(L-x_2)\cos ({\theta }_U+{\theta }_N+\frac{{\theta }_U}{2})\hfill \\ \multicolumn{1}{c}{y_S=y_R+(L-x_2)\sin ({\theta }_U+{\theta }_N+\frac{{\theta }_U}{2})}\end{array}$$

(24)

being

$$\begin{array}{c}{x}_R=x_1\cos \frac{{\theta }_U}{2}+L_n\cos ({\theta }_U+\frac{{\theta }_N}{2})\hfill \\ \multicolumn{1}{c}{y_R=x_1\sin \frac{{\theta }_U}{2}+L_n\sin ({\theta }_U+\frac{{\theta }_N}{2})}\end{array}$$

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

Fig. 10.

View on the $(x,y)$ plane of the nodes in the undeformed (blue) and deformed (red) configuration of the notched body with $n=2$ and $h_m/h_M=0.5$ as computed by FEA. The nodes on the neutral plane are highlighted in black.

Fig. 11.

Results obtained by the four techniques for the body with $n=2$ and $h_m/h_M=0.5$ . The blue solid line corresponds to the black line in Figure 10.

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 $\mathbf{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 ${\mathbf{p}}_{\mathrm{FEA}}^{(\mathrm{ij})}$ and the position ${\mathbf{p}}_{\iota }^{(\mathrm{ij})}$ computed by the technique ${\mathcal{T}}_{\iota }$ . 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_{\mathrm{FEA}}^{(\mathrm{ij})}=||{\mathbf{p}}_{\mathrm{FEA}}^{(\mathrm{ij})}-{\mathbf{p}}_0||$$

(26)

where ${\mathbf{p}}_0=\{L0{\}}^T$ is the position vector in the undeformed configuration. We compute the relative error for the technique ${\mathcal{T}}_{\iota }$ as

$$e_{\iota }^{(\mathrm{ij})}=\frac{||{\mathbf{p}}_{\iota }^{(\mathrm{ij})}-{\mathbf{p}}_{\mathrm{FEA}}^{(\mathrm{ij})}||}{d_{\mathrm{FEA}}^{(\mathrm{ij})}}$$

(27)

and $\iota$ 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 $\chi$ and $\psi$ .

Fig. 12.

Error $e_{\iota }$ for $h_m/h_M$ equal to 0.25, 0.50, and 0.75 computed using EA (blue solid line), NMSA (red dashed line), and LPM (yellow dash-dotted line).

Fig. 13.

Rectangular region of $\mathbb{C}$ , covered by colored patches that represent different modeling techniques. The fully colored areas correspond to an error up to 10% if the denoted technique (shown in the legend) is employed; the dotted or hatched regions correspond to an error up to 20%. The lines that separate adjacent regions are dashed, meaning that they have been sketched without analysing all the infinite points on the plane (corresponding to an infinite number of bodies), and therefore they cannot be considered sharp boundaries between the regions. Numbers from 1 to 9 refer to the SDI for the nine bodies considered in Figure 8.

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_{\mathrm{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 $\psi$ (<0.15) are better described by EA compared with bodies with greater $\psi$ . 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{array}{cc}\multicolumn{1}{c}{h(x)=}& h_M+f_1(x)·H(\overline{x}-x)·H(x-\overline{x}+a{\prime }_1)\hfill \\ \multicolumn{1}{c}{}& +f_2(x)·H(x-\overline{x})·H(-x+\overline{x}+a{\prime }_2)\hfill \end{array}$$

(28)

in which

$$f_i(x)=c·h_M-2\overline{a}\sqrt{1-\frac{{(x-\overline{x})}^2}{a_i^2}}$$

(29)

$$a_i\prime =a_i\sqrt{1-{\left(\frac{c·h_M}{2\overline{a}}\right)}^2}$$

(30)

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

$$\overline{a}=\frac{1}{2}[(1+c)h_M-h_m]$$

(31)

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

Fig. 14.

Body with elliptical notch. The notch is obtained from two different ellipses, sharing the semiaxis $\overline{a}$ . Point $C$ is the center of the ellipses. The coordinates of $C$ are $(\overline{x},\frac{c·h_M}{2})$ . Below the profile, the function $\alpha (x)$ is shown, for comparison with that in Figure 6. Dimensions of the body are reported in Table 2 and in Table 3, where the body is denoted by ${\mathcal{B}}_A$ .

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
$\overline{x}$	Abscissa at which it is $h=h_m$	75 mm
$a_2$	Semiaxis of the ellipse at the right of $\overline{x}$	5.2 mm

Table 3.

Minimum height $h_m$ , semiaxis $a_1$ , applied moment $\overline{M}$ , and SDI for the seven bodies considered.

Body, $\mathcal{B}$	$h_m$	$a_1$	$\overline{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 ${\mathcal{B}}_A$ and ${\mathcal{B}}_B$ we should expect an error between 10% and 20%, whereas for body ${\mathcal{B}}_D$ we should expect an error close to 10%. Results confirm our expectations. This holds for all the other bodies as well: for ${\mathcal{B}}_F$ and ${\mathcal{B}}_G$ the error is less than 10%; for ${\mathcal{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 ${\mathcal{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_{\iota }$ for the seven bodies listed in Table 3, obtained by EA, NMSA, and LPM.

Body, $\mathcal{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_{\mathrm{LPM}}$ for body ${\mathcal{B}}_C$ to be reduced from 11.98% to 4.09%. Moreover, it is worth noting that for ${\mathcal{B}}_G$ the error is 10%: in some cases, LPM may perform with such accuracy also for bodies with $\psi ~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 $\chi =0$ , $0<\chi <1$ , and $\chi <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\ne 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 $\chi =0$

A subset of bodies is characterized by $\chi =0$ . This condition occurs if it is $\alpha (x)=\alpha =\mathrm{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\left(x\right)=\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{\mathrm{Eb}}{12}(h_{1}x+h_2)$$

(34)

and, therefore,

$$\alpha (x)=\frac{\mathrm{Eb}{h}_1}{12}=\alpha$$

(35)

Under this condition, the first factor in the definition of $\chi$ is equal to zero; however, based on the provided definition, also the term $x_M-x_m$ is null. For this reason, $\chi$ 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 $\psi$ 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 ${\overline{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 $\mathrm{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=\pi r_a^2$ , and $J_i=\pi 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_{\mathrm{EA}}=0.003$ . NMSA is affected by a greater error ( $e_{\mathrm{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.

Fig. 15.

Continuum arm subject to end moment load. The blue and red dots represent the nodes in FEA in the undeformed and the deformed configuration, respectively. The solid black line denotes the solution computed by EA and the dashed black line the solution obtained by NMSA.

For this class of bodies, $\alpha (x)=0$ because $k_f(x)=\mathrm{const}.$ Therefore, $\chi =0$ and $\psi =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<\chi <1$

Here, we consider the classes of bodies characterized by $0<\chi <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\left(x\right)=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\left(x\right)=E\frac{\pi }{4}{r}^4\left(x\right)$$

(37)

The first derivative of the stiffness is

$$\alpha \left(x\right)=-E\pi \frac{r_0-r_L}{L}{(r_0-\frac{r_0-r_L}{L}x)}^3$$

(38)

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

$$\frac{d\alpha (x)}{\mathrm{dx}}=3E\pi {\left(\frac{r_0-r_L}{L}\right)}^2{(r_0-\frac{r_0-r_L}{L}x)}^2$$

(39)

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

• for $r_L<r_0$ , $\alpha (x)<0$ $\forall x$ ; therefore, $\alpha (L)={\alpha }_M$ , $\alpha (0)={\alpha }_m$ , and $|{\alpha }_m|>|{\alpha }_M|$ , thus $\frac{{\alpha }_M-{\alpha }_m}{|{\alpha }_m|}<1$ ;

• for $r_L>r_0$ , $\alpha (x)>0$ $\forall x$ ; therefore, $\alpha (L)={\alpha }_M$ , $\alpha (0)={\alpha }_m$ , and $|{\alpha }_M|>|{\alpha }_m|$ , thus $\frac{{\alpha }_M-{\alpha }_m}{|{\alpha }_M|}<1$ .

In both cases, $|x_M-x_m|=L$ . Hence, we conclude that it is always $0<\chi <1$ . Figure 16 shows $\chi$ and $\psi$ 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 $\mathrm{SDI}=0+i0$ , consistently with what was reported in the previous section for a soft continuum arm with uniform cross section.

Fig. 16.

Plots of $\chi$ and $\psi$ for conical bodies versus $r_L/r_0$ . The curves are computed for bodies with $E=10$ MPa, $r_0=13$ mm, and $r_L$ ranging in the interval $[1,26]$ mm.

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 $\overline{M}$ applied on conical bodies for varying $r_L$ and $r_0$ . The conical bodies are denoted by ${\mathcal{C}}_{\gamma }$ with $\gamma$ ranging from 1 to 12.

$\gamma$	$r_0$	$r_L$	$\overline{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 ${\mathcal{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 ${\mathcal{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.

Fig. 17.

Deformation of a conical body having $E=10$ MPa, $r_0=9$ mm, and $r_L=6$ mm under an applied moment at the free end section $\overline{M}=280$ Nmm. The undeformed conical profile is represented in blue, whereas the deformed as computed by FEA is shown in red. The solid and dashed black lines represent the deformed longitudinal axis as computed by EA and NMSA, respectively.

Table 7.

SDI, $e_{\mathrm{EA}}$ , and $e_{\mathrm{NMSA}}$ for conical bodies ${\mathcal{C}}_{\gamma }$ .

$\gamma$	SDI	$e_{\mathrm{EA}}$	$e_{\mathrm{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<\chi <1$ . Another example is represented by the class of bodies with circular cross section whose radius varies exponentially:

$$r(x)=R_0{e}^{b·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{\pi }{4}{R}_0^4{e}^{4b·x/L}$$

(41)

$$\alpha (x)=\frac{b}{L}E\pi R_0^4{e}^{4b·x/L}$$

(42)

$$\frac{d\alpha (x)}{\mathrm{dx}}=4\frac{b^2}{L^2}E\pi R_0^4{e}^{4b·x/L}$$

(43)

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

• for $b<0$ , $\alpha (x)<0$ everywhere; therefore, $\alpha (L)={\alpha }_M$ , $\alpha (0)={\alpha }_m$ , and $|{\alpha }_m|>|{\alpha }_M|$ (because $e^{-4|b|}<1$ );

• for $b>0$ , $\alpha (x)>0$ everywhere; therefore, $\alpha (L)={\alpha }_M$ , $\alpha (0)={\alpha }_m$ , and $|{\alpha }_M|>|{\alpha }_m|$ (because $e^{4b}>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 $\mathrm{SDI}=0.982+i0.127$ . For an applied moment $\overline{M}=150$ Nmm, we obtain $e_{\mathrm{EA}}$ =1.9% and $e_{\mathrm{NMSA}}$ = 2.1%. Results are displayed in Figure 18.

Fig. 18.

Results obtained by FEA, EA (solid line), and NMSA (dashed line) for a body with circular cross section the radius of which decreases exponentially along the length.

Fig. 19.

SDI for bodies with exponential profile (solid blue line) and conical bodies (dashed red line). Numerical values of material and geometrical parameters for which the curves are computed are reported in the text.

6.5. Bodies with $\chi <0$

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

$$\begin{array}{cc}\multicolumn{1}{c}{h(x)}& =h_M-(h_M-h_m)H(x-\frac{L}{3})\hfill \\ \multicolumn{1}{c}{}& +\frac{1}{2}(h_M-h_m)(\cos \frac{n\prime \pi (x-\frac{L}{3}-\frac{L}{n\prime })}{L}-1)·\hfill \\ \multicolumn{1}{c}{}& H(\frac{L}{3}-x)H(x-\frac{L}{3}+\frac{L}{n\prime })\hfill \end{array}$$

(44)

Again, $H(·)$ is the Heaviside function and the parameter $n\prime$ 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\prime =4$ , $h_M=16$ mm, $h_m=6$ mm, $b=10$ mm, and $L=100$ mm.

Fig. 20.

Body with variable cross section (top), characterized by the negative real part of the SDI. The first derivative of the flexural stiffness is plotted (bottom) and the minimum and the closest maximum are highlighted by red dots.

In this case, $\chi =-12.21$ and $\psi =0.15$ . The deformation of bodies of this class can be effectively modeled by EA, same as the notched bodies characterized by $\mathrm{SDI}=-\chi +i\psi$ ; 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 $\chi >0$ , instead, NMSA is suggested as a convenient technique for bodies having $\psi$ 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 $\mathrm{SDI}=-9.47+i0.09$ .

Fig. 21.

Body with decreasing height of the cross section (top) similar to the body in the previous figure but with different ratio $h_m/h_M$ , which leads to a different value of the minimum of the function $\alpha (x)$ (bottom) compared with the previous case.

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.

Fig. 22.

The body is clamped at the origin of the system of coordinates and its free end $P$ must be displaced to point $Q$ , under the loads $F_v$ and $M_a$ .

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

$${\mathbf{v}}_{\mathrm{LP}}=\{D,l,d\}$$

(45)

Fig. 23.

The body, initially straight, deforms under the applied moment $M_a$ and the force $F_v$ , owing to the deflection of the central segment between points $T$ and $S$ , having diameter $d$ (here not denoted).

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 {\theta }_{S_{s-1}}+f_{v}l)\frac{l^2}{2\mathrm{EJ}}-\frac{f_v{l}^3}{6\mathrm{EJ}}$$

(46)

with $a=L-D-L_p$ , $J=\pi d^4/64$ , $f_v=F_v/n_s$ , and $m_a=M_a/n_s$ . The rotation at point $S$ is computed as

$${\theta }_{S_s}={\theta }_{S_{s-1}}+(m_a+f_{v}a\cos {\theta }_{S_{s-1}}+f_{v}l)\frac{l}{\mathrm{EJ}}-\frac{f_v{l}^2}{2\mathrm{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

$${\mathbf{p}}_s=\{x_{S_s},y_{S_s}{\}}^T+\{a\cos {\theta }_{S_s},a\sin {\theta }_{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

$$\varphi \left(\mathbf{v}\right)=||{\mathbf{p}}_{n_s}-\mathbf{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

$${\mathbf{v}}_{0_{\mathrm{LP}}}=\{\frac{2L}{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\le d\le 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 ${\mathcal{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\prime$	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\prime$	Prescribed coordinates	$(187,-60)$ mm
$n_s$	Number of steps	40

Fig. 24.

Deformation of body ${\mathcal{B}}_1$ computed by FEA (blue dots represent the nodes of the mesh) and by LPM (black full line). FEA simulations reported in this section have been performed with a tetrahedral mesh.

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:

$${\mathbf{v}}_{\mathrm{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.

Fig. 25.

Schematic representation of the model based on NMSA. Nodes are numbered from 1 (clamped) to 6 (on which the loads are applied). The diameter and the length of elements between nodes 2 and 5 are treated as design variables, together with the distance $D$ between nodes 1 and 2.

The optimization is performed under the following constraints:

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

(55)

$$0<l_i<\frac{L}{3}\mathrm{for}i=1,2,3$$

(56)

$$d_m\le d_i\le d_M\mathrm{for}i=1,2,3$$

(57)

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

$${\mathbf{v}}_{0_{\mathrm{NMSA}}}=\{D^{(1)},\frac{l^{(1)}}{3},\frac{l^{(1)}}{3},\frac{l^{(1)}}{3},d_{\mathrm{ave}}^{(1)},d^{(1)},d_{\mathrm{ave}}^{(1)}\}$$

(58)

where $d_{\mathrm{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 ${\mathcal{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%.

Fig. 26.

Deformation of body ${\mathcal{B}}_2$ computed by NMSA (full black line) in comparison with the result obtained by FEA (blue dots represent the nodes of the mesh).

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\prime$ (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{7L}{100}$$

(60)

$$d_m\prime \le d\le d_M$$

(61)

and starting the search from the initial guess

$${\mathbf{v}}_0=\{\frac{L}{3},\frac{L}{20},d{\prime }_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\prime$ and $(185.0,-63.4)$ with the two methods, respectively, resulting in an error equal to 6.0%.

Fig. 27.

Deformation of body ${\mathcal{B}}_3$ computed by FEA (blue dots represent the nodes of the mesh) and by LPM (black full line).

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 ${\mathcal{B}}_2$ , we have performed a search in a seven-dimensional space, but for ${\mathcal{B}}_1$ and ${\mathcal{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 ${\mathcal{B}}_1$ and ${\mathcal{B}}_3$ have been found in less than 0.05 seconds, whereas the search for ${\mathcal{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 ${\mathcal{B}}_1$ , ${\mathcal{B}}_2$ , and ${\mathcal{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 $\alpha \left(x\right)$ recurring to finite differences. In the following, we discretize the interval $[0,L]$ with 1,000 evenly spaced points and we approximate

$$\alpha (x_i)=\frac{\mathrm{\Delta }{k}_{f_i}}{\mathrm{\Delta }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 $\alpha (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 $\mathrm{SDI}=7.6+i0.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 $\mathrm{SDI}=31.4+i0.0995$ : for bodies having such high values of the real part of the SDI and $\psi ~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.

Fig. 28.

Profiles of the bodies ${\mathcal{B}}_1$ , ${\mathcal{B}}_2$ , and ${\mathcal{B}}_3$ , from top to bottom. Below each body, the corresponding function $\alpha$ is shown.

Table 9.

Summary: SDI and error computed for each body.

Body	SDI	Error
${\mathcal{B}}_1$	$7.6+i0.0992$	37.0% by LPM
${\mathcal{B}}_2$	$4.5+i0.0992$	21.9% by NMSA
${\mathcal{B}}_3$	$31.4+i0.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 $\overline{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

$$\mathbf{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 $\overline{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

$$\phi (\mathbf{w})=y_f-y_P$$

(65)

and we perform the optimization imposing the following constraints:

$$\frac{h_M}{10}<h<\frac{9h_M}{10}$$

(66)

$$2<n<14$$

(67)

starting the search from

$${\mathbf{w}}_0^{(1)}=\{\frac{3h_M}{20},8\}$$

(68)

The solver succeeds to find a local minimum while satisfying the assigned constraints: at the point ${\mathbf{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:

$${\mathbf{w}}_0^{(2)}=\{\frac{h_M}{4},11.5\}$$

(69)

Fig. 29.

Comparison between FEA and EA for the three bodies obtained by optimization: (a) $n=7.68$ and $h_m=4.84$ mm; (b) $n=11.79$ and $h_m=4.06$ mm; (c) $n=4.00$ and $h_m=6.31$ mm. Blue and red dots denote the nodes of the mesh for the undeformed and deformed body, respectively. The black line denotes the deformed axis of the body as computed by EA.

We find that the value of the objective function for ${\mathbf{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 ${\mathbf{w}}_0^{(2)}$ is a more unfortunate initial guess than ${\mathbf{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_{\mathrm{EA}}$ . Figure 30 discourages us from repeating the optimization with a different initial guess: the function $\phi$ as computed by EA is characterized by a wavy trend and its intersection with the plane $z=0$ consists of several curves (we show $\phi$ , 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:

$$\psi \left(\mathbf{w}\right)<0.15·H(\chi \left(\mathbf{w}\right)-6)$$

(70)

Fig. 30.

Function in Equation (65). The plot is found by means of implementation of EA as described in Section 5. The black transparent patch denotes the plane $z=0$ .

and we perform again the optimization starting the search from ${\mathbf{w}}_0^{(3)}={\mathbf{w}}_0^{(1)}$ . This time, we find that the objective function takes the value $1.2\times {10}^{-7}$ mm at ${\mathbf{w}}^{(3)}=\{6.31,4.00\}$ while satisfying all the assigned constraints. The body is characterized by $\chi =5.85$ and $\psi =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+i0.1555$ and $16.58+i0.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 ${\mathcal{B}}_2$ , NMSA with five elements performs with $e_{\mathrm{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.