Stiffness evaluation of continuum robots based on the energy method and castigliano’s second theorem

Introduction: Continuum robots are studied and applied in neurosurgery due to their high flexibility and adaptability. The basic performance of continuum is mainly evaluated by stiffness, but there is no systematic and universal evaluation system.

Methods: In this paper, a general experimental platform for continuum robots is designed, based on which the fundamental performance of the notched continuum robot used in neurosurgery is evaluated. The continuum stiffness evaluation method based on energy method and Castigliano’s second theorem is proposed. By solving the internal force and energy of the notched continuum in sections, the stiffness model of single-segment and double-segment series continuum is established. The relationship between the stiffness of the continuum and the bending angle is obtained.

Results: The simulation and experimental results show that under the condition of small deformation angle, the spatial stiffness model obtained by strain energy basically conforms to the actual model, which verifies the correctness and rationality of the stiffness calculation method proposed in this paper.

Discussion: This paper is of significant importance to promote the performance evaluation and optimization of continuum.

1 Introduction

Continuum robots are characterized by elastic structures and infinite degrees of freedom, lacking discrete joints and rigid links typical of traditional rigid robots (Webster III and Jones, 2010). They exhibit a high degree of dexterity not found in traditional robots, coupled with strong adaptability to workspace constraints. Consequently, they find extensive applications in specialized fields such as medical devices, search and rescue, demonstrating exceptional performance (He et al., 2018). Continuum robots have shown great potential in neurosurgery, such as cerebral hemorrhage aspiration (Burgner et al., 2013), transnasal skull base tumor resection (Zhang et al., 2022), and other operations. Currently, continuum robot designs include hinge-based, segmented, concentric tube, and notched configurations among others (Li and Du, 2013).

Hinge-based continuum robots connect their joints through ball sockets or hinges. Professor Zheng Li from The Chinese University of Hong Kong has developed corresponding kinematic and workspace models based on the assumption that joint curvatures along the continuum bending curve are equal (Li et al., 2016). While these robots exhibit out-standing flexibility, friction between ball sockets can cause noticeable hysteresis and uneven bending during motion. Segmented continuum robots, powered by drive cables with central elastic rods providing compensatory force, approximate constant curvature behavior during bending but may occlude the central channel. In 2006, researchers in-cluding Nabil Simaan from Columbia University and Kai Xu from Shanghai Jiao Tong University proposed using Nitinol alloy tubes instead of drive cables to increase robot stiffness (Xu and Simaan, 2006; Xu and Simaan, 2008; Simaan et al., 2009). Similarly, researchers such as Bin He from Tongji University utilized superelastic Nitinol alloys to design a three-backbone continuum robot, establishing kinematic and dynamic models (He et al., 2013). Concentric tube continuum robots, introduced by Robert James Webster III from Vanderbilt University, achieve bending motion through pre-curved elastic sheaths and inner tube feeding, applied in medical fields for complex surgical operations. However, inherent elastic forces between the inner and outer tubes contribute to significant errors in the robot’s kinematic model, challenging precise end-effector control (Webster et al., 2008). Notched continuum robots, initially proposed by Johns Hopkins University (Badescu and Mavroidis, 2004), use linear drive mechanisms to form joints at incisions, compensating for uneven bending in the continuum, providing improved stiffness. Researchers from Harbin Institute of Technology (Wilkening et al., 2017; Gao et al., 2016a) investigated triangular and square notched continuum robots, establishing their mechanical and kinematic models. Zhang et al. (2024b), Zhang et al. (2024a) propose to mix the rectangular cut continuum with the concentric tube continuum and establish the kinematic models.

In current research, the focus on continuum robots primarily centers around structural innovations, precise control, and performance evaluation. Due to their unique operating environment, continuum robots are required not only high flexibility but also a certain degree of rigidity, as they must balance their load capacity and motion precision while ensuring adequate adaptability and safety (Gu and Ren, 2023; Lin et al., 2024).

Stiffness analysis is a crucial aspect of the design and control of continuum robots, determining the relationship between deformation and forces in these systems. Numerous experts have conducted in-depth studies on stiffness metrics. Selig and Ding utilized screw theory (Dai, 2012) to analyze the flexibility and stiffness matrices of beams. Pei et al. (2009a) investigated the flexibility of wheel flexible joints. Ding and Dai (Pei et al., 2009b) explored the spatial continuum flexibility of serial and parallel mechanisms based on screw theory and Lie group theory, employing eigenvectors and eigenvalues to identify principal screws within mechanisms. Awtar and Sen (2010) proposed a generalized constraint model for analyzing the flexural flexibility and stiffness of 2D beams. Tunay introduced the concept of equivalent bending stiffness. Gao et al. (2016b) developed a mathematical model to predict the load postures of single cross-section continuum robotic arms. Qi et al. (2015) analyzed the flexibility characteristics of a novel planar spring continuum robot. Gravagne et al. (2003) discussed the dynamics of planar continuum backbone sections using a large deflection dynamic model. Trivedi et al. (2007) introduced a novel modeling approach for flexible robotic arms that incorporates material nonlinearity and the effects of distributed weight and payload. Camarillo et al. (2008) proposed a new linear model for converting desired beam configurations to tendon displacements and vice versa. Fraś et al. (2014) described the design and implementation of a static model used for position estimation of modular medical robotic arms equipped with fiber optic sensors. Sadati et al. (2017) presented a series solution method for static and Lagrangian dynamic analysis of a new variable curvature Cosserat rod. However, a comprehensive and universally applicable evaluation framework is currently lacking.

This paper designs a universal continuum experimental platform, which can realize the driving of multiple segments of cable-driven continuum of different sizes. Based on the energy method and Castigliano’s second theorem, the stiffness of the notched continuum introduced above is modeled. The contributions of this paper as follows:

• This paper designs a universal continuum platform, based on which the stiffness model of the double-segment notched continuum is verified and the stiffness of the serial continuum under different driving strategies is tested.

• Utilizing energy methods and Castigliano’s second theorem, the internal forces and energies of incision-type continuum bodies are sequentially resolved, establishing stiffness models for both single-section and dual-section series-connected continua. The relationship between continuum stiffness and bending angle is derived.

• Simulation and experimental results indicate that, under small deformation conditions, the spatial stiffness models derived from strain energy align closely with practical models. The stiffness relationships of dual-section series-connected notched under varying numbers of driving cables are established as follows: External 4 + Internal 4 $>$ External 2 + Internal 4 $>$ External 4 + Internal 2 $>$ External 2 + Internal 2.

The remaining parts of this paper are as follows: Section 2 describes the structure of the universal experimental platform for continuum robot. Section 3 establishes stiffness models for single-section and dual-section series-connected notched continuum. Section 4 conducts simulations of the proposed stiffness analysis models and validates their effectiveness through experimentation. Section 5 provides a summary of the entire paper.

2 The structure of the universal platform

To realize the driving and testing of continuum of various sizes and different numbers of drive cables, this paper designs a universal continuum experimental plat-form, as shown in Figure 1A. The overall structure of the continuum experimental platform can be divided into two parts, namely, the main driving mechanism and the linear module. The main driving mechanism consists of 4 sets of lead screw slider mechanisms, a front base, a rear base and a front wire guide mechanism. The pressure sensor is installed on the slide. One end of the pressure sensor is installed together with the slide, and the other end is installed with a V-groove guide wheel to transmit the drive cable and measure the pressure of the drive cables. The structure of the double-segment series continuum is shown in Figure 1B.

Figure 1

Figure 1. The structure of the universal platform. (A) The universal continuum experimental platform. The platform consists of three parts: continuum structure, external continuum drive mechanism and internal continuum drive mechanism. (B) The specific structure of the dual-segment notched continuum robot.

The structure of the double-segment series continuum is shown in Figure 1B. The structure of the double-segment series continuum is shown in Figure 1B. The parameters of the external and internal continuum are shown in Table 1.

Table 1

Table 1. The parameters of the external and internal continuum.

3 Stiffness evaluation in series-connected continuum

This section models the stiffness of the notched continuum based on the energy method (Straughan, 2013) and the Castigliano’s second theorem (Eastwood et al., 2016). Additionally, the relationship between the spatial stiffness of serial continuum and bending angles is obtained. Finally, the stiffness of serial continuum under four different driving strategies are compared.

3.1 Energy method and stiffness

When calculating the end displacement of the continuum, this chapter adopts the energy method. Assuming the elastic body has no rigid displacements under the constraints of the supports, and is subjected to n external forces, the strain energy stored in the system due to these external forces is denoted as $V_{\epsilon }$. If an incremental force $\mathrm{\Delta }{F}_i$ is applied to one of these forces $F_i$, the total strain energy of the system becomes $V_{\epsilon }+\mathrm{\Delta }{V}_{\epsilon }$. Altering the sequence of force application involves first applying $\mathrm{\Delta }{F}_i$ to the elastic body before applying the external force. Applying $\mathrm{\Delta }{F}_i$ results in a displacement $\mathrm{\Delta }\left(\mathrm{\Delta }{F}_i\right)$ in the direction of this increment. Therefore, the work done by $\mathrm{\Delta }{F}_i$ is represented as W.

$V_{\epsilon }+\mathrm{\Delta }{V}_{\epsilon }=V_{\epsilon }+\mathrm{\Delta }{F}_i{\mathrm{\Delta }}_i+\frac{1}{2}\mathrm{\Delta }{F}_i\cdot \mathrm{\Delta }\left({\mathrm{\Delta }}_i\right)(1)$

Neglecting higher-order terms, it can be obtained that

${\mathrm{\Delta }}_i=\frac{\mathrm{\Delta }{V}_{\epsilon }}{\mathrm{\Delta }{F}_{\text{i}}}(2)$

From the Equation 1 and Equation 2, it can be seen that the partial derivative of strain energy with respect to any load $F_i$ equals the displacement of the point of action of $F_i$ in the direction of $F_i$. This is known as Castigliano’s second theorem, applicable exclusively to linear elastic structures. Additionally, formulas are provided below for the incremental strain energy $\mathrm{\Delta }{V}_{\epsilon }$ induced by forces of various natures. The incremental strain energy induced by the normal support force $F_N$ is shown in Equation 3.

$\mathrm{\Delta }{V_{\epsilon }}^{F_N}=\frac{{F_N}^2}{2EA}ds(3)$

The incremental strain energy induced by the shear force $F_s$ is shown in Equation 4.

$\mathrm{\Delta }{V_{\epsilon }}^{F_s}=\frac{k{F_s}^2}{2GA}ds(4)$

The incremental strain energy induced by the bending moment M is shown in Equation 5.

$\mathrm{\Delta }{V_{\epsilon }}^M=\frac{M^2}{2EI}ds(5)$

The incremental strain energy induced by the torque T is shown in Equation 6.

$\mathrm{\Delta }{V_{\epsilon }}^T=\frac{T^2}{2G{I}_P}ds(6)$

where E is the modulus of elasticity, A is the cross-sectional area of the continuum, K is the shear shape coefficient, I is the moment of inertia, G is the shear modulus of elasticity, $I_p$ is the polar moment of inertia.

Due to the typical bending state of continua in operational conditions, this section is based on the assumption of constant curvature for analyzing the stiffness of continuum under bending. The continuum is typically subjected to three types of end forces: axial force $F_a$, principal normal force $F_n$, and deputy normal force $F_b$, as shown in Figure 2. These forces exert varying effects on continuum, necessitating distinct analytical approaches.

Figure 2

Figure 2. Force conditions at the end of the continuum joint.

3.2 Continuum robot stiffness model

The continuum can be envisioned as a smooth circular tube. As shown in Figure 3, consider the overall bending angle $\mathrm{\Theta }$ of the continuum robot, with a total length L. To analyze a small segment of the continuum robot near its end, with a bending angle $\theta$ and length s. We note the following forces: axial force $F_a$, normal force $F_N$, shear force $F_s$, and bending moment M acting on this segment of the continuum under stress.

Figure 3

Figure 3. The force on the continuum. (A) Axial force situation of continuum element. (B) Principal normal force situation of element. (C) Tangential normal force situation of element.

As shown in Figure 3A, according to the force balance, the following equation can be obtained:

$\left\{\begin{array}{c}{F}_N=F_{\text{a}}\cos \theta \hfill \\ F_s=F_{\text{a}}\sin \theta \hfill \\ M=F_{a}l\left(1-\cos \theta \right)/\mathrm{\Theta }\hfill \end{array}\right.(7)$

According to Castigliano’s second theorem, the deformation at the end is related to the strain energy. Let the deformation variable at the end be $dx$, then it satisfies

$dx=\frac{\partial \left(\underset{0}{\overset{l}{\int }}\frac{{F_N}^2}{2EA}ds+\underset{0}{\overset{l}{\int }}\frac{k{F_s}^2}{2GA}ds+\underset{0}{\overset{l}{\int }}\frac{M^2}{2EI}ds\right)}{\partial F_a}(8)$

Substituting Equation 7 into Equation 8 yields Equation 9:

$\begin{array}{c}dx=\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{F_{a}l{\cos }^2\theta }{EA\mathrm{\Theta }}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{k{F}_{a}l{\sin }^2\theta }{GA\mathrm{\Theta }}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{F_a{\left(1-\cos \theta \right)}^2}{EI}\frac{l^3}{{\mathrm{\Theta }}^3}d\theta \hfill \\ =\frac{F_{a}l}{EA\mathrm{\Theta }}\frac{2\mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{k{F}_{a}l}{GA\mathrm{\Theta }}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}+\frac{F_a}{EI}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}\hfill \end{array}(9)$

In the above equation, the first to third terms on the right-hand side represent the deformations at the end caused by axial force, shear force, and bending moment, respectively. E and G are constants obtainable through experimentation. Axial stiffness is defined as the ratio of axial force to deformation, thus K_a can be described in Equation 10.

$K_a={\left(\frac{l}{EA\mathrm{\Theta }}\frac{2\mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{kl}{GA\mathrm{\Theta }}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}+\frac{1}{EI}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}\right)}^{-1}(10)$

Furthermore, since the continuum is a slender rod with a large length-to-diameter ratio, shear strain is significantly smaller compared to other strain energies. Therefore, it can be neglected, and $K_a$ can be simplified to Equation 11.

$K_a={\left(\frac{l}{EA\mathrm{\Theta }}\frac{2\mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{1}{EI}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}\right)}^{-1}(11)$

Similarly, as shown in Figure 3B, the following can be derived in Equation 12.

$K_n={\left(\frac{l}{EA\mathrm{\Theta }}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}+\frac{1}{EI}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}\right)}^{-1}(12)$

When the tangential direction is subjected to an external force $F_b$, not only shear force and bending moment occur, but also torsion is generated, as shown in Figure 3C.

$T=F_b\frac{l}{\mathrm{\Theta }}\left(1-\cos \theta \right)(13)$

Similarly, the expression of $dx$ and $K_b$ are as follows.

$\left\{\begin{array}{c}dx=\frac{\partial \left(\underset{0}{\overset{l}{\int }}\frac{T^2}{2G{I}_P}ds+\underset{0}{\overset{l}{\int }}\frac{k{F_s}^2}{2GA}ds+\underset{0}{\overset{l}{\int }}\frac{M^2}{2EI}ds\right)}{\partial F_b}\hfill \\ K_b={\left(\frac{1}{G{I}_P}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{2}{EI}\frac{l^3}{{\mathrm{\Theta }}^3}\left(\mathrm{\Theta }-\sin \mathrm{\Theta }\right)\right)}^{-1}\hfill \end{array}\right.(14)$

The stiffness of the continuum in three different directions is compared under different lengths l and outer diameters R, as shown in Figures 4A–F. The distribution of stiffness $K_a$, $K_n$ and $K_b$ across different orientations shows that within the bending range of (0, 1.5) radians, the stiffness in the principal normal direction significantly exceeds that in the axial and tangential directions. However, at extremely small bending angles, the tangential direction stiffness can surpass the axial direction stiffness.

Figure 4

Figure 4. Comparison of stiffness in different directions. (A) $R=5mm$, $l=40mm$. (B) $R=5mm$, $l=60mm$. (C) $R=5mm$, $l=80mm$. (D) $R=3mm$, $l=60mm$. (E) $R=3mm$, $l=60mm$. (F) $R=8mm$, $l=60mm$.

Analyzing the more generalized spatial stiffness model, as shown in Figure 2, where the continuum subjected to a load F at its end experiences a total deformation $dx$. The force F can be decomposed into components $F_a$, $F_n$ and $F_b$. Given that the plane containing F makes an angle $\beta$ with the plane defined by directions a and b, and within that plane, F forms an angle $\phi$ with direction b. Therefore, the magnitudes of the individual force components are shown in Equation 15.

$\left\{\begin{array}{c}{F}_n=F\sin \varphi \sin \beta \hfill \\ F_{\text{a}}=F\sin \varphi \cos \beta \hfill \\ F_b=F\cos \varphi \hfill \end{array}\right.(15)$

Similarly, the deformation of F at the end can be represented as Equation 16.

$\overrightarrow{dx}=\overrightarrow{d{x}^{\text{a}}}+\overrightarrow{d{x}^{\text{n}}}+\overrightarrow{d{x}^{\text{b}}}(16)$

If the loads F act individually in three directions, resulting in deformations $d{x}^{F_a}$, $d{x}^{F_n}$ and $d{x}^{F_b}$ respectively at the ends, the relationship can be expressed as

$dx=d{x}^{Fa}{\cos }^2{\phi }^a+d{x}^{Fn}{\cos }^2{\phi }^n+d{x}^{Fb}{\cos }^2{\phi }^b(17)$

The angles ${\phi }_a$, ${\phi }_n$ and ${\phi }_b$ represent the angles between the spatial force F and the directions of the three component forces. Further derivation yields:

$\frac{1}{K}=\frac{dx}{F}=\frac{{\cos }^2{\phi }^a}{K_a}+\frac{{\cos }^2{\phi }^n}{K_n}+\frac{{\cos }^2{\phi }^b}{K_b}(18)$

3.3 Single-section notched continuum stiffness model

The skeleton diagram of the notched continuum is shown in Figure 5. The constant curvature assumption is adopted to describe the bending of the continuum for computational convenience. The bending angle ${\theta }_i$ at section i of the continuum is assumed to be uniform. The total bending angle $\theta$ of the continuum satisfies the Equation 19.

$\mathrm{\Theta }={\displaystyle \sum _0^{\text{i}}}{\theta }_{\text{i}}(19)$

Figure 5

Figure 5. The bending condition of three consecutive joints of the notched continuum robot.

The cross-section of the continuum consists of two symmetric annular structures. The cross-section is simplified into two symmetrically distributed rectangles. Each rectangle has the length a, width b, and the distance from the rectangle’s center to the centroid c. Therefore, the moments of inertia of the continuum skeleton section about the x and y-axes, as well as the polar moment of inertia, are given by Equation 20.

$I_x=\frac{a{b}^3}{6},I_y=\frac{b{a}^3}{6}+2abc,I_{p1}=\frac{a{b}^3}{6}+\frac{b{a}^3}{6}+2abc(20)$

The continuum drive cable studied in this paper is made of nickel-titanium alloy and can be approximated to obey Hooke’s law. Its loading conditions are treated as equivalent to a slender rod. When the drive cable is fully tensioned, friction between the cable and the continuum skeleton is neglected. For computational convenience, this study considers four cables as a unified whole for analysis, with each cable having a cross-section of a uniformly distributed circle.

Therefore, the moments of inertia of the drive wires around the x/y-axis and the equivalent polar moment of inertia are shown in Equation 21.

$I_2=\frac{\pi d^4}{16}+\frac{\pi d^2{e}^2}{2},I_{p2}=\frac{\pi d^4}{8}+\pi d^2{e}^2(21)$

To calculate the axial stiffness, strain energy possessed by the continuum should be determined. The upper and lower annular discs of a joint can be regarded as rigid body, hence they possess no strain energy. The strain energy of the joint is concentrated in the elastic rod, as shown in Figure 6, which shows a schematic of joint i near the end of the continuum. $F_a$ represents the axial load.

Figure 6

Figure 6. Three consecutive joint bending situations. (A) Bending condition of joint i. (B) Bending condition of joint i-1. (C) Bending condition of joint i-2.

Similar to axial stiffness, the force equilibrium within the internal elemental units of the continuum can be inferred from the previous section. The force distribution at the base of each unit resembles that of the principal normal stiffness. Hence, the strain energy within this joint is shown in Equation 22.

${V_{\epsilon }}^{\text{i}}=\underset{0}{\overset{l}{\int }}\frac{{\left({F_{dN}}^{\text{i}}\right)}^2}{2E_{1}A}ds+\underset{0}{\overset{l}{\int }}\frac{k{\left({F_{ds}}^{\text{i}}\right)}^2}{2G_{1}A}ds+\underset{0}{\overset{l}{\int }}\frac{{\left({M_d}^{\text{i}}\right)}^2}{2E_1{I}_{\text{x}}}ds(22)$

Considering the case of the i-1 joint. As shown in Figure 6A, based on its geometric features, the equations of static equilibrium can be formulated in Equation 23.

$\left\{\begin{array}{c}{F_s}^{\text{i\hspace{0.17em}}-1}={F_s}^{\text{i}}\cos \theta +{F_N}^{\text{i}}\sin \theta \hfill \\ {F_N}^{\text{i\hspace{0.17em}}-1}={F_s}^{\text{i}}\sin \theta -{F_N}^{\text{i}}\cos \theta \hfill \\ M^{\text{i}}-M^{\text{i\hspace{0.17em}}-1}={F_s}^{i-1}l\sin \theta /{\theta }_{\text{i\hspace{0.17em}}-1}+{F_N}^{i-1}l\left(1-\cos \theta \right)/{\theta }_{\text{i\hspace{0.17em}}-1}\hfill \end{array}\right.(23)$

The strain energy within the i-1 joint shown in Equation 24.

${V_{\epsilon }}^{\text{i\hspace{0.17em}}-1}=\underset{0}{\overset{l}{\int }}\frac{{\left({F_{dN}}^{\text{i\hspace{0.17em}}-1}\right)}^2}{2E_{1}A}ds+\underset{0}{\overset{l}{\int }}\frac{k{\left({F_{ds}}^{\text{i\hspace{0.17em}}-1}\right)}^2}{2G_{1}A}ds+\underset{0}{\overset{l}{\int }}\frac{{\left({M_d}^{\text{i\hspace{0.17em}}-1}\right)}^2}{2E_1{I}_y}ds(24)$

The strain energy within the i-2 joint is computed similarly to that within the i-1 joint. Following this method, for any joint n (where n = 1, 2 … i), the total strain energy of the continuum skeleton can be expressed, then the total strain energy of the continuum skeleton can be expressed as Equation 25.

$V_{\epsilon }={\displaystyle \sum _{n=1}^i}{V_{\epsilon }}^n(25)$

To express the axial stiffness of the notched continuum further in Equation 26.

$\left\{\begin{array}{c}dx=\frac{\partial \left(V_{\epsilon }\right)}{\partial F_a}\hfill \\ K_a=\frac{F_a}{dx}\hfill \end{array}\right.(26)$

To determine the principal normal stiffness, considering the i-1st joint scenario, the input force at this joint is influenced by the end force at the first joint. As shown in Figure 7, based on geometric characteristics, the static equilibrium equations can be derived.

$\left\{\begin{array}{c}{F_s}^{\text{i\hspace{0.17em}}-1}={F_s}^{\text{i}}\cos \theta -{F_N}^{\text{i}}\sin \theta \hfill \\ {F_N}^{\text{i\hspace{0.17em}}-1}={F_s}^{\text{i}}\sin \theta +{F_N}^{\text{i}}\cos \theta \hfill \\ M^{\text{i}}-M^{\text{i\hspace{0.17em}}-1}={F_s}^{i-1}l\sin \theta /{\theta }_{\text{i\hspace{0.17em}}-1}+{F_N}^{i-1}l\left(1-\cos \theta \right)/{\theta }_{\text{i\hspace{0.17em}}-1}\hfill \end{array}\right.(27)$

Figure 7

Figure 7. Three consecutive joint bending situations. (A) Bending condition of joint i. (B) Bending condition of joint i-1. (C) Bending condition of joint i-2.

Through the above equation, the forces at the bottom of the second joint and subsequently determine the strain energy within the second joint can be solved. The strain energy of the i-2nd joint aligns with Equation 27. Similarly, can be derived (n = 1, 2 … i), thus obtaining the total strain energy of the continuum. This enables the stiffness of the notched continuum in the principal normal direction can be expressed.

To solve for the deputy normal stiffness of the notched continuum, as shown in Figure 7A, when subjected to an external force $F_b$, a torque is generated. The force situation is illustrated in Figure 7C. Similar to the static analysis of the i-2nd joint and Equations 13, 14, the strain energy of the continuum joints can be determined. In the i-1st joint, the condition satisfied is shown in Equation 28.

$\left\{\begin{array}{c}{F_s}^i={F_s}^{i-1}\hfill \\ M^{i-1}=M^i-2{F_s}^i\frac{l}{{\theta }_{i-1}}\sin \frac{\theta }{2}\hfill \\ T^{i-1}=T^i+{F_s}^{i-1}l\left(1-\cos \theta \right)/{\theta }_{i-1}\hfill \end{array}\right.(28)$

By analogously applying the static equilibrium relationships in the i-2nd joint, we can similarly derive (n = 1, two … i) the total strain energy of the continuum. The stiffness of the notched continuum in the transverse direction can be expressed.

Analyzing the axial stiffness of the drive cable, when the axial load $F_a$ is applied to the end of the continuum, as indicated by earlier discussions and considering static equilibrium and deformation coordination relationships, dx and K_a can be shown in Equation 29 and Equation 30.

$\begin{array}{c}dx=\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{F_{a}l{\cos }^2\theta }{E_{2}A\mathrm{\Theta }}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{k{F}_{a}l{\sin }^2\theta }{G_{2}A\mathrm{\Theta }}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{F_a{\left(1-\cos \theta \right)}^2}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}d\theta \hfill \\ =\frac{F_{a}l}{E_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{k{F}_{a}l}{G_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}\hfill \\ +\frac{F_a}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}\hfill \end{array}(29)$

$K_a={\left(\frac{l}{E_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{kl}{G_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}+\frac{1}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}\right)}^{-1}(30)$

When analyzing the primary normal stiffness of the driving cable, and considering that the end of the continuum is subjected to a load $F_n$ in the primary normal direction, as previously discussed. Thus dx and K_n can be described in Equation 31 and Equation 32.

$\begin{array}{c}dx=\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{F_{n}l{\sin }^2\theta }{E_{2}A\mathrm{\Theta }}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{k{F}_{n}l{\cos }^2\theta }{G_{2}A\mathrm{\Theta }}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{F_n{\sin }^2\theta }{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}d\theta \hfill \\ =\frac{F_{n}l}{E_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}+\frac{k{F}_{n}l}{G_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{F_n}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}\hfill \end{array}(31)$

$K_n={\left(\frac{l}{E_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}+\frac{kl}{G_{2}A\mathrm{\Theta }}\frac{2\mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{1}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{2\mathrm{\Theta }-\sin 2\mathrm{\Theta }}{4}\right)}^{-1}(32)$

When analyzing the secondary normal stiffness, upon the continuum’s end being subjected to a load $F_b$ in the deputy normal direction, resulting in reactions and moments at the end of the cable. Thus dx and K_b can be described in Equation 35 and Equation 34.

$\begin{array}{c}dx=\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{F_b{\left(1-\cos \theta \right)}^2}{G_2{I}_P}\frac{l^3}{{\mathrm{\Theta }}^3}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{k{F}_{b}l}{G_{2}A\mathrm{\Theta }}d\theta +\underset{0}{\overset{\mathrm{\Theta }}{\int }}\frac{4F_b{\sin }^2\frac{\theta }{2}}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}d\theta \hfill \\ =\frac{F_b}{G_2{I}_{P2}}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{k{F}_{b}l}{G_{2}A}+\frac{2F_b}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}\left(\mathrm{\Theta }-\sin \mathrm{\Theta }\right)\hfill \end{array}(33)$

$K_b={\left(\frac{1}{G_2{I}_{P2}}\frac{l^3}{{\mathrm{\Theta }}^3}\frac{6\mathrm{\Theta }-8\sin \mathrm{\Theta }+\sin 2\mathrm{\Theta }}{4}+\frac{kl}{G_{2}A}+\frac{2}{E_2{I}_2}\frac{l^3}{{\mathrm{\Theta }}^3}\left(\mathrm{\Theta }-\sin \mathrm{\Theta }\right)\right)}^{-1}(34)$

Comparing the coupled stiffness between the continuum and the cables in different directions, the expressions for their stiffness along the axial have been derived, principal normal, and deputy normal directions in preceding sections. Thus, the combined stiffness in these three orientations can be determined, as outlined in Hong et al. (2022).

$\left\{\begin{array}{c}{K}_a=K_{a-C}+K_{a-L}\hfill \\ K_n=K_{n-C}+K_{n-L}\hfill \\ K_b=K_{b-C}+K_{b-L}\hfill \end{array}\right.(35)$

In this context, subscript C denotes the continuum, while subscript L represents the drive cable. The iterative computations are performed on the directional stiffness of the notched continuum and drive cable coupled model using MATLAB. The parameter definitions and values used for modeling are detailed in Table 1, outside the continuum section.

Taking the total bending angle as a variable, the computed results are depicted in Figure 8. It is evident that the trends align closely with those of an ideal continuum; as the bending angle increases, $K_a$ gradually decreases, $K_n$ exhibits increased fluctuations, and $K_b$ demonstrates the U shaped trend, albeit with minimal magnitude of variation.

Figure 8

Figure 8. The continuum drive cable’s coupled model.

Analyzing the stiffness model of the continuum and the cables, as concluded earlier, the deformation can be expressed as Equation 17. The stiffness model of the continuous system’s end under spatial loads can be represented by Equation 18.

3.4 Dual-section notched continuum robot stiffness model

As shown in Figure 9, the dual-segment configuration studied in this paper consists of nested outer and inner continuum with different diameters and lengths. Both inner and outer continuum are purely sectional. The geometric parameters of the single-segment continuous body are provided in Section 3.2. In this section, to distinguish between the inner and outer continuum, the inner tube is designated as Tube 2 and the outer tube as Tube 1. The following analysis considers the combined stiffness when both tubes are bent in arbitrary configurations. It is assumed that the bending planes of the two segments do not overlap, and the angle between these bending planes, denoted as $\gamma$, is between $Y_{1-bottom}$ and $Y_{2-end}$.

Figure 9

Figure 9. Schematic diagram of establishing the coordinate system of a dual-segment notched continuum robot.

Before proceeding with the subsequent derivations, this paper makes several assumptions to simplify the calculations regarding the model: 1). Frictional effects between the drive cable and the continuum are neglected. 2). The notches at the coupling interfaces of the two segments of the continuum overlap and have identical shapes. 3). The internal forces within both Segment two and Segment one at their junction interface are assumed to be exactly equal. 4). It is assumed that the bending directions of the two segments of the continuum lie within the same plane, and rotational effects of the continuum are disregarded. When the two segments of the continuum are coupling, the overall axial stiffness is the sum of the overall stiffness of the wire and the overall stiffness of the continuum.

$K=K_C+K_L(36)$

In the setup, continuum 2 comprises a total of i segments of joints, while Continuum 1 comprises a total of 2i segments of joints. Aside from differences in length, diameter, and cross-sectional dimensions, all other parameters between Continuum 1 and Continuum 2 are identical. The overall bending angles at the coupled part of the continuum are represented as ${\mathrm{\Theta }}_1$, and the overall bending angles where the continuum is not involved in the coupling are denoted as ${\mathrm{\Theta }}_2$. $E_{\mathit{continuum}}$ and $E_{\mathit{cable}}$ refer to the elastic moduli of the continuum and the cable respectively, while GC and GL denote the shear moduli of the continuum and the driving cables. $I{X}_{\mathit{C21}}$, $I{X}_{\mathit{C22}}$, $I{X}_{C1}$, $I{Y}_{\mathit{C21}}$, $I{Y}_{\mathit{C22}}$ and $I{Y}_{C1}$ respectively represent the moments of inertia of continuum two in the uncoupled part, continuum two in the coupled part, and continuum 1 in the X and Y directions.

After concatenation, the total energy of the continuum equals the sum of the total strain energies of continuum 1 and continuum 2.

$V_{\epsilon -C}={\displaystyle \sum _{m=1}^n}{V_{\epsilon }}^m+{\displaystyle \sum _{m=1}^i}{V_{\epsilon }}^m(37)$

The total strain energy in the continuum section:

$\begin{array}{c}\hfill V_{\epsilon -C}={\displaystyle \sum _1^i}{V}_{\epsilon -C2}^n+{\displaystyle \sum _1^i}{V}_{\epsilon -C1}^n\hfill \\ \hfill ={\displaystyle \sum _{n=j}^i}\left(\underset{0}{\overset{l}{\int }}\frac{{\left({F_N}^n\right)}^2}{2E_{C}A2}ds+\underset{0}{\overset{l}{\int }}\frac{k{\left({F_s}^n\right)}^2}{2G_{C}A2}ds+\underset{0}{\overset{l}{\int }}\frac{{\left(M^n\right)}^2}{2E_C{I}_{\mathit{C21}}}ds\right)\hfill \\ \hfill +{\displaystyle \sum _{n=1}^j}\left(\underset{0}{\overset{l}{\int }}\frac{{\left({F_N}^n\right)}^2}{2E_{C}A2}ds+\underset{0}{\overset{l}{\int }}\frac{k{\left({F_s}^n\right)}^2}{2G_{C}A2}ds+\underset{0}{\overset{l}{\int }}\frac{{\left(M^n\right)}^2}{2E_C{I}_{\mathit{C22}}}ds\right)\hfill \\ \hfill +{\displaystyle \sum _{n=1}^m}\left(\underset{0}{\overset{l}{\int }}\frac{{\left({F_N}^n\right)}^2}{2E_{C}A1}ds+\underset{0}{\overset{l}{\int }}\frac{k{\left({F_s}^n\right)}^2}{2G_{C}A1}ds+\underset{0}{\overset{l}{\int }}\frac{{\left(M^n\right)}^2}{2E_C{I}_{C1}}ds\right)\hfill \end{array}(38)$

where $F_{Nn}$, $F_{Sn}$, and $M_n$ are the stress conditions at the bottom of the nth joint of the continuum, with stress values derived recursively from the preceding section. I is the moment of inertia, and the expressions for different joints are distinct.

The strain energy of the cable is the sum of the total strain energy of the two segments of cable2 and the strain energy of cable1, namely,:

$\begin{array}{c}\hfill V_{\epsilon -L}={\displaystyle \sum _1^2}{V}_{\epsilon -L2}^n+V_{\epsilon -L1}\hfill \\ \hfill =\left(\underset{0}{\overset{l1}{\int }}\frac{{\left({F_N}^n\right)}^2}{2E_{L}A}ds+\underset{0}{\overset{l1}{\int }}\frac{{\left(k{F_s}^n\right)}^2}{2G_{L}A}ds+\underset{0}{\overset{l1}{\int }}\frac{{\left(M^n\right)}^2}{2E_L{I}_{\mathit{L21}}}ds\right)\hfill \\ \hfill +\left(\underset{l1}{\overset{l2}{\int }}\frac{{\left({F_N}^n\right)}^2}{2E_{L}A}ds+\underset{l1}{\overset{l2}{\int }}\frac{{\left(k{F_s}^n\right)}^2}{2G_{L}A}ds+\underset{l1}{\overset{l2}{\int }}\frac{{\left(M^n\right)}^2}{2E_L{I}_{\mathit{L22}}}ds\right)\hfill \\ \hfill +\left(\underset{0}{\overset{l2}{\int }}\frac{{\left({F_N}^n\right)}^2}{2E_{L}A}ds+\underset{0}{\overset{l2}{\int }}\frac{{\left(k{F_s}^n\right)}^2}{2G_{L}A}ds+\underset{0}{\overset{l2}{\int }}\frac{{\left(M^n\right)}^2}{2E_L{I}_{L1}}ds\right)\hfill \end{array}(39)$

Given the known bending and deflection angles of the two-segment continuum, the stiffness of the end under axial force $F_a$ can be computed using the aforementioned method. Similarly, the total strain energy in the principal and secondary directions of the two-segment notched continuum can be calculated. Drawing an analogy to Equations 35–39, the stiffness models for the principal and secondary directions can be derived.

The stiffness model of the drive cable is established in Equation 40. When the continuum’s end is subjected to axial load $F_a$, the stiffness of the drive screw can be calculated in two parts: uncoupled and coupled.

$K_{a-L}=K_{a-LU}+K_{a-LUU}(40)$

where U is the coupled part, UU is the uncoupled part. The input force for the uncoupled part of the drive cable is simply $F_a$, whereas for the coupled part, it should be the reaction force exerted on the lower end of the uncoupled part. Therefore, this study only needs to refer to Section 3.4 to separately calculate the stiffness of each part, and then linearly combine them. Similarly, when the continuum’s end is subjected to loads in the principal normal direction $F_n$ and secondary normal direction $F_b$, the stiffness is defined in the same manner.

Compare the coupled stiffness in various directions between the continuum and the drive cable. Referring to Section 3.3, iterative calculations of the coupled stiffness models in each direction for the notched continuum and drive cable are conducted in MATLAB, with parameters provided in Section 2.

Using the total bending angle ${\mathrm{\Theta }}_1$ and of the two-segment continuum as variables, the three-dimensional graph of stiffness distribution is plot, as shown in Figure 10.

Figure 10

Figure 10. Variation of stiffness with angle. (A) Axial stiffness distribution of the continuum. (B) Principal normal stiffness distribution of the continuum. (C) Secondary normal stiffness distribution of the continuum.

It can be observed that the axial stiffness $F_a$ of the continuum is greatly influenced by ${\mathrm{\Theta }}_1$ and ${\mathrm{\Theta }}_2$ from Figure 10A. When both angles approach zero, the axial stiffness reaches its maximum value, and the stiffness value sharply decreases as the angles in-crease. Furthermore, it is noted that for smaller values of $\theta$, the stiffness curve shows significant variations, closely resembling the axial stiffness curve of a single-segment continuum. In contrast, for larger values of $\theta$, the magnitude of stiffness variation diminishes, and the trend of change also alters., It can be seen that the principal normal stiffness $K_n$ is significantly influenced by ${\mathrm{\Theta }}_2$ from Figure 10B. Specifically, as ${\mathrm{\Theta }}_1$ increases, $K_n$ also increases, while the influence of ${\mathrm{\Theta }}_1$ on is relatively small. When both ${\mathrm{\Theta }}_1$ and ${\mathrm{\Theta }}_2$ approach zero, $K_n$ reaches its minimum value, which closely resembles the trend of the principal normal stiffness of a single-segment continuum. It can be observed that the secondary normal stiffness $K_b$ is similarly greatly affected by ${\mathrm{\Theta }}_2$ from Figure 10C, As ${\mathrm{\Theta }}_2$ increases, $K_b$ initially increases rapidly, then its rate of in-crease slows down, and finally it increases rapidly again. Conversely, the impact on $K_n$ is relatively minor, and as ${\mathrm{\Theta }}_1$ increases, $K_b$ initially increases rapidly followed by a slower rate of increase. Similarly, $K_b$ reaches its minimum value as ${\mathrm{\Theta }}_1$ and ${\mathrm{\Theta }}_2$ approach 0. This trend closely aligns with the behavior of secondary normal stiffness observed in a single continuum segment.

Establish the spatial stiffness model for the two-segment continuum. From Section 3.3, let ${\phi }_a$, ${\phi }_n$, and ${\phi }_b$ denote the angles between the spatial load F and the three force directions, respectively. The stiffness model of the two-segment continuum under the action of spatial loads can be expressed as Equation 18. To compare the stiffness under four different driving modes, the coupled stiffness under different conditions by varying the bending angles and angles of forces on the continuum are simulated and calculated in MATLAB. In the calculations, there are a total of five angle variables, denoted as ${\mathrm{\Theta }}_1$, ${\mathrm{\Theta }}_2$, ${\phi }_a$, ${\phi }_n$, and ${\phi }_b$. Since this study only considers the case where the bending angles of the two continua are in the same plane, ${\phi }_b=90^{◦}$. The simulation results are shown in Figure 11.

Figure 11

Figure 11. The stiffness in different angles. (A) Spatial stiffness distribution under different ${\mathrm{\Theta }}_1$ conditions. (B) Spatial stiffness distribution under various ${\mathrm{\Theta }}_2$ conditions.

Comparison of spatial stiffness of the two-segment continuum under different driving strategies. The drive cable quantity is varied, thus changing the inertia moment and total cross-sectional area of the drive cables. Several sets of ${\mathrm{\Theta }}_1$ and ${\mathrm{\Theta }}_2$ are selected to compare the spatial stiffness of the continuum. Partial simulation results are shown in Figure 12, where Figures 12A, B vary with ${\mathrm{\Theta }}_1$, Figures 12C, D vary with ${\mathrm{\Theta }}_2$. From the figures, it can be observed that in practice, the stiffness of the continuum is nearly the same for the External 4+Internal 4 and External 2+Internal 4 strategies, and similarly for the External 4+Internal 2 and External 2+Internal 2 strategies. Overall, the relationship in magnitude is External 4+Internal 4 $>$ External 2+Internal 4 $>$ External 4+Internal 2 $>$ External 2+Internal 2.

Figure 12

Figure 12. The spatial stiffness distribution for different numbers of driving strategies. (A) The spatial stiffness distribution under different ${\mathrm{\Theta }}_1$ conditions. (B) The spatial stiffness distribution under various ${\mathrm{\Theta }}_1$ conditions. (C) The spatial stiffness distribution under different ${\mathrm{\Theta }}_2$ conditions. (D) The spatial stiffness distribution under various ${\mathrm{\Theta }}_2$ conditions.

4 Experiment and result analysis

In this section, the continuum experimental platform is built and the basic performance is tested. The stiffness model verification experiment based on the experimental platform is carried out. Finally, the stiffness of the series continuum under different driving strategies is compared experimentally.

4.1 Experimental platform testing

The serial continuum experimental platform constructed in this study is depicted as shown in Figure 13. Two segments of continuum robots are respectively mounted at the front end. Additionally, the bending angle of the continuum robot is measured using a single electromagnetic tracking sensor installed at the end-effector of the robot. The sensor provides real-time position and orientation data of the end-effector. The bending angle $\theta$ is calculated based on the direction vector $\mathbf{D}$ obtained from the sensor and the initial direction vector ${\mathbf{D}}_0$. The calculation formula is given by Equation 41.

$\theta =\arccos \left(\frac{\mathbf{D}\cdot {\mathbf{D}}_0}{\Vert \mathbf{D}\Vert \Vert {\mathbf{D}}_0\Vert }\right)(41)$

where $\mathbf{D}$ is the direction vector measured by the sensor, and ${\mathbf{D}}_0$ is the initial direction vector of the robot.

Figure 13

Figure 13. The prototype system of continuum robot experimental platform.

Motion tests on the inner and outer continuum bodies at various angles are conducted. Due to the material characteristics of the continuum robot, the bending angle between the two segments to be less than $60^{◦}$ during the experiments are controlled. This setup is sufficient for small deformation angles. For larger deformations, the use of an array of sensors may be considered in future work to enhance measurement accuracy.

4.2 Stiffness model verification

To validate the spatial stiffness model proposed in Section 3.4 a series of experiments are designed in this paper. The continuum robot is bent to a certain extent, and electromagnetic tracking sensors are used to detect its initial position. The weight of 20 g is attached to the end of the continuum robot to provide a constant force, and the change in end position is monitored simultaneously. This allowed for the calculation of the overall stiffness of the continuum robot. Initially, ${\mathrm{\Theta }}_2$ is kept constant while gradually increasing ${\mathrm{\Theta }}_1$. Displacement data for two sets of continuum robots are measured to obtain stiffness curves, as shown in Table 2 and Figure 14. Subsequently, ${\mathrm{\Theta }}_1$ is kept constant while gradually increasing ${\mathrm{\Theta }}_2$ and displacement data are collected to derive stiffness curves, as shown in Figure 15, with corresponding data shown in Table 3 and Figure 15. To ensure clarity and readability, the corresponding figures display the average values with error bars representing the standard deviations.

Table 2

Table 2. Experimental data of stiffness with invariable ${\mathrm{\Theta }}_1$.

Figure 14

Figure 14. Comparison of theoretical and actual stiffness of the continuum under constant ${\mathrm{\Theta }}_1$. (A) ${\mathrm{\Theta }}_2=20^{◦}$. (B) ${\mathrm{\Theta }}_2=40^{◦}$.

Figure 15

Figure 15. Comparison of theoretical and actual stiffness of the continuum under constant ${\mathrm{\Theta }}_2$. (A) ${\mathrm{\Theta }}_1=15^{◦}$. (B) ${\mathrm{\Theta }}_1=35^{◦}$.

Table 3

Table 3. Experimental data of stiffness with invariable ${\mathrm{\Theta }}_2$.

From Figures 14, 15, it can be seen that under small deformation angles, the actual stiffness of the continuum follows a trend closely aligned with the theoretical stiffness curve. This confirms that the spatial stiffness model proposed in Section 3.4 largely conforms to the practical model. However, there are still errors. Preliminary analysis shows that the sources of errors are as follows: (1) The influence of friction is not considered in the theoretical model. (2) The constant curvature model is used for continuum modeling, which has errors. (3) There are systematic errors and accidental errors in the experiment. (4) There are errors in the selection of actual material characteristic parameters, such as Young’s modulus, shear modulus and other parameters.

4.3 Continuum stiffness comparison

The verification work on the stiffness relationships of serial continuum under different driving strategies proposed in Section 3.4 is conducted. The experimental setup is identical to the previous subsection. Due to the conclusion that under the External 4+Internal four and External 2+Internal 4 strategies, the stiffness of continuum is nearly identical, as is the case under the External 4+Internal two and External 2+Internal 2 strategies, with an overall hierarchy of External 4+Internal 4 $>$ External 2+Internal 4 $>$ External 4+Internal 2 $>$ External 2+Internal 2, serial continuum satisfy this regardless of angle. Thus, experiments are conducted with randomly assigned angle values for two-stage continuum, only changing the number of drive cables. The experimental results are shown in Table 4. Table 4 presents the experimental data for stiffness comparison under different driving strategies. The input variables are the driving strategies (e.g., External 4+Internal 4, External 2+Internal 4), while the output variables are the measured displacements under these strategies.

Table 4

Table 4. Stiffness experimental data with fixed ${\mathrm{\Theta }}_1$.

It can be observed that the stiffness order under the four driving strategies is generally consistent with the theoretical results, but there are numerical errors that may cause the conclusion that under the External 4+Internal 4 and External 2+Internal 4 strategies, the stiffness of continuum is nearly identical, as is the case under the External 4+Internal 2 and External 2+Internal 2 strategies to be less obvious.

5 Conclusion

Basic performance evaluation of continuum robots is beneficial to their operation optimization and precise control in neurosurgery. The basic performance of continua is mainly evaluated by stiffness, but there is no systematic and universal evaluation system. In order to realize the driving and testing of continuum with different configurations and different driving strategies, this paper designs a universal experimental platform for continuum robot, and the continuum stiffness evaluation method based on the energy method and Castigliano’s second theorem is proposed. By solving the internal force and energy of the cut continuum in sections, the stiffness model of single-segment and double-segment series continuum is established. The relationship between the individual stiffness and bending angle of the single-segment continuum and the driving cables in the axial direction, main normal direction, and secondary normal direction is obtained. The simulation and experimental results show that under the condition of small deformation angle, the spatial stiffness model obtained by strain energy basically conforms to the actual model, which verifies the correctness and rationality of the stiffness calculation method proposed in this paper. Future work will explore the incorporation of nonlinear material properties and geometric nonlinearities to extend the model’s applicability to larger deformations.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

MY: Conceptualization, Data curation, Methodology, Software, Validation, Writing–original draft. ZA: Formal Analysis, Software, Validation, Writing–original draft. ZL: Investigation, Data curation, Writing–original draft. YW: Writing–review and editing, Formal Analysis, Investigation. TP: Project administration, Writing–review and editing, Visualization. FD: Funding acquisition, Project administration, Resources, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Key Research and Development Program of China (2023YFB4705800), the National Natural Science Foundation of China (52375020), the National Key Research and Development Program of China (2022YFB4703000), the Key R&D Program of Shandong Province, China (2022CXGC010503), Shandong Provincial Postdoctoral Innovative Talents Funded Scheme (238226), the Fundamental Research Funds for the Central Universities and the Young Scholars Program of Shandong University.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Awtar, S., and Sen, S. (2010). A generalized constraint model for two-dimensional beam flexures: nonlinear load-displacement formulation. J. Mech. Des. N. Y. 132. doi:10.1115/1.4002005

CrossRef Full Text | Google Scholar

Badescu, M., and Mavroidis, C. (2004). New performance indices and workspace analysis of reconfigurable hyper-redundant robotic arms. Int. J. Robotics Res. 23, 643–659. doi:10.1177/0278364904044406

CrossRef Full Text | Google Scholar

Burgner, J., Swaney, P. J., Lathrop, R. A., Weaver, K. D., and Webster, R. J. (2013). Debulking from within: a robotic steerable cannula for intracerebral hemorrhage evacuation. IEEE Trans. Biomed. Eng. 60, 2567–2575. doi:10.1109/tbme.2013.2260860

PubMed Abstract | CrossRef Full Text | Google Scholar

Camarillo, D. B., Milne, C. F., Carlson, C. R., Zinn, M. R., and Salisbury, J. K. (2008). Mechanics modeling of tendon-driven continuum manipulators. IEEE Trans. robotics 24, 1262–1273. doi:10.1109/tro.2008.2002311

CrossRef Full Text | Google Scholar

Dai, J. S. (2012). Finite displacement screw operators with embedded chasles’ motion. J. Mech. Robot. 4. doi:10.1115/1.4006951

CrossRef Full Text | Google Scholar

Eastwood, K. W., Azimian, H., Carrillo, B., Looi, T., Naguib, H. E., and Drake, J. M. (2016). “Kinetostatic design of asymmetric notch joints for surgical robots,” in 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Daejeon, Korea (South), 09-14 October 2016 (IEEE), 2381–2387.

CrossRef Full Text | Google Scholar

Fraś, J., Czarnowski, J., Maciaś, M., and Główka, J. (2014). “Static modeling of multisection soft continuum manipulator for stiff-flop project,” in Recent advances in automation, robotics and measuring techniques (Springer), 365–375.

Google Scholar

Gao, A., Murphy, R. J., Liu, H., Iordachita, I. I., and Armand, M. (2016a). Mechanical model of dexterous continuum manipulators with compliant joints and tendon/external force interactions. IEEE/ASME Trans. Mechatronics 22, 465–475. doi:10.1109/tmech.2016.2612833

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, G., Wang, H., Xia, Q., Song, M., and Ren, H. (2016b). Study on the load capacity of a single-section continuum manipulator. Mech. Mach. Theory 104, 313–326. doi:10.1016/j.mechmachtheory.2016.06.010

CrossRef Full Text | Google Scholar

Gravagne, I. A., Rahn, C. D., and Walker, I. D. (2003). Large deflection dynamics and control for planar continuum robots. IEEE/ASME Trans. mechatronics 8, 299–307. doi:10.1109/tmech.2003.812829

CrossRef Full Text | Google Scholar

Gu, X., and Ren, H. (2023). A survey of transoral robotic mechanisms: distal dexterity, variable stiffness, and triangulation. Cyborg Bionic Syst. 4, 0007. doi:10.34133/cbsystems.0007

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B., Wang, Z., Li, Q., Xie, H., and Shen, R. (2013). An analytic method for the kinematics and dynamics of a multiple-backbone continuum robot. Int. J. Adv. Robotic Syst. 10, 84. doi:10.5772/54051

CrossRef Full Text | Google Scholar

He, B., Xu, S., and Wang, Z. (2018). Research on stiffness of multibackbone continuum robot based on screw theory and euler-Bernoulli beam. Math. Problems Eng. 2018, 1–16. doi:10.1155/2018/6910468

CrossRef Full Text | Google Scholar

Hong, W., Feng, F., Xie, L., and Yang, G.-Z. (2022). A two-segment continuum robot with piecewise stiffness for maxillary sinus surgery and its decoupling method. IEEE/ASME Trans. Mechatronics 27, 4440–4450. doi:10.1109/tmech.2022.3157041

CrossRef Full Text | Google Scholar

Li, Z., and Du, R. (2013). Design and analysis of a bio-inspired wire-driven multi-section flexible robot. Int. J. Adv. Robotic Syst. 10, 209. doi:10.5772/56025

CrossRef Full Text | Google Scholar

Li, Z., Ren, H., Chiu, P. W. Y., Du, R., and Yu, H. (2016). A novel constrained wire-driven flexible mechanism and its kinematic analysis. Mech. Mach. Theory 95, 59–75. doi:10.1016/j.mechmachtheory.2015.08.019

CrossRef Full Text | Google Scholar

Lin, B., Song, S., and Wang, J. (2024). Variable stiffness methods of flexible robots for minimally invasive surgery: a review. Biomim. Intell. Robotics 4, 100168. doi:10.1016/j.birob.2024.100168

CrossRef Full Text | Google Scholar

Pei, X., Yu, J., Zong, G., Bi, S., and Su, H. (2009a). The modeling of cartwheel flexural hinges. Mech. Mach. Theory 44, 1900–1909. doi:10.1016/j.mechmachtheory.2009.04.006

CrossRef Full Text | Google Scholar

Pei, X., Yu, J., Zong, G., Bi, S., and Su, H. (2009b). The modeling of cartwheel flexural hinges. Mech. Mach. Theory 44, 1900–1909. doi:10.1016/j.mechmachtheory.2009.04.006

CrossRef Full Text | Google Scholar

Qi, P., Qiu, C., Liu, H., Dai, J. S., Seneviratne, L. D., and Althoefer, K. (2015). A novel continuum manipulator design using serially connected double-layer planar springs. IEEE/ASME Trans. Mechatronics 21, 1281–1292. doi:10.1109/tmech.2015.2498738

CrossRef Full Text | Google Scholar

Sadati, S. H., Naghibi, S. E., Walker, I. D., Althoefer, K., and Nanayakkara, T. (2017). Control space reduction and real-time accurate modeling of continuum manipulators using ritz and ritz–galerkin methods. IEEE Robotics Automation Lett. 3, 328–335. doi:10.1109/lra.2017.2743100

CrossRef Full Text | Google Scholar

Simaan, N., Xu, K., Wei, W., Kapoor, A., Kazanzides, P., Taylor, R., et al. (2009). Design and integration of a telerobotic system for minimally invasive surgery of the throat. Int. J. robotics Res. 28, 1134–1153. doi:10.1177/0278364908104278

PubMed Abstract | CrossRef Full Text | Google Scholar

Straughan, B. (2013). The energy method, stability, and nonlinear convection. Springer Sci. and Bus. Media 91.

Google Scholar

Trivedi, D., Lotfi, A., and Rahn, C. D. (2007). “Geometrically exact dynamic models for soft robotic manipulators,” in 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems, San Diego, CA, USA, 29 October 2007 - 02 November 2007 (IEEE), 1497–1502.

CrossRef Full Text | Google Scholar

Webster, R. J., Romano, J. M., and Cowan, N. J. (2008). Mechanics of precurved-tube continuum robots. IEEE Trans. robotics 25, 67–78. doi:10.1109/tro.2008.2006868

CrossRef Full Text | Google Scholar

Webster III, R. J., and Jones, B. A. (2010). Design and kinematic modeling of constant curvature continuum robots: a review. Int. J. Robotics Res. 29, 1661–1683. doi:10.1177/0278364910368147

CrossRef Full Text | Google Scholar

Wilkening, P., Alambeigi, F., Murphy, R. J., Taylor, R. H., and Armand, M. (2017). Development and experimental evaluation of concurrent control of a robotic arm and continuum manipulator for osteolytic lesion treatment. IEEE robotics automation Lett. 2, 1625–1631. doi:10.1109/lra.2017.2678543

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, K., and Simaan, N. (2006). “Actuation compensation for flexible surgical snake-like robots with redundant remote actuation,” in Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006, Orlando, FL, 15-19 May 2006 (IEEE), 4148–4154.

Google Scholar

Xu, K., and Simaan, N. (2008). An investigation of the intrinsic force sensing capabilities of continuum robots. IEEE Trans. Robotics 24, 576–587. doi:10.1109/tro.2008.924266

CrossRef Full Text | Google Scholar

Zhang, G., Du, F., Xue, S., Cheng, H., Zhang, X., Song, R., et al. (2022). Design and modeling of a bio-inspired compound continuum robot for minimally invasive surgery. Machines 10, 468. doi:10.3390/machines10060468

CrossRef Full Text | Google Scholar

Zhang, G., Du, F., Zhang, X., Wu, K., Zheng, G., Li, Y., et al. (2024a). Continuum robots: a real-time model-based data-driven nonlinear controller. IEEE Trans. Industrial Electron. 71, 16176–16186. doi:10.1109/tie.2024.3398678

CrossRef Full Text | Google Scholar

Zhang, G., Su, J., Du, F., Zhang, X., Li, Y., and Song, R. (2024b). Composite continuum robots: accurate modeling and model reduction. Int. J. Mech. Sci. 276, 109342. doi:10.1016/j.ijmecsci.2024.109342

CrossRef Full Text | Google Scholar