Position analysis in planar parallel continuum mechanisms

Mechanism and Machine Theory 132 (2019) 13–29

Contents lists available at ScienceDirect

Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmachtheory

Research paper

Position analysis in planar parallel continuum mechanisms Oscar Altuzarra∗, Diego Caballero, Francisco J. Campa, Charles Pinto Department of Mechanical Engineering, University of the Basque Country UPV/EHU, Bilbao, Spain

a r t i c l e

i n f o

Article history: Received 23 May 2018 Revised 27 September 2018 Accepted 16 October 2018

Keywords: Parallel continuum mechanism Kirchhoff rod Inverse Forward

a b s t r a c t Parallel continuum mechanisms are closed, kinematic, chained compliant mechanisms that are mobile due to flexible and slender elements. Large deformations and the force-moment equilibrium require a more complex modeling of the kinematic analysis than in their rigidlink counterparts. Although real-time solutions to the position problem using numerical iterative methods are available, a comprehensive kinematic analysis of such mechanisms has not yet been conducted. The first step to this end is a procedure that solves the full forward and inverse kinematic problems, and is proposed in this paper for planar systems to allow for an easier interpretation of their kinematic characteristics. This procedure provides the means to obtain multiple solutions arising from the closed-loop structure and the multiplicity of the buckling mode of slender links. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Compliant mechanisms are ones that are mobile through the relative flexibility of some of their parts. Their advantages include reduced wear, no backlash, and high precision [1], which renders them appealing for such applications as space exploration [2], micromanipulation [3], and surgical tools [4]. Mechanisms where motion occurs owing to slender elements that exhibit large nonlinear deformations, rather than due to notches in elements with lumped flexibility, are called continuum mechanisms [5]. Continuum mechanisms are largely motivated by the problem of manipulation in confined, hard-to-reach work spaces [6], and many are inspired by biological systems such as trunks, tentacles, and snakes. This makes them suitable for minimally invasive procedures due to their dexterity and ease of miniaturization [7]. Within continuum mechanisms, the analysis of parallel architectures is new [5]. The so-called parallel continuum mechanisms can provide higher precision and stiffness than slender, single-member continuum mechanisms [8]. Their kinematic analysis has some similarities with their rigid-link counterparts and cable-driven manipulators alike. This analysis requires a model of the rod’s nonlinear deformation, posing a problem that is not uniquely of a geometrical kind. This means that the pose of this kind of mechanisms is influenced by both actuation and external loading [9]. Moreover, the closed-loop assembly of rods introduces a coupling of equations not only through position and orientation constraints, but also due to the equilibrium between the deformation and the load. This has implications when analyzing problems in forward and inverse kinematics apart from classical singularities, as changes to buckling modes or instabilities. Hence, a comprehensive kinematic analysis of the full position problem will provide a tool for the detailed analysis of the work space and improve the performance of the mechanism considered.

∗

Corresponding authors. E-mail addresses: [email protected] (O. Altuzarra), [email protected] (D. Caballero).

https://doi.org/10.1016/j.mechmachtheory.2018.10.014 0094-114X/© 2018 Elsevier Ltd. All rights reserved.

14

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

An appropriate model for the nonlinear deformation of a rod should be accurate, computationally efficient for real-time implementation, and allow for the control of the buckling mode. The importance of this last condition lies in the fact that the position problem of a slender rod is not unique due to its ability to deform according to different buckling modes. In the literature there are various methods that have been proposed to solve large deformation in slender elements, like circle-arc method [10], pseudo-rigid-body model (PRBM) method [11], finite element models [12], among others. The Cosserat rod model is a commonly used one. It is accurate as it maintains the exact nonlinear deformation, and allows a computationally effective algorithm to solve real-time position problems. It has been used for modeling continuum robots driven by tendons, pre-shaped tubes, and pneumatic actuators [13], or in computer graphics to simulate realistic movements of the hair [14]. It requires the solution to a known pose to solve iteratively close positions, and thus does not allow for the solving of all possible solutions. The Cosserat rod model is a conceptual model that can thus be mathematically expressed in different ways, with a nonlinear system of differential equations the most accurate one. Another characteristic of the Cosserat rod model is that its dynamic expression can be obtained and implemented within a real-time algorithm [8,15]. This is beyond the scope of this paper but highlights its ability as a general approach. When dealing with planar mechanisms—and slender elements with no shear and extension are straight rods in their freestress reference configuration—the Cosserat nonlinear differential system of equations leads to a simplified model called the Kirchhoff model. This can be manipulated conveniently to arrive at a mathematic model in which the rod deformation is expressed through elliptic integrals. The elliptic integral solution method allows us to obtain the deformation of a rod from elliptic integral parameters in a univocal manner, and these are defined within a limited range. The use of this approach is now classical in the literature on Compliant mechanisms. There are various methods of elliptic integral solution, of which the one developed in [1] stands out for its simplicity. It has been used extensively for the analysis of compliant mechanisms where the flexible rod was easily assembled with other rigid elements so that its nonlinear analysis was uncoupled. For example, in [16] this approach is used for the analysis of a microactuator based in a bistable mechanism and the result is compared with finite element analysis. Other similar models require a division of the rod into as many parts as there are inflection points [17]. Our research in this paper follows from earlier work on the finding of solutions for some simple parallel continuum mechanisms appeared at some conference papers. In [18,19] we did the quasistatic analysis of planar deformable systems where, as in [16,17], the solution of the deformation of each rod was uncoupled and boundary conditions were constant and known, hence there was not need for a comprehensive analysis of multiple solutions. However, we found that such multiplicity could appear in similar parallel mechanisms. In [20,21] we did some preliminary work on this direction, first making extensive use of numerical direct integration from a diversity of initial guesses looking for possible solutions, and then making use of particular conditions on a certain mechanism to find systematically all solutions. The natural subsequent step was to establish a comprehensive and systematic method for non-redundant planar parallel continuum mechanisms, that it is the present contribution. To the best of our knowledge, no research to date has been conducted on the full kinematic analysis of a planar parallel continuum mechanism, mainly due to the difficulty of obtaining the number of possible solutions of the position problem. Two features render the elliptic integral solution method suitable to solve the multiplicity problem of this kind of mechanisms. First, it allows control over the choice of the buckling mode of each rod involved. Second, the modular angle and the elliptic modulus, i.e., the two parameters of elliptic integrals, uniquely define rod deformation. Our goal here is to show a unified method to solve both the forward and the inverse kinematic position problems. The remainder of the paper is structured as follows: In Section 2, we start with a review of the elliptic integral method where we recall the fundamentals of the nonlinear analysis of flexible rods using the Cosserat rod model [22], and simplify it to the Kirchhoff rod model so that planar systems are solved using elliptic integrals [1]. In Section 3, we show a method to assemble these sets of equations for closed-loop planar mechanisms using boundary conditions. In Section 4, we describe an iterative procedure to search all possible solutions of the Forward Kinematic (FK) position problem for planar closed-loop flexible mechanisms of 2 and 3 degrees of freedom (DOF). In Section 5, the aforementioned procedure is adapted to search all possible solutions of the Inverse Kinematic (IK) problem. In Section 6, a brief discussion is provided of the singularities of both FK and IK position problems, the stiffness that characterizes each position, and work space analysis. Finally, in Section 7, the conclusions of this paper are drawn.

2. Fundamentals of nonlinear analysis of deformation of rods The purpose of this section is to recall the fundamentals on the solution of the nonlinear analysis of the deformation of elastic planar rods with the elliptic integrals method, a detailed treatment of the contents in this section can be found in [1,22]. The shape of a deformed rod is defined by a pair of mathematical structures: a parametric Cartesian curve in space p(s ) ∈ R3 that depicts the centroid of each cross-section, and the orthonormal rotation matrix R(s) ∈ SO(3) expressing the orientation of a local frame attached to the section. The x- and y-axes are usually placed on the principal cross-section axis, while the z-axis is on the tangent to the deformed rod, see Fig. 1. Both the position and orientation are functions of a scalar reference arc-length parameter s over some finite interval, s ∈ [0, L], where L is the length of the rod in the reference configuration. As it is assumed that the sections remain planar, a mapping from s to a homogeneous rigid-body transformation, T(s) ∈ SE(3), describes the entire deformation of the rod.

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

15

Fig. 1. A rod deformed in space with frames used.

Now, let us consider the derivatives of p(s) and R(s) with respect to s. The position and orientation evolve along the arc length according to linear, v(s ) ∈ R3 , and angular, u(s ) ∈ R3 , rates of change expressed in the local frame:

dp(s )  = R ( s )T p ( s ) ds  0 −uz (s ) T  0 U ( s ) = R ( s ) R ( s ) = uz ( s ) −uy (s ) ux ( s ) v ( s ) = R ( s )T



uy ( s ) −ux (s ) 0

(1)



where the superscript denotes a derivative with respect to s. Changes in these magnitudes with respect to the values of the reference configuration of the rod, namely v(s ) = v(s ) − v◦ (s ) and u(s ) = u(s ) − u◦ (s ), are related to material strain. On the one hand, vz represents the local stretching of the rod; with vz > 1 for extension, vz < 1 for compression, and vz = 1 when the length of the curve remains constant; vx and vy depict the shear deformation along the local axes. On the other hand, ux and uy are related to flexure along the principal cross-section axes, and uz to torsion. The constitutive relationship between strain and stress can be defined depending on the mechanical behavior of the material. In this sense, effects like hyperelasticity, plasticity, or viscoelasticity can be considered. Nevertheless, the classical Cosserat rod model explained here is valid as long as the assumption that there is no cross-sectional deformation is valid. Because the stress field is obtained due to the strain field, the former can be expressed as a function of v(s) and u(s). And because the internal force n(s) and moment m(s) of the entire cross-section are due to the stress field, they can also be expressed as a function of v(s) and u(s). In the following we use a linear constitutive law that relates linear deformation v to internal forces, and angular deformation u to internal moment, expressed in the local frame through the use of some diagonal stiffness matrices, i.e., KSE (s) for shear and extension, and KBT (s) for bending and torsion. The internal forces and moments expressed in the fixed global frame are obtained through the rotation matrix as

n(s ) = R(s )KSE (s )v(s ) m(s ) = R(s )KBT (s )u(s )

(2)

where the stiffness matrices are



KSE =

0 GA(s ) 0

0 0 EA(s )

EIx (s ) 0 0

0

0 0

 KBT =



GA(s ) 0 0

EIy (s ) 0



GJz (s )

(3)

16

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

with E and G as the elastic and shear modules, respectively, A as the area of the cross-section, Ix and Iy as the second moments of area about the local x- and y-axes, respectively, and Jz as the polar moment about the local axis z. To establish the equilibrium of forces and moments for every portion of the deformed rod, we consider an infinitesimal piece of the rod between sections located at s and s + ds. By convention, we denote the internal force by −n(s ) and the moment by −m(s ) that material [0, s) exerts on portion [s, s + ds], and by n(s + ds ), and m(s + ds ) the force exerted by material (s + ds, L]. The portion considered may also be subjected to distributed forces f(s) and moments l(s) defined within the infinitesimal interval ds. The equilibrium of forces and moments yields:

−n(s ) + n(s + ds ) + f(s )ds = 0 −m(s ) + m(s + ds ) + l(s )ds + p(s ) × (−n(s ))+ +p(s + ds ) × n(s + ds ) + p(s + μds ) × f(s )ds = 0

(4)

where 0 < μ < 1 denotes the position where the resultant of the distributed load f(s) is applied. Upon simplification, the classical form of the equilibrium nonlinear differential equations for a Cosserat rod describing the evolution of the internal force n(s) and moment m(s) along arc length arise:

dn(s ) + f (s ) = 0 ds

dm(s ) dp(s ) + × n (s ) + l (s ) = 0 ds ds

(5)

Once the kinematic, constitutive, and equilibrium laws have been stated, they must be combined to obtain a set of equations that determine the shape of the deformed rod. In case of a simplified model of the Kirchhoff rod, shear and extension are neglected because the deformed shape is considered to be caused by only bending and torsion (which is more likely to happen as slenderness increases). Thus, v = v◦ . Considering that there are usually no shear and extension deformations in the reference deformation state, i.e., v = v◦ = e3 where e3 is a unit vector tangent to the deformed shape, the system of differential equations becomes 

p = R e3 

R = RU 

n = −f 



m = −p × n − l

(6)

Such a system of differential equations is suitable for the general case of a rod that deforms in space. When the movement of all points on the rod is within a plane, further simplifications apply. Consider the rod in Fig. 2, the deformation of which occurs within the plane OXY. The orientation of an arbitrary cross-section can be expressed through the rotation matrix as a function of angle θ = θ (s ):



cos θ sin θ 0

R (s ) =

− sin θ cos θ 0

0 0 1

 (7)

the derivative of which with respect to arc length s is

 

R (s ) =

− sin θ cos θ 0



− cos θ − sin θ 0

0 0 0

θ



(8)

and allows finding vector u in:

 

U = RT R =

0

θ



0

−θ 0 0



0 0 0



 =

0 uz 0

−uz 0 0



0 0 0

(9)

Thus, only the component uz of vector u is not null. This makes sense, as there is no bending moment with respect to the local y-axis and no torsion (expressed through ux ).

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

17

Fig. 2. Deformed shape of a clamped-hinged rod within a plane under load R at extreme.

For the Kirchhoff rod model in the plane, the expression of the system of differential equations Eq. (6), considering Eq. (2), is:

⎧ dx ⎪ ds ⎪ dy ⎪ ⎪ ⎪ ⎨ ddsθ

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

⎧ ⎪ ⎪ ⎪ ⎪ ⎨

⎫

cos θ ⎪ ⎪ ⎪ sin θ ⎪ ⎬ mz + u ◦,z ds EI = dnx − fx ⎪ ⎪ ds ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dny ⎪ − fy ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ ds ⎭ ⎩ dm nx sin θ − ny cos θ − lz z

(10)

ds

For the stress-free reference state of a straight rod, i.e., u◦,z = 0, the result is consistent with the well-known Bernoulli– Euler law, which states that the bending moment mz at a point is proportional to the curvature κ :

mz dθ 1 = = =κ EI ds ρ

(11)

where ρ is the radius of the curvature. If no distributed force and moment, f and l, are applied along the rod, we get

⎧ ⎫ ⎨ dndsx ⎬ ⎩

dny ds dmz ds

⎭

 =

nx sin θ

0 0 − ny cos θ

(12)

Hence, internal force n is constant. If only end-point load is applied (expressed in terms of magnitude R and orientation

ψ as in Fig. 2), for constant E and I, with a stress-free reference straight rod, considering Eq. (11), and upon substitution into Eq. (12), we get:

d2 θ 1 R = (R cos ψ sin θ − R sin ψ cos θ ) = sin (θ − ψ ) EI EI d2 s

(13)

Therefore, the system of differential equations Eq. (10) has been simplified to:

⎧ ⎨ ⎩

dx ds dy ds d2 θ d2 s

⎫ ⎬ ⎭

 =



R EI

cos θ sin θ sin (θ − ψ )

(14)

18

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

A first integration in θ of the last equation in Eq. (14) provides:

dθ = ds

2C −

2R cos (θ − ψ ) EI

(15)

where C is a constant of integration. A second integration requires a more complex mathematical manipulation, where the first step is a sort of change of parameters from θ , C to φ , k:

C=

R ( 2k2 − 1 ) EI

 cos

(16)



ψ −θ

= k sin φ

2

(17)

Note that k remains a constant, whereas φ varies continuously along the rod as in the case of θ . Upon differentiation of Eq. (17) we get:

dθ =



2k cos φ 2

1 − k2 sin

dφ

φ

(18)

and simple trigonometric manipulation of Eq. (17) provides:

cos (θ − ψ ) = 2k2 sin

2

φ−1

(19)

Upon substitution of Eqs. (16), (18), and (19) into Eq. (15), and some simplifications, we get:

 0

L

 φ2 R 1 ds = dφ  2 EI φ1 1 − k2 sin φ √ √ EI R= [F (k, φ2 ) − F (k, φ1 )] L

(20)

where F(k, φ ) is the incomplete elliptic integral of the first kind. Angle φ varies continuously along the length of the rod between φ 1 and φ 2 , which are defined by the boundary conditions of the rod at extremes. For example, if we consider a rod clamped horizontally at the first extreme, see Fig. 2, θ (0 ) = 0, and Eq. (17) provides the corresponding value of φ , namely



φ1 = arcsin

1 cos k



ψ



(21)

2

whereas at the other extreme, upon substitution of Eqs. (16) and (17) into Eq. (15) we get:

dθ = 2k ds

R cos φ EI

(22)

For a hinged extreme, this is a point of inflection of the curve, ddsθ is null, and φ at this position is

φ2 = ( 2 q − 1 ) π / 2

for q = 1, 2, . . . ,

(23)

where q identifies the so-called buckling modes. Each buckling mode has a corresponding number of points of inflection, i.e., mode 1 has one inflection point precisely at the hinged end, mode 2 has two, and so on. To integrate Eq. (14) to get the x coordinate of a point on the curve, we can manipulate its expression as follows:

dx = cos θ ds

(24)

where we substitute ds from Eq. (20), and cos θ is a function of k and ψ . As a result of integration, we get the following:

x ( φi ) = −





 φi  φi  EI 2 cos ψ 2 1 − k2 sin φ dφ − R φ1 φ1



1 1−

k2

2

sin

φ

dφ

−

EI 2k sin ψ [cos φi − cos φ1 ] R

(25)

where the first integral is the incomplete elliptic integral of the second kind, E(k, φ ). Thus, we can obtain the coordinate x of the deformed shape of the rod for position φ i along the arc length as:

x ( φi ) = −

EI cos ψ [2E (k, φi ) − 2E (k, φ1 ) − F (k, φi ) + F (k, φ1 )] − R

To find x at the extreme of the rod we use φi = φ2 .

EI 2k sin ψ [cos φi − cos φ1 ] R

(26)

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

19

0.8

0.8

0.6

0.7

0.4

0.6

0.2

0.5

y [m]

k [-]

1

0 -0.2

R

0.4 0.3

-0.4

0.2

-0.6

Y

R

0.1 -0.8

X

0 -1 0

50

100

150

200

250

300

350

-0.2

0

0.2

0.4

0.6

x [m]

[deg]

Fig. 3. Deformed shapes the end-point force R for several buckling modes of a clamped-hinged rod for given values of k and ψ .

To get the y coordinate of a point on the curve, an analogous deduction to the one for x provides:

y ( φi ) = −

EI sin ψ [2E (k, φi ) − 2E (k, φ1 ) − F (k, φi ) + F (k, φ1 )] + R

EI 2k cos ψ [cos φi − cos φ1 ] R

(27)

Hence, Eq. (14) has been integrated to yield Eqs. (26), (27), and (20), i.e., a parametric system that uniquely defines a deformed rod in terms of parameters k and ψ for known boundary conditions and a given buckling mode. Note that feasible values of k depend on ψ , in order to make φ 1 in Eq. (21) real. This condition is indicated with the gridded area in Fig. 3 left. As shown in Fig. 3, if we set a duple of values to k and ψ , inside the feasible range (figure on the left), we will get the deformed shape of the rod for every buckling mode chosen (see figure on the right).

3. Analysis of closed-loop mechanisms of flexible rods The planar parallel continuum mechanisms considered here are those formed by closed-loop assemblies of two flexible rods at a punctual end-effector point P, or three rods on a rigid-body platform. See Fig. 4. An actuator drives each rod, and attachments to the rods include the clamped and hinged conditions. Each flexible rod i is ruled by Eqs. (26), (27) and (20), and uniquely defined by parameters ki and ψ i (the superscript identifies the rod). The quasi-static position problem of the deformed mechanism is solved by finding the set of parameters ki and ψ i for each rod subject to the corresponding boundary conditions. In this section, we define the set of boundary conditions formed in three types: first, slope-related conditions upon the attachment of rods to actuators and the end effector; second, static equilibrium at the end effector; and third, geometric attachment of rods to the end effector. To start with, we use local frames for each rod placed at attachments actuated by the motors, as in Fig. 5. In this way, equations obtained in the previous section can be used directly in these local frames and transformed into a fixed frame using the corresponding transformation matrices that are dependent on the inputs. For example, if rod i is clamped to rotational actuators, as shown in Fig. 4a, the local frame is oriented according to it as in Fig. 5. Hence, we can use a first boundary condition for rod given by Eq. (21), namely



φ = arcsin i 1

1 cos ki



ψi



2

(28)

As the rod ends at the end effector are hinged joints, Eq. (23) applies:

φ2i = (2qi − 1 )

π 2

(29)

where qi is defined by the buckling mode desired for each rod. The different combinations of buckling modes considered are one source of multiple solutions to the position problem. Then, there must be a static equilibrium of forces and moments at the end effector, including not only the end-forces Ri and end moments mzi due to the deformation of the rods, but those due to the load wrench, i.e., FX , FY , and MZ . In the global fixed frame, using the corresponding rotation matrices for each rod RFi that depend on the inputs, and taking

20

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

0.3 0.2

0 -0.1 -0.2 -0.3 -0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

x [m]

1 0.8 0.6 0.4

y [m]

y [m]

0.1

0.2 0 -0.2 -0.4 -0.6 -1

-0.5

0

0.5

x [m]

Fig. 4. Flexible mechanisms of 2 and 3 DOF: 2RFR and 3PFR.

0.6

0.7

0.8

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

21

0.6 0.5

R2

0.4

Fext R

0.3

1

y [m]

0.2 0.1

Y 1

0 -0.1

2

X

O 1 y1 O x

x2

1

O2 y

2

-0.2 -0.3

-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

x [m] Fig. 5. A 2RFR flexible mechanism assembled, indicating the local frames and the equilibrium of end-point forces Ri when no external load is applied.

moments about the loading point on the end effector, we get: n  i=1 n 



RFi

Ri cos ψ i · Ri sin ψ i



ci × RFi ·

i=1



 

F + X FY

Ri cos ψ i Ri sin ψ i



+

 

= n 

0 0

mzi k + MZ k = 0

(30)

i=1

where n is the number of flexible rods attached to the end effector. The moment exerted by the flexible rods on the end effector mzi is given by Eq. (11):

mzi = EI

 dθ i = 2ki EIRi cos φ2i ds

(31)

ci is the location of the end-point attachments on the end effector with respect to the loading point. In case of flexible mechanisms with 2 DOF that have a punctual end effector P, n = 2 and ci = 0. If the rods are hinged to P, mzi = 0 and, hence, no external moment can be applied, the second condition in Eq. (30) is null. Finally, the assembly of rods at the end effector implies that the end coordinates of the rod are related. The local frames for each rod allows us to directly use expressions for the end coordinates given by Eqs. (26) and (27). Transformations of i the homogeneous coordinates of end points p (Li ) = [xi , yi 1]T from local frames to a fixed frame require transformation F matrices Ti that depend on the inputs. If we have a punctual end effector, P as in Fig. 5, the end coordinates for both rods must be equal:

TF1 p (L1 ) = TF2 p (L2 ) 1

2

(32)

For a rigid-body triangular platform, the three end points of the rods must fulfill rigid-body conditions, i.e., that the distances dij between the attachment points and angles β ijk between the sides are known constants, i.e.:

|p1 (L1 ) − p2 (L2 )| = d12 |p1 (L ) − p3 (L3 )| = d13  1   1 1  2 | p (L1 ) − p (L2 ) × p (L1 ) − p3 (L3 ) | = d12 d13 sin β213

(33)

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

k r [-]

k [-]

22

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

-1

-1 0

60

120

180

240

300

360

0

60

120

180

240

300

360

Fig. 6. Plots of feasible values of ki against ψ i versus kirel against ψ i .

4. Solving the full FK problem The forward kinematic problem is defined with the given values of the inputs at the motors, with the objective of finding the pose of the end effector. To this end, we need to find the parameters that define the deformed shape of each flexible rod i in the mechanism, i.e., ki and ψ i , so that static equilibrium and the requirement of a geometric assembly are fulfilled. For example, in 2 DOF flexible mechanisms with a punctual end effector as shown in Fig. 5, for the given poses of local frames, we search (ψ i , ki ) sets corresponding to forces Ri that equilibrate the system and end-point coordinates that assemble. The foregoing objective is approached with an iterative method in three steps: define the grid of the sampled parameters, apply static equilibrium conditions to reduce the set of feasible parameters, and verify the geometric conditions of assembly to obtain i-duples of parameters that provide solutions. During this process, the buckling mode of each rod must be chosen, so all possible combination of buckling mode are evaluated. First, we define the full grid of parameters that are considered. We note that the conditions given in Eq. (28) imply a restriction on the possible values of ki with respect to ψ i . There is a kimin such that real values are obtained for φ1i in Eq. (28), defined by:



kimin

= cos

ψi



(34)

2









As shown in the left side of Fig. 6, the values of ki are in the restricted ranges −1, −kimin and kimin , 1 .

To sample the values of ki and ψ i such that a variation in the conditions is evenly represented, preliminary tests on simple examples showed that the sampled values were better defined with the parameter kirel that provides a rectangular range as shown in the right side of Fig. 6:

kirel = Sign(ki ) ·

|ki | − |kimin | 1 − |kimin |

(35)

Moreover, grid density could be adjusted with a higher density in areas where we observed that the results changed rapidly, as shown in Fig. 6-right. Second, we apply static equilibrium conditions at the end effector to restrict the set of feasible parameters. For each discretized value of ψ i in its range [0, 2π ], we take values of kirel in the range [−1, 1], respectively ki . From each pair of ki ,

ψ i , we obtain the corresponding value of end force Ri at rod i in Eq. (20) using the corresponding boundary conditions of

the rod, e.g., for clamped- hinged Eqs. (28) and (29). Note that such values depend only on the material of the rod, and its length and end conditions. This makes it possible to define the buckling mode sought for each rod. Then, in Eqs. (30), we can obtain the combinations of Ri and ψ i (eventually ki and ψ i ) that fulfill these conditions. For example, for 2 DOF flexible mechanisms with punctual end effectors where two rods are assembled, each value considered for ψ 1 and k1 (corresponding to a value of R1 ) defines a pair of ψ 2 and R2 . With the value of the end-point force R2 and ψ 2 , it is possible to obtain a corresponding k2 with Eq. (20). This reduces the feasible parameter space to be searched for solutions to the problem to the grid (ψ 1 , k1 ) because each of these feasible points has a corresponding (ψ 2 , k2 ). However, multiple solutions are possible for k2 , and must be found with a root-finding method at a significant computational cost. Finally, we obtain the values for xi , yi in Eqs. (26) and (27) with the foregoing sets of feasible parameters. Upon substitution into geometrical assembly conditions, Eqs. (32) or (33), we can evaluate a residue of these conditions where, when the residues of these conditions are null, we obtain the solutions.

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

23

1

k r2

0.5

0

-0.5

-1 360 300 -1

240 180

-0.5 120

0 60

1

[rad]

0.5 0

1

k 1r

Fig. 7. Plots of fulfillment of conditions for R2 − X and R2 − Y, indicating solutions.

0.6 1

0.4 0.8

y [m]

y [m]

0.2

0

0.6

0.4

-0.2 0.2

-0.4 0

-0.6

-0.4

-0.2

0

x [m]

0.2

0.4

0.6

-0.6

-0.4

-0.2

0

0.2

0.4

x [m]

Fig. 8. The three solutions of the FK for given inputs and mode 1 in both rods in 2RFR and 2PFR mechanisms.

Alternatively, to speed up the solution by avoiding the root-finding of k2 , we propose a feasible parameter space (ψ 1 , k1 , k2 ). For each point in this space, we get a unique value of ψ 2 , and three residues of the conditions on R2 , and the coordinates X and Y of assembly at P. The simultaneous fulfillment of the three conditions provides multiple solutions. It is possible to plot the null residues against the grid of sampled values of (ψ 1 , k1 , k2 ) that were computed. To get a clearer view, we define curves of the simultaneous fulfillment of conditions for (R2 , X) and those for (R2 , Y), as shown in Fig. 7. At the intersections of the curves, we obtain multiple solutions of the position problem in terms of parameters (ψ 1 , k1rel , k2rel ), respectively (k1 , k2 ), as ψ 2 uniquely determined by them. For a 2 DOF mechanism with rotational inputs at clamped ends, as in Fig. 4a, the assembled solutions obtained are shown in Fig. 8-left. The FK problem of 2 DOF mechanisms with different types of inputs is solved in the same way, with

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

0.2

0.2

0

0

0

-0.2

-0.2

-0.2

-0.4

-0.4

-0.4

-0.6

y [m]

0.2

y [m]

y [m]

24

-0.6

-0.6

-0.8

-0.8

-0.8

-1

-1

-1

-1.2

-1.2

-1.2

-1.4

-1.4

-1.4

-0.2

0

0.2

0.4

x [m]

-0.2

0

0.2

0.4

-0.2

x [m]

0

0.2

0.4

x [m]

Fig. 9. The three solutions to the IK for a given output and mode 1 in both rods in a 2PFR mechanism.

the transformation matrices for each rod being defined accordingly, as shown in Fig. 8-right, for a mechanism with linear actuators. In case of a rigid-body end effector, the process is analogous. This time, we need to consider the initial sets of four parameters, (ψ 1 , k1 ) and (ψ 2 , ψ 3 ). Applying the end-effector force balance, Eq. (30), we can obtain R2 and R3 . As before, we can find corresponding multiple values of k2 and k3 with a root-finding method on Eq. (20). For each of these potential solution sets, we evaluate the residues for the conditions of the balance of moments, and the assembly conditions of the end-effector platform. When we obtain sets that approximate to conditions, it is convenient to apply numerical integration to refine the solution. The n-tuples of parameters that fulfill all conditions are the multiple solutions. 5. Solving the full IK problem The inverse kinematics problem is defined with given values of the pose of the end effector, with the objective of solving for the input values of the actuators. The boundary conditions to be used are, as before, the conditions on the slope of each rod at the extremes, as in Eqs. (28) and (29), the conditions for static equilibrium in Eqs. (30), and the assembly conditions of the end points of each rod for the end effector. In the IK problem, as the output coordinates are known, the assembly conditions used instead of Eqs. (32) and (33) are, for the punctual end effector P, given by:



TFi

xi yi 1

= pP

for i = 1, 2

(36)

and for rigid body platforms,



TFi

xi yi 1



=

TFEE

xEE i yEE i 1



for i = 1, 2, 3

(37)



T

where TFEE is the known pose of the end-effector platform, and xEE , yEE 1 are homogeneous coordinates in the endi i effector frame of the points of assembly with flexible rods. The unknown inputs are inside the transformation matrices TFi that relate the coordinates of the local frames of each rod to the fixed frame. Hence, the unknowns are coupled in the conditions to be fulfilled. As a consequence, the IK problem is more complex. However, depending on the type of inputs to the flexible mechanism, this difficulty varies. For linear inputs, as λi in the 2PFR, the local frames have a constant orientation, and the rotation matrices RFi are constant and known. Hence, the static equilibrium condition is independent of the inputs, and can be used as in the FK to reduce the

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

25

Fig. 10. Fulfillment of geometric conditions of IK for both rods in mode 1 in the parameter space of a 2RFR mechanism.

7 R1 R

6

2

R [N]

5

4

3

2

1

0 -180

-120

-60

0 1

,(

2

60

120

180

-180º)

Fig. 11. Fulfillment of equilibrium conditions of IK for both rods in mode 1 of a 2RFR mechanism.

set of feasible parameters (ψ i , ki ). Moreover, when applying the known output coordinates in Eq. (36) or (37), it is simple to decouple some conditions for assembly from equations that produce the input values λi . On checking residues in the assembly conditions, we can obtain the solution sets (ψ i , ki ) as in the FK, and find the inputs λi in the remaining equations (see Fig. 9). For rotational inputs, as α i in the 2RFR, as unknown inputs are inside the rotation matrices RFi , the problem is strongly coupled, and we must resort to a more complex approach. The idea is to reduce the feasible parameters so that we can easily search for multiple solutions. We first find the sets of parameters (ψ i , ki ) and inputs α i of each rod that produce null residues in the geometric conditions of assembly given by Eqs. (36) and Eqs. (37), resp. For example, in the aforementioned 2RFR mechanism, for a given output position P, i.e., pP , the fulfillment of each of Eqs. (36) can be plotted as distinct surfaces in the parameter space (ψ i , ki , α i ) for each rod (see Figs. 10). Curves of the simultaneous fulfillment of residues correspond to the reduced set of feasible parameters once the assembly conditions are applied.

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

y [m]

y [m]

26

0.1

0.1

0

0

-0.1

-0.1

-0.2

-0.2 -0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.4

-0.3

-0.2

-0.1

0

x [m]

0.1

0.2

0.3

0.4

0.5

x [m]

Fig. 12. Multiple solutions of IK for both rods in mode 1 of a 2RFR mechanism.

10 -6 2

0.5

0.4

UP [J]

0

y [m]

0.3

-2

0.2

-4

0.1

-6 1

0

-3

0.5 10

1 0.5

0

0

-0.5 -0.3

-0.2

-0.1

0

0.1

0.2

0.3

y P [m]

10 -3

-0.5 -1

-1

x [m]

x P [m]

Fig. 13. Unstable pose of a 2RFR mechanism and its corresponding potential function of elastic energy.

Of the feasible parameters, we need to check the ones that fulfill static equilibrium conditions, i.e., Eqs. (30). To illustrate this point, we plot the corresponding magnitudes of force at the end of the rod, obtained for feasible parameters in both rods and against an input angle adjusted to provide the correct direction along which end forces are parallel, i.e., the direction along which to equilibrate (see Fig. 11). At the curve crossings, both forces are balanced, and we obtain multiple solutions of the IK problem shown in Fig. 12. 6. Discussion The aforementioned procedures provide multiple solutions for the FK and IK problems. This is the first step in a complete analysis of the position problem, where the work space and the joint space are evaluated in detail. The solutions need to be evaluated regarding singularities and instabilities, which requires finding the Jacobian and stiffness matrices. In rigid-link parallel manipulators, as the Jacobian is defined analytically, we can obtain the geometric conditions for singular poses. In continuum parallel manipulators, this is not the case, as Jacobians at a pose are found numerically upon the application of position analysis with small displacements around:

pP = JIK · θ θ = JF K · pP

(38)

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

y [m]

y [m]

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

-1 -1

-0.5

0

0.5

-1 -1

1

-0.5

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

-0.5

0

0.5

-1 -1

1

-0.5

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

y [m]

y [m]

1

0

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

0

0.5

-1 -1

1

-0.5

x [m] 1 0.8

0.6

0.6

0.4

0.4

0.2

0.2

y [m]

1

0

1

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

0

0

0.5

1

0

-0.2

-0.5

0.5

x [m]

0.8

-1 -1

0

0

-0.2

-0.5

1

x [m]

x [m]

-1 -1

0.5

0

-0.2

-1 -1

0

x [m] 1

y [m]

y [m]

x [m]

y [m]

27

0.5

1

-1 -1

-0.5

x [m]

Fig. 14. Workspace for a 2RFR mechanism.

0

x [m]

0.5

1

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

1

1

0.8

0.8

0.6

0.6

0.4

0.4

y [m]

y [m]

28

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4 -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

x [m]

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

x [m] Fig. 15. IK and FK singularities.

Similarly, the stiffness matrix K of the mechanism is found numerically by evaluating deformation under infinitesimal loads:

FP = K · pP

(39)

Fig. 13 shows the mechanisms in an unstable pose and the corresponding potential function of elastic energy when inputs remain fixed and the position of the end effector varies in the vicinity. The Hessian of this function is precisely the value of the stiffness matrix in Eq. (39). This solution is unstable since the obtained function represents a saddle point, i.e., there is a direction towards which the function is a maximum. However, instabilities of the system can also occur at any other point of the rods in the system [23]. Therefore, a variational approach should be consider for a comprehensive analysis of this phenomenon as in [24]. Further work in this direction has to be undertaken. Preliminary work space analysis on a 2RFR mechanism showed some interesting results. Fig. 14 contains plots of the different workspaces obtained for the mechanism from the different multiple solutions as plotted on each figure. Outer and inner limits are identified by the IK singularities and curves for the FK singularities. They are found numerically for the time being, and thus numerical issues need to be resolved so that complete singularity curves can be defined. In Fig. 15 on the left, we plot the mechanism in an IK singularity pose at the outer limit of the work space, and on the right, we plot an FK singularity pose. An important issue is that there are continuous transitions between solutions with different buckling modes. This is shown in Fig. 14 with differently colored areas. There is a principal area where the mechanism is deformed with buckling mode 1 in both rods, connected smoothly with the other areas with different buckling modes. Therefore, this cannot be a characteristic to discriminate among solutions. A deeper analysis of this and other related issues will be the subject of future research. 7. Conclusions In this paper, a method was proposed to solve the full position problem of planar parallel continuum mechanisms with the following conditions: The flexible elements are straight rods in free-stress reference configurations, and no distributed forces or moments are considered to act along the length of the rod. The number of solutions, both for forward and inverse kinematics, was higher than in the case of their rigid-link counterparts. Buckling modes contributed to such multiplicity, and can be addressed through the use of a suitable model, in this case, the elliptic integral solution method. Another important aspect is that to some extent, the method needs a numerical approach. This means that although all solutions can be calculated ultimately, their number, even existence, is not known a priori. Future work should extend an analysis of the full position problem to space mechanisms where the deformation in each rod is not necessarily within a plane. To this end, a space rod needs to be modeled using parameters that describe the deformation of a rod in space in a univocal manner. The approach used to address this problem would be similar to the one proposed here, with the only significant difference being the number of parameters involved. Acknowledgments The authors received financial support from the Spanish Government (DPI2015-64450-R), an MEC PhD. grant (FPU15/04422), and the Regional Government of the Basque Country (Project IT949-16).

O. Altuzarra et al. / Mechanism and Machine Theory 132 (2019) 13–29

29

Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechmachtheory. 2018.10.014. References [1] L.L. Howell, Compliant Mechanisms, John Wiley & Sons, 2002. [2] E.G. Merriam, J.E. Jones, S.P. Magleby, L.L. Howell, Monolithic 2 dof fully compliant space pointing mechanism, Mech. Sci. 4 (2) (2013) 381–390, doi:10.5194/ms- 4- 381- 2013. [3] Y. Li, Q. Xu, A novel design and analysis of a 2-dof compliant parallel micromanipulator for nanomanipulation, IEEE Trans. Autom. Sci. Eng. 3 (3) (2006) 247–254, doi:10.1109/TASE.2006.875533. [4] S. Kota, K.J. Lu, Z. Kreiner, B. Trease, J. Arenas, J. Geiger, Design and application of compliant mechanisms for surgical tools, J. Biomech. Eng. 127 (6) (2005) 981–989. doi: 10.1115/1.2056561. [5] C.E. Bryson, D.C. Rucker, Toward parallel continuum manipulators, in: Robotics and Automation (ICRA), 2014 IEEE International Conference on, IEEE, 2014, pp. 778–785. [6] P.L. Anderson, A.W. Mahoney, R.J. Webster, Continuum reconfigurable parallel robots for surgery: shape sensing and state estimation with uncertainty, IEEE Rob. Autom. Lett. 2 (3) (2017) 1617–1624, doi:10.1109/LRA.2017.2678606. [7] J. Burgner-Kahrs, D.C. Rucker, H. Choset, Continuum robots for medical applications: a survey, IEEE Trans. Rob. 31 (6) (2015) 1261–1280, doi:10.1109/ TRO.2015.2489500. [8] J. Till, D.C. Rucker, Elastic rod dynamics: validation of a real-time implicit approach, in: 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2017, pp. 3013–3019, doi:10.1109/IROS.2017.8206139. [9] D.C. Rucker, R.J.W. III, Statics and dynamics of continuum robots with general tendon routing and external loading, IEEE Trans. Rob. 27 (6) (2011) 1033–1044, doi:10.1109/TRO.2011.2160469. [10] L.F. Campanile, A. Hasse, A simple and effective solution of the elastica problem, Proc. Inst. Mech.Eng. Part C 222 (12) (2008) 2513–2516, doi:10.1243/ 09544062JMES1244. [11] S.M. Lyon, L.L. Howell, A simplified pseudo-rigid-body model for fixed-fixed flexible segments, in: Volume 5: 27th Biennial Mechanisms and Robotics Conference, 2002, doi:10.1115/detc2002/mech-34203. [12] J.W. Wittwer, M.S. Baker, L.L. Howell, Simulation, measurement, and asymmetric buckling of thermal microactuators, Sens. Actuators, A 128 (2) (2006) 395–401, doi:10.1016/j.sna.2006.02.014. [13] S. Hasanzadeh, F. Janabi-Sharifi, An efficient static analysis of continuum robots, J. Mech. Robot. 6 (3) (2014). 031011–031011, doi: 10.1115/1.4027305. [14] T. Kugelstadt, E. Schmer, Position and orientation based cosserat rods, in: L. Kavan, C. Wojtan (Eds.), Eurographics/ ACM SIGGRAPH Symposium on Computer Animation, The Eurographics Association, 2016, doi:10.2312/sca.20161234. [15] H. Lang, J. Linn, M. Arnold, Multi-body dynamics simulation of geometrically exact cosserat rods, Multibody Syst. Dyn. 25 (3) (2011) 285–312, doi:10. 1007/s11044-010-9223-x. [16] G.L. Holst, G.H. Teichert, B.D. Jensen, Modeling and experiments of buckling modes and deflection of fixed-guided beams in compliant mechanisms, J. Mech. Des. 133 (5) (2011) 051002, doi:10.1115/1.4003922. [17] A. Zhang, G. Chen, A comprehensive elliptic integral solution to the large deflection problems of thin beams in compliant mechanisms, J. Mech. Robot. 5 (2) (2013). 021006–021006, doi: 10.1115/1.4023558. [18] O. Altuzarra, M. Diez, J. Corral, G. Teoli, M. Ceccarelli, Kinematic analysis of a continuum parallel robot, in: P. Wenger, P. Flores (Eds.), New Trends in Mechanism and Machine Science, Springer International Publishing, Cham, 2017, pp. 173–180. [19] O. Altuzarra, M. Diez, J. Corral, F.J. Campa, Kinematic analysis of a flexible tensegrity robot, in: B. Corves, E.-C. Lovasz, M. Hüsing, I. Maniu, C. Gruescu (Eds.), New Advances in Mechanisms, Mechanical Transmissions and Robotics, Springer International Publishing, Cham, 2017, pp. 457–464. [20] O. Altuzarra, D. Caballero, Q. Zhang, F.J. Campa, Kinematic characteristics of parallel continuum mechanisms, in: J. Lenarcic, V. Parenti-Castelli (Eds.), Advances in Robot Kinematics 2018, Springer International Publishing, Cham, 2019, pp. 293–301. [21] O. Altuzarra, D. Caballero, J. Campa, C. Pinto, Forward and inverse kinematics in 2-dof planar parallel continuum manipulators, (To appear shortly in EuCoMeS 2018), 2018. [22] S.S. Antman, Nonlinear Problems of Elasticity, Springer, 2005. [23] T. Bretl, Z. McCarthy, Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations, Int. J. Rob. Res. 33 (1) (2014) 48–68. [24] J. Till, D.C. Rucker, Elastic stability of cosserat rods and parallel continuum robots, IEEE Trans. Rob. 33 (3) (2017) 718–733, doi:10.1109/TRO.2017. 2664879.