MAT-DiSMech: A Discrete Differential Geometry-Based Computational Tool for Simulation of Rods,Shells,and Soft Robots

Abstract

Accurate and efficient simulation tools are essential in robotics, enabling the visualization of system dynamics and the validation of control laws before committing resources to physical experimentation. Developing physically accurate simulation tools is particularly challenging in soft robotics, largely due to the prevalence of geometrically nonlinear deformation. A variety of robot simulators tackle this challenge by using simplified modeling techniques—such as lumped mass models—which lead to physical inaccuracies in real-world applications. In contrast, high-fidelity simulation methods for soft structures, such as finite element analysis, offer increased accuracy but lead to higher computational costs. In light of this, we present a discrete differential geometry-based simulator that provides a balance between physical accuracy and computational efficiency. Building on an extensive body of research on rod and shell-based representations of soft robots, our tool provides a pathway to accurately model soft robots in a computationally tractable manner. Our open-source MATLAB-based framework is capable of simulating the deformations of rods, shells, and their combinations, primarily utilizing implicit integration techniques. The software design is modular for the user to customize the code—for example, add new external forces and impose boundary conditions—to suit their requirements. The implementations for prevalent forces encountered in robotics, including gravity, contact, kinetic and viscous damping, and hydrodynamic and aerodynamic drag, have been provided. We provide several illustrative examples that showcase the capabilities and validate the physical accuracy of the simulator. The open-source code is available at https://github.com/StructuresComp/dismech-matlab. We anticipate that the proposed simulator can serve as an effective digital twin tool, enhancing the Sim2Real pathway in soft robotics research.

Keywords

soft robotics discrete differential geometry rods shells

Introduction

Soft robots are robotic systems designed to accomplish desired objectives through the controlled deformation of their flexible structures.¹ These robots are largely inspired by biological creatures who appear to effectively leverage their bodily softness for efficient movement in complex terrains. Soft robots aim to replicate the adaptability, agility, and resilience found in nature. In addition to this, such robots also lead to enhanced safety in human–robot interaction owing to their inherent compliance. To fully harness the potential of soft robots, robust and accurate computational models are essential for understanding and predicting their complex behaviors.

Modeling enables simulation, allowing researchers to create a digital twin for validating designs and concepts in a safe and controlled virtual environment before extensive physical experimentation.² Moreover, simulations offer the ability to explore physical parameters beyond the limitations of current hardware, such as miniaturizing designs—a capability particularly valuable in biomedical applications.³

However, accurately modeling the dynamics of soft robots—which, in principle, possess infinite degrees of freedom (DOF)—presents significant challenges⁴ due to their highly nonlinear behavior. Unlike rigid body dynamics, where nonlinearity arises primarily from rotational kinematics while the body itself remains undeformed, soft body dynamics exhibit nonlinearities stemming from both rotations and large deformations of the continuum. These nonlinearities originate from two sources: geometry and material properties. Geometric nonlinearity arises when slender structures undergo significant bending, twisting, or stretching, leading to a nonlinear relationship between displacements and strains. In contrast, material nonlinearity refers to the nonlinear relationship between stress and strain inherent to the material’s constitutive behavior. In this article, we focus on addressing geometric nonlinearity.

A spectrum of methods ranging from simplistic lumped-mass approaches to high-fidelity tools such as finite element analysis (FEA) have been employed to tackle this challenge.⁵ These techniques exhibit a visible trade-off between accuracy and efficiency, that is, more physically accurate models such as FEA tend to lack computational efficiency, and simpler models tend to compromise on physical accuracy.⁶

A relatively newer approach in this area has been the use of discrete differential geometry (DDG), which is the study of discrete counterparts to smooth geometric surfaces, focusing on the discrete versions of concepts such as curvature and geodesics.⁷ This mathematical discipline has found extensive applications in computer graphics research for visually realistic simulation of hair and clothes. Unlike FEA, in DDG-based methods, the deformation is modeled by nonlinear ordinary differential equations (ODEs) instead of partial differential equations in weak form. These methods use simpler discretized elements, for example, when modeling rods, the individual discrete segments do not undergo bending deformation, effectively reducing the degrees DOF of each segment. This reduction in complexity, along with a focus on local interactions between discrete segments, makes it computationally more efficient than FEA, which requires solving large global stiffness matrices.

Alongside DDG and classical FEA, the absolute nodal coordinate formulation (ANCF) has also been widely explored for simulating structures undergoing large geometric deformation.⁸ ANCF is a finite element-based approach that represents nodal configuration using both position and slope vectors, allowing large rotations to be captured implicitly without introducing angular or director-based variables.⁹ This formulation offers high physical fidelity for continuum modeling, especially in scenarios involving extreme bending or twisting. At the same time, the expanded nodal representation leads to more DOF per node and, consequently, higher computational cost, especially for slender structures. As a result, ANCF and DDG occupy complementary positions in the soft-robotics modeling landscape: ANCF prioritizes generality and continuum accuracy, while DDG offers greater scalability through computational efficiency and geometric exactness for rod and shell structures.

DDG-based approaches have been successfully applied to the mechanics of slender structures, demonstrating remarkable physical accuracy across a variety of scenarios.^10,11 Building on this foundation, Huang et al.¹² introduced DDG-based two-dimensional (2D) beam elements to capture the planar deformations of soft robotic structures. Extending these ideas to fully three-dimensional (3D) settings offers a promising direction for advancing soft robot modeling and motivates the DDG-based rod and shell formulations employed in this work.

A wide variety of soft robotic systems, including pneumatic network (PneuNet) actuators,^13–15 mobile robots,^16–18 and soft robotic grippers,¹⁹ can be modeled as a combination of one or more slender rod-like structures. The discrete elastic rod (DER) algorithm proposed by Bergou et al.²⁰ is a method that has proven to be a promising tool for physically accurate and fast simulations of elastic rods. This technique works especially well for efficient soft contact handling even in complex scenarios such as knot tying,^21,22 hence making it a great choice for rod-like soft robot simulations.

In addition to slender rods, thin plate and shell structures are highly prevalent in soft robotics, and a large variety of soft robots can be modeled as a combination of rods and shells.²³ Modeling soft shells has garnered the interest of the computer graphics community, especially due to its application to realistic cloth simulations. The initial attempts to model deformation in thin plates used the notion of bending at hinges, which are edges shared by two mesh triangles.²⁴ This technique was modified by Grinspun et al.²⁵ to inculcate 3D shells, that is, shells with nonzero natural curvature. While highly effective for generating visually realistic simulations in computer graphics, this method exhibits a significant dependence on the mesh’s shape and orientation,²⁶ which hinders its ability to achieve the physical accuracy essential in robotics. Grinspun et al.²⁶ and Chen et al.²⁷ presented advanced modeling techniques for elastic shells that remedy this mesh-dependent behavior, making them a promising tool for simulating soft structures in robotics.

Our work is motivated by the potential of DDG-based techniques to model flexible body dynamics, along with the growing demand for effective simulators in soft robotics, as evidenced by the recent development of various tools. The widely used robot simulation environments such as Gazebo²⁸ and Bullet/PyBullet²⁹ employ rigid body dynamics to model robot structures and, hence, are not suitable for soft body simulation. SoMo,³⁰ a wrapper around PyBullet, facilitates the simulation of soft links by approximating them as rigid links connected with springs. While computationally efficient, such methods fail to deliver adequate physical accuracy especially when dealing with complex soft structures. Another class of prominent simulation tools includes SOFA,³¹ utilized in the widely recognized Soft Robotics Toolkit developed by Holland et al.,³² and CHRONO,³³ which use FEA to model soft body dynamics. While accurate in dynamic modeling, these simulators come at the cost of a higher computational expense. Yet, another class of more recent soft robot simulators includes Elastica³⁴ and SoroSim,³⁵ which use geometric modeling techniques such as Cosserat rod theory to model soft rod dynamics. Although these simulators achieve a good balance of accuracy and efficiency, their applicability is limited to modeling rod-like flexible structures. Furthermore, Elastica employs an explicit time integration scheme, which requires very small timesteps to converge, especially when solving stiff systems of ODEs. This highlights the need for a simulation framework that combines accuracy, efficiency, and versatility to accommodate a broader range of soft robotic systems.

In our previous work, we introduced a framework based on DDG for simulating multiple connected rodlike structures.³⁶ The simulator coded in C++ leveraged the DER method to model rods and used an implicit integration scheme, showing an improvement in computational efficiency in comparison to the existing tools. Building on this foundation, we now present an enhanced generalizable simulation framework, which extends support to a wide range of soft robots composed of slender rods and thin shells. This new tool is developed fully in MATLAB and enables the simulation of a variety of soft robot configurations, such as one or more rods, shells, and their combinations, making it suitable for complex systems not easily addressed by existing methods. The MATLAB implementation, while somewhat slower than compiled languages, offers significant advantages in accessibility and usability by eliminating dependencies, reducing code complexity, and simplifying installation. Its open-source and dependency-free design makes it an excellent pedagogical tool for teaching and learning simulation methods. For researchers seeking greater computational efficiency, the code provides a clear and modular structure that can be adapted and reimplemented in a compiled language or integrated into existing architectures. Figure 1 shows the key features of MAT-DiSMech. The main contributions of our article are as follows:

FIG. 1.

Examples of robots or soft structures: snake, manta ray, and parachute, modeled using MAT-DiSMech. The key features within MAT-DiSMech. (a) Elasticity: The framework accounts for stretching, bending, and twisting deformation modes for elastic rods and stretching, bending deformations for shells. (b) External forces: Common external forces encountered in soft robotics include ① buoyancy, ② hydrodynamic drag, ③ gravity, ④ contact, and ⑤ friction, illustrated here for an octopus-inspired soft robot locomoting underwater. These forces, among others, are implemented in MAT-DiSMech, and the framework further allows users to define custom external force models. (c) Self-contact: Self-contact is modeled using IMC, which treats contact implicitly by defining a penalty energy as a function of $Δ$ , the minimum distance between two edges that can potentially come in contact. (d) Actuation: Actuation is primarily done by modifying the natural curvature or twist of the elastic elements. The parameters, $\bar{κ}, \bar{τ}$ , and $\bar{ϕ}$ , denote the natural curvature of the rod, the natural twist of the rod, and the natural hinge angle of the shell, respectively. IMC, implicit contact method.

To our knowledge, MAT-DiSMech is the first DDG-based robot simulator capable of handling elastic rods, shells, and their combinations.

We provide an open-source implementation of the discrete elastic shell bending energy that uses the discrete shape operator proposed by Grinspun et al.²⁶ that helps achieve higher mesh independence than hinge-based shell models.

We provide an open-source MATLAB implementation of implicit contact method (IMC),³⁷ an implicit contact modeling technique. This fully implicit method has proven highly effective in simulating elastic contact in contact-rich scenarios such as knot tying and flagella bundling.^22,37

Our article is structured as follows: The section “Dynamical Modeling” presents the kinematics and the equations of motion (EOMs). The section “Elastic energy” describes the formulation of the elastic energy of the entire structure. The “Self-Contact and Friction” section covers the methods to model self-contact and Coulomb friction implicitly. The “External Forces” section elaborates on the modeling of commonly encountered forces in robotics, in particular gravity, contact, coulomb and viscous friction, and aerodynamic drag. In the “Overall Framework” section, the overall framework of the simulator is described, and the procedure to go from a robot to its digital counterpart by providing appropriate inputs to the simulator is detailed. Finally, the section “Results and Discussion” showcases simulations of various soft robots and structures and results from physical validation of our solver.

Dynamical Modeling

A structure is discretized into $N$ nodes, which can represent a rod, a shell, or their connection. These nodes are used to define $E$ edges that form slender rods and $T$ triangles that constitute the shell structure. Each edge is a vector connecting two nodes, and each triangle is defined by three nodes connected together. Figure 2a shows a schematic of the discrete representation of a structure with $N = 16$ , $E = 7$ , and $T = 10$ . The user defines the geometry of the structure by providing the following information: (i) the nodal position coordinates in 3D in the form of an $N \times 3$ array; (ii) the array of size ( $E \times 2$ ) of node indices $m$ and $n$ ( $1 \leq m, n \leq N$ ) for all edges, $e^{i} = x_{n} - x_{m}$ , where $1 \leq i \leq E$ , and (iii) the array of size ( $T \times 3$ ) of node indices $l, m,$ and $n$ ( $1 \leq l, m, n \leq N$ ), for each of the $T$ triangles.

FIG. 2.

Schematics of physical modeling are shown here. (a) Discrete representation of a structure with both rod and shell elements. For this structure, the number of nodes, $N = 16$ , number of edges, $E = 7,$ and the number of triangles, $T = 10$ . (b) Stencil of the bending–twisting spring for an elastic rod, comprising of nodes $x_{m}$ , $x_{n}$ , and $x_{o}$ and edges $e^{i}$ and $e^{j}$ with attached reference frame ( $a_{1}, a_{2}$ ) and material frame ( $m_{1}, m_{2}$ ), used for the DER algorithm. (c) Schematic of a hinge spring for the hinge-based bending model for a discrete elastic shell. Hinge angle $ϕ$ is used to calculate the bending energy. (d) The discretization stencil for the mid-edge normal-based bending energy model. For each shell-edge, $ξ^{i}$ is a scalar DOF that represents the rotation of the mid-edge normal about the edge. (e) Definition of the mid-edge normal. The blue curve labeled “surface” denotes the actual surface of the shell being modeled, and “mesh edge” $e$ denotes the edge of a triangle in the mesh that approximates the surface. The mid-edge normal $n^{m}$ for $e$ is normal to the surface, which, when extrapolated, intersects the triangle edge at its midpoint. (f) The schematic showcasing the edge attached reference frame ${n^{avg, i}, τ^{i}, {\hat{e}}^{i}$ } and other vectors used in the mid-edge normal bending method. (g) An elastic joint between two rods, and a rod and a shell. The node $x_{m}$ shared between the rod-edge $e^{k}$ and the shell triangle is denoted as a “joint node,” and edges in blue color have the properties of both rod-edge and shell-edge. Bending–twisting springs are depicted by the spirals. DER, discrete elastic rod; DOF, degree of freedom.

DOF and kinematics

To model the dynamics of an elastic rod or a network of rods, we use the DER method,²⁰ which builds upon the Kirchhoff rod theory.³⁸ In this framework, rods can undergo three primary deformation modes: stretching, bending, and twisting. While the original Kirchhoff rod model describes inextensible rods, the DER method extends this by introducing a stretching energy, allowing the rods to deform axially. If strict inextensibility is required, it is also possible to impose edge-length preserving constraints, for instance, using the fast projection method.³⁹

The DOF vector for rods comprises the nodal positions, denoted by x, and the twist angles of the edges, denoted by θ, which capture their rotational state. Unlike rods, shells can undergo only two types of deformation: stretching or membrane deformation and bending. We implement two distinct discrete elastic shell theories to model shell bending, allowing the user to select their preferred method: (1) hinge-based²⁵ and (2) mid-edge normal-based.²⁶ For hinge-based bending, the DOF vector for the shell only consists of the nodal positions. Therefore, for a structure with $N$ nodes and $E$ edges, the DOF vector is of size $3 N + E$ and is expressed as

\begin{matrix} q = {[x_{1}, x_{2}, \dots, x_{N}, θ^{1}, θ^{2}, \dots, θ^{E}]}^{⊤}, \end{matrix}

where superscript

^{⊤}

represents the transpose operator. When using the mid-edge normal-based theory, it is necessary to track the edges on the shells—similar to the edges on the rods—by introducing an additional scalar DOF, denoted by

ξ

, for each of the shell-edges. The definition of

ξ

is explained in the subsection “Mid-edge normal-based bending energy”. Given a shell with

T \times 3

triangle array, the shell-edges can be obtained from the triangle array. Hence, if using the mid-edge normal-based formulation, assuming

Z

to be the number of shell-edges obtained as above, the DOF vector is of size

3 N + E + Z

and is expressed as

\begin{matrix} q = {[x_{1}, x_{2}, \dots, x_{N}, θ^{1}, θ^{2}, \dots, θ^{E}, ξ^{1}, ξ^{2}, \dots, ξ^{Z}]}^{⊤} . \end{matrix}

Equations of Motion

The EOMs of the system are

\begin{matrix} M \ddot{q} - F^{elastic} - F^{ext} - F^{IMC} = 0, \end{matrix}

where

M

is a square diagonal lumped mass matrix of size

3 N + E + Z

, dot ( ˙) represents derivative with respect to time, the vector

F^{elastic} = - \frac{\partial E^{elastic}}{\partial q}

represents the elastic force (where

E^{elastic}

is the elastic energy),

F^{ext}

is the external force vector (such as gravity), and

F^{IMC}

represents the contact and friction forces following the IMC.^22,37

The EOMs are solved numerically using the implicit Euler method. To step from time $t = t_{k}$ to $t = t_{k + 1}$ with a time step size of $Δ t = t_{k + 1} - t_{k}$ , we rewrite the EOMs as a system of first-order equations by introducing the velocity $u = \dot{q}$ . The implicit Euler method then updates the state variables (positions and velocities) according to

f \equiv M \frac{1}{Δ t} (\frac{q_{k + 1} - q_{k}}{Δ t} - u_{k}) + \frac{\partial E^{elastic}}{\partial q_{k+1}} - F_{k + 1}^{ext} - F_{k + 1}^{IMC} = 0,

\begin{matrix} u_{k + 1} = \frac{q_{k + 1} - q_{k}}{Δ t} . \end{matrix}

where the subscript

k

and

k + 1

denotes the evaluation of a quantity at time

t_{k}

and

t_{k + 1}

. Note that these equations are implicit in nature because the forces

F_{k + 1}^{elastic}

F_{k + 1}^{ext}

, and

F_{k + 1}^{IMC}

depend on the unknown state at

t_{k + 1}

To solve Eq. 4, we employ the Newton–Raphson method and iteratively solve for $q_{k + 1}$ . A line search method, similar to the one used in Tong el al.,³⁷ is used to decide the step magnitude $α$ in the direction of the gradient to help aid convergence, especially in the case of contact-rich scenarios. The running value of $q_{k + 1}$ is updated at every iteration as follows:

\begin{matrix} \begin{array}{r} q_{k + 1}^{(i + 1)} = q_{k + 1}^{(i)} - α linearSolve (J, f), \end{array} \end{matrix}

where

i

is the iteration index,

linearSolve (J, f)

solves the linear system

J Δ q = f

for

Δ q

f

is the residual of EOMs evaluated at

q_{k + 1}^{(i)}

from Eq. (4), and

J

is the Jacobian matrix of the EOMs defined as

\begin{matrix} \begin{array}{r} J = M \frac{1}{Δ t^{2}} + \frac{\partial^{2} E^{elastic}}{\partial q_{k + 1} \partial q_{k + 1}} - \frac{\partial F^{ext}}{\partial q_{k + 1}} - \frac{\partial F^{IMC}}{\partial q_{k + 1}} . \end{array} \end{matrix}

where

M \frac{1}{Δ t^{2}}

represents the inertial contribution and is straightforward to compute;

\frac{\partial^{2} E^{elastic}}{\partial q_{k + 1} \partial q_{k + 1}}

is the Hessian of the elastic energy;

\frac{\partial F^{IMC}}{\partial q_{k + 1}}

accounts for the contribution of contact and friction forces; and

\frac{\partial F^{ext}}{\partial q_{k + 1}}

is the gradient of the external force. The last term can be ignored if their analytical expressions are not available, although doing so may reduce the convergence rate of the Newton–Raphson method. Once the positions

q_{k + 1}

converge within a specified tolerance, the velocities

u_{k + 1}

are updated using Eq. 5. The forces in Eq. 4—arising from elasticity, contact/friction, and external effects—and their gradients with respect to the DOFs in Eq. 7 are described in the subsequent sections.

Elastic Energy

MAT-DiSMech decomposes the structure into a network of spring-analog elements, where each spring is associated with an elastic energy. These spring-like elements allow us to compute the total elastic energy of the structure by summing the contributions from all individual springs. A structure can have four types of springs, namely stretching spring, bending–twisting spring, hinge spring, and triangle. The total elastic energy for S stretching springs, B bending–twisting springs, H hinge springs, and T triangles is

\begin{matrix} \begin{array}{r} E^{elastic} = \sum_{i = 1}^{S} E_{i}^{stretch} + \sum_{i = 1}^{B} (E_{i}^{bend} + E_{i}^{twist}) + E^{shell}, \end{array} \end{matrix}

where

E^{shell} = \{\begin{matrix} \sum_{i = 1}^{H} E_{i}^{hinge} if hinge bending, \\ \sum_{i = 1}^{T} E_{i}^{mid-edge} if mid - edge bending . \end{matrix}

Stretching springs ( $E_{i}^{stretch}$ ) can represent deformation in either rods or shells, while bending–twisting springs ( $E_{i}^{bend}$ and $E_{i}^{twist}$ ) are exclusive to rods, capturing their bending and twisting behavior. Hinge springs and triangles are specific to shells, with the choice of hinge-based ( $E_{i}^{hinge}$ ) or mid-edge normal-based ( $E_{i}^{mid-edge}$ ) models determining the method used to model shell bending.

Stretching energy

Stretching energy is stored in stretching springs associated with edges in rods and shells. For example, in Figure 2b, the $i$ -th stretching spring lies between nodes $x_{m}$ and $x_{n}$ on the edge $e^{i} = x_{n} - x_{m}$ . The axial strain in the $i$ -th edge is given by

\begin{matrix} \begin{array}{r} ϵ^{i} = \frac{∥ e^{i} ∥}{∥ e^{¯^{i}} ∥} - 1, \end{array} \end{matrix}

where overbar represents a quantity in underformed or stress-free state, and

∥ e ¯^{i} ∥

is the undeformed length of the

i

-th edge. The stretching energy along edge

e^{i}

\begin{matrix} \begin{array}{r} E_{i}^{stretch} = \frac{1}{2} k_{s} {(ϵ^{i})}^{2} ∥ e ¯^{i} ∥, \end{array} \end{matrix}

where

k_{s}

is the stretching stiffness. For edges on rods,

k_{s} = E A

, where

E

is the Young’s modulus, and

A

is the area of cross-section. The default value for

A

for rod-edges used in the simulator is

π r_{0}^{2}

assuming a circular cross-section of radius

r_{0}

. The default value of

k_{s}

for shell-edges is

\frac{\sqrt{3}}{4} E h ||e ¯^{i}||

, where

E

is the Young’s modulus, and

h

is the thickness of the shell.⁴⁰ If desired, users can modify the stiffness values of individual stretching springs in the simulator. Furthermore, the stiffness of each spring does not have to be the same.

Bending and twisting energy

Bending and twisting energy for slender rods is stored in bending–twisting springs. A bending–twisting spring comprises three nodes—the center node at which the bending and twisting strains are measured and its two adjacent nodes, and two rod-edges—the edge entering the center node and the edge exiting the center node. In Figure 2b, three nodes (x_m, x_n, and x _o ) and two edges (e ⁱ and e ^j ) form a bending–twisting spring. The edge e ⁱ (and e ^j ) is the vector from x _m to x _n (and from x _n to x _o ). Each rod-edge $e^{i}$ has two sets of orthonormal frames, a reference frame ${d_{1}^{i}, d_{2}^{i}, t^{i}}$ and a material frame ${m_{1}^{i}, m_{2}^{i}, t^{i}}$ . Both of these frames share the tangent vector along the edge $t^{i} = e^{i} / ∥ e^{i} ∥$ as one of the directors. The reference frame is initialized at time $t = 0$ and updated at each time step by parallel transporting it from the previous time step. The use of so-called time-parallel transport is a key feature of DER, leading to high computational efficiency.⁴¹ The material frame for the edge can be obtained from the reference frame by applying the twist angle for the edge, $θ^{i}$ , about their shared director $t^{i}$ along the edge.

Bending strain is measured at the middle node $x_{n}$ through the curvature binormal vector

{(κ b)}_{k} = \frac{2 e^{i} \times e^{j}}{∥ e^{i} ∥ ∥ e^{j} ∥ + e^{i} \cdot e^{j}} .

The scalar curvatures along the first and second material directors using the curvature binormal are

\begin{matrix} κ_{k}^{(1)} = \frac{1}{2 (m_{2}^{i} + m_{2}^{j})} \cdot {(κ b)}_{k} . \end{matrix}

\begin{matrix} κ_{k}^{(2)} = - \frac{1}{2 (m_{1}^{i} + m_{1}^{j})} \cdot {(κ b)}_{k} . \end{matrix}

The bending energy for the $k$ -th bending–twisting spring is

\begin{matrix} \begin{array}{r} E_{k}^{bend} = \frac{1}{2} \frac{E I}{Δ l_{k}} [{(κ_{k}^{(1)} - {κ ¯}_{k}^{(1)})}^{2} + {(κ_{k}^{(2)} - {κ ¯}_{k}^{(2)})}^{2}], \end{array} \end{matrix}

where

Δ l_{k} = \frac{1}{2} (∥ e ¯^{i} ∥ + ∥ e ¯^{j} ∥)

is the Voronoi length for

x_{k}

{κ ¯}_{k}^{(1)}

and

{κ ¯}_{k}^{(2)}

are the natural curvatures in undeformed state, and

E I

is the bending stiffness. The default value of

I

π r_{0}^{4} / 4

, assuming a circular cross-section of radius

r_{0}

The twist between the two edges $i$ and $j$ corresponding to the $k$ -th bending–twisting spring is

\begin{matrix} \begin{array}{r} τ_{k} = θ^{j} - θ^{i} + Δ m_{k, ref}, \end{array} \end{matrix}

where

Δ m_{k, ref}

is the reference twist, which is the twist of the reference frame as it moves from the

i

-th edge to the

j

-th edge. The twisting energy for the

k

-th bending–twisting spring is

\begin{matrix} \begin{array}{r} E_{k}^{twist} = \frac{1}{2} \frac{G J}{Δ l_{k}} (τ_{k} - {\bar{τ}}_{k})^{2}, \end{array} \end{matrix}

where

{\bar{τ}}_{k}

is the twist along the centerline in the undeformed state,

Δ l_{k}

is the Voronoi length for the center node,

G J

is the twisting stiffness,

G

is the Shear modulus of the rod’s material given by

G = E / (2 (1 + ν))

ν

is the Poisson’s ratio of the rod’s material, and

J

is the polar moment of inertia. The default value of

J

π r_{0}^{2} / 2

, assuming a circular cross-section. As in the case of stretching springs, if desired, users can modify the bending or twisting stiffness values of individual bending–twisting springs in the simulator.

An edge may be shared by multiple bending–twisting springs, especially while modeling a network of rods. The above formulation assumes a specific orientation for the edges ( $x_{m} \to x_{n}$ for the first edge and $x_{n} \to x_{o}$ for the second edge). As shown in Figure 2g, for bending–twisting spring ${x_{a}, x_{b}, x_{c} and e^{i}, e^{j}}$ , since both rod-edges $e^{i}$ and $e^{j}$ are pointing toward the node $x_{b}$ , the negative of $e^{j}$ is used for calculations so that it points away from and $e^{i}$ points toward the node $x_{b}$ . The reference frame vector $d_{1}^{j}$ , material frame vector $m_{1}^{j}$ , and edge rotation angle $θ^{j}$ are multiplied by $- 1$ for computing the bending and twisting forces. Additionally, the resulting force vector and the Jacobian are reoriented to their original states by multiplying the affected terms by $- 1$ if the edge and its related quantities were adjusted during the computation. Note that we could have equivalently flipped the edge $e^{i}$ instead of $e^{j}$ and proceeded as above, and the result would be mathematically equivalent. In the simulator, bending–twisting springs are initialized during the problem setup, and information about whether an edge needs to be negated is stored at the beginning of the simulation. This straightforward bookkeeping approach enables efficient handling of rod-like network structures.

Shell bending energy

Two different methods are available in MAT-DiSMech to formulate the bending energy for an elastic shell. The hinge-based formulation,²⁵ while simple, leads to mesh-dependent behavior, meaning that the mesh orientation significantly influences the deflection. In contrast, the mid-edge normal-based method²⁶ yields more mesh-independent behavior, making it more appropriate for applications that require higher accuracy.

Hinge-based bending energy

This method was initially proposed by Baraff et al.²⁴ to simulate plates and was further modified by Grinspun et al.²⁵ to include shells with nonzero natural curvatures. Referring to Figure 2c, a hinge is defined as any shell-edge that is shared by two adjacent shell triangles about which the two faces can rotate relative to each other. The hinge spring comprises four nodes ( $x_{l}, x_{m}, x_{n},$ and $x_{o}$ ), and two of these nodes define the hinge, which in this case is the edge vector from $x_{m}$ to $x_{n}$ . The hinge angle $ϕ$ is defined as the angle between the normal vectors on these two shell triangles. The bending energy corresponding to the $i$ -th hinge is

\begin{matrix} E_{i}^{hinge} = \frac{1}{2} k_{b} {(ϕ_{i} - \bar{ϕ_{i}})}^{2}, \end{matrix}

where

k_{b} = \frac{1}{\sqrt{3}} \frac{E h^{3}}{12}

is the bending stiffness,

E

is the Young’s modulus, and

h

is the thickness of the shell.

Mid-edge normal-based bending energy

In this method, the notion of the shape operator is used to describe curvature. The shape operator, denoted as $Λ_{p}$ , at a point $p$ on a smooth surface, is a linear map that relates the tangent vector $t$ at $p$ to the rate of change of the surface normal vector $n$ at $p$ in the direction of $t$ (see Grinspun et al.²⁶). Mathematically, this relationship is expressed as

\begin{matrix} Λ_{p} t = - \partial_{t} n, \end{matrix}

where

\partial_{t} n

denotes the directional derivative of

n

with respect to

t

that quantifies how the surface normal

n

changes in the direction of

t

. The eigenvectors and eigenvalues of the shape operator at each point determine the directions in which the surface bends at that point.

To compute the discrete shape operator for a meshed surface, the concept of mid-edge normal is introduced. As shown in Figure 2e, for a particular $e^{i}$ on the mesh, the mid-edge normal $n^{m, i}$ for that edge is defined as the normal to the smooth surface (approximated by the mesh), which, when extrapolated, intersects the edge at its mid-point. In Figure 2f, a reference frame { $n^{avg, i}, τ^{i}, {\hat{e}}^{i}$ } is attached to edge $e^{i}$ . Here, $n^{avg, i}$ is defined as the unit vector along $\frac{n_{1} + n_{2}}{2},$ where $n_{1}$ and $n_{2}$ are the unit normal vectors of the two faces that share the edge $e^{i}$ , ${\hat{e}}^{i}$ is the unit vector along $e^{i}$ and $τ^{i} = n^{avg, i} \times {\hat{e}}^{i}$ . At the start of a time step, $τ^{i}$ is computed for each shell-edge, and this quantity is denoted as $τ^{i, 0}$ . For the $i$ -th edge, a scalar $ξ^{i}$ is defined as the projection of the mid-edge normal on $τ^{i, 0}$ , that is, $ξ^{i} = n^{m, i} \cdot τ^{i, 0}$ in Figure 2f. This scalar $ξ^{i}$ is considered as a DOF associated with each shell-edge that represents the rotation of the mid-edge normal about the edge. Then, if a triangle has three nodes with indices ${l, m, n}$ and three edges with indices ${p, q, r}$ as shown in Figure 2d, we have a DOF vector of size $(12 \times 1)$ , $q_{mid-edge} = {[x_{l}^{⊤}, x_{m}^{⊤}, x_{n}^{⊤}, ξ^{p}, ξ^{q}, ξ^{r}]}^{⊤}$ . We use the following per-triangle discrete shape operator presented by Grinspun et al. to compute the bending energy associated with the $i$ -th triangle,

Λ_{i} = \sum_{k} \frac{s^{k} ξ^{k} - (n_{i} \cdot τ^{k, 0})}{{\bar{A}}_{i} ∥ {\bar{e}}^{k} ∥ (t^{k} \cdot τ^{k, 0})} t^{k} \otimes t^{k}

where

k \in {p, q, r}

denotes the

k

-th edge of the triangle,

n_{i}

denotes the unit normal to the triangle,

{\bar{A}}_{i}

denotes the undeformed area of the triangle,

∥ e ¯^{k} ∥

denotes the undeformed length of the

k

-th edge,

t^{k}

is the tangent to the surface perpendicular to the edge (note that the

\hat{()}

operator indicates unit vector),

s^{k}

is a variable that takes values

1

- 1

depending on the ownership of the mid-edge normal, and the operator

\otimes

denotes the outer product. Figure 2e shows all the above-defined vectors. Note that the “shape operator” defined in Eq. (19) is not operating on a point. Instead, it refers to the (3 × 3) matrix that represents the linear map from the shape operator definition in Eq. (18). For the i-th triangle, the discrete shell bending energy, in terms of the discrete shape operator

Λ_{i}

, is

E_{i}^{mid-edge} = k_{b} {\bar{A}}_{i} [(1 - ν) Tr ({(Λ_{i} - {\bar{Λ}}_{i})}^{2}) + ν {(Tr (Λ_{i}) - Tr ({\bar{Λ}}_{i}))}^{2}],

where

k_{b} = \frac{E h^{3}}{24 (1 - ν^{2})}

E

is the Young’s modulus,

h

is the thickness,

ν

is the Poisson’s ratio of the shell,

{\bar{A}}_{i}

and

{\bar{Λ}}_{i}

denote the area and the shape operator of the i-th triangle in the undeformed or natural configuration, and

Tr ()

denotes the trace of a matrix. The complete analytical expressions for the gradient and Hessian of the energy in Eq. (20) are not available in previous publications but are needed in Eqs. (4) and (7); hence, we provide the complete expressions in Appendix A.

Bending and twisting at the rod-shell joint

A rod-shell joint is formed when a node is shared between an edge of a rod and one or more shell triangles, which we refer to as the “joint node.” In Figure 2f, $x_{m}$ is such a node. To model its elastic energy, all stretching, bending, and twisting deformations involving the joint node must be considered. For this purpose, we treat the edges ( $e^{p} = x_{n} - x_{m}$ , $e^{q} = x_{l} - x_{n},$ and $e^{r} = x_{m} - x_{l}$ ) that form the triangles sharing the joint node as if they are also rod-edges, augmenting them with material and reference frames and associated twist angles. As illustrated in Figure 2g, bending and twisting springs are then associated with every possible three-node, two-edge combination resulting from this assumption, that is, ${x_{b}, x_{m}, x_{n} and e^{k}, e^{p}}$ , ${x_{m}, x_{n}, x_{l} and e^{p}, e^{q}}$ , ${x_{n}, x_{l}, x_{m} and e^{q}, e^{r}}$ , ${x_{l}, x_{m}, x_{n} and e^{r}, e^{p}}$ , and ${x_{l}, x_{m}, x_{b} and e^{r}, e^{k}}$ . Consequently, if the mid-edge normal-based bending method is used, the edges $e^{p}, e^{q}$ , and $e^{r}$ will possess both $θ$ and $ξ$ DOF. The joint node, in turn, will experience elastic forces from bending and twisting of rods, bending of shells, and stretching of edges.

Self-Contact and Friction

In addition to modeling inherent elasticity, we account for contact and friction forces resulting from the structure’s self-interactions for a realistic representation of its behavior during deformations and collisions. We do so by using a fully implicit penalty energy method, IMC,^22,37 that integrates seamlessly into our framework. We divide the force into two parts: contact force to ensure nonpenetration and friction forces due to Coulomb friction.

The contact force is evaluated for each pair of nonadjacent edges that are in close proximity, denoted as a “contact pair.” Figure 1c shows a contact pair comprising two edges ${e^{i}, e^{j}}$ or four nodes ${x_{a}$ , $x_{b}$ , $x_{c}$ , $x_{d}}$ , since $e^{i} = x^{b} - x^{a}$ and $e^{j} = x^{d} - x^{c}$ , and $Δ (\geq 0)$ is the shortest distance between $e^{i}$ and $e^{j}$ measured using Lumelsky’s algorithm.⁴² Penetration occurs if $Δ < 2 r_{0}$ . To compute the contact force for each pair, IMC defines a contact energy that penalizes penetration as follows:

E^{con} = \{\begin{array}{l} (2 r_{0} - Δ)^{2} & Δ \leq 2 r_{0} - δ, \\ 0 & Δ \geq 2 r_{0} + δ, \\ {(\frac{1}{K_{1}} l og (1 + e^{K_{1} (2 r_{0} - Δ)}))}^{2} & otherwise \end{array}

where

K_{1} = 15 / δ

is a stiffness parameter and

δ

is a user-defined contact distance tolerance. Similar to elastic energy, for a contact pair, the contact force on the

i

-th node of the contact pair is derived as the gradient of the contact energy with respect to the DOF of the node as,

\begin{matrix} F_{i}^{con} = - \frac{\partial E^{con}}{\partial Δ} \frac{\partial Δ}{\partial q_{i}} . \end{matrix}

If friction in nonzero, Coulomb friction force on the i-th node of the contact pair arising due to the contact force $F_{i}^{con}$ is

\begin{matrix} \begin{matrix} F_{i}^{fr} = - μ γ \hat{u} ∥ F_{i}^{con} ∥, \\ γ = \frac{2}{1 + e^{- K_{2} ∥ u ∥}} - 1, \end{matrix} \end{matrix}

where

u

is the tangential relative velocity of the two edges in contact, which can be computed from the velocities of the nodes in the contact pair;

μ

is the friction coefficient;

K_{2} = 15 / ν

is a stiffness parameter used to calculate

γ \in [0, 1]

that is a smooth scaling factor that models the transition from sticking and sliding modes between the contacting bodies; and

ν

is a user-defined slipping tolerance, that is, the velocity beyond which the bodies are considered to undergo slipping. The total force due to self-contact and friction is obtained by summing over the contribution from all contact pairs.

F^{IMC} = \sum_{a, b, c, d \in contact pair} (\sum_{i \in a, b, c, d} F_{i}^{con} + F_{i}^{fr}) .

The expressions for the Jacobian of contact and friction forces are omitted here but can be found in Ref.³⁷ The above contact model is designed for soft contact between slender rods. However, the same code can be used by approximating shell triangles using their edges. For greater accuracy, the minimum distance (Δ) can be defined as the shortest distance between two triangles instead of two edges, which can then be used to formulate the contact energy. The extension to accurately model shell contact is left for future work.

External Forces

In addition to the elastic forces and self-contact (and friction), the algorithm models the commonly encountered environmental forces in robotics such as gravity, buoyancy, contact, coulomb and viscous friction, and aerodynamic drag. Moreover, since external forces are usually problem-dependent, the framework allows users to add custom external forces to the system as explained in the “Custom external force” subsection. The external force formulations are provided per node, and the total external force is calculated by summing over all the nodes in the structure.

Gravity and buoyancy

If the mass associated with the $i$ -th node is $M_{i}$ and $g$ is the acceleration vector due to gravity, the gravity force on that node is $M_{i} g$ . To model the effect of the buoyancy force for a body submerged in a fluid medium, we replace the value of $g$ by $\frac{ρ - ρ_{med}}{ρ} g$ , where $ρ$ is the density of the structure and $ρ_{med}$ is the density of the fluid medium.

Contact and friction with ground

When a node $x_{i}$ comes in contact with the ground surface, that is, the distance between the node and ground ( $Δ_{i}$ ) is less than some tolerance value ( $δ$ ), it will experience a repulsive force along the normal $\hat{n}$ to the ground at the point of contact. This normal force encountered and its Jacobian is calculated using the following equations, respectively.

F_{i}^{c} = k_{c} \frac{- 2 e^{- K_{1} Δ_{i}} \log (e^{- K_{1} Δ_{i}} + 1)}{K_{1} (e^{- K_{1} Δ_{i}} + 1)} \hat{n},

J_{i i}^{c} = k_{c} \frac{2 e^{- K_{1} Δ_{i}} \log (e^{- K_{1} Δ_{i}} + 1) + 2 e^{- 2 K_{1} Δ_{i}}}{(e^{- K_{1} Δ_{i}} + 1)^{2}} \hat{n} {\hat{n}}^{⊤},

where

k_{c}

is the contact stiffness, which depends on the firmness of the ground one wants to model, and

K_{1}

is a tunable parameter which we have set to

15 / δ

Coulomb friction force due to the ground is computed using the above contact force similar to friction force due to self-contact (see Eq. 23), except the velocity ( $\hat{u}$ ) in this case is the velocity of the node in contact with the ground along the tangent to the ground. By default, ground is assumed to be parallel to the $x - y$ plane; hence, $\hat{n}$ is $[0, 0, 1]^{⊤}$ . Note that the above penalty energy-based method for ground contact force is especially useful to model contact when friction is nonzero. For frictionless contact, a simpler predictor–corrector approach can be taken similar to Ref.¹¹ that ensures no penetration into the ground.

Viscous damping

Viscous friction is modeled by applying a dissipative force to each node, proportional to the nodal velocity. The viscous force on the i-th node at the n-th time instant is given by

\begin{matrix} F_{i}^{visc} = - η \cdot (\frac{x_{i} (t_{n}) - x_{i} (t_{n - 1})}{Δ}) \cdot Δ l_{i} \end{matrix}

where

η

(Pa

\cdot

s) is the viscosity coefficient of the medium, and

Δ l_{i}

is the Voronoi length of the node. The Jacobian of this force is

J_{i i}^{visc} = - \frac{η}{Δ t} \cdot Δ l_{i} \cdot I,

where

I

is the square identity matrix of size

(3 \times 3)

Hydrodynamic force using RFT

Resistive force theory (RFT)⁴³ is among the most widely used hydrodynamic models for describing interactions between low-Reynolds-number flows and slender structures. For the i-th node, the hydrodynamic force at n-th time instant is calculated by summing contributions from all edges connected to it, as follows,

F_{i}^{RFT} = \sum_{k \in E} - \frac{l^{k}}{2} [(C_{t} - C_{n}) {\hat{e}}^{k} ({\hat{e}}^{k})^{⊤} + C_{n} I] u_{i},

where

E

is the set of indices of edges that are connected to the

i

-th node,

u_{i} = (\frac{x_{i} (t_{n}) - x_{i} (t_{n - 1})}{Δ t})

is the velocity of the

i

-th node at

n

-th time instant,

Δ t

is the simulation time step,

\hat{}

denotes unit vector,

C_{t}

and

C_{n}

are the tangential and normal resistive coefficients, respectively, and

l_{k}

is the length of the

k

-th edge in undeformed state. Since each edge is shared by two nodes,

l^{k}

is divided by

2

to get the length corresponding to one of those. The negative sign denotes that the force acts in the direction opposite to the velocity component. The Jacobian of this force is

\begin{matrix} J_{i j}^{RFT} = \frac{\partial F_{i}^{RFT}}{\partial x_{j}} \\ = \sum_{k \in E} - \frac{l^{k}}{2} [(C_{t} - C_{n}) (u_{i}^{⊤} {\hat{e}}^{k}) \frac{\partial {\hat{e}}^{k}}{\partial x_{j}} + (C_{t} - C_{n}) {\hat{e}}^{k} (({\hat{e}}^{k})^{⊤} \frac{\partial u_{i}}{\partial x_{j}} + u_{i}^{⊤} \frac{\partial {\hat{e}}^{k}}{\partial x_{j}}) + C_{n} \frac{\partial u_{i}}{\partial x_{j}}], \end{matrix}

where

I

is the identity matrix of size

(3 \times 3)

. Note that

J_{i j}^{RFT}

will be nonzero only for

j = i

and the indices

j

that correspond to the nodes adjacent to the

i

-th node. The expressions for

\frac{\partial {\hat{e}}^{k}}{\partial x_{j}}

and

\frac{\partial u_{i}}{\partial x_{j}}

are provided in Appendix B.

Aerodynamic drag

Aerodynamic drag is a dissipative force proportional to the square of the nodal velocity and acts on the nodes contained in the shell structure only. The drag force on the i-th node is obtained by summing over all shell triangles that share the node and adding a force proportional to the square of the component of the node velocity ( $u_{i}$ ) along the normal of the triangle in the direction opposite to the velocity component. Thus, the drag force acting on the i-th node at the n-th time instant is given by,

F_{i}^{drag} = \sum_{k \in T} - \frac{ρ_{med} C_{D}}{2} \frac{A_{k}}{3} sgn (u_{i} \cdot {\hat{n}}_{k}) {(u_{i} \cdot {\hat{n}}_{k})}^{2} {\hat{n}}_{k}

where

ρ_{med}

denotes the density of the fluid medium in which the robot is operating,

n_{k}

denotes the normal vector of the

k

-th triangle,

T

is the set of indices of triangles which share the

i

-th node,

u_{i} = (\frac{x_{i} (t_{n}) - x_{i} (t_{n - 1})}{Δ t})

is the velocity of the

i

-th node at

n

-th time instant,

Δ t

is the simulation time step,

C_{D}

is the dimensionless drag coefficient,

A_{k}

is the area of the

k

-th triangle in its undeformed state,

sgn

denotes the signum function, and the negative sign implies that the force acts in the direction opposite to the velocity component. Since each triangle is shared by three nodes,

A_{k}

is divided by

3

to get the area corresponding to one of those. The Jacobian of this force is given by,

\begin{matrix} J_{i j}^{drag} = \frac{\partial F_{i}^{drag}}{\partial x_{j}} = \sum_{k \in T} - \frac{ρ_{med} C_{D} A_{k}}{6} s g n (u_{i} \cdot {\hat{n}}_{k}) \\ [2 {\hat{n}}_{k} (u_{i} \cdot {\hat{n}}_{k}) ({\hat{n}}_{k}^{⊤} \frac{\partial u_{i}}{\partial x_{j}} + u_{i}^{⊤} \frac{\partial {\hat{n}}_{k}}{\partial x_{j}}) + (u_{i} \cdot {\hat{n}}_{k})^{2} \frac{\partial {\hat{n}}_{k}}{\partial x_{j}}] \end{matrix}

Note that $J_{i j}^{drag}$ will be nonzero only for $j = i$ and the indices $j$ that correspond to the nodes that share a triangle with the $i$ -th node. The expressions for $\frac{\partial {\hat{n}}_{k}}{\partial x_{j}}$ and $\frac{\partial u_{i}}{\partial x_{j}}$ are provided in Appendix C.

Custom external force

If a user desires to apply an external force of their choice other than the ones implemented in our simulator, they can do so by creating a custom function to calculate the external force and its Jacobian similar to the functions for the provided external forces. The inputs to the function can be the DOF vector, the velocity vector, and any other system parameters as necessary, and the output would be the force, a vector of the same length as the DOF vector and the Jacobian, a square matrix of the same size as the DOF vector. Note that even though our method uses the force and Jacobian of the force to propagate the dynamics forward, the Jacobian is used to provide a direction of search for the solution to the equation of motion at every time step and is not a requirement. Hence, if the external force to be applied does not have an analytical gradient, one can still apply just the force and simulate the dynamics. Having a Jacobian aids in the convergence of the Newton–Raphson iterations, and one can use automatic differentiation or numeric differentiation tools as well to find the Jacobian if needed.

Overall Framework

The overall framework of our simulator is shown in Figure 3.

FIG. 3.

Overall framework of MAT-DiSMech. The first (leftmost) column shows how a physical structure can be modeled as a discrete shell structure. The second column details the inputs to be given to the simulator. The third column depicts the step-by-step simulation procedure. The fourth column represents the results in the form of simulation videos or plots of variation of useful physical quantities, such as displacements with time.

Input

The input to the simulator consists of

Geometry

A text (.txt) file specifying the geometric configuration of the soft structure. The file should contain 1.

the nodal position coordinates $(x, y, z)$ ,

the edge array, that is, the node indices corresponding to each rod-edge, in case there are rod-like structures in the robot body. For example, in Figure 2a, the edge array would be of size $(7 \times 2)$ as follows:

[\begin{matrix} 1 & 2 \\ 2 & 3 \\ 3 & 4 \\ 7 & 8 \\ 8 & 9 \\ 9 & 10 \\ 2 & 8 \end{matrix}]

the triangle array for the shell, that is, the node indices corresponding to each triangle, in case there is a shell structure in the robot body. This can be generated using MATLAB or any other external tool. For example, in Figure 2a, the triangle array would be of size $(10 \times 3)$ as follows:

[\begin{matrix} 10 & 4 & 11 \\ 4 & 5 & 11 \\ 11 & 5 & 12 \\ 12 & 5 & 6 \\ 12 & 6 & 13 \\ 10 & 11 & 14 \\ 11 & 15 & 14 \\ 11 & 12 & 15 \\ 15 & 12 & 16 \\ 12 & 13 & 16 \end{matrix}] .

Other than the above, the dimensions of the cross-section of the rod and the thickness of the shell are to be provided. Some example input text files have been provided in the GitHub repository for MAT-DiSMech.

Material properties

The material densities $ρ_{rod}, ρ_{shell}$ , Young’s moduli $Y_{rod}, Y_{shell}$ , and Poisson’s ratios $ν_{rod}, ν_{shell}$ are inputs to the simulator.

Boundary conditions

The user will specify the fixed nodes and/or edges. Additionally, one can also fix specific DOF for nodes, for instance, at a roller support, the motion is constrained in the $z$ -direction only; hence, to model this, the $x, y$ DOFs of the node can be left free, and $z$ fixed.

Simulation parameters

The user specifies the time step $(Δ t)$ , the total time of simulation, and tolerance on the acceptable error while solving the Newton–Raphson iterations. Additionally, there are flags to decide if the simulation is static or dynamic, if the simulation is 2D or 3D, if line search is enabled to aid in convergence, and if the shell simulation uses the hinge-based bending or the mid-edge-based bending model. Logging and visualization methods are also provided, and the user can decide the frequency at which the data are saved, and the visualization is plotted in MATLAB.

Actuation

Actuation is achieved by modifying the natural curvature/twist. To do so, the user can provide an Excel (.xls) or a text (.txt) file with time-varying values of natural curvature/twist and use the helper functions provided in the framework to apply them at every time step. We have provided a few examples of simulating actuated rod and shell structures in the GitHub repository.

Environment and contact variables

The environment variables are used to apply external forces to the system. The user provides the list of external forces that are to be used in the simulation. The different variables required for the forces implemented within the simulator are 1.

Gravity: acceleration due to gravity ( $g$ ),

Buoyancy: density of the medium $ρ_{med}$

Contact: contact stiffness ( $k_{c}$ ) and contact distance tolerance ( $δ$ ), for self-contact and contact with the ground,

Friction: friction coefficient ( $μ$ ) for the structure’s own material for self-contact as well as friction coefficient with the floor,

Viscous friction: viscosity coefficient $η$ ,

Aerodynamic drag: drag coefficient $C_{D}$ and density of the medium $ρ_{med}$ .

Initial values

Initial values can be specified for twist angle for rod-edges ( $θ$ ), velocity for nodes ( $v$ ), angular velocity for rod-edges, or shell hinge angles ( $ω$ ).

Time stepping

We employ the backward Euler method, an implicit time integration scheme that enables the use of larger time steps compared with explicit methods. However, backward Euler introduces significant artificial numerical damping, which can overly suppress vibrations in the structure. To address this, we also offer an implicit midpoint method, which helps reduce numerical damping and is preferable when accurate vibration behavior is important for a specific application.

Implicit integration provides substantial computational advantages over explicit schemes, particularly for systems with higher Young’s moduli. For context, Choi et al.³⁶ present a comparison of computational times between simulators employing explicit and implicit time-stepping methods. Nevertheless, for very simple systems with low Young’s moduli, explicit methods such as forward Euler can be faster, as they avoid the computational overhead of computing and inverting Jacobians. Therefore, we also include an implementation of the forward Euler method for users who prefer an explicit scheme.

Our framework allows users to easily replace the integration scheme to suit their requirements. It is worth noting that in contact-rich simulations involving soft robotic structures, the implicit solver may sometimes struggle to converge. To improve convergence in such scenarios, we provide an optional line search strategy that determines an appropriate step size, ensuring the error decreases with every iteration.

The complete pseudocode for MAT-DiSMech is given in Algorithm 1.

Algorithm 1: MAT-DiSMech 1: Input: q_k, u_k , Δt, tolerance 2: Output: q_k ₊₁, u_k ₊₁ 3: Require: boundary conditions → free 4: Function timeStepper: 5: Guess: q⁽⁰⁾ ← q_k 6: i ← 0, ϵ ← ∞ 7: while ϵ > tolerance do 8: //prepare for iteration 9: //calculate F^elastic and J^elastic 10: //calculate F^ext and J^ext 11: //calculate F^IMC and J^IMC 12: f ← Mq¨⁽ⁱ⁾ − F^ext − F^elastic − F^IMC 13: J ← M/Δt² − J^ext − J^elastic − J^IMC 14: f_free ← f (free) 15: J_free ← J(free) 16: Δq_free ← linearSolve(J_free, f_free) 17: α ← lineSearch(·) 18: q⁽ⁱ ^{+ 1)}(free) ← q⁽ⁱ ⁾(free) − αΔq_free 19: ϵ ← ∥f_free∥ 20: i ← i + 1 21: end 22: qk + 1 ← q⁽ⁱ ⁾ 23: u_k ₊₁ ← (q_k ₊₁ − q_k )/Δt 24: return q_k ₊₁, u_k ₊₁

Results and Discussion

In this section, we present simulations of various rod-like and shell-like soft structures. We also compare the deflections obtained using our simulator with the analytical values for a cantilever beam.

Simulation experiments

PneuNet

We model a PneuNet actuator as a DER. The change in curvature upon application of pressure is modeled by modifying the natural curvature $\bar{κ}$ of the rod. Figure 4 shows three PneuNet actuators with natural curvatures of $15.70$ , $31.45$ , and $47.15 m^{- 1}$ fixed at one end. Note that discrete curvature is dimensionless. The standard definition of curvature of the dimension 1/length relates to the discrete curvature for the $k$ -th node as $κ_{k} / Δ l_{k}$ . Upon applying the same force at the end effector, these actuators take different shapes owing to their different natural curvatures. The deformed shapes of the actuators under a point force of $[175, 0, 7]^{⊤}$ N at the free end and gravity with $g = [0, 0, - 9.8]^{⊤}$ m/s² are shown in Figure 4. We use $ρ = 1200$ kg/m³ $, E = 20$ GPa, and $ν = 0.5$ as the material parameters for the system. Each rod is $0.1 m$ long and has a cross-sectional radius of $1$ mm. This example is inspired by the work of Payrebrune et al.¹³

FIG. 4.

PneuNet actuator modeled as a discrete elastic rod. (a) Snapshots of the PneuNet actuator with natural curvature $15.7 m^{- 1}$ with only gravity (top) and with gravity and point force on the free end (bottom). (b) Snapshots of the PneuNet actuator with natural curvature $31.45 m^{- 1}$ with only gravity (top) and with gravity and point force on the free end (bottom). (c) Snapshots of the PneuNet actuator with natural curvature $47.15 m^{- 1}$ with only gravity (top) and with gravity and point force on the free end (bottom). Note that the cross-sectional radius of the rod is $1$ mm; it has been enlarged in the figure for better visualization. (d) The plot of the natural shapes of the rods in (a), (b), and (c) under gravity. (e) Plot of the deformed shape of the rods under gravity and point force at the free end.

Earthworm

As our second example, we simulate rectilinear locomotion in earthworms driven by axial peristaltic waves. We represent the worm as a single elastic rod composed of three nodes (forming two edges). Actuation is achieved by varying the natural length $(| | \bar{e} | |)$ of the edges through alternating cycles of contraction and expansion. The contractions and expansions of the two edges are out of phase—while one edge contracts, the other remains unchanged, and vice versa. This alternating pattern generates longitudinal asymmetry, allowing the worm to advance forward by exploiting frictional forces from the ground. Indeed, if there is no friction between the rod and the ground, the earthworm does not move forward and keeps oscillating in the longitudinal direction at its original position. This example is inspired by the work of Rafsanjani et al.⁴⁴ We use $E = 200$ MPa, $ν = 0.5$ and $ρ = 1200$ kg/m³ and a cross-sectional radius of $0.001 m$ for the rod. The length of the rod is $0.1 m,$ and each edge contracts/expands by $0.001 m$ . The external forces acting on the earthworm are gravity, with $g = [0, 0, - 9.8]^{⊤}$ m/s² and frictional contact with the ground. The contact stiffness and coefficient of friction of the ground are 1000 and 0.25, respectively. Figure 5 shows the intermediate poses of the earthworm during one actuation cycle and the trajectory of the front end when friction is present and when it is absent.

FIG. 5.

An earthworm modeled as an elastic rod in MAT-DiSMech. Snapshots of the earthworm’s motion at different time stamps (a1–e1) with friction ( $μ = 0.25$ ) and (a2–e2) without friction. Here, time is given in seconds. (f) Displacement of the $x$ -coordinate of the front end of the earthworm with time. Notice how, when there is no friction, the earthworm does not move forward, and the contraction and expansion of the edges are also less pronounced compared with when there is friction. The time stamps corresponding to the snapshots (a–e) are marked in the plot using vertical dashed lines.

Manta ray

For our third example, we simulate the flapping and undulating locomotion in a manta ray. We model the manta ray body as a discrete elastic shell by discretizing the system into an equilateral triangle mesh. The wingspan and the length of the ray are both $1 m$ . The flapping and undulation are achieved through actuation at the center line. We modify the natural curvature angle ( $\bar{ϕ}$ ) of the hinges in a sinusoidal manner with a frequency of $3.4$ Hz and an amplitude of $2 π / 3$ to achieve flapping. A phase lag of $π$ is introduced between the sine signals at the leading and trailing edge hinges, by varying the phase for the hinges on the center line linearly. We use $E = 6$ GPa, $ν = 0.3, ρ = 1057$ kg/m³ for the manta ray body. The external forces, in this case, are gravity, buoyancy, and aerodynamic drag, and we use $g = [0, 0, - 9.8]^{⊤}$ m/s² $, ρ_{med} = 1000$ kg/m³ and $C_{D} = 0.5$ . Figure 6 shows the intermediate poses from the simulation of a manta ray in motion and the trajectory of the midpoint of its leading edge. It can be seen that the y-coordinate of the leading-edge increases with time, indicating motion in the forward direction, the z-coordinate varies periodically, almost maintaining the altitude, and the x-coordinate does not change since the structure and the loads experienced are symmetric in the x-direction about the centerline.

FIG. 6.

Manta ray modeled as an elastic shell in MAT-DiSMech. (a–d) Snapshots of the Manta ray’s motion at different time stamps. Here, the time is given in seconds. (e) Displacement of the mid-point of the leading edge with time. The time stamps corresponding to the snapshots (a–d) are marked in the plot using vertical dashed lines.

Snake

As our fourth example, we simulate the serpentine locomotion of a snake swimming in a fluid environment, subject to hydrodynamic forces. The snake is modeled as an elastic rod with a length of $0.1 m$ , a cross-sectional radius of $1$ mm, Young’s modulus $E = 2$ MPa, Poisson’s ratio $ν = 0.5$ , and density $ρ = 1200$ kg/m³. The surrounding fluid exerts hydrodynamic forces computed using RFT (Gray et al.⁴³), with tangential and normal drag coefficients set to $0.01$ and $0.1$ , respectively. This anisotropy in drag coefficients enables forward propulsion when serpentine undulations are applied. Serpentine motion is generated by modulating the natural curvature, $\bar{κ}$ , of the rod’s bending–twisting springs to produce a sinusoidal wave along the snake’s body, with an amplitude of $0.05 m$ and a frequency of $1$ Hz. Figure 7 presents snapshots of the snake in motion and the trajectory of its leading end. The plots show that the $x$ -coordinate steadily increases over time as the snake propels forward, and the $y$ -coordinate oscillates sinusoidally reflecting the serpentine path, while the $z$ -coordinate remains constant at 0.

FIG. 7.

Snake modeled as an elastic rod in MAT-DiSMech. (a–e) Snapshots of the snake’s motion at different time stamps. Here, the time is given in seconds. (f) Displacement of the leading end of the snake with time. The time stamps corresponding to the snapshots (a–e) are marked in the plot using vertical dashed lines.

Parachute

For our fifth showcase, we simulate a parachute that is modeled as a combination of elastic rods and shell. The canopy is modeled as an initially flat plate in the shape of a regular hexagon of side length $1 m$ and thickness $1$ mm. The ropes are modeled as elastic rods. A rigid payload mass of $0.13$ kg is attached at the end of the ropes, which is modeled simply by increasing the mass at the node shared between all rods. The material parameters used in this simulation are $E_{rod} = 10$ MPa, $E_{shell} = 1$ GPa, $ν_{rod} = 0.5, ν_{shell} = 0.3,$ and $ρ = 1500$ kg/m³ for both rod and shell. The external forces, in this case, are gravity and aerodynamic drag. We use $g = [0, 0, - 9.8]^{⊤}$ m/s², and, $ρ_{med} = 1$ kg/m³ and $C_{D} = 10$ to calculate the gravity and drag force, respectively. Figure 8 shows the snapshots of the parachute as it falls to the ground. It can be seen that as the simulation progresses, the initially flat canopy takes a 3D-inverted bowl-like shape as expected.

FIG. 8.

A point mass falling with a parachute modeled as a combination of rods and shell within MAT-DiSMech. (a–e) Snapshots of the parachute at different time stamps. Here, the time is given in seconds. (f) Displacement of the hanging mass with time. The time stamps corresponding to the snapshots (a–e) are marked in the plot using vertical dashed lines.

Rod dropped on ground

For our sixth example, we simulate a rod falling onto the ground from a height of $5$ cm. The density of the rod’s material and its Poisson ratio are taken to be $1500$ kg/m³ and $0.3$ , respectively. We simulate the rod for two different values of Young’s modulus $Y$ : (a) $Y = 2$ MPa and (b) $Y = 2$ GPa. The external forces, in this case, are gravity, ground contact, and friction. The values of $[0, 0, - 9.8]^{⊤}$ m/s², $0.25$ , and $20$ are used, respectively, for the acceleration due to gravity $g$ , the friction coefficient $μ$ , and contact stiffness $k_{c}$ between the ground and the rod. Figure 9 shows snapshots of the rod as it hits the ground. As Young’s modulus for the second case is much higher than the first case, we observe significant deformation of the rod in the first case, while almost no deformation in the second.

FIG. 9.

Elastic rods with different Young’s moduli dropped on the ground. (a1–f1) Snapshots of the elastic rod with Young’s modulus 2 MPa, and (a2–f2) snapshots of the elastic rod with Young’s modulus 2 GPa. Here, time is given in seconds. Note that the cross-sectional radius of the rod is 1 mm; it has been enlarged in the figure for better visualization. (g) The trajectory of the z-coordinate of the lowermost point on the rod with time. The time stamps corresponding to the snapshots (a–f) are marked in the plot using vertical dashed lines.

Rod gripping a sphere

For our seventh example, we demonstrate a rod gripping a fixed spherical object. The rod has material properties of $E = 200$ MPa, $ν = 0.5$ , and $ρ = 1200$ kg/m³. It measures $0.1 m$ in length with a cross-sectional radius of $0.001 m$ . Initially positioned along the $x$ -axis, the rod is clamped at one end at the origin, while the sphere is located at the point $[0.1, 0.0, 0.05]^{⊤}$ . Actuation is performed by gradually varying the natural curvature of the bending–twisting springs along the final three-quarters of the rod’s length, transitioning from $[0, 0]^{⊤}$ to $[0.1886, 0]^{⊤}$ . The external forces in this case arise from contact and friction between the rod and the sphere. The simulations use a contact stiffness of 20 and a coefficient of friction of 0.25. Figure 10 illustrates intermediate configurations of a clamped slender rod undergoing actuation via changes in its natural curvature. In case (a), with contact disabled, the rod adopts its natural shape without regard for the spherical obstacle, passing through it if necessary. In contrast, case (b) demonstrates that when contact forces are included through our penalty energy formulation, the rod successfully wraps around the sphere.

FIG. 10.

Intermediate configurations of a clamped slender rod actuated by varying its natural curvature to bend inward. The rod’s color transitions from yellow to red to indicate progression over time. (a) Without contact modeling, the rod passes through the spherical object to achieve its natural curvature-driven shape. (b) With contact and friction forces modeled via our penalty energy method, the rod interacts with the sphere and successfully wraps around it.

Overhand knot

As our final example, we simulate the tightening of an overhand knot to demonstrate the capability of the simulator to accurately capture self-contact and frictional interactions. The rod is initialized in a looped configuration, and tightening is induced through a moving-boundary condition in which the two nodes at the free ends are pulled apart and subsequently released over a total duration of 62 s (60 s of tightening followed by 2 s of release). Figure 11 shows representative snapshots of the rod configuration during the knot-tightening process. We use a friction coefficient of $μ = 0.3,$ and the material parameters are Young’s modulus $0.18 MPa$ , Poisson’s ratio $0.5$ , density $1180 {kg/m}^{3}$ , and cross-sectional radius $1.6 mm$ .

FIG. 11.

Tightening of an overhand knot. (a–d) Rod configurations at different time instances ( $μ =$ 0.3). Here, time is shown in seconds.

In addition to these examples, prior studies have demonstrated that the modeling techniques utilized in our simulator can capture significantly more complex phenomena such as flagella buckling,¹⁰ tying a variety of knots,²² and shear-induced pitchfork bifurcations in ribbons.⁴⁵ While our examples in this article are not intended to precisely replicate specific biological trajectories, they are deliberately chosen to showcase the simulator’s ability to model sophisticated motions and handle challenging interactions with the environment—including contact dynamics, frictional forces, and hydrodynamic resistance.

Validation for cantilever beam case

In Figure 12, we show the plots of the static displacement of the free end of a cantilever beam under gravity. The system is simulated as a rod, a hinge-based shell, and mid-edge normal-based shell. The density of the system is $1200$ kg/m³, and the Poisson’s ratio is $0.5$ . We perform a static simulation in MAT-DiSMech for various Young’s moduli, specifically $20$ GPa, $2$ GPa, $200$ MPa and $20$ MPa, and compare the obtained displacement with the theoretical values obtained using the Euler–Bernoulli beam theory.⁴⁶ From Figure 12a, it is clear that the algorithm is able to capture the deflections in the cantilever beam accurately. Further, from Figure 12b, it can be seen that the values from the simulator match the values obtained from the Euler–Bernoulli theory more closely at larger values of Young’s modulus and move away from the theoretical values at lower values of Young’s modulus. This is expected since Euler–Bernoulli’s theory is based on approximations that work well for low-deformation cases that correspond to the higher Young’s modulus situation.

FIG. 12.

Comparison between the cantilever free end deflection obtained from our simulator and the theoretical value from Euler–Bernoulli beam theory. (a) The log-log plot of the deflection with Young’s modulus for DER, hinge-based, midge-based shell, and Euler–Bernoulli beam. (b) A plot of the difference between the value of the deflection obtained from static simulation in MAT-DiSMech and using the Euler–Bernoulli beam theory with Young’s modulus in log scale for DER, hinge-based shell, and midge-based shell.

To examine a limitation of the hinge-based bending model for shells, we perform simulations assessing its sensitivity to mesh quality. It is well known that the hinge model performs reliably only on meshes composed of nearly equilateral triangles with good aspect ratios,²⁶ whereas ill-conditioned meshes with skewed elements or extreme aspect ratios can result in significant errors or convergence issues. Figure 13 presents the normalized deflection error for a cantilever beam under gravity, comparing the hinge and mid-edge bending models across five mesh configurations, from highly regular (equilateral) to severely distorted (nonuniform). The results indicate that the mid-edge bending model maintains consistently low errors across all mesh types, demonstrating robustness against mesh irregularities. In contrast, the hinge model shows significant errors for the right-isosceles and nonuniform meshes. The right-isosceles mesh features many hinges oriented obliquely to the main bending direction, while the nonuniform mesh contains triangles with poor aspect ratios and skewness. Overall, these results highlight the hinge method’s mesh dependence, whereas the mid-edge method’s edge-associated DOFs substantially improve robustness to mesh quality.

FIG. 13.

Normalized deflection at the free end of a cantilever beam (length 0.1 m, width 0.02 m, thickness 0.001 m, $E = 2$ GPa, $ν = 0.5$ ) under gravity, modeled with hinge and mid-edge bending across different mesh types: equilateral, random, right-isosceles, equilateral-aligned, and nonuniform.

Conclusions

The MAT-DiSMech simulation framework presented in this article offers a comprehensive and versatile approach to simulating soft robotic systems, specifically those composed of rods, shells, or their combinations. By leveraging DDG methods, our MATLAB-based simulator achieves a balance between accuracy and efficiency, effectively addressing the challenges posed by geometrically nonlinear dynamics in soft robots. The modular design supports a broad range of soft robotic configurations and dynamic forces, including gravity, contact, Coulomb, and viscous friction, as well as hydrodynamic and aerodynamic drag, allowing researchers to tailor simulations to specific system requirements. We present a diverse set of showcase examples involving complex environmental interactions, such as frictional ground contact and viscous fluid forces, for both rod and shell-like structures. We also demonstrate numerical validation against theoretical solutions, highlighting the strengths and limitations of different modeling techniques.

Since controlled actuation plays a central role in soft robotics, MAT-DiSMech supports time-dependent adjustments of natural stretch, curvature, and twist. This enables the simulation of dynamic behaviors such as shape morphing, as demonstrated in examples including an earthworm, manta ray, snake, PneuNet actuator, and gripping manipulator. In conclusion, MAT-DiSMech provides an accurate and efficient platform for simulating soft robot dynamics. It facilitates rapid prototyping and validation of diverse soft robotic designs, thereby advancing research along Sim2Real pathways. In the future, we plan on expanding the scope of MAT-DiSMech by integrating data-driven modeling techniques to model material nonlinearity and black-box external forces.

Footnotes

Author Disclosure Statement

No competing financial interests exist.

Funding Information

The authors acknowledge financial support from the National Science Foundation (grant numbers 20476630, 2209782) and the National Institute of Neurological Disorders and Stroke (1R01NS141171-01).

Supplemental Material

Appendix

References

Wang

, Iida

. Deformation in soft-matter robotics: A categorization and quantitative characterization. IEEE Robot Automat Mag, 2015; 22(3):125–139.

Choi

, Crump

, Duriez

, et al. On the use of simulation in robotics: Opportunities, challenges, and suggestions for moving forward. Proc Natl Acad Sci U S A, 2021; 118(1):e1907856118.

Bergeles

, Yang

G-Z

. From passive tool holders to microsurgeons: Safer, smaller, smarter surgical robots. IEEE Trans Biomed Eng, 2014; 61(5):1565–1576.

Lipson

. Challenges and opportunities for design, simulation, and fabrication of soft robots. Soft Robot, 2014; 1(1):21–27.

Armanini

, Boyer

, Mathew

, et al. Soft robots modeling: A structured overview. IEEE Trans Robot, 2023; 39(3):1728–1748.

Qin

, Peng

, Huang

, et al. Modeling and simulation of dynamics in soft robotics: A review of numerical approaches. Curr Robot Rep, 2023; 5(1):1–13.

Crane

, Wardetzky

. A glimpse into discrete differential geometry. Notices Amer Math Soc, 2017; 64(10):1153–1159.

Shabana

, Hussien

, Escalona

. Application of the absolute nodal coordinate formulation to large rotation and large deformation problems. Journal of Mechanical Design, 1998; 120(2):188–195.

Shabana

. An overview of the ancf approach, justifications for its use, implementation issues, and future research directions. Multibody Syst Dyn, 2023; 58(3–4):433–477.

10.

Jawed

, Khouri

, Da

, et al. Propulsion and instability of a flexible helical rod rotating in a viscous fluid. Phys Rev Lett, 2015; 115(16):168101.

11.

Jawed

, Da

, Joo

, et al. Coiling of elastic rods on rigid substrates. Proc Natl Acad Sci USA, 2014; 111(41):14663–14668.

12.

Huang

, Huang

, Majidi

, et al. Dynamic simulation of articulated soft robots. Nat Commun, 2020; 11(1):2233.

13.

de Payrebrune

, O’Reilly

. On constitutive relations for a rod-based model of a pneunet bending actuator. Extreme Mech Lett, 2016; 8:38–46.

14.

, Wang

, Ge

, et al. Analytical modeling and design of generalized pneu-net soft actuators with three-dimensional deformations. Soft Robot, 2021; 8(4):462–477.

15.

Lotfiani

, Yi

, Shao

, et al. Analytical modeling and optimization of a corrugated soft pneumatic finger considering the performance of pinch and power grasps. Extreme Mech Lett, 2021; 44:101215.

16.

Liao

, Zang

, Chen

, et al. Soft rod-climbing robot inspired by winding locomotion of snake. Soft Robot, 2020; 7(4):500–511.

17.

Wan

, Sun

, Qin

, et al. Design, analysis, and real-time simulation of a 3d soft robotic snake. Soft Robot, 2023; 10(2):258–268.

18.

Shepherd

, Ilievski

, Choi

, et al. Multigait soft robot. Proc Natl Acad Sci USA, 2011; 108(51):20400–20403.

19.

Shintake

, Cacucciolo

, Floreano

, et al. Soft robotic grippers. Advanced Materials, 2018; 30(29):1707035.

20.

Bergou

, Wardetzky

, Robinson

, et al. Discrete elastic rods. ACM Trans Graph, 2008; 27(3):1–12.

21.

Jawed

, Dieleman

, Audoly

, et al. Untangling the mechanics and topology in the frictional response of long overhand elastic knots. Phys Rev Lett, 2015; 115(11):118302.

22.

Choi

, Tong

, Jawed

, et al. Implicit contact model for discrete elastic rods in knot tying. Journal of Applied Mechanics, 2021; 88(5).

23.

Rus

, Tolley

. Design, fabrication and control of soft robots. Nature, 2015; 521(7553):467–475.

24.

Baraff

, Witkin

. Large Steps in Cloth Simulation , 1st ed.. Association for Computing Machinery: New York, NY; 2023.

25.

Grinspun

, Hirani

, Desbrun

, et al. Discrete shells. In: Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ‘03 , Eurographics Association: Goslar, DEU; 2003; pp. 62–67.

26.

Grinspun

, Gingold

, Reisman

, et al. Computing discrete shape operators on general meshes. Computer Graphics Forum (Eurographics), 2006; 25(3):547–556.

27.

Chen

H-Y

, Sastry

, van Rees

, et al. Physical simulation of environmentally induced thin shell deformation. ACM Trans Graph, 2018; 37(4):1–13.

28.

Koenig

, Howard

. “Design and use paradigms for gazebo, an open-source multi-robot simulator,” In: 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), vol. 3. IEEE; 2004; pp. 2149–2154.

29.

Coumans

, Bai

. “Pybullet, a python module for physics simulation for games, robotics and machine learning, 2016–2021.” Available from: http://pybullet.org

30.

Graule

, Teeple

, McCarthy

, et al. Somo: Fast and accurate simulations of continuum robots in complex environments. In: 2021 IEEE International Conference on Intelligent Robots and Systems (IROS) , IEEE; 2021, pp. 3934–3941.

31.

Coevoet

, Morales-Bieze

, Largilliere

, et al. Software toolkit for modeling, simulation, and control of soft robots. Advanced Robotics, 2017; 31(22):1208–1224.

32.

Holland

, Park

, Polygerinos

, et al. The soft robotics toolkit: Shared resources for research and design. Soft Robot, 2014; 1(3):224–230.

33.

Tasora

, Serban

, Mazhar

, et al. “Chrono: An open source multiphysics dynamics engine”. In: High Performance Computing in Science and Engineering – Lecture Notes in Computer Science. Springer; 2016; pp. 19–49.

34.

Naughton

, Sun

, Tekinalp

, et al. Elastica: A compliant mechanics environment for soft robotic control. IEEE Robot Autom Lett, 2021; 6(2):3389–3396.

35.

Mathew

, Hmida

IMB

, Armanini

, et al. Sorosim: A matlab toolbox for hybrid rigid-soft robots based on the geometric variable-strain approach. IEEE Robot Automat Mag, 2022; 30:2–18.

36.

Choi

, Jing

, Sabelhaus

, et al. Dismech: A discrete differential geometry-based physical simulator for soft robots and structures. IEEE Robot Autom Lett, 2024; 9(4):3483–3490.

37.

Tong

, Choi

, Joo

, et al. A fully implicit method for robust frictional contact handling in elastic rods. Extreme Mech Lett, 2023; 58:101924.

38.

Kirchhoff

. Ueber das gleichgewicht und die bewegung eines unendlich dünnen elastischen stabes. J Für Die Reine Und Angewandte Mathematik, 1859; 56:285–313.

39.

Goldenthal

, Harmon

, Fattal

, et al. Efficient simulation of inextensible cloth. ACM Trans Graph, 2007; 26(3):49.es.

40.

Savin

, Kurpios

, Shyer

, et al. On the growth and form of the gut. Nature, 2011; 476(7358):57–62.

41.

Bergou

, Audoly

, Vouga

, et al. Discrete viscous threads. ACM Trans Graph, 2010; 29(4):1–10.

42.

Lumelsky

. On fast computation of distance between line segments. Inf Process Lett, 1985; 21(2):55–61.

43.

Gray

, Hancock

. The propulsion of seaurchin spermatozoa. J Exp Biol, 1955; 32:802–814.

44.

Rafsanjani

, Zhang

, Liu

, et al. Kirigami skins make a simple soft actuator crawl. Sci Robot, 2018; 3(15):eaar7555.

45.

Huang

, Wang

, Li

, et al. Shear induced supercritical pitchfork bifurcation of prebuckled bands, from narrow strips to wide plates. J Mech Phys Solids, 2020; 145:104168.

46.

Bauchau

, Craig

. Euler-Bernoulli beam theory . Springer Netherlands: Dordrecht; 2009, pp. 173–221.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

37.32 MB

0.00 MB

37.32 MB