The longest soft robotic arm

Journal Pre-proof The longest soft robotic arm András A. Sipos, Péter L. Várkonyi

PII: DOI: Reference:

S0020-7462(19)30443-3 https://doi.org/10.1016/j.ijnonlinmec.2019.103354 NLM 103354

To appear in:

International Journal of Non-Linear Mechanics

Received date : 20 June 2019 Revised date : 8 November 2019 Accepted date : 8 November 2019 Please cite this article as: A.A. Sipos and P.L. Várkonyi, The longest soft robotic arm, International Journal of Non-Linear Mechanics (2019), doi: https://doi.org/10.1016/j.ijnonlinmec.2019.103354. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Journal Pre-proof

The longest soft robotic arm Department of Mechanics, Materials and Structures, Budapest University of Technology and Economics, M˝ uegyetem rkp. 1-3. K261, Budapest 1111, Hungary b MTA-BME Morphodynamics Research Group, Budapest, Hungary

pro

a

of

Andr´as A. Siposa,b , P´eter L. V´arkonyia

Abstract

urn al P

re-

We investigate the maximal outreach of soft, elastic robotic arms with controllable intrinsic curvature by using the three-dimensional geometrically exact Cosserat rod theory and a linear elastica model. We compare results with previous predictions for rods without intrinsic curvature based on a planar rod theory. We point out that planar models often overestimate the maximal outreach achievable by arms not constrained to a plane per se, by neglecting instabilities associated with spatial configurations. At the same time, intrinsic curvature can be utilized to stabilize the rod. Here, we study the interplay between spatial buckling and intrinsic curvature on a clamped, prismatic rod carrying a vertical tip load. The numerical solution of the associated boundary value problem and discrete stability tests on the obtained equilibrium solutions reveal significant improvement in the maximal outreach, especially for arms with wide and flat cross-sections. Keywords: buckling, instability, optimization, soft robotic arm, intrinsic curvature 1. Introduction

Jo

Structural optimization [1, 2] often results in soft and slender structures exhibiting compliance-induced instability. In 1966, Keller and Niordson’s seminal paper The Tallest Column investigated the optimal profile of a tapered column subject to a combination of its own weight and a tip load. The height of the column was found to be limited by Euler buckling, i.e. the most common form of elastic instability in structural mechanics. Their investigations initiated lively mathematical debates [3, 4, 5, 6] and also inspired other researchers to solve a series of similar optimization problems Preprint submitted to International Journal of Non-linear Mechanics

November 8, 2019

Journal Pre-proof

re-

pro

of

in structural mechanics. Some of these works address the optimal profiles of columns [7, 8]. Others investigate beams subject to bending with focus on minimized tip deflection [9, 10], or on maximum safety against lateral torsional buckling [11]. The present paper investigates linear elastic, prismatic cantilever beams subject to a static vertical tip load F acting at the center-line of the rod. In this case, Euler buckling and lateral torsional buckling are the most common forms of destabilization. We use a geometrically exact elastica beam model assuming small strains but allowing large rotations, in which the rod is assumed to be inextensible and unshearable. The effect of warping is also neglected, which is reasonable for close-to-circular cross sections. For a horizontal beam with constant lateral bending stiffness B1 , in-plane bending stiffness B2 , torsional stiffness B3 , and arc length L. S. Timoshenko was the first to determine the condition of stability in the asymptotic limit B2 >> 1: L ≤ Rltb := χ(B1 B3 )1/4 /F 1/2 ,

(1)

urn al P

where χ ≈ 2.003 [12]. The value of χ, and hence of the critical length, becomes larger as B2 is decreased [12]. In particular, lateral torsional buckling is not possible at all when B2 < B1 (Fig. 1a). Nevertheless, lateral torsional buckling is not the only compliance-induced limitation for the extension of a soft cantilever. Several recent works investigated the in-plane deformations of soft cantilever beams (B2 small) under the assumption that B1 , B3 are large enough to eliminate lateral torsional buckling (Fig. 2a) The aim of these works was to determine the maximum ’reachout’ (i.e. distance between the endpoints) for a fixed value of B2 [13, 14, 15, 16]. Because of the in-plane deformations, the distance between the endpoints is significantly smaller than L in general. More specifically, it has been found that the height −x of the free end above the fixed end exhibits an upper bound −x ≤ Rb := π/2(B2 /F )1/2 ,

(2)

Jo

due to in-plane buckling even if L is arbitrary. Additionally the horizontal distance z between the endpoints also exhibits a nontrivial finite upper bound |z| ≤ Rdef := 2(B2 /F )1/2 ,

(3)

and a slightly more restrictive limit applies when we seek a purely horizontal reachout, i.e. x = 0 is required: ∗ |z| ≤ Rdef := µ(B2 /F )1/2 .

2

(4)

urn al P

re-

pro

of

Journal Pre-proof

Jo

Figure 1: (a): illustration of classical lateral torsional buckling under tip load (thick arrow) in the case of a high and narrow cross-section (i.e. B2 > B1 ). The unbuckled shape is in the vertical x − z plane marked by semi-transparent shading, and the strong axis with stiffness B2 (thin arrow at free end, red online) is parallel to the y axis. Torsional deformation of the beam rotates the strong axis and enables larger deflection at the tip. At the same time, the beam moves laterally in the y direction. (b): lateral torsional buckling of a beam with nonzero intrinsic curvature compensating for load-induced deflection. The unbuckled shape is in the vertical x − z plane, and the intrinsic curvature vector (thin arrow at free end, red online) is parallel to the y axis. Torsional deformation of the beam rotates the intrinsic curvature vector, and allows larger deflection at the tip due to the reduced y component of the intrinsic curvature. This is combined with sideways deformation in the y direction due to the emerging x component of the intrinsic curvature. The latter mechanism is active regardless of the ratio B2 /B1 .

3

Journal Pre-proof

450

- 2.5

1

-

0.5

x

-

0

-0.5 -1 -1.5

stable

360

un s

270

180

ta

stable

-2 -1

0

z (a)

1

2

90 90

bl

e

180

270

360

450

(L)at free end (°) tangent angle

(b)

re-

-2.5 -2

of

- 1.5

pro

tangent angle (0) at fixed end (°)

2

-

urn al P

Figure 2: (a): numerically found stable (thick curves, blue online) and unstable (thin curves, red online) equilibria of a planar cantilever of length L = 2.4m and stiffness B2 = 1Nm2 with a clamped end at point (0, 0) and a free end subject to a vertical, concentrated force F = 1N. Note that each value of the tangent angle φ(L) at the free end corresponds to a unique solution. (b): φ(L) as function of tangent angle φ(0) at the clamped end. Note that multiple solutions occur for some values of φ(0). This phenomenon occurs whenever L exceeds the first critical length Rb = π/2 of Euler buckling under pure compression. Instability may occur only when multiple solutions are present. Simulation methods are presented in Sec. 3.

Jo

where µ ≈ 1.81. Figure 2 depicts some deformed shapes and their in-plane stability by using our own numerical simulations. Figure 3 illustrates the known theoretical limits of reachout. Our aim here is to discuss the extension of previous results for cantilevers with nonzero intrinsic curvature, which has natural applications in the field of soft robotics [17, 18]. Soft arms or tentacles or various animals such as the octopus [19] are characterized by curved shapes controlled by internal muscles and subject to large deformations under external loads. Robotic arms of this type have become increasingly popular during the past decade [20]. Their main advantages include low mass, the potential of adaptation to a geometrically complex environment by passive compliance (especially useful in the case of invasive surgery [21] or rescue robots), safe operation in a space shared with humans and adaptability to motion of the soft human body (as

4

z=Rdef

z=-Rdef

-

* z=-Rdef

-x=Rb

-

urn al P

x

* z=Rdef

re-

-

pro

of

Journal Pre-proof

z

Jo

Figure 3: Possible positions of the free end of a cantilever with stiffness B2 = 1Nm2 subject to a vertical, concentrated tip load F = 1N. Dark (blue online) dots denote stable configurations for any length L, and light gray dots correspond to unstable configurations ∗ (for any L ≤ 6m). Dashed lines at x = −Rb and at z = ±Rdef , ±Rdef indicate the theoretical limits of reachout.

5

Journal Pre-proof

Jo

urn al P

re-

pro

of

in the case of exo-suits and medical rehabilitation robots). Shape control is most often realized through pneumatic actuators, dielectric elastomers or shape memory alloys. Unlike traditional robotic arms, whose performance is typically limited by mechanical strength or actuator capacity, soft arms (just like tall columns) may also fail to work properly due to large deformations or by elastic instability. Such phenomena can be examined by using elastic rod theory, but special care should be taken because of the presence of large deformations and geometric nonlinearity. Here, we adopt an idealized model of a soft robotic arm including controllable length, tangent angle at the fixed endpoint, and intrinsic curvature function along the bar. Our aim is to find the optimal performance of the arm by tuning the control parameters. Performance means largest possible value of the load for a given position of the free endpoint, or equivalently, the maximum reachout of the cantilever in a given direction for a given value of the load. We have seen that the performance of an elastica cantilever without intrinsic curvature (Figure 3) is limited primarily by bending-induced elastic deformations. The same mechanism is however not limiting the performance of an ideal robotic arm subject to curvature control. In this case, one can prescribe an arbitrary load along with an arbitrary smooth deformed shape (i.e. curvature). The corresponding bending moment function of the deformed arm can be determined uniquely. Then, the known stiffness B2 allows one to calculate a unique intrinsic curvature function along the cantilever, which deforms into the prescribed geometry under the effect of the prescribed load (Fig. 4). Based on this argument, one could naively conclude that perfect control of intrinsic curvature eliminates the compliance-induced performance limits of cantilevers. This is however false, since intrinsic curvatures compensating bending-induced deformations make the system susceptible to lateral torsional buckling, even when B1 and B3 are large (Fig. 1b). This observation was our motivation to systematically explore instability in the presence of intrinsic curvatures, and to identify the compliance-induced theoretical performance limits of elastica arms. In Sec. 2, we present the mechanical model of the arm. Sec. 3 is devoted to describe the applied numerical method. In Sec. 4, equilibrium configurations of the arm are found by solving the Euler-Lagrange equations of the problem numerically. The second variation of the potential energy of a 6

urn al P

re-

pro

of

Journal Pre-proof

Jo

Figure 4: Compensation of bending-induced deformations by intrinsic curvature. An arbitrary preferred loaded physical shape of the cantilever (curve a in left panel) can be achieved by considering the corresponding curvature function (curve b in right panel) and the bending moment function dictated by the loaded shape in order to determine the corresponding intrinsic curvature function (curve d). From that, the unloaded physical shape of the cantilever can be determined (curve e). For linear elastic beams, the intrinsic curvature is obtained as the difference between the actual curvature and the curvature induced by bending moments (curve c). The dimensions of length and curtvature are [m] and [1/m], respectivlely.

7

Journal Pre-proof

of

hyperelastic Cosserat rod is used for stability analysis and numerical optimization methods are used to tune the control parameters towards maximum performance. The paper is then closed by a Conclusion section. 2. Model development

re-

pro

2.1. General setting In what follows, we derive the geometrically exact governing equations of a rod with intrinsic curvature. Let {e1 , e2 , e3 } denote a fixed, right-handed orthonormal basis in R3 , with e1 called the ’vertical direction’ in the sequel. Note that the x, y and z axes of Figure 1 are associated with the basis vectors e1 , e2 and e3 , respectively. Let s ∈ [0, L] be the arc length parameter of the unloaded shape in the reference configuration with (.)0 denoting derivation with respect to s. Let r(s) denote the position vector of the material point s. The rotation tensor of a cross section R(s) defines the directors of the special Cosserat Rod theory [22] via i = 1, 2, 3.

urn al P

di = R(s)ei ,

(5)

One can show [23], that there is a unique vector field κ that fulfills d0i = κ × di , κ = κi di , r0 = νi di ,

i = 1, 2, 3,

(6) (7) (8)

Jo

where the Einstein summation is applied. Here νi and κi are the strains and the curvatures in the theory. In [23] a system of 13 first order ODEs are derived to express the variation of the internal forces and the shape along the rod under external loads. In this work we assume the rod to be unshearable and inextensible made of a hyperelastic material, hence ν1 = ν2 = 0 and ν3 = 1. Equivalently, the existence of a scalar valued, non-negative C 2 function W : R3 −→ [0, ∞) function called the stored energy density that depends on the curvature fields is postulated. Note that a formulation of the potential energy with unknown fields κ1 ,κ2 ,κ3 and the director fields would require extra Lagrange-multiplier fields to enforce (6). Instead of entering this cumbersome setting, based on (6) and

8

Journal Pre-proof

κ1 (s) = d02 · d3 = −d03 · d2 , κ2 (s) = d03 · d1 = −d01 · d3 , κ3 (s) = d01 · d2 = −d02 · d1 .

of

(7) the curvatures are expressed as: (9) (10) (11)

pro

With these in hands and employing the orthogonality of the director fields (d3 = d1 × d2 ) the stored energy density might be represented as ˜ (d1 , d2 ) = W (κ1 (d1 , d2 ), κ2 (d1 , d2 ), κ3 (d1 , d2 )). Now, the potential energy W functional L can be formulated as Z L ˜ (d1 , d2 ) − F (d2 × e1 · d1 − d0 · e1 )}ds, {W (12) L(d1 , d2 ) = 3 0

urn al P

re-

where the second term is the external work of the force F with d03 being the tangent vector-field associated with the initial (stress-free shape). Note that in this setting the unknown fields of the problem are d1 and d2 only. Furthermore, a fixed end at s = 0 suggests the initial values d1 (0) = p1 e1 + p2 e3 with p21 + p22 = 1 and d2 (0) = e2 . Let η1 and η2 denote admissible variations of d1 and d2 , respectively. Considering the unit lengths of d1 and d2 and their orthogonality it follows that η1 (s) = χ(s)d2 + ζ(s)d3 , η2 (s) = −χ(s)d1 + ψ(s)d3 ,

(13) (14)

where χ(s), ζ(s), ψ(s) are arbitrary C 1 ([0, L] → R) functions with χ(0) = ζ(0) = ψ(0) = 0. Let ω : R → R3 with ω(s) := [χ(s), ζ(s), ψ(s)]T .

(15)

Jo

Using these expressions the first variation of the energy is readily computed as Z L ∂W 0 δL(d1 , d2 )[ω] = (d2 · (η1 × d2 ) + η20 · (d1 × d2 ) + d02 · (d1 × η2 )) + ∂κ 1 0 ∂W ((d1 × d2 )0 · η1 + (η1 × d2 )0 · d1 + (d1 × η2 )0 · d1 ) + ∂κ2  ∂W 0 0 (η · d2 + d1 · η2 ) − F d2 × e1 · η1 − F η2 × e1 · d1 ds. ∂κ3 1 (16) 9

Journal Pre-proof

pro

of

Substitution of η1 (s) and η2 (s) yields Z L ∂W ∂W δL = (−ζd02 · d1 − χd01 · d3 + ψ 0 ) + (χd03 · d2 − ζ 0 − ψd02 · d1 ) + ∂κ1 ∂κ2 0  ∂W 0 0 0 (χ + ζd3 · d2 + ψd1 · d3 ) − F d2 × e1 · ζd3 − F ψd3 × e1 · d1 ds. ∂κ3 (17)

re-

Now the curvatures from eqs. (9-11) can be substituted. At an equilibrium solution the first variation vanishes, hence Z L ∂W ∂W (ζκ3 + χκ2 + ψ 0 ) + (−χκ1 − ζ 0 + ψκ3 ) + δL = ∂κ ∂κ 1 2 0  ∂W 0 (χ − ζκ1 − ψκ2 ) + ζF d1 · e1 + ψF e1 · d2 ds = 0. (18) ∂κ3

urn al P

Integration by parts and the fact, that the admissible variations ζ(s), χ(s) and ψ(s) are arbitrary deliver the Euler-Lagrange equations: 0  ∂W ∂W ∂W =F e1 · d2 + κ3 − κ2 , (19) ∂κ1 ∂κ2 ∂κ3  0 ∂W ∂W ∂W + κ1 , (20) = − F e1 · d1 − κ3 ∂κ2 ∂κ1 ∂κ3  0 ∂W ∂W ∂W = − κ1 + κ2 . (21) ∂κ3 ∂κ2 ∂κ1 Integration by parts in the previous step also delivers the following natural boundary conditions: ∂W ∂W ∂W (L) = (L) = (L) = 0. ∂κ1 ∂κ2 ∂κ3

(22)

Jo

Any solution in equilibrium fulfills eqs. (19-21), the boundary conditions in (22) and eqs.(6). Generally speaking, this system of (altogether 9) first order ODEs determines the shape of the spatial elastica. The second variation of the potential energy, evaluated at an equilibrium solution characterized by the curvature functions κ1 (s), κ2 (s), κ3 (s) and the director fields

10

Journal Pre-proof

d1 (s), d2 (s), d3 (s), reads

2

0

∂ W ∂κ23 ∂W ∂κ1 ∂W ∂κ2 ∂W ∂κ3

L



∂ 2W ∂ 2W 2 0 2 (ζκ + χκ + ψ ) + (−χκ1 − ζ 0 + ψκ3 ) + 3 2 2 2 ∂κ1 ∂κ2

of

δ L(d1 , d2 )[ω, ω] =

Z

2

(χ0 − ζκ1 − ψκ2 ) + F e1 · (ζχd2 − ψχd1 + ζ 2 d3 + ψ 2 d3 )+
ζχ0 − χζ 0 − ψζκ2 − ζ 2 κ1 − χ2 κ1 + χψκ3 +

pro

2


ψχ0 − χψ 0 − χζκ3 − ψ 2 κ2 + χ2 κ2 − ψζκ1 + 
0 0 2 2 ψζ − ζψ + ψχκ1 − ψ κ3 − ζ κ3 − χζκ2 ds. (23)

urn al P

re-

2.2. Governing equations for linear elastic rods with intrinsic curvature Let κ0i (s), i = 1, 2, 3 denote the intrinsic curvature functions for the stressfree (i.e. reference) shape. In this case the stored energy density W of a simple, linearly elastic material with bending and torsional rigidities B1 , B2 , B3 , respectively, can be written as W (κ1 , κ2 , κ3 ) =

3 X 1 i=1

2

Bi (κi − κ0i )2 .

(24)

Now the governing equations from (6), (19), (20) and (21) read d01 = κ2 (d2 × d1 ) + κ3 d2 , d02 = κ1 (d1 × d2 ) − κ3 d1 ,


0 B1 κ01 = F e1 · d2 + B2 κ3 (κ2 − κ02 ) − B3 κ2 (κ3 − κ03 ) + B1 κ01 ,
0 B2 κ02 = −F e1 · d1 − B1 κ3 (κ1 − κ01 ) + B3 κ1 (κ3 − κ03 ) + B2 κ02 ,
0 B3 κ03 = −B2 κ1 (κ2 − κ02 ) + B1 κ2 (κ1 − κ01 ) + B3 κ03 .

(25) (26) (27) (28) (29)

Jo

2.3. Equations for planar shapes Without loss of generality, planar solutions might be characterized by d2 ≡ e2 along with κ1 = κ3 = 0. It implies that κ01 ≡ 0 and κ03 ≡ 0 also hold. Regarding the shape itself, in the triple (x(s), y(s), z(s)) the component y(s) ≡ 0 as long as the shape is planar (Fig. 1). Without out-of-plane 11

Journal Pre-proof

d1 = [cos(φ), 0, − sin(φ)]T ,

of

deformations the rhs. of equations (26), (27) and (29) are identically zero. To simplify notation let the director d1 be represented via (30)

pro

where φ(s) is the turning angle between the local tangent d3 and axis e3 . Obviously, φ(s) uniquely determines x(s) and z(s). These equations and substitution of (30) into (25) and (28) yield the following system of ODEs:

re-

x0 = sin(φ), z 0 = cos(φ), φ0 = κ2 ,
0 F cos φ + κ02 . κ02 = − B2

(31) (32) (33) (34)

These four equations are equipped with the following boundary conditions: x(0) = 0,

z(0) = 0,

angle(x(L), z(L)) = α,

κ2 (L) = κ02 (L),

(35)

urn al P

where the angle(ξ, ζ) function stands for the angle between vectors (ξ, ζ), and (0, 1), and α is a prescribed value. The first two conditions prescribe the position of the fixed end, the third one determines the direction of the reachout, and the last one is the natural boundary condition from eq. (22). Note that the angle φ(0) at the fixed end has not been specified as it is determined implicitly by the third boundary condition. 2.4. Non-dimensional equations In order to simplify computation the following normalization is applied: ¯ := 1, ... s¯ := sL−1 , x¯ := xL−1 , z¯ := zL−1 , L ¯ := φ(s), κ φ(s) ¯ i := κi L, δ L¯ := δLL, δ 2 L¯ := δ 2 LL, B¯1 := B1 /B2 , B¯3 := B3 /B2 .

(36) (37) (38) (39)

Jo

Now s¯ ∈ [0, 1] and equations (31)-(34) take the dimensionless form ¯ x¯0 = sin(φ), ¯ z¯0 = cos(φ), φ¯0 = κ ¯2,


0 κ ¯ 02 = −γ cos φ¯ + κ ¯ 02 , 12

(40) (41) (42) (43)

Journal Pre-proof

where the non-dimensional parameter is defined via F L2 . B2

(44)

of

γ :=

pro

The second variation associated with the non-dimensional form of the governing equations reads Z 1n 2 2 2 2 ¯ ¯1 (χ¯ ¯3 (χ0 − ψ¯ δ L= B κ2 + ψ 0 ) + (ζ 0 ) + B κ2 ) + (45) 0
¯ + (¯ ¯ 2 + χ2 κ ¯ 2 Ld¯ s. γ(−ψχ cos φ¯ + ζ 2 sin φ¯ + ψ 2 sin φ) κ2 − κ ¯ 02 ) ψχ0 − χψ 0 − ψ 2 κ

x¯(0) = 0,

z¯(0) = 0,

re-

Finally, the boundary conditions in (35) can be expressed in dimensionless variables as angle(¯ x(1), z¯(1)) = α,

κ ¯ 2 (1) = κ ¯ 02 (1).

(46)

urn al P

2.5. Stability of planar shapes Stability of a loaded shape hinges on the positivity of the second variation ¯ of the potential energy in eq. (45). Note that φ(s) = φ(s) uniquely deter0 mines both κ ¯ 2 (s) and κ ¯ 2 (s). First, we demonstrate that there exists a unique γ = γcrit at which the second variation is positive semidefinite. Integrating eq. (43) and considering the natural boundary condition κ ¯ 2 (1) = κ ¯ 02 (1) yields Z s Z 1 0 ¯ ¯ κ ¯ 2 (s) − κ ¯ 2 (s) = −γ cos φ(t)dt + γ cos φ(t)dt = γ(−¯ z (s) + z¯(1)). 0

0

(47)

Substitution of (47) into (45) shows, that the second variation of the potential energy is a linear function of the parameter γ. It is straightforward ¯1 and B ¯3 to demonstrate, that due to the strict positivity of B lim δ 2 L¯ > 0

(48)

lim δ 2 L¯ < 0

(49)

γ→0

Jo

for all admissible variations and

γ→∞

for some admissible variations. Linearity in γ implies the existence of γ = γcrit at which there exist an admissible variation that makes δ 2 L¯ vanish. Next we 13

Journal Pre-proof

pro

of

show, that stability results with a planar model provides an upper bound for the maximal outreach as long as the loaded shape is fixed. Let Z 1 2 2 ¯ ¯ d s¯, {(ζ 0 ) + γζ 2 sin φ}L (50) δ Lp := 0 Z 1n 2 2 2 ¯ ¯1 (χ¯ ¯3 (χ0 − ψ¯ δ Ls := B κ2 + ψ 0 ) + B κ2 ) + (51) 0

γ −ψχ cos φ¯ + ψ 2 sin φ¯ + (¯ z (1) − z¯(s)) ψχ0 − χψ 0 − ψ 2 κ ¯ 2 + χ2 κ ¯2 Ld¯ s.

urn al P

re-

Obviously δ 2 L¯ = δ 2 L¯p + δ 2 L¯s . While δ 2 L¯p depends only on the test function ζ and δ 2 L¯s is independent of it, we conclude that the discretization of δ 2 L¯ should be a block-diagonal matrix. As ζ is an admissible variation of the shape in the [xz] plane, it follows that the block associated with δ 2 L¯p refers to planar perturbations. Similarly, the block that results from the discretization of δ 2 L¯s belongs to out-of-plane perturbations. As eigenvalues of a block diagonal matrix are the union of eigenvalues of the blocks, it follows, that a planar model provides an upper estimate on stability of an equilibrium shape. Furthermore, δ 2 L¯p is independent of the intrinsic curvature κ ¯ 02 . 3. Computational methods

All computations are carried out using the non-dimensional version of the governing equations (40)-(43). We exploit two important observations: ¯ • For a fixed dimensionless loaded shape (i.e. φ(s)) and parameter value 0 γ, the intrinsic curvature function κ ¯ 2 is unique as eq. (43) is a first order ODE in κ ¯ 02 (see also Fig. 4 for an intuitive explanation). This fact lead us to parametrize solutions by using the truncated Taylor expansion of the dimensionless loaded shape.

Jo

¯ • We have seen that for fixed φ(s) there exists a unique critical value of γ = γcrit at which the positivity of the second variation is lost as γ is ¯ increased. Accordingly, for any given choice of φ(s), we can efficiently determine the maximal reachout of the rod . The details of the computations are given below.

14

Journal Pre-proof

M X

c¯i s¯i−1 .

i=1

pro

¯ s) ≈ φ(¯

of

3.1. Parametrization of solutions ¯1 , and B ¯3 of As first step of the solution, we specify the parameters α, B the optimization problem. Then we fix a positive integer M and take the following Taylor polynomial to approximate the loaded shape by the turning angle field: (52)

Let the set C denote all, but the first Taylor coefficient associated with the shape, i.e. (53)

re-

M −1 C := {¯ ci } M . i=2 ≡ R

urn al P

For any fixed choice of C, the remaining coefficient, c¯1 (i.e. the slope at s¯ = 0) is determined by using the third condition in (46). There is always one unique solution (modulo 2π) as variations of c¯1 simply establishes a rigid body rotation of the the loaded configuration without changing its shape, hence angle(¯ x(1), z¯(1)) is a linear function of c¯1 . Therefore, the maximal reachout is sought via optimization in the the (M − 1) dimensional real space C.

Jo

3.2. The objective function ¯ s) has been specified, one can determine κ Once φ(¯ ¯ 2 by simple derivation according to (42) and the functions x¯ and z¯ by integration of the relevant differential equations and boundary conditions. Then, we have all ingredients for computing the discrete second variation in eq. (45) as a function of γ. We employ a uniform spatial discretization of the unit interval [0, 1] with N = 100 vertices. Derivatives of the test functions are approximated by finite difference approximations as in [24]. A positive definite discrete second variation of the potential energy implies a stable configuration, which is tested by a try of Cholesky decomposition. In practice, the eigenvalue with the smallest magnitude is used to obtain γcrit , which can be readily computed with a simple root finder algorithm (such as the function fzero in MATLAB). At this point, we can use (36), and (44) determine the maximal reachout corresponding to fixed values of C as R(C) = (γcrit B2 /F )1/2 |(¯ x(1) − x¯(0), z¯(1) − z¯(0))|. 15

(54)

Journal Pre-proof

ˆ R(C) = (F/B2 )1/2 R(C)

of

which will serve as objective function of the optimization problem. The same optimization problem can also be solved by using the rescaled reachout (55)

which is independent of the dimensional parameters F , B2 . Accordingly, it is convenient to represent rod shapes in rescaled coordinates. (56)

zˆ(s) := z(s)(F/B2 )1/2

(57)

pro

xˆ(s) := x(s)(F/B2 )1/2

to eliminate the need of specifying any dimensional parameters.

re-

3.3. Numerical optimization As the critical value γcrit always exists and is unique, the outlined algorithm assigning a maximal reachout R(C) to fixed values of C is robust. This algorithm can be combined with any numerical optimization algorithm to obtain ˆ R∗ := max R(C)

urn al P

C

(58)

i.e. the maximum of the (scaled) distance between the endpoints for given ¯1 , B ¯3 , and α. Note that R(C) ˆ values of the parameters B is a nonlinear function, hence several local maxima may exist. Once R∗ and the corresponding loaded shape have been found, one can also recover the corresponding κ ¯ 02 initial curvature function by simple integration of (43). Now, the numerical fit of a Taylor polynomial κ ¯ 02 (s)

≈

M X

¯bi s¯i−1

(59)

i=1

Jo

is straightforward. The coefficients ¯bi will be used to illustrate the beneficial effect of choosing nonzero intrinsic curvature. So far we have formulated an unconstrained optimization problem to obtain the maximal outreach. In our work we are also interested in the behavior of initially straight rods. In this case we need to meet κ ¯ 02 (¯ s) ≡ 0, which is implemented as a different unconstrained problem with a leastP optimization ¯b2 instead of R(C). ˆ square type objective function B(C) = M The numeri=1 i ical optimization procedure aims to minimize B(C), and the corresponding reachout is denoted by R0 . 16

Journal Pre-proof

4. Results

c¯2 ∈ [−4, 4],

c¯3 ∈ [−2, 2],

pro

of

4.1. The beneficial effect of the intrinsic curvature First we examine the horizontal outreach of the cantilever, i.e. we specify α = 0 (cf. (46)). Here we consider rods with possible spatial buckling modes (such as lateral torsional buckling), therefore one would expect smaller values for the maximal outreach, than predicted by eq. (3). On the other hand, the presence of intrinsic curvature may potentially improve the performance of the arm. For this investigation we take M = 4 with the following limits in the Taylor coefficients: c¯4 ∈ [−1/2, 1/2].

(60)

urn al P

re-

Recall that c¯1 is determined to meet α = 0. In the case of a rod with circular cross section, dependence of the reachout R∗ on the non-dimensional shape parameters c¯2 and c¯3 is depicted in panel a) of Figure 5. The very same results are also plotted in the (¯b1 ,¯b2 ) space spanned by the Taylor-coefficients of the non-dimensional intrinsic curvature κ ¯ 02 (s) (c.f. eq. (59)). Observe, that the maximum in the reachout belongs to an initially curved configuration (plotted in panel c) of the figure). In other words, R∗ exceeds the R0 value corresponding to initially straight cantilevers (bi = 0 ∀i, plotted in panel d)) by approximately 40%. These findings support the intuitive picture outlined in the introduction: lateral torsional buckling may reduce the critical lengths of robotic arms, whereas initial curvature might indeed significantly improve the performance of a robotic arm even if stability against lateral torsional buckling is required.

Jo

4.2. The optimal curvatures for horizontal outreach In the previous subsection we focused solely on a rod with a circular cross section and limited the dimension of the space C for the sake of illustration. Now we systematically investigate the optimal intrinsic curvatures yielding R∗ for various values of the rigidities. We keep the fixed value of α = 0. To that end, the numerical optimization process has been repeated many ¯1 and B ¯3 were systematically varied within times at fixed M = 8, while B the [1/4, 7] and [1/6, 14/3] intervals, respectively. The maximal outreach obtained this way is presented as a pair of contour plots in Fig. 6. Interestingly, ¯1 and we have identified two local optimum values of C for each value of B ¯3 . The corresponding maxima have been depicted in panels a) and b) of B 17

urn al P

re-

pro

of

Journal Pre-proof

Jo

Figure 5: Dependence of the maximal horizontal outreach of cantilevers with intrinsic curvature on the shape parameters defining the loaded shape (panel a) and the intrinsic, ¯1 = 1, B ¯3 = 2/3. unloaded shape (panel b). The cross section of the rod is circular with B The physical shape plotted in panel c) belongs to the bar with the maximal horizontal outreach, denoted by ’P’ in panels a) and b). The shape in panel d) belongs to the solution in data set closest to an intrinsically straight rod, denoted by ’Q’. Observe the significant increase in the reachhout due to the presence of intrinsic curvature.

18

Journal Pre-proof

pro

of

the figure, respectively. The first set of optima is characterized by high curvatures of the loaded shape near the tip, whereas the second one involves nearly circular loaded shapes. Clearly the first one is superior to the second one throughout the whole parameter regime considered. Systematic investigations of all local optima, and potential improvements in the case of higher M values are beyond the scope of our work. ¯1 and It is worth noting that the outreach is an increasing function of B ¯ B3 . This is not surprising, since the critical length is determined by lateral torsional buckling. In the classical case of an initially straight bar under small deflections, eq. (1) also predicts similar dependence on B1 and B3 . Contours of constant outreach hints towards a close to inverse proportion in ¯1 , B ¯3 ) space. In specific, the (B

re-

¯3 = ΩB ¯ −ι , B 1

(61)

urn al P

where ι ≈ 0.9. In panel a) of Fig. 6 some of these fitted curves (with Ω = 0.6, Ω = 1.6 and Ω = 5.0, respectively) are plotted with dashed lines for comparison. The numerical values of Fig. 6 can be compared directly with the value µ ≈ 1.81 predicted by (4). For relatively small values of B¯1 and B¯3 , we obtain similar or even slightly lower values due to sensitivity of the system to lateral torsional buckling. At the same time huge improvements (exceeding 300%) are possible for larger values of the rigidity parameters. In practice, this regime of parameters corresponds to wide and flat cross-sections, where adding intrinsic curvature dramatically increases to performance of an arm.

Jo

4.3. Maximal outreach in other directions The optimization procedure has been repeated once more. This time the angle α was varied systematically. We also considered several fixed values of stiffness. Figure 7 shows envelope diagrams of the positions of the free endpoint if curvature is optimized (panel a) and in the lack of intrinsic curvature (panel b). Equivalently the two panels can also be interpreted as a polar plot of R∗ (panel a) and the reachout R0 corresponding to no intrinsic curvatures (panel b) as functions of α. Both panels contain raw datapoints along with smooth curves fitted to these points. The only difference between panel b) and Fig. 3 is that now we consider the risk of lateral torsional buckling and thus recover a smaller set of stable solutions. In the view of Fig. 3, this computation recovers the stability 19

urn al P

re-

pro

of

Journal Pre-proof

Jo

Figure 6: Maximal horizontal outreach of cantilevers with intrinsic curvature as the rigidity parameters are varied. The panels on the left depict contours of the distance (F/B2 )1/2 z ¯1 and B ¯3 . Panels a) reachable by a cantilever characterized by non-dimensional stiffness B and b) show numerically identified, separate local optima, with panel a) representing the global optimum. The dashed curves are fitted graphs of the form eq. (61) to the contours. ¯1 , B ¯3 ) = (7.5, 5.0) (denoted The insets c) and d) on the right depict the optimal shape at (B to ’P’ in panel a)) and the optimal shape for a circular cross section (denoted to ’Q’ in ¯1 , B ¯3 ) = (7.5, 4.3) (denoted to panel a)), respectively. Similarly, the optimal shape at (B ’R’ in panel b) and the optimal shape for a circular cross section (denoted to ’S’ in panel b)) are plotted in insets e) and f), respectively. Observe the close to straight loaded shapes as optimal solutions. Also note that very high intrinsic curvature appears as optimal. Accordingly, spirals are formed by the unloaded (reference) configuration at large values ¯1 and B ¯3 . of B 20

Journal Pre-proof

urn al P

re-

pro

of

boundary in the scaled (ˆ z , xˆ) space. Nevertheless, the stability boundaries of ¯1 → ∞ Fig. 6.b approach the envelope of stable points in Fig. 3 in the B limit. Comparison of the two panels reveals that optimized curvature offers significant improvement in maximal reachout in almost all direction. For example the maximal horizontal reachout for arbitrary x(L) takes values well above 2 in panel a), i.e. it exceeds the value given by (3), whereas the values of panel b) are strictly below the limit given by (3). As a notable exception, the maximal reachout vertically upwards does not exceed (2). Here the maximal reachout is determined by Euler buckling. In fact, the optimal shape appears to be curvature-free in accordance with the common practice of using straight columns for holding loads above the ground. This result is plausible based on the analysis of the second variation of the potential energy. We have seen, that positivity of δ 2 L¯p cannot be modified by introducing intrinsic curvature in the plane. Improvement in (a nearly) vertical direction might perhaps be possible for spatial intrinsic curvature, i.e. non-vanishing κ ¯ 01 and/or κ ¯ 03 . This case is not covered by our investigation focusing on rods with planar unloaded shapes.

Jo

4.4. Relevant regimes of model parameters In order to illustrate the relevance of the findings summarized in Figures ¯1 , and B ¯3 for a simple one param5, 6 and 7, we calculate the values of B eter family of cross-sections. Consider a soft robotic arm with a steady elliptical tube section along its center-line. Let the principal axes of the ellipse (measured in meters) be denoted by a and b in directions x and y, respectively. Thepcircumference C of the ellipse is approximated [25] as C ≈ π(3(a + b) − 10ab + 3(a2 + b2 )). Let the thickness t of the section be sufficiently small with respect to a such that the interior and exterior curves of the cross section might be approximated as ellipses with a marginal error. Without loss of generality, fix the modulus of elasticity E = 1 [N/m2 ] and the Poisson ratio ν = 0.25. Let J stand for the torsional constant of the cross

21

urn al P

re-

pro

of

Journal Pre-proof

Jo

¯3 = Figure 7: Stability boundaries (dots) with numerically fitted smooth curves at fixed B ¯ 2/3 as B1 is varied with (a) and without (b) intrinsic curvature. Note the significantly larger stable region (i.e. the region below the curves) in the case of intrinsic curvature. In panel (b) the B1 → ∞ limit recovers the envelope of stable solutions presented in Figure 3. In panel (b) the predicted limits of eqs. (2) and (3) are plotted by dotted lines and the maximal horizontal outreach in eq. (4) is represented by the red dot.

22

Journal Pre-proof

pro

of

section. Now, the rigidities in the Cosserat rod model can be approximated:  3    3   π t t π t t B1 ≈ a+ b+ − a− b− , (62) 4 2 2 4 2 2   3   3 π t t π t t B2 ≈ a+ b+ − a− b− , (63) 4 2 2 4 2 2 E 8 a2 b 2 t 8 a2 b 2 t p B3 = J ≈ π2 ≈ π . (64) 2(1 + ν) 5 C 5 3(a + b) − 10ab + 3(a2 + b2 ) Let us fix a = 1 and B2 = 1, hence b is found to follow 4πt + 3πt3 − 16 πt(12 + t2 )

re-

b=−

(65)

¯1 and B ¯3 as functions of with limt→0+ b = ∞. Now, as B2 = 1 we obtain B ¯1 (t) and B ¯3 (t) are the parameter t. It is straightforward to show, that both B monotonous decreasing in t with the following limits ¯1 (t) = ∞, lim B

urn al P

t→0+

¯3 (t) ≈ 1.6825. lim B

t→0+

(66) (67)

This example demonstrates, that the effect of intrinsic curvature is indeed significant in practical situations. For instance, a ’flat’ elliptical cross section with t = 0.15, a = 1.0 and b ≈ 2.50 has B1 ≈ 4.00, B2 = 1.0 and B3 = 1.27. ∗ According to Fig. 6 R∗ ≈ 3.70, which exceeds the double of the Rdef ≈ 1.81 value predicted by (4). 5. Conclusion

Jo

By numerical investigation of the spatial instability of in-plane configurations of soft elastic cantilevers subject to a tip load, we demonstrated the previously unexplored importance of out-of-plane modes of instability. At the same time, we demonstrated the capability of controlled intrinsic curvatures to stabilize the rod in configurations with large reachout. We found that extremely large improvements occur in the case of flattish cross-sections, when the rigidities associated with out-of-plane bending and torsion are significantly larger than the in-plane stiffness. At the same time, the performance of arms with high and narrow cross-sections was improved slightly. 23

Journal Pre-proof

urn al P

re-

pro

of

Our findings offer a new paradigm to design efficient expandable robotic arms. A flat and wide arm can be stored efficiently in a rolled configuration, and it offers high resistance against lateral torsional buckling. Nevertheless, it suffers from large in-plane deformations, which seriously limit its performance. Our new results imply that compensation of deformations by controlled intrinsic curvatures, may result in a dramatically improved performance. Note that rods with significant improvement in their stable length are highly curved when unloaded. Direct application of the tip load on such a shape would not results in the predicted, outreached shape. In terms of soft robotic arms this difficulty can be overcome by either enforcing a straight shape on the rod until the load is applied or applying curvature control in such a way, that intrinsic curvature is increased as the load is applied. Our optimization technique might be readily applied for such a control problem. Our study assumed a simple mechanical model neglecting the effect of warping, compression and shear deformations as well as material nonlinearities or non-elastic effects. While all these factors could be implemented in our modeling approach, we believe that they do not play crucial roles in the studied phenomenon, hence we chose to use a simple model with low number of material-related parameters. We also used the a-priori assumption of planar unloaded shapes. Extension to spatial unloaded shapes is also a possible further direction, which might potentially uncover some non-intuitive optimal solutions. At the same time, we assumed that intrinsic curvatures are arbitrary while the intrinsic curvature profiles of existing pneunet and shape memory alloy actuators are limited [26, 27]. Such a situation can be addressed by adding constraints to the optimization problem. We focused on cantilevers carrying a vertical tip load, which is highly relevant as an idealized model of a light arm carrying a heavy manipulator or sensor device. Our study can be extended in the future for other types of loads including for example the weight of the arm. Another possible extension addresses the optimal performance of non-prismatic arms, which is also beyond the scope of the present work.

Jo

Acknowledgment

Research was supported by Grant 124002 of the National Research, Development, and Innovation Office, Hungary.

24

Journal Pre-proof

Author contribution

of

A.A.S and P.L.V developed the theoretical part, developed the numerical code and wrote the paper. References

pro

[1] A. Gajewski, M. Zyczkowski, Optimal structural design under stability constraints, volume 13, Springer Science & Business Media, 2012. [2] N. V. Banichuk, Problems and methods of optimal structural design, volume 26, Springer Science & Business Media, 2013.

re-

[3] S. J. Cox, C. McCarthy, The shape of the tallest column, SIAM Journal on Mathematical Analysis 29 (1998) 547–554. [4] C. M. McCarthy, The tallest column—optimality revisited, Journal of computational and applied mathematics 101 (1999) 27–37.

urn al P

[5] J. D´ıaz, M. Sauvageot, Euler’s tallest column revisited, Nonlinear Analysis: Real World Applications 11 (2010) 2731–2747. [6] Y. V. Egorov, On the tallest column, Comptes Rendus M´ecanique 338 (2010) 266–270. [7] D. Naicu, C. Williams, The use of dynamic relaxation to solve the differential equation describing the shape of the tallest possible building, in: VII International Conference on Textile Composites and Inflatable Structures, University of Bath, 2015, pp. 34–45. [8] I. U. Cagdas, S. Adali, Optimal shapes of clamped–simply supported columns under distributed axial load and stress constraint, Engineering Optimization 45 (2013) 123–139.

Jo

[9] R. H. Plaut, L. N. Virgin, Optimal design of cantilevered elastica for minimum tip deflection under self-weight, Structural and Multidisciplinary Optimization 43 (2011) 657–664. [10] A. Bratus, The optimum shape of a bending beam, Journal of Applied Mathematics and Mechanics 64 (2000) 993–1004.

25

Journal Pre-proof

of

[11] R. Drazumeric, F. Kosel, Shape optimization of beam due to lateral buckling problem, International Journal of Non-Linear Mechanics 47 (2012) 65–74. [12] S. P. Timoshenko, J. M. Gere, Theory of elastic stability, 2nd ed., McGraw-Hill International Book Company, NY, 1963.

pro

[13] C. Wang, Longest reach of a cantilever with a tip load, European Journal of Physics 37 (2015) 012001. [14] M. Batista, Comment on ’longest reach of a cantilever with a tip load’, European Journal of Physics 37 (2016) 058004.

re-

[15] R. H. Plaut, L. N. Virgin, Furthest reach of a uniform cantilevered elastica, Mechanics Research Communications 83 (2017) 18–21. [16] C. Armanini, F. Dal Corso, D. Misseroni, D. Bigoni, From the elastica compass to the elastica catapult: an essay on the mechanics of soft robot arm, Proc. R. Soc. A 473 (2017) 20160870.

urn al P

[17] R. J. Webster III, B. A. Jones, Design and kinematic modeling of constant curvature continuum robots: A review, The International Journal of Robotics Research 29 (2010) 1661–1683. [18] C. Majidi, Soft robotics: a perspective—current trends and prospects for the future, Soft Robotics 1 (2014) 5–11. [19] Y. Yekutieli, R. Sagiv-Zohar, R. Aharonov, Y. Engel, B. Hochner, T. Flash, Dynamic model of the octopus arm. i. biomechanics of the octopus reaching movement, Journal of neurophysiology 94 (2005) 1443– 1458. [20] C. Laschi, M. Cianchetti, B. Mazzolai, L. Margheri, M. Follador, P. Dario, Soft robot arm inspired by the octopus, Advanced Robotics 26 (2012) 709–727.

Jo

[21] J. Burgner-Kahrs, D. C. Rucker, H. Choset, Continuum robots for medical applications: A survey, IEEE Transactions on Robotics 31 (2015) 1261–1280. [22] S. S. Antman, Nonlinear Problems in Elasticity, 2nd ed., SpringerVerlag, NY, 2005. 26

Journal Pre-proof

of

[23] T. J. Healey, P. G. Mehta, Straightforward computation of spatial equilibria of geometrically exact cosserat rods, Int. J. Bifur. Chaos 15 (2005) 949–965.

pro

[24] T. J. Healey, A. A. Sipos, Computational stability of phase-tip splitting in the presence of small interfacial energy in a simple two-phase solid, Physica D: Nonlinear Phenomena 261 (2013) 62 – 69. [25] S. Ramanujan, Modular equations and approximations to π, Quart. J. Pure App. Math. 45 (1914) 350–372.

re-

[26] K. M. de Payrebrune, O. M. O’Reilly, On constitutive relations for a rod-based model of a pneu-net bending actuator, Extreme Mechanics Letters 8 (2016) 38–46.

Jo

urn al P

[27] N. N. Goldberg, X. Huang, C. Majidi, A. Novelia, O. M. O’Reilly, D. A. Paley, W. L. Scott, On planar discrete elastic rod models for the locomotion of soft robots, Soft robotics (2019).

27

Journal Pre-proof

Highlights

pro

reurn al P

 

The interplay between spatial buckling and intrinsic curvature and the performance of soft robotic arms is investigated. Optimization of the intrinsic curvature significantly improves the maximal outreach of the arm. Maximal outreach for a tip load at the free end of the arm is numerically determined as the stiffness of the rod is varied.

Jo



of

for the paper: András A. Sipos, Péter L. Várkonyi: The longest soft robotic arm, submitted to the International Journal of Non-linear Mechanics

Journal Pre-proof

Declaration of conflicts of interests for the paper: András A. Sipos, Péter L. Várkonyi:

The longest soft robotic arm, submitted to

of

the International Journal of Non-linear Mechanics

The authors of the paper titled The longest soft robotic arm declare no competing financial interests.

pro

Budapest, November 8, 2019

Jo

urn al P

re-

Andras A. Sipos Peter L. Varkonyi