Real-time Dynamics of Soft Manipulators with Cross-section Inflation: Application to the Octopus Muscular Hydrostat

In plain words

A soft robot arm bends and stretches in ways a jointed robot arm cannot, much like an octopus arm. Engineers usually model such a slender body as a thin rod whose cross-section keeps its shape. That assumption breaks when the arm is squeezed sideways by its own muscles, because the cross-section then changes size, and the simple rod model misses the resulting force.

This paper adds one extra variable that says how much each cross-section grows or shrinks in its own plane. It also gives the material a realistic kind of stiffness and internal damping. The full equations are too heavy for real time, so the authors describe the arm's shape with a small set of smooth building blocks. That turns the problem into a handful of ordinary differential equations that can be solved quickly.

They use the model for two tasks: tuning how stiff a soft rod is by squeezing it, and simulating an octopus-like arm reaching for and fetching a target. The authors report that the new model keeps the arm's volume nearly constant, which the classic model does not, and that its motion looks more like a real octopus.

TL;DR

• Extends Cosserat rod dynamics with an in-plane inflation ratio ρ(s), so the cross-section can change size under lateral (transverse) loads §1.
• Uses a Saint-Venant Kirchhoff elastic law with Kelvin-Voigt viscous damping, reduced by a small-strain approximation and a near-incompressibility choice §2.2.
• Reduces the equations with the Geometric Variable Strain approach, using Legendre bases for strains and Hermite splines for ρ, and integrates them numerically §4, app. D–E.
• Author-reported results: an octopus-like arm elongates by about 20% of the bend-point travel with volume change under 3%, against 18–36% volume change for the classic model §5.2.
• Code is public (MATLAB, MIT). Several config values and table entries differ from the paper; see Paper vs code.

Problem

Classic Cosserat rod theory lets a slender body bend, twist, shear and stretch, but it keeps every cross-section rigid. It therefore cannot represent in-plane deformation during tasks with active lateral traction, nor a hollow tube §1. Full 3D finite-element models capture this but are too expensive for real-time simulation and control §1. Reduced-order GVS models are cheap but inherit the rigid-section assumption §1.

Assumptions that matter for implementation:

• Small strain: the displacement-gradient norm is taken to be much smaller than one, and the quadratic term is dropped §2.1, after eq. 2.
• Axisymmetric cross-section, reflectively symmetric, with the inflation ratio ρ applied uniformly across it §2.1; §2.2 text before eq. 7.
• Homogeneous body force over each cross-section, so gravity and buoyancy give no moment about the centroid §3.1; §3.3.
• Isotropic, nearly incompressible material: Poisson's ratio close to 0.5 §2.2, after eq. 5.
• The model does not cover anisotropic cross-sections or out-of-plane deformation §6.

Method

Kinematics §2.1. A local frame is attached to the rod centreline at arclength s. The position of a material point is the centreline plus ρ(s) times its cross-section offset in the two directors. The deformation gradient is F = R(D + I), where D collects the strain variables (ν1, ν2, ν3) for shear and stretch, (κ1, κ2, κ3) for bending and twist, and the in-plane inflation ratio ρ and its derivative ρ′. The small-strain tensor keeps only the linear part of D eq. 1, eq. 2.

Constitutive law §2.2. The elastic energy is the Saint-Venant Kirchhoff form. Two modelling choices make it usable:

1. For isotropic material with Poisson ratio ν0 tending to 0.5, the normal in-plane stress is forced to vanish, which removes the λ → ∞ blow-up §2.2, after eq. 5. The result is an energy that depends on the shear and stretch strains and the Young's modulus eq. 5.
2. ρ is kept close to one inside a, b and c, so the energy per unit length is quadratic in the strain variables. This gives Φ_E eq. 7. The authors note that this form looks ad hoc but is justified by the internal-force analysis that follows §2.2, text after eq. 7.

Kelvin-Voigt damping adds a dissipation density Φ_D eq. 8. Differentiating the total energy by chain rule defines the stress resultants n, m (force and moment), Q (weighted shear traction) and q (weighted normal traction) eq. 9–14. Explicit forms, with reference strains ν* and κ* and the lateral part n^d, are eq. 15–18. The classic Cosserat resultants are the same with the inflation terms removed eq. 19–20.

Dynamics §2.3. Momentum balance on a segment gives the local equations for linear and angular momentum eq. 24, eq. 25. The material momentum balance for the cross-section gives the inflation equation eq. 26–28. The three sets of equations need 14 boundary-condition sets. The inflation ODE is second order in ρ and needs two boundary conditions, one of the three types in tab. 1.

Lateral (transverse) actuation §3.1. A homogeneous inward pressure f0 on the cross-section boundary gives a lateral traction r = 2 f0 π z² eq. 29. Here z is the reference radius.

Longitudinal and oblique actuation §3.2. A tendon of tension T along a route gives force n_a = T t and moment m_a = T d × t, where t is the unit tangent of the route in the local frame eq. 30.

External loads §3.3. The total external force per unit length is the sum of gravity/buoyancy, added mass, drag and point loads; the external moment is the point moment eq. 31, eq. 32. Added mass eq. 33 and drag eq. 34 are written for a circular cross-section, and point loads are smeared with a Dirac function eq. 35.

Reduced-order solver (GVS) §4. The strain twist and inflation ratio are expanded in basis functions:

\boldsymbol\xi(s,\mathbf q_\xi) = \boldsymbol\Phi_\xi(s)\,\mathbf q_\xi + \boldsymbol\xi^*(s) eq. 42

\rho(s,\mathbf q_\rho) = \boldsymbol\Phi_\rho(s)\,\mathbf q_\rho + \rho^*(s) eq. 48

The configuration is recovered by a recursive Magnus-based update along the rod eq. 40, eq. 43, eq. 44. The geometric Jacobian J and its derivative give the velocities and accelerations eq. 45–47. Projecting the momentum and inflation equations onto the joint spaces gives the reduced equations of motion eq. 53, eq. 54. Static equilibrium uses the same matrices with time derivatives removed eq. 55, eq. 56.

Activation §5.2, §5.3. Octopus actuators are driven by a travelling sigmoid wave eq. 57 whose amplitude and position vary with time; the values are in the Activation table below.

Numerics app. C, app. E. The Magnus expansion uses a fourth-order Zannah quadrature. Integrals along the rod use Gauss-Legendre quadrature. The dynamics are integrated with MATLAB ode15s; statics are solved with MATLAB fsolve with an initial guess.

Pseudocode

Resultants at one cross-section, from the explicit forms eq. 15–18:

# Lateral stiffness K_l and damping B_l are diagonal  [§2.2, after eq. 18]
# Inputs: nu, kappa (3 each), rho, rho_prime, and their time derivatives
def section_resultants(nu, kappa, rho, rho_prime, d_nu, d_kappa, d_rho):
    # Force resultant n  [eq. 15]
    n = R @ (K_l @ (nu - nu_ref) + B_l @ d_nu + n_d)  # (3,)
    # Moment resultant m  [eq. 16]
    m = R @ (K_theta @ (kappa - kappa_ref) + B_theta @ d_kappa)
    # Weighted shear traction Q  [eq. 17]
    Q = mu * I33 * rho_prime + eta * I33 * d_rho_prime
    # Weighted normal traction q  [eq. 18]
    q = (4 * (lam + mu) * A0 * (rho - 1)
         + 2 * lam * A0 * (nu[2] - 1)
         + 4 * eta * A0 * d_rho)
    return n, m, Q, q

Octopus actuation waves, in each actuator's own row of the activation table eq. 57:

# Travelling sigmoid, X = s / L in [0, 1]  [eq. 57]
# alpha_t: amplitude, mu_t: wave position, sigma: ramp steepness  [tab. 3]
def activation(X, alpha_t, mu_t, sigma):
    # Leading minus sign as written in the paper  [eq. 57]
    ramp = 1 - 1 / (1 + exp(-sigma * (X - mu_t)))
    return -alpha_t * ramp

Reduced-order step, from the matrices in app. D and eq. 53–54:

# Basis matrices on the rod grid, polynomial for xi, Hermite for rho  [§4.1]
phi_xi = legendre_basis(quadrature_points)  # (6, n_xi) per point
phi_rho = hermite_basis(quadrature_points)  # (1, n_rho) per point
# Recursive kinematics along the rod, then integrate  [eq. 43–46] [app. E]
J_mats = recurse_jacobian(q_xi, phi_xi)     # (6, n_xi) per point
M_xi = sum(w * J.T @ M_bar @ J for each point)
# Equations of motion  [eq. 53–54]
# Solve for accelerations, then advance with ode15s  [app. E]
q_xi_dd, q_rho_dd = solve_reduced_system(M_xi, M_rho, state)

Hyperparameters

Octopus-arm values, all from the paper tab. 2 and §5.2:

Axial stiffness test §5.1:

Fetching motion §5.3:

Activation functions. Each row is the amplitude α(t), position μ(t) and steepness σ used in eq. 57:

Other numerical settings: ρ boundary conditions are those of tab. 1. The mixed Hermite spline on each piece is defined in app. E. The Magnus expansion uses a fourth-order Zannah quadrature app. C.

Results

Author-reported. Setup: a MATLAB model of a 50 cm soft arm, run on a MacBook Pro §5.2. Baselines are the classic Cosserat model and the extended model without transverse muscle contraction §5.2.

The paper describes the reaching result as consistent with biological profiles of a real octopus arm, citing a published study §5.2. Treat this as the authors' interpretation, not a validation against measured data from this paper.

Paper vs code

Code: MATLAB repository yuchensun97/sorosim_v7.0 at code_commit, MIT licence (LICENSE), linked from the paper's footnote. It is modified from the SoRoSim toolbox README.md. Checked in the code: material and geometry constants in MEG.m, the Lamé relation, added mass, drag, the octopus activation functions (reaching and fetching), and the basis orders in the configs.

Differences:

• Poisson ratio in the configs. The paper uses ν0 = 0.4999 tab. 2; §5.1. Both shipped configs use 0.499 Octopus.json Stiffness.json. Lamé λ = Eν/((1+ν)(1−2ν)) SorosimLink.m scales with 1/(1−2ν), so λ is about ten times smaller in the code than with the paper's value. I did not check which value the published figures used.
• Stiffness-test diameter. §5.1 gives a 1.5 cm diameter. The config sets base and tip radius to 0.015 m, which is 1.5 cm as a radius Stiffness.json. The paper's value and the config differ by a factor of two in diameter; which one generated the figure is not determinable from the code.
• Stiffness-test viscosity. The config sets η = 300 Pa·s Stiffness.json. The paper gives no viscosity for §5.1.
• Torsional damping may lack η. §2.2 gives the torsional damping as η I33 §2.2, after eq. 18, B_θ. The code writes the torsional entry of the damping matrix as the polar moment alone, with no η MEG.m — the first entry of Gs is Ix_p, while the other entries carry 3*Eta or Eta. The inflation damping matrix uses the same entry findD_rho.m, where App. D gives η I33 app. D. I did not find a place where this entry is multiplied by η. If none exists, the torsional damping in the code is smaller than the paper's by a factor of η = 120 Pa·s. Verify before relying on it.
• Extra boundary term in K_ρ. The code subtracts σ times (φ(end)φ′(end) − φ(1)φ′(1)) findK_rho_part.m. The σ is the value left over from the loop at …findK_rho_part.m, applied at both ends. The comment on that line reads "Todo: check this". App. D has no such term app. D.
• Table 3, reaching LM 1. The printed coefficient 0.072 does not match the linear term and the floor. The code uses (Fmax − Fmin)/Tp² = 0.18/6.25 = 0.0288 with Tp = 2.5 s main.m LMrelease.m. With 0.072 the parabola has its minimum 0.128 at t = 1 s, so the 0.02 floor never applies. The code's form is continuous at 2.5 s.
• Table 4, fetching LM 2–4. The printed amplitude is min(−0.048 t + 0.04, 0.16), which decreases. The code uses +0.048 t + 0.04, rising to the 0.16 cap at t = 2.5 s fetch.m LM.m. The cap is reached exactly at Tp, so I take the code to be the intended form. The authors have not confirmed this.
• Gravity. The reaching run sets Gravity = false, while §3.3 describes gravity and buoyancy main.m. The paper does not say whether §5 runs include gravity.
• Time step and end time. Not stated in the paper. The reaching run uses dt = 0.01 s and stops at 3.5 s main.m. The fetching run stops at 9 s fetch.m.
• Transverse pressure units. Table 3 and Table 4 print 800 and 1600 without units. The code states Pa main.m fetch.m.
• Boundary condition of ρ for the octopus. §4.2 defines a mixed Hermite spline. §5.2 does not state the boundary condition used. The octopus config uses neumann Octopus.json.
• Not checked. The stiffness-test driver (Custom/Stiffening/axialStiffness.m), the derivative routine, and the stiffness and damping functions for ρ (EG_rho_bar, EG_xi_bar) were not read. The 0–5 N and 0–8 kPa ranges are unverified.

No discrepancy found for: ρ0, E, η, ρ_a, B1, B2, C_L, C_D, the Legendre and Hermite orders in the octopus config, the reaching sigmoid form and the LM 2–4 and TM functions of tab. 3, the fetching TM and OM functions of tab. 4, and the inertia I_αα used in the inflation mass matrix findM_rho.m.

Pitfalls

• Equation numbers. Counted in order from the LaTeX source; the compiled PDF was not checked. If a number looks wrong, check it against the PDF; the section numbers are safe.
• Notation clash. E is both Young's modulus and the Green-Lagrange strain tensor; ν is both the strain vector and (as ν0) Poisson's ratio; ρ is the inflation ratio, ρ0 the material density and ρ_a the water density. Read the subscript and context.
• Equations 53–56. The printed coupling for the inflation equation applies K_ρ to q_ξ eq. 54, eq. 56. K_ρ is square in q_ρ, which the code's ndof_rho shape confirms findK_rho_part.m. Implement K_ρ q_ρ. App. D also writes K_ρ as a function of (q_ξ, q̇_ξ), which is inconsistent with its own use.
• Near-incompressibility is numerically delicate. With E = 2 kPa and ν0 = 0.4999, λ = Eν/((1+ν)(1−2ν)) is about 3.3 MPa while μ = E/(2(1+ν)) is about 0.67 kPa. The ratio is about 5000 (derived from tab. 2 with the standard Lamé relations; the paper does not state λ). Compute λ from E and ν in double precision with 1−2ν kept exact. Using 0.499 instead of 0.4999 lowers λ tenfold, and 0.5 makes it infinite.
• Equation 7 is an approximation. It keeps only quadratic terms in ν, κ and ρ−1, and the paper calls it ad hoc §2.2, text after eq. 7. The large-strain behaviour is not captured. The paper lists a neo-Hookean model as future work §6.
• Transverse traction assumes homogeneous pressure. r = 2 f0 π z² eq. 29 holds only for a uniform pressure on the boundary and for a circular section.
• Basis choices are case-specific. The stiffness test keeps only stretch and inflation; the octopus runs keep only planar bending κ2, stretch and inflation §5.1; §5.2. Extend the basis set for 3D motion, as the fetching run does §5.3.
• Transverse amplitudes have no printed units. The reaching and fetching TM amplitudes in tab. 3 and tab. 4 carry no unit. Use the code's Pa values as the reference main.m.
• Cost. The paper reports one real-time factor, 0.45, on a laptop §5.2. It gives no timing for other basis sizes or quadrature counts, which set the cost app. E.
• Solver tolerances. Not stated. The paper uses ode15s for dynamics and fsolve for statics app. E.

Lineage

• Builds on: the Geometric Variable Strain approach, cited as renda_geometric_2020 §1; §4; classic Cosserat dynamics, cited as renda_dynamic_2014 §1; Saint-Venant Kirchhoff and Kelvin-Voigt laws, cited as tunay_spatial_2013 and linn_geometrically_2013 §2.2. None of these has a card in this wiki.
• Compared against: the classic Cosserat rod without inflation §2.2, eq. 19–20; §5.2; the extended model without transverse-muscle contraction §5.2; the quadratic-constitutive variant of [angles_viper_2019] §1; static inflation models kumar_geometrically_2011 and tunay_spatial_2013 §1.
• Related card: Soft Robots Modeling: a Structured Overview covers Cosserat and finite-element models of soft robots.
• Superseded by: not known.