Driving torque model of the bionic soft arm’s hyperelastic bellows

ABSTRACT

The technology of soft continuum robots represents an advance in the field of robotics to benefit a wide range of industries such as healthcare, manufacturing or environmental exploration. Soft continuum robots can be driven pneumatically by bellows as the soft continuum manipulator presented in this article. A driving torque model of the bellows considering hyperelastic material properties, friction and restoring torques is derived, whose purpose is suitability for model-based control design and subsequent trajectory generation. Hence, the driving torque model must be real-time capable. This is realized by an iterative algorithm calculating the bellows’ torque transmission by assuming a two-dimensional no-slip membrane contact. Nonlinear strain behaviours and hysteresis effects of the bellows are considered by the Ogden material model. The driving torque model’s performance is validated experimentally by measuring the external torques of the bellows.

KEYWORDS:

• Soft robotics
• hyperelastic bellows
• dynamic actuator model

1. Introduction

Soft continuum robotics is a rapidly developing field that has attracted increasing attention in recent years due to the numerous benefits. Unlike traditional robots that are composed of rigid materials, soft continuum robots are flexible and deformable. These flexible and deformable structures enable soft continuum robots to adapt to complex environments and to operate safely in close proximity to humans or delicate objects. Soft continuum robots can bend, twist, and deform in ways that are difficult to achieve with traditional robots, and they offer novel opportunities for sensing and actuation [1]. These properties enable soft robotics to further develop industries such as healthcare [2], manufacturing [3] or environmental exploration [4].

Soft bellows are the key component of the bionic soft arm, a pneumatically driven soft continuum manipulator (Figure 1). The bionic soft arm is actuated pneumatically by four bellows segments (Figure 2) and three rotary drives. As the modelling of the rotary drives relies on standard engineering mechanics, modelling the bellows is more complicated due to its nonlinear material characteristics. Therefore, the design and characterization of the soft bellows is crucial to the performance of the bellows’ driving torque model (henceforth, driving torque model).

Figure 1. Bionic soft arm, source: Festo SE & Co. KG. The $1.05$m long continuum manipulator is actuated pneumatically by four bellows segments and three rotary drives.

Figure 2. Detailed view of a bellows segment, source: Festo SE & Co. KG. Three aluminium ribs are shown, which divide the bellows chain into four interconnected bellows.

Soft bellows can be made of various materials [5] and their behaviour can be modelled using nonlinear hyperelastic material models, such as the Mooney-Rivlin model [6] or the Ogden model [7]. However, Della Santina et al. [8] points out that most papers neglect the modelling of the actuators since using the kinematic and dynamic models of continuum manipulators for model-based control are already difficult. These material models are used for the materials, rubber [9,10] and knitted fabric [11], the bellows consists of. Feng et al. [12] and Liu et al. [13] studies the contact problem of a membrane with a point or a plane using a frictionless contact assumption based on the Mooney-Rivlin material model. These three-dimensional approaches require expensive finite element simulations. Hence, Srivastava et al. [14] investigates an infinitely long rectangular elastic membrane, which simplifies the membrane contact problem from a three-dimensional to a two-dimensional problem. Müller et al. [15] adapts the approach of Srivastava et al. [14] to the setup presented in this paper (Figures 1,2).

In this article, the approach of Müller [15] et al. is extended by in depth analysis of the actuator’s properties and their effects on the performance of the driving torque model leading to the following modifications. The friction inside the bellows and in the bearings are described by pressure-dependent parameters additionally. This also holds for the occurring restoring torques. Hence, the fitted parameters are improved by having an interpretable physical meaning. Furthermore, the Ogden material model is considered instead of the Mooney-Rivlin material model. So, a larger strain range can be represented, which matches the behaviour of the nonlinear hyperelastic bellows better. The driving torque model’s parameters are identified by various measurements and the resulting driving torque model is validated considering hysteresis and under high dynamics. An iterative algorithm is used to approximate the no-slip membrane contact between the bellows and the ribs. Thus, the computation of the bellows’ torque transmission is real-time capable and suitable for model-based control design and subsequent trajectory generation. Therefore, the driving torque model is compared to the driving torque model from [16], which is used in closed-loop control.

This article is structured as follows. First, the setup of a bellows segment is introduced in Section 2. Then, the bellows segment is modelled in Section 3. The parameter identification is described in Section 4 and in Section 5, the driving torque model performance is validated and discussed. Finally, this article is concluded in Section 6.

This section describes the construction of the bionic soft arm (Figure 1) and its bellows segments (Figure 2). The bionic soft arm is actuated pneumatically by four bellows segments and three rotary drives (Figure 3). The rotary drives enable the bionic soft arm to rotate about their axes. In contrast, the bellows segments allow the bionic soft arm to bend along their axes. Hence, the tool-centre-point positions of the bionic soft arm, which are located at the tip of it, can be achieved by a large redundant configuration space.

Figure 3. Mechanics of the bionic soft arm, where ${\delta }_1$, ${\delta }_2$, ${\delta }_3$ are the rotation angles of the rotary drives and ${\mathrm{\Psi }}_1$, ${\mathrm{\Psi }}_2$, ${\mathrm{\Psi }}_3$, ${\mathrm{\Psi }}_4$ are the bending angles of the bellows segments. The gripper is located at the tool-centre-point of the bionic soft arm.

Each bellows segment consists of two opposing bellows chains, which are attached at the beginning and end of a segment to the frame of the bionic soft arm. In between, the bellows chains are held in shape by three aluminium ribs for each bellows chain. These divide a bellows chain into four interconnected bellows.

The construction of a bellows segment is shown in Figure 4. The $\mathit{j}$^th bellows chain within the $\mathit{i}$^th bellows segment is actuated by the mass flow ${\dot{m}}_{\mathrm{B},\mathit{i},\mathit{j}}$ causing the bellows to bulge outward. The resulting bellows pressure and volume are denoted by $p_{\mathrm{B},i,j}$ and $V_{\mathrm{B},i,j}$. The contact area between the frame and bellows at the end of a bellows chain is given as $A_{\mathrm{B},i,j}$. The minimal lever arm between the centre of gravity of a bellows segment and the centre of gravity of a bellows chain is denoted by $r_{\mathrm{B},\min }$. The theory of operation can be explained as follows. If the pressure of the first bellows chain $p_{\mathrm{B},i,1}$ is larger than the pressure of the second bellows chain $p_{\mathrm{B},i,2}$, the contact area of the first bellows chain $A_{\mathrm{B},i,1}$ is also larger than the contact area of the second bellows chain $A_{\mathrm{B},i,2}$ due to the larger elongation of the first bellows chain and vice versa. Furthermore, the lever arm of the first bellows chain is larger than the lever arm of the second bellows chain due to larger contact area of the first bellows chain $A_{\mathrm{B},i,1}$. So, the resulting torque of the first bellows chain is larger than the resulting torque of the second bellows chain causing the bellows segment to bend to the left according to Figure 4. The bending angle of a bellows segment ${\mathrm{\Psi }}_i$ is described by the sum of the joint angles within a bellows segment.

Figure 4. Construction of a bellows segment.

(1) ${\mathrm{\Psi }}_i={\psi }_{\mathit{i},1}+{\psi }_{\mathit{i},2}+{\psi }_{\mathit{i},3}+{\psi }_{\mathit{i},4}.$(1)

The construction of the bellows chain is shown in Figure 5. A bellows has three layers. The innermost layer consists of an airtight rubber. The second layer is knitted fabric and lies directly on top of the rubber. The third layer, which is the outermost layer, is made of stiff knitted fabric allowing the bellows chain to withstand high pressures being up to $7.5$bar absolute pressure. The knitted fabric of the outermost layer has different knitting patterns at different places. For example, the knitting pattern marked in blue (Figure 5, bottom), is more closely meshed and thus exhibits larger stiffness against tensile stresses. The area $A_{\mathrm{B},\min }$ marked in blue (Figure 5, top) describes the minimum area resulting from the design and is attached to the frame at the end of a bellows segment.

Figure 5. Construction of a bellows chain. The area which connects the bellows to the frame is coloured blue in the upper subfigure.

2. Modelling

The driving torque of a bellows segment to be modelled results from several superimposed effects. These effects are listed in the following:

• Restoring torque of the bellows due to natural curvature as explained in the next section.

• Pressure-dependent friction in the bellows due to different layers of the bellows.

• Pressure-dependent friction due to increase in bearing force.

• Curvature- and pressure-dependent contact area $A_{\mathrm{B},i,j}$ between rib and bellows.

• Material effects due to nonlinear material properties as well as friction effects between the bellows and the ribs.

The equation of motion for the bellows segment is derived by the Lagrangian formalism:

(2) $\underset{{\tau }_{\mathrm{B},i,\mathrm{motion}}}{\underbrace{M({\mathrm{\Psi }}_{\mathrm{i}}){\ddot{\mathrm{\Psi }}}_{\mathit{i}}+\mathit{C}({\mathrm{\Psi }}_i,{\dot{\mathrm{\Psi }}}_{\mathit{i}}){\dot{\mathrm{\Psi }}}_{\mathit{i}}+\mathit{N}({\mathrm{\Psi }}_i)}}={\tau }_{\mathrm{ext}}+\underset{{\tau }_{\mathrm{B},i,\mathrm{drive}}}{\underbrace{{\mathit{\tau }}_{\mathrm{B},i,\mathrm{a}}+{\tau }_{\mathrm{B},i,\mathrm{s}}+{\tau }_{\mathrm{B},i,\mathrm{f}}}},$(2)

where $\mathit{M}({\mathrm{\Psi }}_i)$ is the mass matrix, $\mathit{C}({\mathrm{\Psi }}_i,{\dot{\mathrm{\Psi }}}_{\mathit{i}}){\dot{\mathrm{\Psi }}}_{\mathit{i}}$ is the vector of Coriolis and centrifugal terms and $\mathit{N}({\mathrm{\Psi }}_i)$ is the vector of gravitational terms. The right side of the equation contains the external torque ${\tau }_{\mathrm{ext}}$, and the driving torque ${\tau }_{\mathrm{drive}}$, which consists of the actuator torque ${\tau }_{\mathrm{B},i,\mathrm{a}}$, the restoring torque ${\tau }_{\mathrm{B},i,\mathrm{s}}$ and the friction torque ${\tau }_{\mathrm{B},i,\mathrm{f}}$.

2.1. Restoring torque

Preliminary investigations show that the construction of the bellows chain and the differently meshed knitted fabric result in a natural curvature at different pressure levels, as shown in Figure 6. So, if, at a given bellows pressure $p_{\mathrm{B},i,j}$, the angle of a bellows segment ${\mathrm{\Psi }}_i$ deviates from the natural curvature, forces and torques will result which will push the bellows back to its natural curvature.

Figure 6. A bellows chain at different absolute pressure levels. The subfigures show different pressure levels and the resulting curvatures of the bellows chain. The bellows’ volume increases when the pressure rises.

Figure 6 shows a bellows chain consisting of eight bellows at different pressure levels. The bellows chain is already curved at ambient pressure of $1$ bar due to tighter knitting pattern at the right side of the bellows chain. The curvature increases for rising pressure since the left side elongates more as the right side, as there is a larger area on which the internal pressure acts. These two effects result in a curvature of the bellows causing a restoring torque.

The restoring torque of the bellows segment is given by,

(3a) ${\tau }_{\mathrm{B},i,\mathrm{s}}=\underset{X_{B,s}^T}{\underbrace{\sum _{\mathit{j}=1}^2\left[{\mathrm{\Psi }}_i-{\mathrm{\Psi }}_{\mathit{i},\mathit{j},0}\left(p_{\mathrm{B},i,j,\mathrm{scl}}\right)p_{\mathrm{B},i,j,\mathrm{scl}}\left({\mathrm{\Psi }}_i-{\mathrm{\Psi }}_{\mathit{i},\mathit{j},0}\left(p_{\mathrm{B},i,j,\mathrm{scl}}\right)\right)\right]}}\underset{{\theta }_{\mathrm{B},\mathrm{s}}}{\underbrace{\left[\begin{array}{c}{\mathit{k}}_{\mathrm{s}}\\ k_{\mathrm{s},\mathrm{p}}\end{array}\right]}},$(3a)

(3b) ${\tau }_{\mathrm{B},i,\mathrm{s}}={\mathit{X}}_{\mathrm{B},\mathrm{s}}^{\mathbf{T}}{\theta }_{\mathrm{B},\mathrm{s}},$(3b)

where ${\mathit{X}}_{\mathrm{B},\mathrm{s}}$ is the regressor and ${\theta }_{\mathrm{B},\mathrm{s}}$ is the parameter vector to be identified. The natural curvature of the bellows is described by the following ansatz function

(4) ${\mathrm{\Psi }}_{\mathit{i},\mathit{j},0}\left(p_{\mathrm{B},i,j,\mathrm{scl}}\right)=(-1{)}^{\mathit{j}+1}\left[p_{\mathrm{B},i,j,\mathrm{scl}}tanh\left(p_{\mathrm{B},i,j,\mathrm{scl}}\right)\right]\left[\begin{array}{c}{c}_{{\mathrm{\Psi }}_0}\\ c_{{\mathrm{\Psi }}_0,\mathrm{p}}\end{array}\right].$(4)

The latter is dependent on the relative scaled pressure in the bellows $p_{\mathrm{B},i,j,\mathrm{scl}}$ and is zero at ambient pressure. Therefore, the parameters of the natural curvature $c_{{\mathrm{\Psi }}_0}$ and $c_{{\mathrm{\Psi }}_0,\mathrm{p}}$ can only be identified if the pressure in the bellows is larger than the ambient pressure. So, the restoring torque from (3a) is divided into a pressure-independent part with the parameter $k_{\mathrm{s}}$ and a pressure-dependent part with the parameter $k_{\mathrm{s},\mathrm{p}}$.

2.2. Friction torque

Due to the superimposed layers of the bellows, internal bellows friction occurs as the pressure increases. This can be explained by the fact that the inner layers press against the outer layer at a given pressure $p_{\mathrm{B},i,j}$. At the same time, the layers move over each other at different curvatures, creating frictional forces.

In addition to the pressure-dependent friction within the bellows, the bearing force in the joints also increases as the pressure rises. If the pressure in both bellows chains $p_{\mathrm{B},i,1}$ and $p_{\mathrm{B},i,2}$ within a bellows segment rises, a greater force is exerted on the ribs, which must be absorbed by the plain bearings. At the pressure of $7.5$bar, the additional bearing force corresponds to a weight of more than $160$kg, which is distributed between two bearings. So, the bellows friction is given by,

(5a) ${\tau }_{\mathrm{B},i,\mathrm{f}}=\underset{X_{\mathrm{B},\mathrm{f}}^{\mathbf{T}}}{\underbrace{\left[tanh\left({\mathit{c}}_{\mathrm{B},\mathrm{vel}}{\dot{\mathrm{\Psi }}}_{\mathit{i}}\right)\sum _{\mathit{j}=1}^2{p}_{\mathrm{B},i,j,\mathrm{scl}}tanh\left(c_{\mathrm{B},\mathrm{vel}}{\dot{\mathrm{\Psi }}}_{\mathit{i}}\right)\right]}}\underset{{\theta }_{\mathrm{B},\mathrm{f}}}{\underbrace{\left[\begin{array}{c}{\mathit{c}}_{\mathrm{B},\mathrm{f}}\\ c_{\mathrm{B},\mathrm{f},\mathrm{p}}\end{array}\right]}}$(5a)

(5b) ${\tau }_{\mathrm{B},i,\mathrm{f}}={\mathit{X}}_{\mathrm{B},\mathrm{f}}^{\mathbf{T}}{\theta }_{\mathrm{B},\mathrm{f}},$(5b)

with $c_{\mathrm{B},\mathrm{vel}}=10$s. ${\mathit{X}}_{\mathrm{B},\mathrm{s}}$ is the regressor and ${\mathit{\theta }}_{\mathrm{B},\mathrm{s}}$ is the parameter vector to be identified. The bellows friction is divided into a pressure-independent part with the parameter $c_{\mathrm{B},\mathrm{f}}$ and a pressure-dependent part with the parameter $c_{\mathrm{B},\mathrm{f},\mathrm{p}}$.

2.3. Actuator torque

The hyperelastic material behaviour can be described by the energy density $\mathit{W}$. The energy density depends on the current strain of the material, the material density and material constants. The relationship between the stress and the strain of the material is determined via the derivative of the energy density. This allows the material behaviour to be accurately described with few parameters. Compared to material models generated by machine learning approaches, e.g. using neural networks, the energy density-based models generalize better [17].

Formally, any point of the bellows in the three-dimensional space is described by ${\mathit{x}}_{\mathrm{u}}\in {\mathbb{R}}^3$ in the undeformed state and by ${\mathit{x}}_{\mathrm{v}}\in {\mathbb{R}}^3$ in the deformed state. Thus, the deformation gradient is given by,

(6) $\mathit{F}=\frac{\mathrm{\partial }{\mathit{x}}_{\mathrm{v}}}{\mathrm{\partial }{\mathit{x}}_{\mathrm{u}}}.$(6)

The main strain directions are calculated via the eigenvalues of $\mathrm{Eig}(\mathit{F})={\lambda }_{1,2,3}$. The energy density is defined as a function of the main strain directions. Common material models are often represented using the left Cauchy-Green tensor and the invariants $\mathit{B}$ [18]

(7) $\mathit{B}=\mathit{F}{\mathit{F}}^{\mathrm{T}},$(7)

(8a) $I_1=\mathrm{tr}(\mathit{B})$(8a)

(8b) $={\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2$(8b)

(8c) $I_2=\frac{1}{2}\left(\mathrm{tr}{(\mathit{B})}^2+\mathrm{tr}({\mathit{B}}^T\mathit{B})\right)$(8c)

(8d) $={\lambda }_1^2{\lambda }_2^2+{\lambda }_2^2{\lambda }_3^2+{\lambda }_3^2{\lambda }_1^2$(8d)

(8e) $I_3=\det (\mathit{B})$(8e)

(8f) $={\lambda }_1^2{\lambda }_2^2{\lambda }_3^2.$(8f)

For a better interpretability, the energy density is defined as a function of the compressibility $J_{\mathrm{k}}=\sqrt{I_3}={\lambda }_1{\lambda }_2{\lambda }_3$. In general, the models used in the literature [19], such as the Neo-Hookean model [20], Mooney-Rivlin model [6] or Yeoh model [21], can be described by a polynomial hyperelastic material model with the energy density function

(9) $W_{\mathrm{poly}}=\sum _{k_1,k_2=0}^{N_{1,2}}{C}_{k_1,k_2}{\left({\bar{I}}_1-3\right)}^{k_1}{\left({\bar{I}}_2-3\right)}^{k_2}+\sum _{k_3=1}^{N_3}\frac{1}{C_{k_3}}{\left({\mathit{J}}_{\mathrm{k}}-1\right)}^{2k_3},$(9)

where ${\bar{I}}_1=J_{\mathrm{k}}^{-23}{I}_1,{\bar{I}}_2=J_{\mathrm{k}}^{-43}{I}_2$ and $k_1,k_2,k_3\in {\mathbb{N}}_0$ holds. The material parameters $C_{k_1,k_2}$ and $C_{k_3}$ are determined by tensile and compression tests. The order of the polynomial energy density function is determined by $N_{1,2}$ and $N_3$.

The incompressible Mooney-Rivlin model of the form is suitable for the description of both rubber and knitted fabrics [9–11]. However, large strains are applied. Therefore, the Ogden model, a material model for large deformations, is intended. The Ogden model [7] extends (9) by also allowing rational powers.

(10) $W_{\mathrm{Mooney}}=C_{1,0}\left({\bar{I}}_1-3\right)+C_{0,1}\left({\bar{I}}_2-3\right)$(10)

(11) $\mathit{W}=\sum _{\mathit{k}=1}^N\frac{{\mu }_k}{{\gamma }_k}\left({\lambda }_1^{{\gamma }_k}+{\lambda }_2^{{\gamma }_k}+{\lambda }_3^{{\gamma }_k}-3\right),$(11)

with ${\mu }_k,{\gamma }_k\in \mathbb{R}$. This makes the material model more complex, but thus a greater strain range can be efficiently represented [7,22].

It should be noted that for (9), (10) and (11) the assumption of an isotropic material is made, which is generally not fulfiled for knitted fabrics. If the material is considered to behave in an anisotropic manner, the energy density function has to be extended by the invariants of the fibre direction, which requires a projection of the deformation gradient in fibre direction. In order to obtain a solution for the position of all material points ${\mathit{x}}_{\mathrm{v}}$ over time in the load case, a computationally expensive finite element simulation is necessary in the three-dimensional case. This applies in particular if additional contact forces occur between the frame and the bellows. Since only the effects of the material on the torque of the bellows segment are of interest, the temporal solution of the position of all material points ${\mathit{x}}_{\mathrm{v}}$ is not relevant. Instead, the three-dimensional case is reduced to a two-dimensional problem in the further course, which also eliminates the need of a more complex anisotropic consideration of the energy density function.

For the approximation of the material model several assumptions are made and their effects on (11) are described below. By these assumptions, the relationship between the pressure and the material stress is established. The static solution of the position of all material points, results from the stresses in the material and the geometric relationships. It is assumed that the actuator torque of the bellows is described by the Ogden material model (11). Furthermore, it is assumed that rubber [7] and the knitted fabric [11] are incompressible materials, which means that the relation holds.

(12) ${\lambda }_1{\lambda }_2{\lambda }_3=1$(12)

The wall thickness of the bellows, which is $1$mm, is very small in comparison to the other dimensions, which are in the two-digit centimetre range. Therefore, the bellows is modelled as a membrane. So, the bending stiffness of the bellows is neglected. The contact between the bellows and the ribs is taken into account, which is referred as the membrane contact problem in the literature. This problem cannot be solved analytically even for simple geometric shapes, such as hemispheres, but requires the solution of a boundary value problem in the free case [23] and in the contact case without friction [12], as well as a finite element simulation in the contact case with friction [24].

In order to circumvent these numerically expensive procedures and to fulfil the requirement of real-time capability, a further simplification is made: The material behaviour is approximated by projecting the problem onto the midplane of the bellows segment. An illustration of the resulting two-dimensional model is shown in Figure 7. The figure shows the last joint of a bellows segment in a detailed view. The dashed black line indicates the bisector of the joint and at the same time forms the symmetry axis. The bellows itself is represented above the bisector by the dark blue circle segment and the light blue line. Below the angle bisector the bellows is indicated by the dashed lines. The curve parameter $\bar{a}\in [0,1]$ runs along the bellows material and assigns the corresponding point in the deformed state to a material point in the undeformed state. Thus, regardless of the deformation state, $\bar{a}=0$ describes the point at which the bellows is attached and $\bar{a}=1$ the point at which the bellows intersects the bisector. The variable $\mathit{a}\in [0,1)$ indicates what part of the material is in contact with the rib. Thus, $\bar{a}=\mathit{a}$ is always the transition point between the adjacent part and the free part of the bellows. Due to the small lever arm and the closely meshed knitted fabric on the side of the bellows facing the joint as shown in Figure 5, this side is neglected. In the bellows, there is a pressure of $p_{\mathrm{B},i,1}$, which presses against the membrane and thus expands it. The length $l_0$ indicates the distance between the centre of the joint and the rib. $s_0$ indicates the distance between the attachment point of the bellows $\bar{a}=0$ and the centre of the rib. For a simpler description the angle is introduced. The length of the adjacent line segment of the bellows is described by the variable $s_a$. The goal is to create a calculation rule for $s_a$. In this way, the surface $A_{\mathrm{B},i,j}$ is approximated and the effects of the material on the driving torque are considered.

Figure 7. Sketch of the two-dimensional no-slip membrane contact between bellows and ribs.

(13) $\mathit{\alpha }=-\frac{1}{2}{\psi }_{\mathit{i},4}$(13)

In the following it is assumed that ${\lambda }_1$ describes the strain of the material in the plane of Figure 7 and ${\lambda }_2$ describes the strain perpendicular to this plane. Thus, for the three-dimensional case, ${\lambda }_1$ describes the strain in the longitudinal direction and ${\lambda }_2$ describes the strain in the transverse direction. Due to the incompressibility (12) the strain in the normal direction is given by ${\lambda }_3=\frac{1}{{\lambda }_1{\lambda }_2}$, which is in the plane of Figure 7. In order to be able to evaluate (11), the following assumption is done regarding ${\lambda }_2$.

The transverse strain ${\lambda }_2$ is greatest along the centre of the bellows, i.e. in Figure 7 at $\bar{a}=1$, and decreases towards the attachment point of the bellows. Since the bellows is fixed at $\bar{a}=0$, the deformation gradient there corresponds to the unit matrix $\mathit{F}={\mathit{I}}_{3\times 3}$ and thus for the strain in the transverse direction holds ${\lambda }_2=1$. Given the setup, it is obvious that in the vicinity of $\bar{a}\approx 0$, the transverse strain takes values in the range ${\lambda }_2\approx 1$. For the calculation of $s_a$, the region close to the attachment, i.e. at $\bar{a}=0$, and less the region at $\bar{a}=1$ is crucial. Therefore, the assumption is made that ${\lambda }_2=1$ is fulfiled.

Thus, the energy density is represented as a function of ${\lambda }_1.$

(14) $\mathit{W}=\sum _{\mathit{k}=1}^N\frac{{\mu }_k}{{\gamma }_k}\left({\lambda }_1^{{\gamma }_k}+1+\frac{1}{{\lambda }_1^{{\gamma }_k}}-3\right).$(14)

From the material model, the material stress $\mathit{T}\in {\mathbb{R}}^1$ along the material is given by

(15a) $\mathit{T}=\frac{\mathrm{\partial W}}{\mathrm{\partial }{\lambda }_1}$(15a)

(15b) $=\sum _{\mathit{k}=1}^N{\mu }_k\left({\lambda }_1^{{\gamma }_k-1}-{\lambda }_1^{-{\gamma }_k-1}\right).$(15b)

In general, the material behaviour of rubbers is described with an order of the material model of $\mathit{N}=3$, where usually one parameter ${\gamma }_k<0$ [7,22]. Since (14) is invariant with respect to a negative ${\gamma }_k$, the negative parameter is neglected and the order is set to $\mathit{N}=2$. Thus, the material parameters to be identified are ${\mu }_1,{\mu }_2,{\gamma }_1$ and ${\gamma }_2$. The stresses in the material and the static solution of the free material are determined. The stress–strain curve is shown in Figure 8 for the identified material parameters.

Figure 8. Stress–strain curve of the bellows. The correlation between the material stress and the material strain of the identified ogden model is shown.

The relationship between the pressure $p_{\mathrm{B},i,j}$ and the stress $\mathit{T}$ for the free material results from the observation of an infinitesimally small section of material, illustrated in Figure 9. Using the force equilibrium in the $\mathit{y}$-direction follows

Figure 9. Infinitesimal material stresses arising due to pressure.

where the relative pressure is given by

(16) $0=(\mathit{T}+\mathrm{d}\mathit{T})cos(\mathrm{d}\mathit{\beta })-\mathit{T}cos(0)+{\int }_0^{\mathrm{d}\beta }-p_{\mathrm{B},i,j,\mathrm{rel}}\tilde{R}(\tilde{\beta })sin(\tilde{\beta })\mathrm{d}\tilde{\beta },$(16)

(17) $p_{\mathrm{B},i,j,\mathrm{rel}}=p_{\mathrm{B},i,j}-p_0.$(17)

With the relationship ${\int }_0^{\mathrm{d}x}\mathit{f}(\tilde{x})\mathrm{d}\tilde{x}=\mathit{f}(0)\mathrm{d}\mathit{x}$ and the assumption of small angles in the trigonometric functions, the following is obtained:

(18a) $0=(\mathit{T}+\mathrm{d}\mathit{T})\underset{\approx 1}{\underbrace{cos(\mathrm{d}\beta )}}-\mathit{T}cos(0)-\underset{=0}{\underbrace{{\mathit{p}}_{\mathrm{B},i,j,\mathrm{rel}}\tilde{R}(0)sin(0)\mathrm{d}\mathit{\beta }}}$(18a)

(18b) $0=\mathrm{d}\mathit{T}.$(18b)

From the force equilibrium in $\mathit{x}$-direction results

(19a) $0=\mathit{T}sin(\mathrm{d}\mathit{\beta })-p_{\mathrm{B},i,j,\mathrm{rel}}\tilde{R}(0)cos(0)\mathrm{d}\mathit{\beta }$(19a)

(19b) $\mathit{T}=p_{\mathrm{B},i,j,\mathrm{rel}}\tilde{R}(0).$(19b)

From (18b) follows that the stress in the free part, $\bar{a}>\mathit{a}$, is constant. If the pressure in the bellows $p_{\mathrm{B},i,j}$ is constant too, it follows from (19b) that the radius $\tilde{R}$ in the free part of the bellows must also be constant. This means that the solution of the line segment $\bar{a}>\mathit{a}$ must lie on a segment of a circle. For the transition point at $\bar{a}=\mathit{a}$, the condition of continuity of the material also applies. Thus, the arc is constructed geometrically with centre and radius. From the geometrical view the following results for the radius of the free material and for the angle of the circle section

(20) $\mathit{R}=l_0+\left(s_0+s_a\right)tan(\mathit{\alpha }),$(20)

(21) $\mathit{\beta }=\frac{\pi }{2}+\mathit{\alpha }.$(21)

From (19b) and (15b) results the determining equation for the strain of the free material

(22) $p_{\mathrm{B},i,j,\mathrm{rel}}\mathit{R}=\sum _{\mathit{k}=1}^N{\mu }_k\left({\lambda }_1^{{\gamma }_k-1}-{\lambda }_1^{-{\gamma }_k-1}\right).$(22)

The calculation of the strain is thus reduced to the solution of a static, nonlinear, one-dimensional equation, which can, for example, be efficiently solved by means of a characteristic curve.

The length of the free bellows from $\bar{a}=\mathit{a}$ to $\bar{a}=1$ is given by the strain, where

(23) $\mathit{L}=(1-\mathit{a}){\lambda }_1{l}_0\frac{\pi }{2},$(23)

and $l_0\frac{\pi }{2}$ is the length of the undeformed bellows section, which is calculated from $\mathit{\alpha }=0,s_a=0$ and $p_{\mathrm{B},i,j}=p_0$. At the same time, the length is also relevant for the geometric compatibility condition

(24) $\mathit{L}=\mathit{R\beta },$(24)

which is determined by means of angle and radius. The transition point $\mathit{a}$ is determined from these two equations. For a frictionless contact between rib and bellows, the strain of the free material ${\lambda }_1$ corresponds to the strain of the adjacent material ${\lambda }_{\mathrm{a}}$. However, measurements show that the assumption of frictionless contact is not correct, which is why a no slip is assumed. This means that a material point which is in contact with the frame $(\bar{a}<\mathit{a})$ do not changes its position and therefore the strain ${\lambda }_{\mathrm{a}}$ remains constant at this point if the pressure and the angle of curvature remain the same. At the transition point $\bar{a}=\mathit{a}$, the strain of the adjacent material corresponds to that of the free material due to the continuity of the material.

(25) $\mathit{a}=1-\frac{\mathit{R\beta }}{l_0\frac{\pi }{2}{\lambda }_1}.$(25)

(26) ${\lambda }_{\mathrm{a}}(\mathit{a})={\lambda }_1.$(26)

Since the strain ${\lambda }_1$ varies with the pressure, the result is a strain profile ${\lambda }_{\mathrm{a}}(\bar{a})$ over the adjacent material for $0\le \bar{a}\le \mathit{a}$. This is derived by using the history of (26). The length of the line segment $s_a$ is determined directly by means of the strain profile

(27) $s_a=l_0\frac{\pi }{2}{\int }_0^a{\lambda }_{\mathrm{a}}(\bar{a})\mathrm{d}\bar{a}.$(27)

The diameter of the rib is limited, so, it holds that $s_a\in \left[0,s_{\mathit{a},max}\right]$.

For the calculation of the line segment $s_a$, not only the strain ${\lambda }_{\mathrm{a}}(\bar{a})$ but also the transition point $\mathit{a}$ is necessary. These values represent the state of the bellows, but cannot be calculated explicitly. Therefore, an iterative calculation rule is presented in the following.

Initially at the start time $\mathit{a}=0$ is assumed. Then, iteratively in each time step first (27), (20) and (21) are calculated. Afterwards, the nonlinear equation (22)(22) $p_{\mathrm{B},i,j,\mathrm{rel}}\mathit{R}=\sum _{\mathit{k}=1}^N{\mu }_k\left({\lambda }_1^{{\gamma }_k-1}-{\lambda }_1^{-{\gamma }_k-1}\right).$(22) is solved and the state $\mathit{a}$ is adjusted using (25). Finally, ${\lambda }_{\mathrm{a}}$ is fitted according to (26) and again the length of the line segment (27) is calculated.

Since the solution of the integral (27) requires knowledge about ${\lambda }_{\mathrm{a}}(\bar{a})$, this function is approximated by $N_{\lambda }$ sampling points with ${\lambda }_{\mathrm{s}}\in {\mathbb{R}}^{N_{\lambda }}$. For a given strain profile ${\lambda }_{\mathrm{a}}(\bar{a})$, dependent on $\mathit{a}$, the first $n_{\lambda }=\mathrm{ceil}\left(\mathit{a}{N}_{\lambda }\right)+1$ sampling points of ${\lambda }_{\mathrm{s}}$ are assigned values. The function ceil$(\cdot )$ rounds up the argument to the next integer value. The approximation of the length of the line segment is done by the function.

(28) ${\tilde{s}}_{\mathit{a}}=l_0\frac{\pi }{2}\frac{1}{N_{\lambda }}\sum _{\mathit{k}=1}^{n_{\lambda }}{\lambda }_{\mathrm{s},k},$(28)

where ${\lambda }_{\mathrm{s},k}$ corresponds to the $\mathit{k}$^th entry of ${\lambda }_{\mathrm{s}}$. By this computational rule, it is possible to calculate the material behaviour in real time. Thus, an approximation of the contact length ${\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}$ for each bellows exist. The pseudocode for the calculation of ${\tilde{s}}_{\mathit{a}}$ is presented in algorithm 1 in Appendix A. The convergence of the iteration rule of $s_a$ is shown in Appendix B.

For the mapping from ${\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}$ to ${\tau }_{\mathrm{B},i,\mathrm{a}}$ following considerations are made: From the physical point of view, it is clear that ${\tau }_{\mathrm{B},i,\mathrm{a}}$ is composed of a surface on the rib, an effective lever arm, and the pressure in the bellows. It is obvious that the actual surface is at least quadratically dependent on ${\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}$. Since ${\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}$ is only an approximation, which shows deviations with respect to the curvature angle ${\mathrm{\Psi }}_i$ and the bellows pressure $p_{\mathrm{B},i,j}$, it is recommended to combine these quantities. The scaled relative pressure is described by $p_{\mathrm{B},i,j,\mathrm{scl}}$ and is zero if the pressure in the bellows corresponds to the ambient pressure. In this case, the driving torque, which originates from the $\mathit{j}$^th bellows, disappears. These considerations lead to the following ansatz function.

(29a) ${\mathit{X}}_{\mathrm{B},\mathrm{a}}=\sum _{\mathit{j}=1}^2{p}_{\mathrm{B},i,j,\mathrm{scl}}\left[\begin{array}{c}1\\ {\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}\\ {\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}^2\\ {\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}{p}_{\mathrm{B},i,j,\mathrm{scl}}\\ {\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}{\mathrm{\Psi }}_{\mathit{i},\mathit{j}}\\ {\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}{\mathrm{\Psi }}_{i,j}^2\\ p_{\mathrm{B},i,j,\mathrm{scl}}\\ {\mathrm{\Psi }}_{\mathit{i},\mathit{j}}\\ {\mathrm{\Psi }}_{i,j}^2\end{array}\right],$(29a)

(29b) ${\tau }_{\mathrm{B},i,\mathrm{a}}=X_{\mathrm{B},\mathrm{a}}^{\mathbf{T}}{\theta }_{\mathrm{B},\mathrm{a}},$(29b)

where ${\theta }_{\mathrm{B},\mathrm{a}}$ describes the parameter to be identified. Based on the assumption that the behaviour of the bellows is the same, but the bellows are mirrored within a segment, the angle ${\mathrm{\Psi }}_i$ for the opposing bellows is negative, i.e. ${\mathrm{\Psi }}_{\mathit{i},1}={\mathrm{\Psi }}_i$ and ${\mathrm{\Psi }}_{\mathit{i},2}=-{\mathrm{\Psi }}_i$.

3. Identification of the driving torque model

For the identification of the restoring torque, friction torque and actuator torque parameters the following measurements are done:

• Variation of the angle of curvature with vented bellows.

• Variation of the angle of curvature with pressurized bellows.

• Simultaneous variation of the angle of curvature and pressure in the bellows.

Hereby, the angles and the pressures of the bellows segments are measured by encoders and pressure sensors. The parameters to be identified are ${\theta }_{\mathrm{B}}={\left[{\theta }_{\mathrm{B},\mathrm{a}}^{\mathrm{T}}{\theta }_{\mathrm{B},\mathrm{s}}^{\mathrm{T}}{\theta }_{\mathrm{B},\mathrm{f}}^{\mathrm{T}}\right]}^{\mathrm{T}}$ and ${\mu }_1,{\mu }_2,{\gamma }_1$, ${\gamma }_2$. For the given measurements with a total of $N_{\mathrm{MP}}$ measurement points.

(30a) $Y_{\mathrm{B}}^k={\tau }_{\mathrm{B},i,\mathrm{m}\mathrm{o}\mathrm{tion}}^k-{\tau }_{\mathrm{e}\mathrm{x}\mathrm{t}}^k,$(30a)

(30b) ${\mathit{Y}}_{\mathrm{B}}=\left[\begin{array}{c}{Y}_{\mathrm{B}}^1\\ ⋮\\ Y_{\mathrm{B}}^{{\mathit{N}}_{\mathrm{MP}}}\end{array}\right],$(30b)

where the $\mathit{k}$^th measurement point is marked by the superscript index, and the given regressor consisting of the actuator torque regressor (29b), the restoring torques regressor (3b) and the friction torque regressor (5b), the optimal parameters are computed by

(31) ${\mathit{X}}_{\mathrm{B}}=\left[\begin{array}{ccc}{X}_{\mathrm{B},\mathrm{a}}^1& X_{\mathrm{B},\mathrm{s}}^1& X_{\mathrm{B},\mathrm{f}}^1\\ ⋮& ⋮& ⋮\\ X_{\mathrm{B},\mathrm{a}}^{{\mathit{N}}_{\mathrm{MP}}}& X_{\mathrm{B},\mathrm{s}}^{{\mathit{N}}_{\mathrm{MP}}}& X_{\mathrm{B},\mathrm{f}}^{{\mathit{N}}_{\mathrm{MP}}}\end{array}\right],$(31)

(32) ${\theta }_{\mathrm{B}}={\left({\mathit{X}}_{\mathrm{B}}^{\mathrm{T}}{\mathit{X}}_{\mathrm{B}}\right)}^{-1}{\mathit{X}}_{\mathrm{B}}^{\mathrm{T}}{\mathit{Y}}_{\mathrm{B}}.$(32)

It should be noted that ${\mathit{X}}_{\mathrm{B}}$ depends on ${\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}^{\mathit{k}}$, which requires a forward simulation of the bellows. Since ${\tilde{s}}_{\mathit{a},\mathit{i},\mathit{j}}^{\mathit{k}}$ in turn depends on the material parameters ${\mu }_1,{\mu }_2,{\gamma }_1$ and ${\gamma }_2$, ${\mathit{X}}_{\mathrm{B}}\left({\mu }_1,{\mu }_2,{\gamma }_1,{\gamma }_2\right)$ also depends on these parameters. Since the material parameters are not given, ${\theta }_{\mathrm{B}}$ cannot be determined directly according to (32). However, the parameters ${\theta }_{\mathrm{B}}$ result from the material parameters, which is why they can be solved successively in an optimization problem. The following optimization problem is solved to identify the material parameters ${\gamma }_k>2$ results from the convergence of the iteration rule of ${\tilde{s}}_{\mathit{a}}$, which is derived in the appendix B. The optimization problem is solved for several initial values since its solution usually has a local minimum.

(33a) $\underset{{\mu }_1,{\mu }_2,{\gamma }_1,{\gamma }_2}{min}{\left({\mathit{Y}}_{\mathrm{B}}-{\mathit{X}}_{\mathrm{B}}{\theta }_{\mathrm{B}}\right)}^{\mathrm{T}}\left({\mathit{Y}}_{\mathrm{B}}-{\mathit{X}}_{\mathrm{B}}{\theta }_{\mathrm{B}}\right),$(33a)

(33b) $\mathrm{s}.\mathrm{t}.(32),$(33b)

(33c) ${\gamma }_k>2.$(33c)

4. Validation and discussion

In this section, the driving torque model is validated and discussed. In addition to the measured external torque ${\tau }_{\mathrm{ext}}$ and the torque ${\tilde{\tau }}_{\mathrm{ext}}$ predicted from the driving torque model, the estimated torque of the previous driving torque model from [16] is also shown. The parameters from the previous driving torque model are identified using the same data as the model developed in this paper.

(34) ${\tilde{\tau }}_{\mathrm{ext},\mathrm{prev}}={\bar{\tau }}_{\mathrm{B},\mathit{i},\mathrm{f}}+{\bar{\tau }}_{\mathrm{B},\mathit{i},\mathrm{s}}+\sum _{\mathit{j}=1}^2{p}_{\mathrm{B},i,j,\mathrm{scl}}\left[1p_{\mathrm{B},i,j,\mathrm{scl}}{\mathrm{\Psi }}_{\mathit{i},\mathit{j}}{p}_{\mathrm{B},i,j,\mathrm{s}\mathrm{c}\mathrm{l}}^2\right]{\theta }_{\mathrm{prev}}$(34)

The angles of each joint in a bellows segment ${\mathrm{\Psi }}_{\mathit{i},\mathit{j}}$ and the pressures of each bellows chain $p_{\mathrm{B},i,j}$ are measured by encoders and pressure sensors. Figure 10 illustrates the hysteresis of the bellows, shown for the fourth bellows segment. Before external movement, the first bellows chain is pressurized with a pressure of $4$bar at point in time $\mathit{t}=28$s. The second bellows chain remains at ambient pressure of $1$bar. The angle of the bellows segment is ${\mathrm{\Psi }}_4={74.5}^{\circ }$ and the measured torque by an external force-torque sensor is about $6.7$Nm at this time. After external movement, the bellows segment is returned to same state regarding pressure and angle as at $\mathit{t}=28$s. However, the measured torque is about $9.1$Nm showing the hysteresis of the bellows. The previous driving torque model [16] is not able to consider hysteresis since it neglects the bellows’ material properties. Before and after external movement, the previous driving torque model estimates the same external torque of $9.9$Nm, which corresponds to an error of $3.2$Nm with respect to the measured value. Even in the best case, the previous driving torque model could only estimate the mean value of the hysteresis difference, which would be an error of about $1.2$Nm. In comparison, the presented driving torque model achieves an error of about $0.6$Nm in the worst case. Thus, the error of the previous driving torque model is halved by considering the material properties. For the application considered here, the presented driving torque model reduces the error from $3.2$Nm to about $0.6$Nm with an absolute value of about $6.7$Nm.

Figure 10. High dynamic validation of the bellows’ driving torque model. The angle ${\mathrm{\Psi }}_4$, the pressure $p_{\mathrm{B},4,j}$ and the approximated membrane contact length ${\tilde{s}}_{\mathit{a},4,\mathit{j}}$ are exemplarily shown for the fourth bellows segment. The measured external torque ${\tau }_{\mathrm{ext}}$ and the estimated external torque ${\tilde{\tau }}_{\mathrm{ext}}$ are also illustrated. The pressure in the bellows is set by a pressure controller and the position is applied to the system by an external force.

The hysteresis is explained by Figure 11. The top measurement plot shows ${\lambda }_1$ over time of the first bellows chain of the fourth bellows segment from Figure 10. While the angle ${\mathrm{\Psi }}_4$ decreases during external movement, the angle $\mathit{\alpha }$ increases, stretching the bellows, which results in a larger strain. At the same time, the contact of the bellows with the rib, so, the transition point $\mathit{a}$, decreases. Subsequently, ${\mathrm{\Psi }}_4$ increases again and the strain for ${\lambda }_1$ decreases. As the bellows are compressed, the more stretched material snuggles the rib, increasing the contact area. The strain profile is shown over time in Figure 11. The grey line shows the transition point $\mathit{a}$ in the strain profile. It is shown how, in the interval $40\mathrm{s}\le \mathit{t}\le 50\mathrm{s}$ the strain profile builds up and $\mathit{a}$ increases again. In addition, Figure 11 shows that towards the end of the measurement, the strain profile is significantly larger than at the beginning of the measurement. The additional measured torque is due to the increased contact area on the rib, which is caused by the preceding strain of the bellows.

Figure 11. Material behaviour of the first bellows chain of the fourth bellows segment. For the pressurized bellows from Figure 10, the time course of the strain of the free material ${\lambda }_1$ and the transition point $\mathit{a}$ as well as the strain profile are shown. In addition, $\mathit{a}$ is also drawn in grey in the strain profile. The strain of the free material ${\lambda }_1$ after the ventilation of the bellows at $\mathit{t}=28\mathrm{s}$ as well as the additional strain due to the stretching of the bellows between $\mathit{t}=34\mathrm{s}$ and $\mathit{t}=45\mathrm{s}$, caused by the externally imposed movement of the bellows segment, are illustrated.

In order to make a final evaluation of the driving torque model, another validation measurement is presented in which both the pressure and the angle of the bellows segment vary simultaneously. The results are shown in Figure 12. The estimated external torque curve ${\tilde{\tau }}_{\mathrm{ext}}$ mimics the measured torque curve ${\tau }_{\mathrm{ext}}$. The signals are largely superimposed despite the highly dynamic measurement with abrupt changes in motion and rapidly changing pressures. Deviations mostly occur when the direction of motion changes, such as at $\mathit{t}\approx 6$s, $\mathit{t}\approx 20$s and $\mathit{t}\approx 24$s. It is assumed that these deviations are caused by model inaccuracies for difficult replicable friction effects, which are particularly noticeable in the case of zero crossings in the velocity. The absolute error averaged over time between the measured and the estimated torque of the dynamic validation measurement is $0.85$Nm. Overall, it can be concluded that the courses are mapped well and that the expected error in the torque curve is below $1$Nm.

Figure 12. High dynamic validation of the bellows’ driving torque model. The angle ${\mathrm{\Psi }}_4$, the pressure $p_{\mathrm{B},4,j}$ and the approximated membrane contact length ${\tilde{s}}_{\mathit{a},4,\mathit{j}}$ are exemplarily shown for the fourth bellows segment. The measured external torque ${\tau }_{\mathrm{ext}}$ and the estimated external torque ${\tilde{\tau }}_{\mathrm{ext}}$ are also illustrated. The pressure in the bellows is set by a pressure controller and the position is applied to the system by an external force.

5. Conclusion

In this article, the driving torque model of a bellows segment of the bionic soft arm is derived and validated experimentally regarding real-time capability to enable model-based control design. The driving torque model considers the effects of the bellows causing restoring torques, friction torques and actuator torques. The restoring torque appears due to the bellows construction. Differently meshed knitted fabric on the bellows amplifies the curvatures of a bellows chain and thus the restoring torque. The friction torque is caused by friction between the superimposed layers of the bellows and by bearing force. The actuator torque is based on the Ogden model, a material model which considers the hyperelastic behaviour of the bellows. The bellows’ torque transmission is described by a no-slip membrane contact between bellows and ribs. The bellows are assumed to have membrane characteristics due to their thin walls. In order to make the driving torque model suitable for control design, real-time capability of the model is required. So, the no-slip membrane contact is considered to be two-dimensional instead of three-dimensional by mapping it onto the midplane of the bellows. Hence, a contact line is considered instead of a contact area and its iterative algorithm is derived and presented. Further assumptions are made to simplify the driving torque model. The bellows are assumed to be incompressible and the strain only in longitudinal direction of the bellows is considered. The parameters of the driving torque model are identified by regression using experimental data. The performance of the driving torque model is evaluated experimentally by measuring the external torque with a force-torque sensor for static and dynamic validation including external movements.

Overall, the driving torque model is validated successfully in static and high dynamic experiments by considering friction, restoring torques and hysteresis effects by the bellows’ nonlinear material behaviour. It is suitable for model-based control design and subsequent trajectory generation due to its real-time capability. In future work, the driving torque model is included in closed-loop control to enable autonomous operations of the bionic soft arm.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) under grant no. SA 847/20-2.

Appendix A.

Iterative Calculation Rule of ${\tilde{s}}_{\mathit{a}}$

The pseudocode for the calculation of the approximation of the length of the line segment ${\tilde{s}}_{\mathit{a}}$ is presented in algorithm 1.

Algorithm 1

Iterative calculation rule of ${\tilde{s}}_{\mathit{a}}$

Initialize:

$a^0\leftarrow 0$$n_{\lambda }\leftarrow 1$

${\lambda }_{\mathrm{s}}^0$ is empty${\lambda }_{\mathrm{s},1}^0\leftarrow 1$

end

At each time step $k_{\mathrm{t}}$:

${\tilde{s}}_{\mathit{a}}^{k_{\mathrm{t}}}\leftarrow l_0\frac{\pi }{2}\frac{1}{N_{\lambda }}{\sum }_{\mathit{k}=1}^{n_{\lambda }}{\lambda }_{\mathrm{s},k}$$p_{\mathrm{B},i,j,\mathrm{rel}}\leftarrow p_{\mathrm{B},i,j}^{{\mathit{k}}_{\mathrm{t}}}-p_0$$\mathit{\alpha }\leftarrow -\frac{1}{2}{\psi }_{\mathit{i},\mathit{j}}^{{\mathit{k}}_{\mathrm{t}}}$$\mathit{R}\leftarrow l_0+(s_0+{\tilde{s}}_{\mathit{a}}^{k_{\mathrm{t}}})tan(\mathit{\alpha })$$\mathit{\beta }=\frac{\pi }{2}+\mathit{\alpha }$

${\lambda }_1\leftarrow$ solve (22) for ${\lambda }_1$$a^{{\mathit{k}}_{\mathrm{t}}}\leftarrow 1-\frac{\mathit{R\beta }}{l_0\frac{\pi }{2}{\lambda }_1}$$n_{\lambda }\leftarrow \mathrm{ceil}(a^{{\mathit{k}}_{\mathrm{t}}}{N}_{\lambda })+1$

${\lambda }_{\mathrm{s}}^{{\mathit{k}}_{\mathrm{t}}}\leftarrow$ Clear the last $N_{\lambda }-n_{\lambda }$ fields of ${\lambda }_{\mathrm{s}}$${\lambda }_{\mathrm{s},n_{\lambda }}^{k_{\mathrm{t}}}\leftarrow {\lambda }_1$

$\mathcal{E}\leftarrow$ Find all empty fields of ${\lambda }_{\mathrm{s}}^{{\mathit{k}}_{\mathrm{t}}}$ between 1 and $n_{\lambda }$

if $\mathcal{E}$ is not empty then${\lambda }_{\mathrm{s},i}^{{\mathit{k}}_{\mathrm{t}}}\leftarrow {\lambda }_1\mathrm{\forall }\mathit{i}\in \mathcal{E}$

end if

end

Initially, $\mathit{a}=0$ is assumed. Then (27), (17), (13), (20) and (21) are calculated iteratively at each time step. The nonlinear equation (22)(22) $p_{\mathrm{B},i,j,\mathrm{rel}}\mathit{R}=\sum _{\mathit{k}=1}^N{\mu }_k\left({\lambda }_1^{{\gamma }_k-1}-{\lambda }_1^{-{\gamma }_k-1}\right).$(22) is then solved and the variable $\mathit{a}$ is updated via (25). Finally, ${\lambda }_{\mathrm{a}}$ is updated according to (26) and the length of the line segment (27) is recalculated. The algorithm converges in only a few steps to the desired value, which was determined in preliminary experiments. For this reason, only one iteration of the algorithm is required per time step keeping the computational complexity low.

Appendix B.

Convergence of the Iteration Rule of $s_a$

In the following, it is shown that the calculation rule for the length of the line segment $s_a$ converges to a fixed value at static pressure $p_{\mathrm{B},i,j,\mathrm{rel}}$ and static angle of curvature $\mathit{\alpha }$. If $s_a^{l+1}=\mathit{\phi }\left(s_a^l\right)$ is the iteration rule, it converges according to Banach’s fixed point theorem if a Lipschitz constant $\mathit{L}<1$ can be determined for $\mathit{\phi }\left(s_a^l\right)$. This is equivalent to $\frac{\mathrm{d}\varphi \left(s_a^l\right)}{\mathrm{d}{s}_a^l}\le \mathit{L}<1$. If ${\mathrm{\Lambda }}_{\mathrm{a}}$ is the antiderivative of ${\lambda }_{\mathrm{a}}$, the following results from (27) according to the Leibniz rule for integrals.

(B1a) $\frac{\mathrm{d}\phi \left(s_a^l\right)}{\mathrm{d}{s}_a^l}=l_0\frac{\pi }{2}\frac{\mathrm{d}{\mathrm{\Lambda }}_{\mathrm{a}}}{\mathrm{d}{s}_a^l}$(B1a)

(B1b) $=l_0\frac{\pi }{2}{\lambda }_{\mathrm{a}}(\mathit{a})\frac{\mathrm{d}a}{\mathrm{d}{s}_a^l}$(B1b)

(B1c) $\stackrel{(25)}{=}{l}_0\frac{\pi }{2}{\lambda }_1\frac{\mathrm{d}}{\mathrm{d}{s}_a^l}\left(1-\frac{\mathit{R\beta }}{l_0\frac{\pi }{2}{\lambda }_1}\right)$(B1c)

(B1d) $=l_0\frac{\pi }{2}{\lambda }_1\left(-\frac{\frac{\mathrm{d}R}{\mathrm{d}{s}_a^l}\mathit{\beta }}{l_0\frac{\pi }{2}{\lambda }_1}+\frac{\mathit{R\beta }}{l_0\frac{\pi }{2}}\frac{\mathrm{d}{\lambda }_1^{-1}}{\mathrm{d}{s}_a^l}\right)$(B1d)

(B1e) $\stackrel{(20)}{=}{l}_0\frac{\pi }{2}{\lambda }_1\left(-\frac{tan(\mathit{\alpha })\mathit{\beta }}{l_0\frac{\pi }{2}{\lambda }_1}+\frac{\mathit{R\beta }}{l_0\frac{\pi }{2}{\lambda }_1^2}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}\mathit{R}}\frac{\mathrm{d}R}{\mathrm{d}{s}_a^l}\right)$(B1e)

(B1f) $=-tan(\mathit{\alpha })\mathit{\beta }+\frac{\mathit{R\beta }}{{\lambda }_1}\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}\mathit{R}}tan(\mathit{\alpha }).$(B1f)

Since ${\lambda }_1$ cannot be expressed explicitly as a function of $\mathit{R}$ according to (22), the derivative $\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}\mathit{R}}$ is determined explicitly using the theorem for implicit functions

(B2) $\frac{\mathrm{d}{\lambda }_1}{\mathrm{d}\mathit{R}}=\frac{{\mathit{p}}_{\mathrm{B},i,j,\mathrm{rel}}}{{\sum }_1^{N=2}{\mu }_i\left(\left({\gamma }_k-1\right){\lambda }_1^{{\gamma }_k-2}+\left({\gamma }_k+1\right){\lambda }_1^{-{\gamma }_k-2}\right)}.$(B2)

With (22) follows from (B1f)

(B3) $\frac{\mathrm{d}\phi \left(s_a^l\right)}{\mathrm{d}{s}_a^l}=\mathit{\beta }tan(\mathit{\alpha })\left(\frac{{\sum }_i^{N=2}{\mu }_i\left({\lambda }_1^{{\gamma }_k-1}-{\lambda }_1^{-{\gamma }_k-1}\right)}{{\sum }_i^{N=2}{\mu }_i\left(\left({\gamma }_k-1\right){\lambda }_1^{{\gamma }_k-1}+\left({\gamma }_k+1\right){\lambda }_1^{-{\gamma }_k-1}\right)}-1\right).$(B3)

$\mathit{\beta }tan(\mathit{\alpha })<1$ holds for all permissible $\mathit{\alpha }$. On the relevant range, i.e. ${\lambda }_1>1$, it holds that ${\lambda }_1^{{\gamma }_k-1}>0$ and ${\lambda }_1^{-{\gamma }_k-1}>0$. Whether the numerator or denominator is larger depends on which coefficients are larger. In particular, the fraction for ${\gamma }_k>2$ is positive for the relevant range and also smaller than 1 since $1<\left({\gamma }_k-1\right)\mathrm{\forall }{\gamma }_k>2$. The inequality follows from this argument

(B4) $\sum _{\mathit{i}=1}^{\mathit{N}=2}{\mu }_i\left({\lambda }_1^{{\gamma }_k-1}-{\lambda }_1^{-{\gamma }_k-1}\right)<\sum _{\mathit{i}=1}^{\mathit{N}=2}{\mu }_i\left(\left({\gamma }_k-1\right){\lambda }_1^{{\gamma }_k-1}+\left({\gamma }_k+1\right){\lambda }_1^{-{\gamma }_k-1}\right).$(B4)

Thus, the expression in the parenthesis in (B3) for the relevant range is also less than 1, which means that is fulfilled and the calculation rule converges.

(B5) $\frac{\mathrm{d}\phi \left(s_a^l\right)}{\mathrm{d}{s}_a^l}\le 1$(B5)