Computational mechanics of soft filaments

Computational mechanics of soft filaments

Mattia Gazzola Department of Mechanical Science and Engineering, National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA    Levi H. Dudte John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA    Andrew G. McCormick Google, Mountain View, CA 94043, USA John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA    L. Mahadevan lm@seas.harvard.edu John A. Paulson School of Engineering and Applied Sciences, and Department of Physics, Harvard University, Cambridge, MA 02138, USA
Abstract

Soft slender structures are ubiquitous in natural and artificial systems and can be observed at scales that range from the nanometric to the kilometric, from polymers to space tethers. We present a practical numerical approach to simulate the dynamics of filaments that, at every cross-section, can undergo all six possible modes of deformation, allowing the filament to bend, twist, stretch and shear, while interacting with complex environments via muscular activity, surface contact, friction and hydrodynamics. We examine the accuracy of our method by means of several benchmark problems with known analytic solutions. We then demonstrate the capabilities and robustness of our approach to solve forward problems in physics and mechanics related to solenoid and plectoneme formation in twisted, stretched filaments, and inverse problems related to active biophysics of limbless locomotion on solid surfaces and in bulk liquids. All together, our approach provides a robust computational framework to characterize the mechanical response and design of soft active slender structures.

I Introduction

Quasi one-dimensional objects are characterized by having one dimension, the length , much larger than the others, say the radius , so that . Relative to three-dimensional objects, this measure of slenderness allows for significant mathematical simplification in accurately describing the physical dynamics of strings, filaments and rods. It is thus perhaps not surprising that the physics of strings has been the subject of intense study for centuries Kirchhoff:1859 (); Clebsch:1883 (); Love:1906 (); Cosserat:1909 (); Antman:1973 (); Dill:1992 (); Langer:1996 (); vanDerHeijden:2000 (); Goss:2005 (); Maugin:2010 (), and indeed their investigation substantially pre-dates the birth of three-dimensional elasticity.

Following the pioneering work of Galileo on the bending of cantilevers, one-dimensional analytical models of beams date back to 1761 when Jakob Bernoulli first introduced the use of differential equations to capture the relation between geometry and bending resistance in a planar elastica, that is an elastic curve deforming in a two-dimensional space. This attempt was then progressively refined by Huygens, Leibniz and Johann Bernoulli Goss:2009 (), until Euler presented a full solution of the planar elastica, obtained by minimizing the strain energy and by recognizing the relation between flexural stiffness and material and geometric properties. Euler also showed the existence of bifurcating solutions in a rod subject to compression, identifying its first buckling mode, while Lagrange formalized the corresponding multi-modal solution Antman:1973 (). Non-planar deformations of the elastica were first tackled by Kirchhoff Kirchhoff:1859 (); Dill:1992 () and Clebsch Clebsch:1883 () who envisioned a rod as an assembly of short undeformable straight segments with dynamics determined by contact forces and moments, leading to three-dimensional configurations. Later, Love Love:1906 () approached the problem from a Lagrangian perspective characterizing a filament by contiguous cross sections that can rotate relative to each other, but remain undeformed and perpendicular to the centerline of the rod at all times; in modern parlance this assumption is associated with dynamics on the rotation group SO(3) at every cross-section. The corresponding equations of motion capture the ability of the filament to bend and twist, but not shear or stretch. Eventually the Cosserat brothers Cosserat:1909 () relaxed the assumption of inextensibility and cross-section orthogonality to the centerline, deriving a general mathematical framework that accommodates all six possible degrees of freedom associated with bending, twisting, stretching and shearing, effectively formulating dynamics on the full Euclidean group SE(3).

The availability of these strong mathematical foundations Antman:1973 () prompted a number of discrete computational models Yang:1993 (); Spillmann:2007 (); Bergou:2008 (); Arne:2010 (); Audoly:2013 () that allow for the exploration of a range of physical phenomena. These include, for example, the study of polymers and DNA Yang:1993 (); Wolgemuth:2000 (), elastic and viscous threads Bergou:2010 (); Arne:2010 (); Audoly:2013 (); Brun:2015 (), botanical applications Goriely:1998 (); Gerbode:2012 (), elastic ribbons Bergou:2008 (), woven cloth Kaldor:2008 () and tangled hair Ward:2007 (). Because the scaled ratio of the stretching and shearing stiffness to the bending stiffness for slender filaments is , the assumption of inextensibility and unshearability is usually appropriate, justifying the widespread use of the Kirchhoff model in the aforementioned applications.

However, new technologies such as soft robotics and artificial muscles Kim:2013 (); Haines:2014 (); Park:2016 () are generating applications for deformable and stretchable elastomeric filamentous structures Haines:2014 () for which the assumption of inextensibility and unshearability is no longer valid. Motivated by these advancements, we move away from the Kirchhoff model in favor of the complete Cosserat theory. We present a robust and relatively simple numerical scheme that tracks both the rod centerline and local frame allowing for bend, twist, stretch and shear Cosserat:1909 () consistent with the full Euclidean group SE(3), while retaining the Hamiltonian structure of the system and fast discrete operators. This allows us to substantially increase the spectrum of problems amenable to be treated via this class of rod models.

Moving beyond the passive mechanics of individual filaments, we also account for the interaction between filaments and complex environments with a number of additional biological and physical features, including muscular activity, self-contact and contact with solid boundaries, isotropic and anisotropic surface friction and viscous interaction with a fluid. Finally, we demonstrate the capabilities and the robustness of our solver by embedding it in an inverse design cycle for the identification of optimal terrestrial and aquatic limbless locomotion strategies.

The paper is structured as follows. In Section II we review and introduce the mathematical foundations of the model. In Section III we present the corresponding discrete scheme and validate our framework against a battery of benchmark problems. In Section IV we detail the physical and biological enhancements to the original model, and finally in Section V we showcase the potential of our solver via the study of solenoids and plectonemes as well as limbless biolocomotion. Mathematical derivations and additional validation test cases are presented in the Appendix.

Ii Governing Equations

We consider filaments as slender cylindrical structures deforming in three-dimensions with a characteristic length which is assumed to be much larger than the radius () at any cross section. Then the filament can be geometrically reduced to a one-dimensional representation, and its dynamical behavior may be approximated by averaging all balance laws at every cross section Antman:1973 (). We start with a description of the commonly used Kirchhoff-Love theory that accounts for bend and twist at every cross-section but ignores stretch and shear, before moving on to the Cosserat theory that accounts for these additional degrees of freedom as well.

ii.1 Kirchhoff-Love theory for inextensible, unshearable rods

Figure 1: The Cosserat rod model. A filament deforming in the three-dimensional space is represented by a centerline coordinate and a material frame characterized by the orthonormal triad . The corresponding orthogonal rotation matrix with row entries , , transforms a vector from the laboratory canonical basis to the material frame of reference so that and vice versa . If extension or compression is allowed, the current filament configuration arc-length may no longer coincide with the rest reference arc-length . This is captured via the scalar dilatation field . Moreover, to account for shear we allow the triad to detach from the unit tangent vector so that (we recall that the condition and correspond to the Kirchhoff constraint for unshearable and inextensible rods, and implies that ). The dynamics of centerline and material frame are related through quadratic energy functionals that give rise at each cross section to the internal force and torque resultants, and . External loads are represented via the force and couple line densities.

As illustrated in Fig. 1, a filament in the Cosserat rod theory can be described by a centerline and an oriented frame of reference equivalent to the orthonormal triad of unit vectors . Here, is the centerline arc-length coordinate in its current configuration and is time.

Denoting by any generic vector represented in the Eulerian frame and as the body-convected (Lagrangian) frame of reference allows us to write

laboratory: (1)
body-convected: (2)

where Eq. (1) expresses in the laboratory canonical basis , while Eq. (2) expresses the same vector in the body-convected director basis . Then, the matrix transforms any vector from the laboratory to the body-convected representation via and conversely, , since . In general we need to distinguish between the arc-length coordinate that corresponds to the current filament configuration and the arc-length coordinate associated with the reference configuration of the filament, due to stretching (Fig. 1 (throughout we will use hatted quantities to denote the reference configuration).

We will first start by presenting the equations of motion under the assumption of inextensibility (i.e. ), before generalizing them to the stretchable case in the subsequent sections. Denoting the rod angular velocity as and the generalized curvature as , where denotes the 3-vector associated with the skew-symmetric matrix , the following transport identities hold

(3)

Using the above equations (full derivation in the Appendix) we can express the advection of the rod positions and local frames, as well as the linear and angular momentum balance in a convenient Eulerian-Lagrangian form

(4)
(5)
(6)
(7)

where is the constant material density, is the cross sectional area, is the velocity, and are respectively the internal force and couple resultants, and are external body force and torque line densities, and the tensor is the second area moment of inertia (throughout this study we assume circular cross sections, see Appendix).

To close the above system of equations, Eqs. (4-7) and determine the dynamics of the rod, it is necessary to specify the form of the internal forces and torques generated in response to bend and twist, corresponding to the three degrees of freedom at every cross-section. The strains are defined as the relative local deformations of the rod with respect to its natural strain-free reference configuration. Bending and twisting strains are associated with the spatial derivatives of the material frame director field and are characterized by the generalized curvature. Specifically, the components of the curvature projected along the directors () coincides with bending () and twist () strains in the material frame (Table 1).

Assuming a linear material constitutive law implies linear stress-strain relations. Integration of the torque densities over the cross sectional area yields the bending and twist rigidities (Table 1), so that the resultant torque-curvature relations can be generically expressed in vectorial notation as

(8)

where (, , ) is the bend/twist stiffness matrix with the flexural rigidity about , the flexural rigidity about , the twist rigidity about . Here, the vector characterizes the intrinsic curvatures of a filament that in its stress-free state is not straight. We wish to emphasize here that the constitutive laws are most simply expressed in a local Lagrangian form; hence the use of and not .

The Kirchhoff rod is defined by the additional assumptions that there is no axial extension or compression or shear strain. Then the arc-length coincides with at all times, and the tangent to the centerline is also normal to the cross-section, so that Antman:1973 (). This implies that serves as a Lagrange multiplier, and that the torque-curvature relations of Eq. 8 are linear.

This completes the formulation of the equations of motion for the Kirchhoff rod, and when combined with boundary conditions suffices to have a well-posed initial boundary value problem. For the general stretchable and shearable case, all geometric quantities (A, , , etc.) must be rescaled appropriately, as addressed in the following sections.

ii.2 Cosserat theory of stretchable and shearable filaments

In the general case of soft filaments, at every cross-section we also wish to capture transverse shear and axial strains in addition to bending and twisting. Since we wish to account for all six deformation modes associated with the six degrees of freedom at each cross section along the rod, we must augment the Kirchhoff description in the previous section and add three more constitutive laws to define the local stress resultants . In fact, they are no longer defined as Lagrange multipliers that enforce the condition that normals to the cross-section coincide with tangents to the centerline, i.e. we now must have .

The shear and axial strains are associated with the deviations between the unit vector perpendicular to the cross-section and the tangent to the centerline, and thus may be expressed in terms of the derivatives of the centerline coordinate . In the material frame of reference, we characterize these strains by the vector (Fig. 1) which then takes the form

(9)

Here, the scalar field expresses the local stretching or compression ratio (Fig. 1) relative to the rest reference configuration () and is the unit tangent vector.

Whenever the filament undergoes axial stretching or compression the corresponding infinitesimal elements deform and all related geometric quantities are affected. By assuming that the material is incompressible and that the cross sections retain their circular shapes at all times, we can conveniently express the governing equations with respect to the rest reference configuration of the filament (denoted by a hat) in terms of the local dilatation . Then, the following relations hold

(10)
deformation modes   strains   rigidities   loads
bending about         
bending about         
twist about         
shear along         
shear along         
stretch along         
Table 1: Constitutive laws. The generalized curvature is associated with the bending , about the principal directions (, ) and the twist about the longitudinal one (), while is associated with the shears , along the principal directions (, ) and the axial extensional or compression along the longitudinal one (). The material properties of the rod are captured through the Young’s () and shear (G) moduli, while its geometric properties are accounted for via the cross sectional area , the second moment of inertia and the constant for circular cross sections Gere:2001 (). The diagonal entries of the bending/twist and shear/stretch matrices are, respectively, (, , ) and (, , ). Pre-strains are modeled via the intrinsic curvature/twist and shear/stretch .

As with the Kirchhoff rod, assuming a linear material constitutive law implies linear stress-strain relations. Integration of the stress and couple densities over the cross sectional area yields both the rigidities associated with axial extension and shear (Table 1), so that the resultant load-strain relations can be generically expressed in vectorial notation as

(11)

where (, , ) is the shear/stretch stiffness matrix with the shearing rigidity along , the shearing rigidity along , and the axial rigidity along . Here, as with the Kirchhoff rod, the vector corresponds to the intrinsic shear and stretch, and must be accounted for in the case of stress-free shapes that are non-trivial. Although the intrinsic strains , are implemented in our solver to account for pre-strained configurations, to simplify the notation in the remaining text we will assume that the filament is intrinsically straight in a stress-free state, so that .

The rigidities associated with bending, twisting, stretching and shearing are specified in Table 1, and can be expressed as the product of a material component, represented by the Young’s () and shear () moduli, and a geometric component represented by , and the constant for circular cross sections Gere:2001 (). We note that the rigidity matrices and are assumed to be diagonal throughout this study, although off diagonal entries can be easily accommodated to model anisotropic materials such as composite elements. In general, this mathematical formulation can be extended to tackle a richer set of physical problems including viscous threads Arne:2010 (); Audoly:2013 (), magnetic filaments Landau:1984 (), etc., by simply modifying the entries of and and introducing time-dependent constitutive laws wherein and , as for example in Arne:2010 (). We also emphasize here that in the case of stretchable rods, and are no longer constant, rendering the load-strain relations non-linear (Eqs. 8, 11), even though the stress-strain relations remain linear.

Having generalized the constitutive relations to account for filament stretchiness and shearability, we now generalize the equations of motion for this case. Multiplying both sides of Eqs. (6,7) by and substituting the above identities together with the constitutive laws of Eqs. (11) into Eqs. (4-7), yields the final system

(12)
(13)
(14)
(15)

where is the infinitesimal mass element, and is the infinitesimal mass second moment of inertia. We note that the left hand side and the unsteady dilatation term of Eq. (15) arise from the expansion of the original rescaled angular momentum via chain rule. We also note that the external force and couple are defined as and (with and the force and torque line densities, respectively) so as to account for the dependence on .

Combined with some initial and boundary conditions, Eqs. (12-15) express the dynamics and kinematics of the Cosserat rod with respect to its initial rest configuration, in a form suitable to be discretized as described in Section III.

Iii Numerical Method

To derive the numerical method for the time evolution of a filament in analogy with the continuum model of Section II, we first recall a few useful definitions for effectively implementing rotations. We then present the spatially discretized model of the rod, and the time discretization approach employed to evolve the governing equations.

Figure 2: Discretization model. A discrete filament is represented through a set of vertices and a set of material frames . Two consecutive vertices define an edge of length along the tangent unit vector . The dilatation is defined as , where is the edge rest length. The vector represents the discrete shear and axial strains. The mass of the filament is discretized in pointwise concentrated masses at the locations for the purpose of advecting the vertices in time. For the evolution of in time, we consider instead the mass second moment of inertia associated with the cylindrical elements depicted in blue.

iii.1 Rotations

Bending and twisting deformations of a filament involve rotations of its material frame in space and time. To numerically simulate the rod, it is critical to represent and efficiently compute these nonlinear geometric transformations fast and accurately. A convenient way to express rotations in space or time is the matrix exponential Baillieul:1987 (); Levi:1996 (); Goriely:2000 (). Assuming that the matrix denotes the rotation by the angle about the unit vector axis , then this rotation can be expressed through the exponential matrix , and efficiently computed via the Rodrigues formula Rodrigues:1840 ()

(16)

Here represents the skew-symmetric matrix associated with the unit vector

where the operator allows us to transform a vector into the corresponding skew-symmetric matrix, and vice versa .

Conversely, given a rotation matrix , the corresponding rotation vector can be directly computed via the matrix logarithm operator

It is important to notice that the rotation axis is expressed in the material frame of reference associated with the matrix (or ). With these tools in hand, we now proceed to outline our numerical scheme.

iii.2 Spatial discretization

Drawing from previous studies of unshearable and inextensible rods Spillmann:2007 (); Bergou:2008 (); Sobottka:2008 (), we capture the deformation of a filament in three-dimensional space via the time evolution of a discrete set of vertices and a discrete set of material frames , as illustrated in Fig. 2.

Each vertex is associated with the following discrete quantities

(17)

where is the velocity , is a pointwise concentrated mass, and is the external force given in Eq. (14).

Each material frame is associated with an edge connecting two consecutive vertices, and with the related discrete quantities

(18)

where , , are the edge current length, reference length and dilatation factor, is the discrete tangent vector, is the discrete shear/axial strain vector, is the discrete angular velocity, , , , are the edge reference cross section area, mass second moment of inertia, bend/twist matrix and shear/stretch matrix, and finally is the external couple given in Eq. (15).

Whereas in the continuum setting (Section II) all quantities are defined pointwise, in a discrete setting some quantities, and in particular , are naturally expressed in an integrated form over the domain along the filament Grinspun:2006 (); Bergou:2008 (). Any integrated quantity divided by the corresponding integration domain length is equivalent to its pointwise average. In the context of our discretization the domain becomes the Voronoi region of length

(19)

which is defined only for the interior vertices . Each interior vertex is then also associated with the following discrete quantities

(20)

where is the Voronoi domain length at rest and is Voronoi region dilatation factor. Recalling that the generalized curvature expresses a rotation per unit length about its axis, then the quantity naturally expresses the rotation that transforms a material frame to its neighbour over the segment size along the rod. Therefore, the relation holds, so that . Finally, we introduce the bend/twist stiffness matrix consistent with the Voronoi representation.

Then, we may discretize the governing Eqs. (12-15) so that they read

(21)
(22)
(23)
(24)

where is the discrete difference operator and is the averaging operator to transform integrated quantities over the domain to their point-wise counterparts. We note that and operate on a set of vectors and returns vectors, consistent with Eqs. (21-24) (see the Appendix for further details).

iii.3 Time discretization

While the above equations are conservative and preserve energy, in general when additional physical effects and interactions with complex environments are considered, the equations of motion are not conservative. We choose a symplectic, second-order Verlet scheme to integrate the equations of motion in time so that our numerical scheme is energy-preserving in the case of conservative dynamics. We note that despite the failure of Verlet schemes to integrate rotational equations of motion when represented by quaternions, in our case their use is acceptable as rotations are represented instead by Euler angles Kol:1997 ().

The second order position Verlet time integrator is structured in three blocks: a first half-step updates the linear and angular positions, followed by the evaluation of local linear and angular accelerations, and finally a second half-step updates the linear and angular positions again. Therefore, it entails only one right hand side evaluation of Eqs. (23, 24), the most computationally expensive operation (see Appendix for details).

This algorithm strikes a balance between computing costs, numerical accuracy and implementation modularity: by foregoing an implicit integration scheme we can incorporate a number of additional physical effects and soft constraints, even though this may come at the expense of computational efficiency.

iii.4 Validation

We first validate our proposed methodology against a number of benchmark problems with analytic solutions and examine the convergence properties of our approach. Three case studies serve to characterize the competition between bending and twisting effects in the context of helical buckling, dynamic stretching of a loaded rod under gravity, and the competition between shearing and bending in the context of a Timoshenko beam. Further validations reported in the Appendix include Euler and Mitchell buckling due to compression or twist, and stretching and twisting vibrations.

iii.4.1 Helical buckling instability

We validate our discrete derivative operators beyond the onset of instability (see Euler and Mitchell buckling tests in the Appendix) for a long straight, isotropic, inextensible, and unshearable rod undergoing bending and twisting. The filament is characterized by the length and by the bending and twist stiffnesses and . The clamped ends of the rod are pulled together in the axial direction with a slack and simultaneously twisted by the angle , as illustrated in Fig. 3a. Under these conditions the filament buckles into a localized helical shape (Fig. 3e).

The nonlinear equilibrium configuration of the rod can be analitycally determined Coyne:1990 (); vanDerHeijden:2000 (); Neukirch:2002 (); vanDerHeijden:2003 () in terms of the total applied slack and twist . We denote the magnitude of the twisting torque and tension acting on both ends and projected on by and , respectively. Their normalized counterparts and can be computed via the ‘semi-finite’ correction approach Neukirch:2002 () by solving the system

Then, the analytical form of can be expressed vanDerHeijden:2003 () as

where is the normalized arc-length . Here we make use of Eq. (III.4.1) to investigate the convergence properties of our solver in the limit of refinement. To compare analytical and numerical solutions, a metric invariant to rotations about is necessary. Following Bergou et al.Bergou:2008 (), we rely on the definition of the envelope

(26)

where is the angular deviation of the tangent from the axial direction , and is the corresponding maximum value along the filament. The envelope relative to the analytical solution of Eq. (III.4.1), and relative to a numerical model of discretization elements can be estimated via finite differences. This allows us to determine the convergence order of the solver by means of the norms , and of the error .

We simulate the problem illustrated in Fig. 3 at different space-time resolutions. The straight rod originally at rest is twisted and compressed at a constant rate during the period . Subsequently, the ends of the rod are held in their final configurations for the period to allow the internal energy to dissipate (according to the model of Section IV.0.1) until the steady state is reached. Simulations are carried out progressively refining the spatial discretization by varying and the time discretization is kept proportional to , as reported in Fig. 3.

As can be seen in Fig. (3)b,c,e the numerical solutions converge to the analytical one with second order in time and space, consistent with our spatial and temporal discretization schemes. Moreover, to further validate the energy conserving properties of the solver, we turn off the internal dissipation (Fig. 3d) and observe that the total energy of the filament is constant after and matches its theoretical value .

Figure 3: Time-space convergence study for localized helical buckling. (a) We consider a rod originally straight whose ends are pulled together in the axial direction with a slack and simultaneously twisted by the angle . (b) Comparison between the analytical envelope function and numerical approximations at different levels of time-space resolution. Here, the time discretization is slaved by the spatial discretization according to s. (c) Norms (black), (blue) and (red) are plotted against the number of discretization elements . (d) Time evolution of the total energy of a rod (, here the energy is computed assuming quadratic functionals, a suitable representation for an inextensible rod) simulated assuming no dissipation (red line) versus the theoretical total energy (black dashed line). (e) Equilibrium rod configuration numerically obtained given the discretization , and assuming dissipation. For all studies, unless specified otherwise, we used the following settings: length  m, twist , slack  m, linear mass density  kg/m, bending stiffness  Nm, twisting stiffness  Nm, shear/stretch matrix  N, bend/twist matrix  Nm, dissipation constant  kg/(ms), radius  m, twisting time  s, relaxation time  s.

iii.4.2 Vertical oscillations under gravity

We consider a system in which a rod hanging from one end and subject to gravity oscillates due to a mass suspended at the other end, and due to its own mass , as depicted in Fig. 4a,d. This system is analogous to a mass-spring oscillator. The static solution is then obtained by integrating the infinitesimal elongations along the spring due to the local load Galloni:1979 (), yielding the total equilibrium extension

(27)

where is the spring constant, is a constant factor, and is the equivalent mass. Thus, the final equilibrium length of the rod reads , with being the rest unstretched length.

The dynamic solution is instead characterized by oscillations of period and by the time varying length of the spring

(28)

In this case, unlike the static solution, the factor depends on the ratio . In fact it can be shown Galloni:1979 () that for , and for .

The analytical results rely on the assumption of being constant in space and time, given a fixed ratio . However, this condition is not met here since is a function of space and time, due to dilatation and mass conservation. Nevertheless, as the Young’s modulus , that is as a soft filament becomes stiff, the constant and our rod model must recover the behavior of the mass-spring oscillator. Indeed, Fig. 4b,c,e,f shows how the proposed numerical method converges to the analytical oscillation period and normalized longitudinal displacement as increases.

Figure 4: Vertical oscillation under gravity. (a,d) We consider a vertical rod of mass clamped at the top and with a mass attached to the free end. Assuming that the rod is stiff enough (i.e. ), it oscillates due to gravity around the equilibrium position , where with a period with for , and for . Therefore, the rod oscillates according to . (a-b) Case with  kg and  kg. (b) By increasing the stiffness , 2, 3, 5, ,  Pa, the simulated oscillations (red lines) approach the analytical solution (dashed black line). (c) Convergence to the analytical solution in the norms (black), (blue) and (red) with , where is the length numerically obtained as a function of . (c-d) Case with  kg and  kg. (e) By increasing the stiffness , 2, 3, 5, , 2,  Pa, the simulated oscillations approach the analytical solution. (f) Convergence to the analytical solution in the norms , and as a function of . For all studies, we used the following settings: gravity  m/s, rod density  kg/m, shear modulus  Pa, shear/stretch matrix  N, bend/twist matrix  Nm, rest length  m, rest cross sectional area  m, number of discretization elements , timestep , dissipation constant .

iii.4.3 Cantilever beam

We consider now the effect of bend and shear simultaneously by validating our numerical methods against the Timoshenko cantilever of Fig. 5a. Timoshenko’s model accounts for bending elasticity, rotary inertia and shear deformations, building on classical beam theories by Rayleigh (bending elasticity and rotary inertia) and Euler-Bernoulli (bending elasticity only). The model captures the behavior of short or composite beams in which shear deformations effectively lower the stiffness of the rod Rao:1995 (); Gere:2001 ().

We consider a beam clamped at one end and subject to the downward force at the free end , as illustrated in Fig. 5a. The static solution for the displacement along the vertical direction of the rod can then be analytically expressed as

(29)

where is the length of the rod, is the constant cross sectional area, is the area second moment of inertia about the axis , and are the Young’s and shear moduli, and is the Timoshenko shear factor for circular sections and accounts for the fact that the shear stress varies over the section Gere:2001 (). Furthermore, the Timoshenko (as well as Rayleigh and Euler-Bernoulli) theory relies on the assumption of small deflections, so that the horizontal coordinate along the direction can be approximated by the arc-length (Fig. 5a and Appendix for further details and derivation), hence the use of in the above equation .

If the shear modulus approaches infinity or if the ratio , then the Timoshenko model in the static case reduces to the Euler-Bernoulli approximation, yielding

(30)

as the shear term of Eq. (29) becomes negligible.

We compare our numerical model with these results by carrying out simulations of the cantilever beam of Fig. 5a in the time-space limit of refinement. As can be noticed in Fig. 5b the discrete solution recovers the Timoshenko one. Therefore, the solver correctly captures the role of shear that reduces the effective stiffness relative to the Euler-Bernoulli solution. Moreover, our approach is shown to converge to the analytical solution in all the norms , , of the error , where is the vertical displacement numerically obtained in the refinement limit.

We note that the norms and exhibit first order convergence, while decays with a slope between first and second order. We attribute these results to the fact that while the Timoshenko solution does not consider axial extension or tension, it does rely on the assumption of small deflections (), therefore effectively producing a dilatation of the rod. On the contrary, our solver does not assume small deflections and does not neglect axial extension, since the third entry of the matrix has the finite value (see Fig. 5 for details). This discrepancy is here empirically observed to decrease the convergence order.

These studies, together with the ones reported in the Appendix, complete the validation of the proposed numerical scheme and demonstrate the accuracy of our methodology in simulating soft filaments in simple settings.

Figure 5: Time-space convergence study for a cantilever beam. (a) We consider the static solution of a beam clamped at one end and subject to the downward force at the free end . (b) Comparison between the Timoshenko analytical (black lines) and numerical (with , red dashed lines) vertical displacements with respect to the initial rod configuration. As a reference we report in blue the corresponding Euler-Bernoulli solution. (c) Norms (black), (blue) and (red) of the error at different levels of time-space resolution are plotted against the number of discretization elements . Here, the time discretization is slaved by the spatial discretization according to seconds. For all studies, we used the following settings: rod density  kg/m, Young’s modulus  Pa, shear modulus  Pa, shear/stretch matrix  N, bend/twist matrix  Nm, downward force  N, rest length  m, rest radius  m, dissipation constant  kg/(ms), simulation time = s.

Iv Including interactions and activity: solid and liquid friction, contact and muscular effects

Motivated by the advancements in the field of soft robotics Kim:2013 (); Haines:2014 (); Park:2016 (), we wish to develop a robust and accurate framework for the characterization and computational design of soft slender structures interacting with complex environments. To this end, we expand the range of applications of our formalism by including additional physical effects, from viscous hydrodynamic forces in the slender-body limit and surface solid friction to self-contact and active muscular activity. As a general strategy, all new external physical interactions are accounted for by lumping their contributions into the external forces and couples and on the right hand side of the linear and angular momentum balance Eqs. (14, 15). On the other hand, all new internal physical and biophysical effects are captured by adding their contributions directly to the internal force and torque resultants before integrating Eqs. (14, 15).

iv.0.1 Dissipation

Real materials are subject to internal friction and viscoelastic losses, which can be modeled by modifying the constitutive relations so that the internal toques and forces of Eqs. (11) become functions of both strain and rate of strain, i.e. and . Keeping track of the strain rates increases computational costs and the memory footprint of the solver. However, for the purpose of purely dissipating energy, a simple alternative option is to employ Rayleigh potentials Torby:1984 (); Audoly:2013 (). In this case viscous forces and torques per unit length are directly computed as linear functions of linear and angular velocities through the constant translational and rotational internal friction coefficients, so that

(31)
(32)

This approach does not model the physical nature of viscoelastic phenomena, although it does dissipate energy, effectively mimicking overall material friction effects. In the context of our numerical investigations, we did not observe any appreciable difference between the two outlined methods, so that, for the sake of simplicity and computational efficiency, we opted for the second one. Throughout the remainder of the text we will then employ Eqs. (31, 32) with a single dissipation constant , therefore assuming .

iv.0.2 Muscular activity

To study limbless biolocomotion on solid substrates and in fluids, we allow our soft filaments to be active, by generating internal forces and torques corresponding to coordinated muscular activity driven, for example, by a central pattern generator Altringham:1990 (); Gazzola:2015 ().

Following the approach detailed in Gazzola:2012 (); Rees:2013 (), we express the muscular activity magnitude as a traveling wave propagating head to tail along the filament

(33)

where is the phase, is time, and are, respectively, the activation period and wavelength. The amplitude of the traveling wave is represented by the cubic B-spline characterized by control points with , as illustrated in Fig. 6. The first and last control points are fixed so that and , therefore assuming the ends of the deforming body to be free. One of the main advantages of the proposed parametrization is that it encompasses a large variety of patterns with a relatively small number of parameters Rees:2013 ().

A given activation mode can be achieved by multiplying the scalar amplitude with the appropriate director. For example, if we wish to study earthworm-like locomotion we may employ a wave of longitudinal dilatation and compression forces, so that

(34)

Similarly, if we wish to investigate a slithering snake characterized by a planar kinematic wave, we may consider a torque activation of the form

(35)

assuming and to be coplanar to the ground. These two contributions are directly added to the internal force and torque resultants.

In the most general case, all deformation modes can be excited by enabling force and torque muscular activity along all directors , and , providing great flexibility in terms of possible gaits.

Figure 6: Muscular activity. Example of muscular activity amplitude profile (solid black line) described by cubic B-spline through control points with . The control points are located along the filament at the positions , and are associated with the amplitude values . The first and last control points are fixed so that and , therefore assuming the ends of the deforming body to be free.

iv.0.3 Self-contact

To prevent the filament from passing through itself, we need to account for self-contact. As a general strategy, we avoid enforcing the presence of boundaries via Lagrangian constraints as their formulation may be cumbersome Riewe:1996 (), impairing the modularity of the numerical solver. We instead resort to calculating forces and torques directly and replacing hard constraints with ‘soft’ displacement-force relations.

Our self-contact model introduces additional forces