# Research on Stiffness of Multibackbone Continuum Robot Based on Screw Theory and Euler-Bernoulli Beam

Authors: Bin He; Shoulin Xu; Zhipeng Wang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Continuum robots have become a focus for extensive research, since they can work well in complex and confined environments. The main contribution of this paper is to establish a stiffness model of a single section multibackbone continuum robot and analyze the effect of the structural parameters of continuum robot on the overall rotation and translation stiffness. First, a stiffness model which indicates the end configuration of continuum robot under external load is deduced by the screw theory and Euler-Bernoulli beam. Then, the stiffness elements are fully analyzed, therefore, obtaining the influence of the structural parameters of continuum robot on the stiffness elements. Meanwhile, a numerical analysis of stiffness elements is given. Furthermore, the minimum and maximum rotation/translation stiffness are introduced to analyze the effect of the structural parameters of continuum robot on the overall rotation and translation stiffness. Finally, the experiments are used to validate the proposed stiffness model. The experimental results show that the proposed stiffness model of continuum robot is correct and the errors are less than 7%.


---

## Body

## 1. Introduction

Continuum robot is a new kind of bionic robot, inspired by elephant trunks and octopus tentacles. Continuum robot does not have its own joints, which can produce flexible deformation in any part, so it has a strong ability to avoid obstacles and better adapt to the complex unstructured environments. Continuum robots offer a number of potential advantages over the traditional rigid-link robots in applications involving disaster relief [1], industrial applications [2], and medical aid [3].The main types of continuum robots include the rod-driven continuum robot [4], cable-driven continuum robot [5], pneumatic continuum robot [6], and concentric tube continuum robot [7]. Du et al. [8] developed an optimization notched continuum manipulator based on performance evaluation indices. Tian et al. [9] proposed the kinematic analysis of a continuum bionic robot with three flexible actuation rods. Moreover, Bergeles et al. [10] presented optimization framework based on anatomical and surgical task constraints. Li et al. [11] developed a new constrained wire-driven flexible mechanism which has a larger workspace and is more dexterous compared to the existing surgical arms. He et al. [12] proposed a multibackbone continuum robot which is driven by NiTi alloy wire. Trivedi et al. [13] discussed the novel capabilities of soft robots, described examples from nature that provide biological inspiration, surveyed the state of the art, and outlined existing challenges in soft robot design, modeling, fabrication, and control. Webster III and Jones [14] reviewed several modeling approaches in a common frame and notational convention, illustrating that, for piecewise constant curvature, they produce identical results.The analysis of stiffness is an important step in the design and control of continuum robots, since it determines the relationship between the deformation and the force of continuum robots. Selig and Ding [15] applied the screw theory [16] to analyze compliance and stiffness matrices of a beam. Pei et al. [17] studied the compliance of cartwheel flexural hinges. Ding and Dai [18] investigated spatial continuous compliance on both serial and parallel mechanisms based on screw theory and Lie groups, applying eigenvectors and eigenvalues to identify principal screws in the mechanisms. Awtar and Sen [19] proposed a generalized constraint model for compliance and stiffness analysis of 2D beam flexures. Krishnan et al. [20] studied serial and parallel concatenation of building blocks. Tunay [21] proposed a concept of equivalent bending stiffness to establish the kinematic model of continuum robot. Gao et al. [22] established a mathematical model for predicting the loaded posture of a single section continuum manipulator. Qi et al. [23] analyzed the compliance characteristics of a new planar spring continuum robot. Gravagne et al. [24] discussed the dynamics of a planar continuum backbone section, incorporating a large-deflection dynamic model. Trivedi et al. [25] presented a new approach for modeling soft robotic manipulators that incorporates the effect of material nonlinearities and distributed weight and payload. Camarillo et al. [26] proposed a new linear model for transforming desired beam configuration to tendon displacements and vice versa. Fras et al. [27] described the design and implementation of a static model used for position estimation of a flexible modular medical manipulator equipped with optic-fiber based sensors. Till and Rucker [28] exhibited low numerical damping, handled arbitrarily large time steps, and provided an accurate, high-order representation of the rod shape in steady state. Hadi Sadati et al. [29] introduced a novel series solution for variable-curvature Cosserat rod static and Lagrangian dynamic method.Several different performance indices have been proposed for stiffness evaluation, including determinant of stiffness matrix, average stiffness, and minimum and maximum stiffness [30, 31]. The determinant of stiffness matrix does not distinguish specific stiffness values, and the average stiffness cannot give enough information of stiffness values [32]. The minimum and maximum stiffness are minimum and maximum eigenvalue of stiffness matrix, respectively, which can indicate variation range of stiffness values, and the corresponding eigenvector directions represent the minimum and maximum stiffness directions, respectively. However, it is well known that the conventional minimum and maximum stiffness cannot be applied to a 6 × 6 symmetric stiffness/compliance matrix. This is due to the fact that the eigenvalues of a stiffness/compliance matrix are not consistent unit and no way to compare the sizes of eigenvalues. As a consequence, we define the minimum and maximum rotation/translation stiffness, to evaluate the influence of the structural parameters of continuum robot on the rotation/translation stiffness of continuum robot.In the paper, a general loading is represented by a wrench, while a general deformation is represented by a twist. The main contribution is to present a method to establish the stiffness model of a single section multibackbone continuum robot based on screw theory and Euler-Bernoulli beam, as well as analyzing the effect of the structural parameters of continuum robot on the overall rotation and translation stiffness by the minimum and maximum rotation/translation stiffness. The remainder of the paper is organized as follows: In Section2, we give structure overview of a single section multibackbone continuum robot. In Section 3, based on the screw theory and Euler-Bernoulli beam, a brief stiffness model which illustrates the end configuration of the continuum robot under external load is deduced. In Section 4, the stiffness elements are analyzed to fully indicate relationship between wrench and twist of continuum robot and obtain the influence of the structural parameters of continuum robot on the stiffness elements. In addition, the minimum and maximum rotation/translation stiffness are introduced to evaluate the effect of the structural parameters of continuum robot on the overall rotation and translation stiffness. In Section 5, the proposed stiffness model is validated by the experiments. The experimental results indicate that proposed stiffness model of continuum robot is correct. In Section 6, some discussions are given. In Section 7, conclusions are provided.


## 2. Overview of a Single Section Multibackbone Continuum Robot

The simplified structure of a single section multibackbone continuum robot is shown in Figure1. The continuum robot is composed of a base disk, several spacer disks, an end disk, and four super elastic NiTi wires as its backbones. The central NiTi wire is the primary backbone and the remaining NiTi wires are the secondary backbones. The primary backbone is attached to all disks, and all disks are distributed in equal distance. The secondary backbones are equidistant from primary backbone and from each other. The secondary backbones are only attached to the end disk and slide through holes in base disk and spacer disks, which have double effect: driving robot to achieve two degrees of freedom bending motion and providing auxiliary support to increase the stiffness of continuum robot.Figure 1
The simplified structure of a single section multibackbone continuum robot.


## 3. Stiffness Model of Continuum Robot

In order to facilitate the derivation of the presented theory, the model assumptions of continuum robot are summarized in this paper as follows:(1)
The disks of robot are thin, so the thickness of the disks is ignored.(2)
The friction between the disks and backbones is neglected.(3)
The continuum robot is under static balance.(4)
The gravity of all the disks and the four backbones, which is ignored, is very smaller.(5)
We assume that the deformation of backbones is sufficiently small so that the principle of Euler-Bernoulli beam theory applies.When the continuum robot is not exerted by an external load and the secondary backbone does not pull and push the continuum robot, then the simplified geometric configuration of continuum robot is shown in Figure2(a) and the length of all the backbones are equal. The global coordinate system O-xyz is placed at the end disk, as we are interested in the motion of the end of continuum robot. Point O is geometric center of end disk, the axis x points from point O to the first secondary backbone, the axis z is perpendicular to the end disk, and the axis y is then defined according to the right-hand rule. Moreover, the coordinate system Ob-xbybzb is established at point Ob, the axis xb is parallel to the axis x, the axis yb is parallel to the axis y, and the axis zb is parallel to the axis z. The secondary backbones are the average interval 120 degrees around the primary backbone. Assuming that the continuum robot is subjected to the external load at the end of continuum robot, the configuration of continuum robot is shown in Figure 2(b).Figure 2
The continuum robot is subjected to the external load at the end of continuum robot.For easy comprehension of the deduction of the presented theory, in this paper, symbols are defined as follows:i: the index of secondary backbones i=1,2,3.r: the distance from the primary backbone to each secondary backbone on the disk.L: the length of primary backbone.L i: the length of secondary backbones i=1,2,3.r b: the radius of backbones.In the framework of a screw theory [14] and the global coordinate system O-xyz, a small deformation of continuum robot is defined as a twist ξ=ϕ,δT in axis coordinates. ξ is an element of Lie algebra se(3) of Lie group SE(3), where ϕ=ϕx,ϕy,ϕzT represents the corresponding rotational deflections around their corresponding axes and δ=δx,δy,δzT reveals the three translational deflections along their corresponding axes. W=τ,fT is defined as the wrench at the end of the continuum robot in ray coordinates, which is an element of the dual space of Lie algebra se∗(3), where τ=τx,τy,τzT represents three moments around the axes x, y, and z, respectively; f=fx,fy,fzT reveals three forces along axes x, y, and z, respectively. According to the elastic theory, wrench W and twist ξ are related by(1)W=KPξ,ξ=CPW,KP=CP-1,where KP and CP are the stiffness and compliance matrices of continuum robot in the global coordinate system O-xyz, respectively.In the following, we derive stiffness model of continuum robot. First, the compliance matrix of primary backbone in the global coordinate systemO-xyz can be derived by screw theory [17] and Euler-Bernoulli beam [13], which is expressed as(2)C0=LEIx000-L22EIx00LEIy0L22EIy0000LGJ0000L22EIy0L33EIy00-L22EIx000L33EIx000000LEA,where A=πrb2 is the area of the cross section, Ix=1/4πrb4 and Iy=1/4πrb4 are the area moments, J=1/2πrb4 is the torsional moment of inertia, E denotes the elastic module of the material, and G denotes the shear module of the material.For convenience, by extracting common factor, (2) is simplified as(3)C0=LEIx1000-L20010L2000012μ0000L20L2300-L2000L23000000rb24,where(4)μ=GE=121+νand ν is the Poisson’s ratio.It is assumed that theK0 is stiffness matrix of primary backbone in the global coordinate system O-xyz, and then the stiffness matrix K0 is calculated as(5)K0=C0-1.The local coordinate systemOi-xiyizi is established for the i(i=1,2,3)th secondary backbone at the end disk, as shown in Figure 2(a). Point Oi is the intersection point of the ith secondary backbone and the end disk, the axis xi is parallel to axis x, the axis yi is parallel to axis y, and the axis zi is parallel to axis z. Suppose that Csi is the compliance matrix of the ith secondary backbone in the local coordinate system Oi-xiyizi, since Li=L(i=1,2,3), which can be obtained as(6)Csi=LEIx1000-L20010L2000012μ0000L20L2300-L2000L23000000rb24,where i=1,2,3.Since compliance matrixCsi(i=1,2,3) is represented in the different local coordinate system, in order to establish the stiffness model of continuum robot, we should transform the compliance matrix Csi(i=1,2,3) to the global coordinate system by coordinate transformation. Here, to shift the local coordinate system Oi-xiyizi of each secondary backbone into the global coordinate system O-xyz, an adjoint action of Lie group SE(3) on its Lie algebra is introduced by a 6 × 6 matrix representation [21]:(7)Adg=R0DRR,where R is a 3 × 3 rotation matrix representing the orientation of the local coordinate system Oi-xiyizi, relative to the global coordinate system O-xyz; for instance, the 1st secondary backbone R is a unit matrix. D is a skew-symmetric matrix spanned by the position vector d from the origin of frame O-xyz to the origin of frame Oi-xiyizi; for example, the 1st secondary backbone d=r,0,0.Here, the adjoint transformation matrix of theith(i=1,2,3) secondary backbone is given by (8)Adgi=I0DiII,Adgi∈Se3,where(9)I=100010001,D1=00000-r0r0,D2=0022r0022r-22r-22r0,D3=00-22r0022r22r-22r0.In addition, the transpose of the adjoint transformation matrix is given by(10)AdgiT=I-ITDi0I.It is supposed thatCi is the compliance matrix of the ith secondary backbone in the global coordinate system O-xyz and ξi′ and Wi′ are twist and wrench at the end of the ith secondary backbone in the global coordinate system O-xyz, respectively. ξi and Wi are twist and wrench at the end of the ith secondary backbone in the local coordinate system Oi-xiyizi, yielding(11)ξi′=CiWi′,ξi=CsiWi.Moreover, the coordinates of a twist and a wrench in the global coordinate systemO-xyz are calculated as [16](12)ξi′=Adgiξi,Wi′=Adgi-TWi.To obtain the compliance matrixCi in the global coordinate system O-xyz, we deduct it as follows based on (11) and (12):(13)ξi′=Adgiξi=AdgiCsiWi=AdgiCsiAdgiTWi′.Accordingly, the compliance matrixCi can be obtained according to the relation(14)Ci=AdgiCsiAdgiT,where i=1,2,3.Similarly, we can derive the stiffness matrixKi of the i(i=1,2,3)th secondary backbone in the global coordinate system O-xyz as(15)Ki=Ci-1=AdgiCsiAdgiT-1.Since the continuum robot is composed of four NiTi wires in parallel, if all deformations are presented into the same global coordinate systemO-xyz, then the overall stiffness matrix of the parallel flexible mechanism is the sum of the stiffness matrix of each flexible mechanism [16]; the overall stiffness matrix KP of continuum robot can be obtained as(16)KP=∑i=03Ki,where(17)K0=EIxL40006L0040-6L00002μ0000-6L012L2006L00012L20000004rb2,K1=EIxL406rL06L004rb2-12r2rb2-4r20-6L0-4rrb2-4r26rL02μL2+12r2L2012rL200-6L012L2006L012rL2012L200-4rrb2-4r20004rb2-4r2,K2=EIxL4rb2+2r2rb22r2rb2-32rL06L22rrb22r2rb24rb2+2r2rb232rL-6L022rrb2-32rL32rL2μL2+12r2L2-62rL2-62rL200-6L-62rL212L2006L0-62rL2012L2022rrb222rrb20004rb2,K3=EIxL4rb2+2r2rb2-2r2rb2-32rL06L-22rrb2-2r2rb24rb2+2r2rb2-32rL-6L022rrb2-32rL-32rL2μL2+12r2L262rL2-62rL200-6L62rL212L2006L0-62rL2012L20-22rrb222rrb20004rb2.Substituting (17) into (16) yields(18)KP=kpij6×6=EIxLk110k130k1500k220k240k26k310k330k3500k420k4400k510k530k5500k62000k66,where i=1,2,3,4,5,6; j=1,2,3,4,5,6; and(19)k11=16+4r2rb2,k13=k31=61-2rLk15=k51=24L,k22=16+8r22r2-rb2rb24r2-rb2k24=k42=-24L,k26=k62=42rrb2+4r4r2-rb2k33=8μ+36r2L2,k35=k53=121-2rL2k44=k55=48L2,k66=12rb2-44r2-rb2.It can be seen thatKP in (18) is determined by six independent design parameters, E, Ix, L, r, rb, and μ. Simultaneously, the stiffness matrix KP gives the relationship between the wrench and twist at the end of the continuum robot. Furthermore, the compliance matrix is the inverse of the stiffness matrix which is calculated as(20)CP=LEIxc11000c1500c220c240c2600c330c3500c420c440c46c510c530c5500c620c640c66,where(21)c11=Γ,c15=c51=-LΓ2c22=1.2r2rb2-4rb4σ×105c24=c42=2Lrb23r2-rb2σ×105c26=c62=2c46=2c64=125rb2r313rb2-4452r2σc33=2.5L2g×105,c35=c53=25875L2rg,c44=q48σc55=e48gr2+rb2,c66=sσ,Γ=rb24r2+rb2σ=1.665r2rb2-1.6rb4+1.652r4×106g=2μL2-8.871r2×106q=L216.065r2rb2-6.4rb4+1.652r4×106e=L28μL2rb2+2μL2r2+35.614r2rb2+9r4×106s=2r2rb4-rb6+4r4rb2×105.


## 4. Stiffness Analysis of Continuum Robot

Based on the stiffness model of continuum robot in Section3, we will first analyze the effect of the structural parameters of continuum robot on the stiffness elements in this section. Then, the influence of the structural parameters on the overall rotation and translation stiffness is analyzed. Finally, to more intuitively understand the stiffness elements of continuum robot, the stiffness elements are analyzed by the numerical example.


### 4.1. Analysis of Stiffness Elements

Each element in the stiffness matrix reflects the mapping relationship between each element in the wrench and each element in the twist. Hence, through the analysis of the stiffness elements of stiffness matrixKP in (18), we can obtain the following conclusions.( 1 ) The rotational stiffness around axes x, y, and z is calculated as(22)kP11=τxϕx=EIxLk11=πErb24rb2+r2L(23)kP22=τyϕy=EIxLk22=πErb24r4+14r2rb2-4rb4L4r2-rb2(24)kP33=τzϕz=EIxLk33=πErb42μL2+9r2L3.By (22), (23), and (24), assuming that r and rb are constant values, it can be seen that when length L increases, with a decreasing of the rotational stiffness kP11, kP22, and kP33, obviously, kP11 and kP22 are inversely proportional to L. By (22) and (24), supposing that L and r are constant values, when length rb increases, kP11 and kP33 also increase. Assuming that L and rb are constant values, when length r increases, kP11 and kP33 also increase. By (23), assuming that L and rb are constant values, when length r increases, kP22 also increases, as shown Figure 3.Figure 3
Relationship betweenkP22 and r.( 2 ) The translational stiffness along axes x, y, and z is expressed as(25)kP44=fxδx=EIxLk44=12πErb4L3(26)kP55=fyδy=EIxLk55=12πErb4L3(27)kP66=fzδz=EIxLk66=4πErb23r2-rb2L4r2-rb2.By (25) and (26), changing r has no effect on the kP44 and kP55. Assuming that rb is constant value, when length L increases, with a decreasing of the translational stiffness kP44 and kP55, which are inversely proportional to L3. Assume that rb is constant value, when length L increases, with decreasing of the translational stiffness kP44 and kP55, which are inversely proportional to L3. Supposing that the length L is constant value, when length rb increases, kP44 and kP55 also increase, which are directly proportional to rb4. By (27), assuming that r and rb are constant values, it can be seen that when length L increases, but kP66 decreases, obviously, kP66 is inversely proportional to L. Supposing that L and rb are constant values, when length r increases, kP66 decreases. By (18), the translational stiffness of kP44, kP55, and kP66 is decoupled; that is, the forces fx, fy, and fz only produce translation along axes x, y, and z, respectively.( 3 ) The rotational stiffness which is produced by force fx is calculated as(28)kP42=fxϕy=EIxLk42=-6πErb4L2.By (28), fx only produce rotation around axis y. Dividing (25) by (28) and taking the absolute value, which is defined as ε1, yield(29)ε1=kP44kP42=ϕyδx=2L.By (29), if length L increases, ε1 decreases. The result shows that when the force fx is a constant value and length L increases, the value ϕy/δx decreases.( 4 ) The rotational stiffness which is produced by force fy is calculated as(30)kP51=fyϕx=EIxLk51=6πErb4L2(31)kP53=fyϕz=EIxLk53=31-2πErb4rL3.By (30) and (31), fy produce two rotational stiffness elements. One is the rotation stiffness kP51 around the axis x, which is a function of the variables L and rb, The effects of variables L and rb on kP51 can be easily obtained by (30). The other is the rotational stiffness kP53 around the axis z, which is a function of the variables L, r, and rb, The effects of variables L, r, and rb on kP53 can be easily given by (31). Equation (26) is divided by (30), which is defined as ε2. Equation (26) is divided by (31), taking the absolute value, which is defined as ε3. Equation (31) is divided by (30), taking the absolute value, which is defined as ε4. This yields(32)ε2=kP55kP51=ϕxδy=2L(33)ε3=kP55kP53=ϕzδy=41+2r(34)ε4=kP53kP51=ϕxϕz=2-1r2L.By (32), when length L increases, ε2 decreases. The result explains that when the force fy is a constant value and length L increases, the value ϕx/δy decreases. By (33), when length r increases, ε3 decreases, indicating that when the force fy is a constant value and length r increases, the value ϕz/δy decreases. By (34), assuming that L is constant value, it can be seen that when r increases, ε4 also increases. The result indicates that when fy, L are constant values and r increases, the value ϕx/ϕz increases. Supposing that r is constant value, when L increases, ε5 also decreases. Indicating that when fy, r are a constant values and L increases, the value ϕx/ϕz decreases.( 5 ) The rotational stiffness which is produced by force fz is calculated as(35)kP62=fzϕy=EIxLk62=πErb2r42r2-2-1rb2L4r2-rb2.Dividing (27) by (35), which is defined as ε5, we obtain(36)ε5=kP66kP62=ϕyδz=43r2-rb2r42r2-2-1rb2.Whenrb≪r, (36) is simplified as(37)ε5=kP66kP62≈32r.By (37), when rb≪r and r increases, ε5 decreases. The result explains that when the force fz is a constant value, rb≪r, and r increases, the value ϕy/δz decreases. When r≪rb, (36) is simplified as(38)ε5=kP66kP62⟶+∞.By (38), when r≪rb, ε5 tends to positive infinity. The result indicates that when the force fz is a constant value and r≪rb, the value ϕy/δz tends to positive infinity.


### 4.2. Eigen-Stiffness Analysis of Continuum Robot

Since the eigenvalues of the stiffness matrix in (18) are not consistent unit, there is no way to compare the sizes of eigenvalues; the conventional minimum/maximum stiffness cannot be applied to a 6 × 6 symmetric stiffness/compliance matrix. Here, we define the minimum and maximum rotation/translation stiffness to evaluate the effect of the structural parameters of continuum robot on the rotation/translation stiffness. In order to facilitate the presentation of this definition, we first block the stiffness matrix KP in (18), which can be generally represented in the form of(39)KP=ABBTC,where the symmetric 3 × 3 block matrices A and C denote the pure rotation and translation matrices, B is the coupling one, BT is transpose of the matrix B, and(40)A=EIxLk110k130k220k310k33,B=EIxL0k150k240k260k350,C=EIxLk44000k55000k66.Then, by (1), we obtain(41)τ=Aϕ+Bδf=BTϕ+CδSince the pure translation matrixC is already a diagonal matrix, the eigenvalues of matrix C are EIx/Lk44, EIx/Lk55, and  EIx/Lk66, respectively. Define KTmin=min{EIx/Lk44,EIx/Lk55,EIx/Lk66} and KTmax=maxEIx/Lk44,EIx/Lk55,EIx/Lk66 as the minimum and maximum translation stiffness for the matrix C, respectively. The pure rotation matrix A is symmetric, but not a diagonal matrix; thus, there is an orthogonal matrix Q, making QTAQ a diagonal matrix [33]; that is,(42)QTAQ=λ1000λ2000λ3,where Q=q1,q2,q3∈SO(3) represents an orthogonal matrix whose columns are just the eigenvectors of matrix A and QTAQ=diagλ1,λ2,λ3 is a diagonal matrix consisting of the corresponding eigenvalues of matrix A. We define KRmin=min{λ1,λ2,λ3} and KRmax=max{λ1,λ2,λ3} as the minimum and maximum rotation stiffness for the matrix A, respectively.The matrixC is diagonal matrix, it is easy to obtain the minimum and maximum translation stiffness; therefore, we do not do a lot of analysis here. However, the matrix A is not a diagonal matrix and is coupled, In the following, we analyze the effect of mainly structural parameters L and r on minimum rotation stiffness KRmin and maximum rotation stiffness KRmax. The parameters of the NiTi wire are shown in Table 1.Table 1
The parameters of the NiTi wire.
Parameter Value rb 6 × 10 - 4 m A 1.13 × 10 - 6 m2 Ix 1.017 × 10 - 13 m4 Iy 1.017 × 10 - 13 m4 J 2.034 × 10 - 13 m4 G 2.5 × 10 10 Pa E 6.5 × 10 10 PaFirst, the influence of parameterL on the minimum and maximum rotation stiffness is analyzed. For ease of research, assuming that r is a constant and r=3 mm, the images of minimum and maximum rotation stiffness with the change L are shown in Figure 4.Figure 4
(a) The influence of parameterL on the minimum rotation stiffness. (b) The influence of parameter L on the maximum rotation stiffness.
(a) (b)In Figure4(a), it can be seen that when L increases, KRmin decreases. When L∈50,125 mm, the speed of change of KRmin is bigger. When L∈125,400 mm, the speed of change of KRmin is smaller. In Figure 4(b), when L increases, KRmax decreases. When L∈[50,150] mm, the speed of change of KRmax is bigger. When L∈150,400 mm, the speed of change of KRmax is smaller. Considering the minimum and maximum rotation stiffness comprehensively, if L∈50,125 mm, the influence of L on the minimum and maximum rotation stiffness is bigger. If L∈125,400 mm, the influence of L on the minimum and maximum rotation stiffness is smaller.Next, we analyze the influence of parameterr on the minimum and maximum rotation stiffness. Here, it is supposed that L is a constant and L=200 mm; the minimum and maximum rotation stiffness of the image with the change r are shown in Figure 5.Figure 5
(a) The influence of parameterr on the minimum rotation stiffness. (b) The influence of parameter r on the maximum rotation stiffness.
(a) (b)In Figure5(a), if r increases, KRmin decreases. When r∈1,17 mm, the speed of change of KRmin is bigger. When r∈[17,30] mm, the speed of change of KRmin is smaller. In Figure 5(b), when r increases, KRmax decreases. When r∈1,10 mm, the speed of change of KRmax is bigger. When r∈10,30 mm, the speed of change of KRmin is smaller. In summary, if r∈[1,15] mm, the influence of r on the minimum and maximum rotation stiffness is bigger. If r∈[15,30] mm, the influence of r on the minimum and maximum rotation stiffness is smaller.The block matricesA and C have been analyzed above. Next, the matrix B is researched. In order to facilitate the study, the lemma is first given as follows.Lemma 1.
If a matrixA is an n-order asymmetric matrix, then the matrix A must not be orthogonal to diagonalization.Proof.
Set a matrixA as an n-order asymmetric matrix, and suppose there is an n-order orthogonal matrix P, making(43)PTAP=Λ,where Λ is a diagonal matrix and the diagonal elements are the eigenvalues of the matrix A. Transposing (43) yields(44)PTAPT=ΛT.
Since(45)PTAPT=PTATP,ΛT=Λ,
by (43), (49), and (50), we obtain(46)PTATP=PTAP.
MultiplyingP and its transpose PT of (51) yields(47)PPTATPPT=PPTAPPT.
Since the matrixP is an orthogonal matrix, that is, PPT is a unit matrix, then, this yields(48)AT=A.
Thus, it is obtained that the matrixA is symmetric matrix. It is inconsistent with the hypothetical matrix A which is an n-order asymmetric matrix at the beginning of the proof. Therefore, the matrix A cannot be orthogonal to diagonalization.
The matrixB reflects the relationship between τ and δ, and the units of each element in the matrix B are consistent. However, the matrix B is not a symmetric matrix in (40), by Lemma 1, obtaining that matrix B cannot be orthogonal to diagonalization. Therefore, the matrix B has no minimum eigenvalue and maximum eigenvalue. Here, we cannot define the minimum stiffness and maximum stiffness of matrix B. The transpose matrix BT also cannot define the minimum stiffness and maximum stiffness.


### 4.3. Numerical Analysis of Stiffness Elements

In this section, in order to understand the stiffness elements of continuum robot more intuitively, the stiffness elements of continuum robot are analyzed by the numerical example. Here, we select the structural parameters of continuum robot such asL=210 mm, r=5 mm. The parameters of NiTi wire are shown in Table 1.The results are tabulated in Table2, where the unit for force component is Newton, the unit of moment is N∙m, and the units for translational and rotational displacements are meter and radian, respectively. By analyzing the numerical results, the following conclusions can be obtained.Table 2
Numerical results of stiffness elements.
Stiffnesselement Value τ x / ϕ x 9.26 τ x / ϕ z - 1.86 × 10 - 3 τ x / δ y 3.6 τ y / ϕ y 9.22 τ y / δ x - 3.6 τ y / δ z 2.48 × 10 3 τ z / ϕ x - 1.86 × 10 - 3 τ z / ϕ z 9.76 × 10 - 2 τ z / δ y - 1.77 × 10 - 2 f x / ϕ y - 3.6 f x / δ x 3.43 × 10 f y / ϕ x 3.6 f y / ϕ z - 1.77 × 10 - 2 f y / δ y 3.43 × 10 f z / ϕ y 2.48 × 10 3 f z / δ z 1.05 × 10 6 Major stiffness elements. Elements not listed are zeros; Units:τx,y,z (N⋅m), fx,y,z (N), ϕx,y,z (rad), δx,y,z (m).( 1 ) The rotational stiffness elements τx/ϕx and τy/ϕy are approximately equal. It is indicated that the same moment is, respectively, applied around the axes x and y at the end of continuum robot; then the rotational angles of the end are the approximately equal. Moreover, τx/ϕx is about 94.88 times larger than τy/δz, and τy/ϕy is about 94.47 times larger than τy/δz. It is shown that the rotation around the axes x and y is harder than that around the axis z.( 2 ) The translational stiffness elements fx/δx and fy/δy are equal. It is shown that the same force is, respectively, applied along the axes x and y at the end of continuum robot; then the translational distances of the end are equal. Meanwhile, the translational stiffness element fz/δz=1.05×106 (N/m) is significantly larger than translational stiffness element fx/δx=fy/δy=3.43×10 (N/m). It is indicated that the translation along the axis z is harder than translation along the axes x and y, respectively. It all fits our intuition.


## 4.1. Analysis of Stiffness Elements

Each element in the stiffness matrix reflects the mapping relationship between each element in the wrench and each element in the twist. Hence, through the analysis of the stiffness elements of stiffness matrixKP in (18), we can obtain the following conclusions.( 1 ) The rotational stiffness around axes x, y, and z is calculated as(22)kP11=τxϕx=EIxLk11=πErb24rb2+r2L(23)kP22=τyϕy=EIxLk22=πErb24r4+14r2rb2-4rb4L4r2-rb2(24)kP33=τzϕz=EIxLk33=πErb42μL2+9r2L3.By (22), (23), and (24), assuming that r and rb are constant values, it can be seen that when length L increases, with a decreasing of the rotational stiffness kP11, kP22, and kP33, obviously, kP11 and kP22 are inversely proportional to L. By (22) and (24), supposing that L and r are constant values, when length rb increases, kP11 and kP33 also increase. Assuming that L and rb are constant values, when length r increases, kP11 and kP33 also increase. By (23), assuming that L and rb are constant values, when length r increases, kP22 also increases, as shown Figure 3.Figure 3
Relationship betweenkP22 and r.( 2 ) The translational stiffness along axes x, y, and z is expressed as(25)kP44=fxδx=EIxLk44=12πErb4L3(26)kP55=fyδy=EIxLk55=12πErb4L3(27)kP66=fzδz=EIxLk66=4πErb23r2-rb2L4r2-rb2.By (25) and (26), changing r has no effect on the kP44 and kP55. Assuming that rb is constant value, when length L increases, with a decreasing of the translational stiffness kP44 and kP55, which are inversely proportional to L3. Assume that rb is constant value, when length L increases, with decreasing of the translational stiffness kP44 and kP55, which are inversely proportional to L3. Supposing that the length L is constant value, when length rb increases, kP44 and kP55 also increase, which are directly proportional to rb4. By (27), assuming that r and rb are constant values, it can be seen that when length L increases, but kP66 decreases, obviously, kP66 is inversely proportional to L. Supposing that L and rb are constant values, when length r increases, kP66 decreases. By (18), the translational stiffness of kP44, kP55, and kP66 is decoupled; that is, the forces fx, fy, and fz only produce translation along axes x, y, and z, respectively.( 3 ) The rotational stiffness which is produced by force fx is calculated as(28)kP42=fxϕy=EIxLk42=-6πErb4L2.By (28), fx only produce rotation around axis y. Dividing (25) by (28) and taking the absolute value, which is defined as ε1, yield(29)ε1=kP44kP42=ϕyδx=2L.By (29), if length L increases, ε1 decreases. The result shows that when the force fx is a constant value and length L increases, the value ϕy/δx decreases.( 4 ) The rotational stiffness which is produced by force fy is calculated as(30)kP51=fyϕx=EIxLk51=6πErb4L2(31)kP53=fyϕz=EIxLk53=31-2πErb4rL3.By (30) and (31), fy produce two rotational stiffness elements. One is the rotation stiffness kP51 around the axis x, which is a function of the variables L and rb, The effects of variables L and rb on kP51 can be easily obtained by (30). The other is the rotational stiffness kP53 around the axis z, which is a function of the variables L, r, and rb, The effects of variables L, r, and rb on kP53 can be easily given by (31). Equation (26) is divided by (30), which is defined as ε2. Equation (26) is divided by (31), taking the absolute value, which is defined as ε3. Equation (31) is divided by (30), taking the absolute value, which is defined as ε4. This yields(32)ε2=kP55kP51=ϕxδy=2L(33)ε3=kP55kP53=ϕzδy=41+2r(34)ε4=kP53kP51=ϕxϕz=2-1r2L.By (32), when length L increases, ε2 decreases. The result explains that when the force fy is a constant value and length L increases, the value ϕx/δy decreases. By (33), when length r increases, ε3 decreases, indicating that when the force fy is a constant value and length r increases, the value ϕz/δy decreases. By (34), assuming that L is constant value, it can be seen that when r increases, ε4 also increases. The result indicates that when fy, L are constant values and r increases, the value ϕx/ϕz increases. Supposing that r is constant value, when L increases, ε5 also decreases. Indicating that when fy, r are a constant values and L increases, the value ϕx/ϕz decreases.( 5 ) The rotational stiffness which is produced by force fz is calculated as(35)kP62=fzϕy=EIxLk62=πErb2r42r2-2-1rb2L4r2-rb2.Dividing (27) by (35), which is defined as ε5, we obtain(36)ε5=kP66kP62=ϕyδz=43r2-rb2r42r2-2-1rb2.Whenrb≪r, (36) is simplified as(37)ε5=kP66kP62≈32r.By (37), when rb≪r and r increases, ε5 decreases. The result explains that when the force fz is a constant value, rb≪r, and r increases, the value ϕy/δz decreases. When r≪rb, (36) is simplified as(38)ε5=kP66kP62⟶+∞.By (38), when r≪rb, ε5 tends to positive infinity. The result indicates that when the force fz is a constant value and r≪rb, the value ϕy/δz tends to positive infinity.


## 4.2. Eigen-Stiffness Analysis of Continuum Robot

Since the eigenvalues of the stiffness matrix in (18) are not consistent unit, there is no way to compare the sizes of eigenvalues; the conventional minimum/maximum stiffness cannot be applied to a 6 × 6 symmetric stiffness/compliance matrix. Here, we define the minimum and maximum rotation/translation stiffness to evaluate the effect of the structural parameters of continuum robot on the rotation/translation stiffness. In order to facilitate the presentation of this definition, we first block the stiffness matrix KP in (18), which can be generally represented in the form of(39)KP=ABBTC,where the symmetric 3 × 3 block matrices A and C denote the pure rotation and translation matrices, B is the coupling one, BT is transpose of the matrix B, and(40)A=EIxLk110k130k220k310k33,B=EIxL0k150k240k260k350,C=EIxLk44000k55000k66.Then, by (1), we obtain(41)τ=Aϕ+Bδf=BTϕ+CδSince the pure translation matrixC is already a diagonal matrix, the eigenvalues of matrix C are EIx/Lk44, EIx/Lk55, and  EIx/Lk66, respectively. Define KTmin=min{EIx/Lk44,EIx/Lk55,EIx/Lk66} and KTmax=maxEIx/Lk44,EIx/Lk55,EIx/Lk66 as the minimum and maximum translation stiffness for the matrix C, respectively. The pure rotation matrix A is symmetric, but not a diagonal matrix; thus, there is an orthogonal matrix Q, making QTAQ a diagonal matrix [33]; that is,(42)QTAQ=λ1000λ2000λ3,where Q=q1,q2,q3∈SO(3) represents an orthogonal matrix whose columns are just the eigenvectors of matrix A and QTAQ=diagλ1,λ2,λ3 is a diagonal matrix consisting of the corresponding eigenvalues of matrix A. We define KRmin=min{λ1,λ2,λ3} and KRmax=max{λ1,λ2,λ3} as the minimum and maximum rotation stiffness for the matrix A, respectively.The matrixC is diagonal matrix, it is easy to obtain the minimum and maximum translation stiffness; therefore, we do not do a lot of analysis here. However, the matrix A is not a diagonal matrix and is coupled, In the following, we analyze the effect of mainly structural parameters L and r on minimum rotation stiffness KRmin and maximum rotation stiffness KRmax. The parameters of the NiTi wire are shown in Table 1.Table 1
The parameters of the NiTi wire.
Parameter Value rb 6 × 10 - 4 m A 1.13 × 10 - 6 m2 Ix 1.017 × 10 - 13 m4 Iy 1.017 × 10 - 13 m4 J 2.034 × 10 - 13 m4 G 2.5 × 10 10 Pa E 6.5 × 10 10 PaFirst, the influence of parameterL on the minimum and maximum rotation stiffness is analyzed. For ease of research, assuming that r is a constant and r=3 mm, the images of minimum and maximum rotation stiffness with the change L are shown in Figure 4.Figure 4
(a) The influence of parameterL on the minimum rotation stiffness. (b) The influence of parameter L on the maximum rotation stiffness.
(a) (b)In Figure4(a), it can be seen that when L increases, KRmin decreases. When L∈50,125 mm, the speed of change of KRmin is bigger. When L∈125,400 mm, the speed of change of KRmin is smaller. In Figure 4(b), when L increases, KRmax decreases. When L∈[50,150] mm, the speed of change of KRmax is bigger. When L∈150,400 mm, the speed of change of KRmax is smaller. Considering the minimum and maximum rotation stiffness comprehensively, if L∈50,125 mm, the influence of L on the minimum and maximum rotation stiffness is bigger. If L∈125,400 mm, the influence of L on the minimum and maximum rotation stiffness is smaller.Next, we analyze the influence of parameterr on the minimum and maximum rotation stiffness. Here, it is supposed that L is a constant and L=200 mm; the minimum and maximum rotation stiffness of the image with the change r are shown in Figure 5.Figure 5
(a) The influence of parameterr on the minimum rotation stiffness. (b) The influence of parameter r on the maximum rotation stiffness.
(a) (b)In Figure5(a), if r increases, KRmin decreases. When r∈1,17 mm, the speed of change of KRmin is bigger. When r∈[17,30] mm, the speed of change of KRmin is smaller. In Figure 5(b), when r increases, KRmax decreases. When r∈1,10 mm, the speed of change of KRmax is bigger. When r∈10,30 mm, the speed of change of KRmin is smaller. In summary, if r∈[1,15] mm, the influence of r on the minimum and maximum rotation stiffness is bigger. If r∈[15,30] mm, the influence of r on the minimum and maximum rotation stiffness is smaller.The block matricesA and C have been analyzed above. Next, the matrix B is researched. In order to facilitate the study, the lemma is first given as follows.Lemma 1.
If a matrixA is an n-order asymmetric matrix, then the matrix A must not be orthogonal to diagonalization.Proof.
Set a matrixA as an n-order asymmetric matrix, and suppose there is an n-order orthogonal matrix P, making(43)PTAP=Λ,where Λ is a diagonal matrix and the diagonal elements are the eigenvalues of the matrix A. Transposing (43) yields(44)PTAPT=ΛT.
Since(45)PTAPT=PTATP,ΛT=Λ,
by (43), (49), and (50), we obtain(46)PTATP=PTAP.
MultiplyingP and its transpose PT of (51) yields(47)PPTATPPT=PPTAPPT.
Since the matrixP is an orthogonal matrix, that is, PPT is a unit matrix, then, this yields(48)AT=A.
Thus, it is obtained that the matrixA is symmetric matrix. It is inconsistent with the hypothetical matrix A which is an n-order asymmetric matrix at the beginning of the proof. Therefore, the matrix A cannot be orthogonal to diagonalization.
The matrixB reflects the relationship between τ and δ, and the units of each element in the matrix B are consistent. However, the matrix B is not a symmetric matrix in (40), by Lemma 1, obtaining that matrix B cannot be orthogonal to diagonalization. Therefore, the matrix B has no minimum eigenvalue and maximum eigenvalue. Here, we cannot define the minimum stiffness and maximum stiffness of matrix B. The transpose matrix BT also cannot define the minimum stiffness and maximum stiffness.


## 4.3. Numerical Analysis of Stiffness Elements

In this section, in order to understand the stiffness elements of continuum robot more intuitively, the stiffness elements of continuum robot are analyzed by the numerical example. Here, we select the structural parameters of continuum robot such asL=210 mm, r=5 mm. The parameters of NiTi wire are shown in Table 1.The results are tabulated in Table2, where the unit for force component is Newton, the unit of moment is N∙m, and the units for translational and rotational displacements are meter and radian, respectively. By analyzing the numerical results, the following conclusions can be obtained.Table 2
Numerical results of stiffness elements.
Stiffnesselement Value τ x / ϕ x 9.26 τ x / ϕ z - 1.86 × 10 - 3 τ x / δ y 3.6 τ y / ϕ y 9.22 τ y / δ x - 3.6 τ y / δ z 2.48 × 10 3 τ z / ϕ x - 1.86 × 10 - 3 τ z / ϕ z 9.76 × 10 - 2 τ z / δ y - 1.77 × 10 - 2 f x / ϕ y - 3.6 f x / δ x 3.43 × 10 f y / ϕ x 3.6 f y / ϕ z - 1.77 × 10 - 2 f y / δ y 3.43 × 10 f z / ϕ y 2.48 × 10 3 f z / δ z 1.05 × 10 6 Major stiffness elements. Elements not listed are zeros; Units:τx,y,z (N⋅m), fx,y,z (N), ϕx,y,z (rad), δx,y,z (m).( 1 ) The rotational stiffness elements τx/ϕx and τy/ϕy are approximately equal. It is indicated that the same moment is, respectively, applied around the axes x and y at the end of continuum robot; then the rotational angles of the end are the approximately equal. Moreover, τx/ϕx is about 94.88 times larger than τy/δz, and τy/ϕy is about 94.47 times larger than τy/δz. It is shown that the rotation around the axes x and y is harder than that around the axis z.( 2 ) The translational stiffness elements fx/δx and fy/δy are equal. It is shown that the same force is, respectively, applied along the axes x and y at the end of continuum robot; then the translational distances of the end are equal. Meanwhile, the translational stiffness element fz/δz=1.05×106 (N/m) is significantly larger than translational stiffness element fx/δx=fy/δy=3.43×10 (N/m). It is indicated that the translation along the axis z is harder than translation along the axes x and y, respectively. It all fits our intuition.


## 5. Experiments

The proposed stiffness model is validated by a single section NiTi wire driving continuum robot, as shown in Figure6. The structural parameters of continuum robot are L=210 mm, r=5 mm. The parameters of NiTi wire are shown in Table 1. The disks are made of 3D printed machine (EDEN-260V) with hard material (VeroClear RGD810), and the disks are very hard. The secondary backbones are connected with three stepping motors, which produce the driving force to bend the continuum robot. The tension sensor is used to record the driving force which is produced by motors.Figure 6
A single section NiTi wire continuum robot.The experimental setup is shown as Figure7, which consists of a continuum robot, a control cabinet, a teaching apparatus, a robotic arm, a computer, and a force/torque sensor. The robotic arm (STEP-SD500) is used to provide an external load at the end of continuum robot, and the 6 DOF force/torque commercial sensor (OPTOFORCE-HEX-70-XE-450) is mounted at the end of robotic arm to record the external load force.Figure 7
Experimental setup for stiffness experiments.Here, we facilitate the experimental operations by verifying the compliance matrix to indicate that the stiffness matrix is correct. Meanwhile, we pay more attention to the relations between forces and deformations at the end of continuum robot along the three translational directions. Therefore, we operate the experiment to verify the relations between forces and deformations in three translational directions. It is assumed that the external forces applied at the end of continuum robot are twistsξ1, ξ2, and ξ3, where ξ1=0,0,0,fx,0,0T, ξ2=0,0,0,0,fy,0T, and ξ3=0,0,0,0,0,fzT. Substituting twists ξ1, ξ2, and ξ3 into (20), respectively, we obtain external load force and translational distance that meet the following relationship:(49)δx=2.785fx(50)δy=2.624fy(51)δz=9.5×10-4fz,where fx (N), fy (N), fz (N), δx (cm), δy (cm), and δz (mm).Next, through the experiments we verify that(49), (50), and (51) are correct. The experiments were done at low moving speed to reduce the effect of vibration. First, the experiments along the axis x are shown in Figure 8. The force and translational distance at the end of continuum robot is recorded by 6 DOF force/torque commercial sensor and the operation of robotic arm, respectively.Figure 8
The experiments along the axisx.
(a) (b) (c) (d) (e) (f)The experimental values and the functional image of(49) are shown in Figure 9. The experimental values are distributed around the theoretical image, and the errors between the theoretical values and the experimental values are no more than 4%, thus, indicating that (49) is correct.Figure 9
Theoretical and experimental values along the axisx.Then, the experiments along the axisy are shown in Figure 10.Figure 10
The experiments along the axisy.
(a) (b) (c) (d) (e) (f)The experimental values and the functional image of(50) are shown in Figure 11. The errors between the theoretical values and the experimental values are no more than 5%, therefore, demonstrating that (50) is correct.Figure 11
Theoretical and experimental values along the axisy.Finally, the experiments along the axisz are shown in Figure 12.Figure 12
The experiments along the axisz.
(a) (b) (c) (d) (e) (f)The experimental values and the functional image of(51) are shown in Figure 13. The errors between the theoretical values and the experimental values are no more than 7%, thence, showing that (51) is correct.Figure 13
Theoretical and experimental values along the axisz.Through the above experiments, we verify that the translational stiffness of continuum robot along the axesx, y, and z, respectively, is correct. By comparing the translational experimental values along the axes x and y, respectively, if the external force fx=fy, the translational distances δx and δy are approximately equal. By comparing the translational experimental values along the axes x and z, respectively, if the external force fx=fz, the translational distance δx is much larger than translational distance δz; it also fits our intuitive perception.


## 6. Discussions

In this paper, we mainly establish the stiffness model of a single section multibackbone continuum robot, with the initial configuration of continuum robot being a straight line andLi=L(i=1,2,3). If the continuum robot is multisection and the initial configuration of continuum robot is a straight line, then the stiffness model of a multisection continuum robot can be studied by each lower section as a payload of its upper sections. In addition, the stiffness model of continuum robot is based on the initial configuration which is a straight line. If the initial configuration of continuum robot is a bending shape, then the length of the four NiTi wires of continuum robot is not equal. Thus, the stiffness model in (18) is not suitable for analyzing the load model of continuum robot. Simultaneously, it should be pointed out that the deformations of all backbones should be small; otherwise, the elastic properties of the backbone are damaged; then the stiffness model in (18) is also not suitable for analyzing the load model of continuum robot.In our model, the6×6 stiffness matrix in (18) indicates relationship between force/moment and motion/rotation at the end of continuum robot. If the end of continuum robot is subjected to force fx along the forward direction of the axis x, the shape of continuum robot at this time is defined as A. Next, on the basis of shape A, the force fy is applied to the end of continuum robot along the forward direction of the axis y; then the final bending shape of continuum robot is defined as B. The shape B is equivalent to the resultant force of forces fx and fy, which is applied to the end of continuum robot; the end position of continuum robot with the shape B can be obtained by substituting fx and fy into the stiffness matrix of (18). Furthermore, on the basis of the shape A, if the force fz is applied to the end of continuum robot along the opposite direction of the axis z, then the final bending shape of continuum robot is defined as C. However, the end position of continuum robot with the shape C cannot be obtained by substituting fx and fz into the stiffness matrix of (18). Therefore the stiffness element kp66=fz/δz is only applicable to the case in which the continuum robot does not have any deformations. Through the same analytical method, the stiffness element kp33=τz/ϕz is only applicable to the case in which the continuum robot does not have any deformations. However, the other main stiffness elements kp11, kp22, kp44, and kp55 are applicable to the case in which the continuum robot is deformed and nondeformed.


## 7. Conclusions

In this paper, the main contribution is to propose a simplified method to establish the stiffness model of a single section multibackbone continuum robot, so as to obtain the relationship between force and deformation at the end of continuum. First, the stiffness model of continuum robot is established by means of the screw theory and Euler-Bernoulli beam. Second, the stiffness elements of the stiffness model are analyzed, obtaining the influence of the structural parameters of continuum robot on the stiffness elements. Third, the minimum and maximum rotation/translation stiffness are defined to analyze the effect of the structural parameters of continuum robot on the overall rotation and translation stiffness. In addition, the stiffness elements of continuum robot are analyzed by the numerical example. Finally, the stiffness model is verified by experiments; simultaneously, the experimental results indicate that the proposed stiffness model of continuum robot is correct. The research of this paper fully indicates the relationship between the force and deformation at the end of continuum robot and, thus, provides an effective theoretical basis for the real-time control of continuum robot.



---

*Source: 6910468-2018-05-16.xml*