# A Discrete Element Method model for frictional fibers

Jérôme Crassous  Univ Rennes, CNRS, IPR (Institut de Physique de Rennes) - UMR
6251, F-35000 Rennes, France<EMAIL_ADDRESS>
###### Abstract

We present a Discrete Element Method algorithm for the simulation of elastic
fibers in frictional contacts. The fibers are modeled as chains of cylindrical
segments connected to each other by springs taking into account elongation,
bending and torsion forces. The frictional contacts between the cylinders are
modeled using a Cundall and Strack model routinely used in granular material
simulations. The physical scales for simulations, the determination and the
tracking of contacts, and the algorithm are discussed. Tests on different
situations involving few or many contact points are presented and compared to
experiments or to theoretical predictions.

##  I introduction

The use of natural or artificial fibers allows to design materials with
original mechanical properties. At the nanometric or micrometric scales,
carbon nanotubes [1] or polymer fibers [2] can be assembled into threads or
networks. At the micrometer and millimeter scales, the frictional forces act
with the elasticity of the fibers to produce a wide variety of materials. The
fibers can just be deposited without any special preparation to form highly
elastic media [3] such as cushions or non-woven fabrics [4]. Textile fibers
can be twisted to produce yarns [5, 6, 7], which are then assembled into cords
[8], woven [9] or knitted fabrics [9, 10]. Cyclic mechanical stresses can form
very compact natural structures [11], and birds also assemble fibers to build
their nests [12, 13]. The contacts between fibers play a fundamental role in
describing the physics of knots, which is a subtle competition between tension
and friction [14, 15], as well as eventual bending of the fibers [16, 17, 18,
19].

Several approaches have been proposed to numerically simulate these
structures. One approach is to use finite element algorithms to discretize the
fibers [20]. This approach allows a complete solution of the elasticity
equations in complex geometries such as nodes [18], but is only possible for
systems with small numbers of contacts. Another approach is to model the
fibers as connected spheres [21] or sphero-cylinders [22] and to use the
discrete element method algorithm widely used for the study of granular
materials. However, the periodic variations of diameter of such fiber may
induce very specific physical properties as interlocked granular chains
stiffening [23].

More realistic approaches are the simulations of fibers as discrete [24] or
continuous [25, 26] cylindrical elastic chains of circular cross-sections. The
non-interpenetration condition between fibers and surfaces, or between fibers,
is then treated as constraints on the displacements. The introduction of
frictional tangential forces in such model has been proposed using methods for
finding forces that match the Coulomb conditions [25, 26, 27]. In those
algorithms, the fibers are moved in order to find the positions of the
surfaces that match the non-penetration of fibers, with forces verifying the
Coulomb condition. Those positions are found using an iterative procedure with
proper regularization of Coulomb law to ensure the convergence towards one
solution verifying the force balance. In the case where many frictional
contacts are present, the problem becomes hyperstatic, and the solution is
expected not to be unique. This is a well known situation in granular material
[28] simulations, and the solutions selected by iterative algorithms are not
well controlled [28], and presumably depend on the algorithm itself. Those
drawbacks are of course of minimal importance in situations where the
indeterminacy in contact forces is absent (hypo- or iso-static problem) such
as in knots with few contacts [29], or if qualitative simulations are needed
as in computer graphic community [30]. In explicit methods, the forces are
obtained directly from the kinematic of the body in contacts. The selection of
one solution of the Coulomb friction forces among many ones is then ensured by
the dynamics of the system. In counterpart, explicit algorithm are usually
slower.

Chains of cylinders with frictional contacts have been first introduced in
Discrete Element Method by Chareyre et al.. These authors used them for the
study of the mechanical properties of granular materials reinforced with
fibers [31, 32], with geotextiles [33], and for the behavior of suspensions of
frictional fibers in viscous flow [34].

This bibliography shows that the modeling of fibers in chains of discrete
elements has been the subject of many studies, but scattered in different
fields. Moreover, the ability of these different models to quantitatively
reproduce the behavior of fibers systems with many frictional contacts has
never been shown. Systems of fibers in frictional interactions are the object
of a growing interest of physicists and mechanics. The object of this study is
to propose to the community a simple Discrete Element Method, easily
reproducible, and founded on the Discrete Element Rod model which include
frictional contacts, and whose capacity to reproduce the behavior of various
frictional fibers is clearly demonstrated.

We will base the model on the theory of elastic chains as proposed by Bergou
et al. [24]. We will keep a formulation with independent elastic constants of
torsion, bending, and torsion, i.e. not linked by a cylindrical beam
elasticity. This will allow to simulate various systems, such as arbitrarily
flexible wires. The contacts will be treated following an approach proposed by
Chareyre et al. [31]. The ingredients of the modeling, as well as the
calculations, will be presented in the simplest possible way so that this
simulation can easily be reproduced by physicists from various fields.

The manuscript is organized in the following way. In the section II, we first
describe the mechanical model of our fibers, including internal elastic forces
and contact forces. The numerical resolution is then detailed in section III,
where we insist on points that are specific compared to DEM simulations of
frictional beads, i.e. the numerical scales that are used, the integration of
displacement, and the search of neighbors. In section IV, we illustrate this
algorithm on various situations including static and dynamics, with few and
many contacts.

##  II Mechanical model of fibers in contact

###  II.1 Description of the fiber

Figure 1: (a) Ensemble of connected point forming the skeleton of the fiber.
(b) Cylinders and spheres forming the shell of the fibers.

Following [24], we model a fiber as an ensemble of $N$ connected points (see
figure 1(a)). Let ${\bf r}_{i}$, with $0\leq i\leq N-1$ be the position of the
point, and ${\bf e}_{i}=({\bf r}_{i+1}-{\bf r}_{i})/\|{\bf r}_{i+1}-{\bf
r}_{i}\|$, with $0\leq i\leq N-2$ the unit vector joining two successive
points. We note $l_{i}=\|{\bf r}_{i+1}-{\bf r}_{i}\|$. The segment joining two
successive points is the generatrix of a cylinder of circular basis of
diameter $d$. In addition, each point ${\bf r}_{i}$ is the center of a sphere
of diameter $d$. So each fiber is a set of $N$ spheres connected by $N-1$
cylindrical segments. A mass $m_{0}$ is assigned to each node of the string,
and a moment of inertia $J$ is assigned to each cylinder.

The kinematic of the deformation is the following. The different nodes of one
fibre may translate, allowing the bending and the stretching of the fibre. The
cylinders joining the different nodes stay straight cylinders and are not bent
when the fiber is deformed. The rectilinear shape allows to determine the
contacts between fibers as contacts between cylinders.The cylinders may rotate
around their axis, allowing the twist of the fibers. The kinematic of the
chain is then determined by the set of $N$ node positions ${\bf r}_{i}$, and
$N-1$ cylinder rotations $\theta_{i}$. To those degrees of freedom, we
associate forces that act on nodes, and torques along the axis of cylinders.
Any system of forces or torques, such as contact forces or elastic forces,
acting on a cylinder will be decomposed as an axial torque and forces on
nodes. This decomposition will be detailed below for elastic twist torques and
contact forces.

###  II.2 Internal elastic forces

The internal elastic forces that we consider in the following are elongation,
flexion and twist forces. The elongational forces are modeled using springs of
stiffness $k_{0}$ with dashpots of damping $\lambda$. The equilibrium length
of the spring is $l_{0}$, and the elongation force exerted by point $i+1$ on
the mass located at ${\bf r}_{i}$ is:

${\bf
f}_{i+1;i}^{(e)}=\bigl{[}k_{0}~{}(l_{i}-l_{0})+\lambda~{}\dot{l_{i}}\bigr{]}{\bf
e}_{i}$ (1)

Each point $i$ is submitted to forces from points $i-1$ and $i+1$ so that
${\bf f}_{i}^{(e)}={\bf f}_{i+1;i}^{(e)}+{\bf f}_{i-1;i}^{(e)}$, excepted the
first $i=0$ and last $i=N-1$ points.

The flexion forces acting on the point $i$ is obtained from the elastic
bending energy $E^{(b)}=\int_{s}(B/2)~{}\kappa^{2}~{}ds$ with $B$ the bending
stiffness of the fiber, and $\kappa$ the curvature. The bending energy of the
discrete fiber is:

$E^{(b)}=\frac{B~{}l_{0}}{2}\sum_{i=1}^{i=N-2}\kappa_{i}^{2}$ (2)

where $\kappa_{i}$ is the curvature at node $i$, and the summation is extended
to all nodes except ending ones. Writing the curvatures $\kappa_{i}$ as
function of nodes positions ${\bf r}_{i}$, the flexion force ${\bf
f}_{i}^{(b)}=-(\partial E^{(b)}/\partial{\bf r}_{i})$ acting on nodes $i$ is
(see Appendix A):

${\bf f}_{i}^{(b)}=-\frac{B}{l_{0}^{3}}\bigl{[}{\bf r}_{i-2}-4{\bf
r}_{i-1}+6{\bf r}_{i}-4{\bf r}_{i+1}+{\bf r}_{i+2}\bigr{]}$ (3)

for $(N-3)\geq i\geq 2$. Expressions of the forces ${\bf f}_{i}^{(b)}$ for
$i<2$ and $i>(N-3)$ are given in Appendix A. The calculation supposes that the
fibers are weakly extended and bent (see Appendix A).

The internal elastic torque is obtained from the twisting energy [24]:
$E^{(t)}=\int_{s}(C/2)~{}\tau^{2}~{}ds$ with $C$ the torsion modulus of the
fiber, and $\tau$ the twist of the fiber. The twist may be written as [35, 36,
37]: $\tau=\tau_{int}+\tau_{s}$. The internal twist $\tau_{int}$ is the twist
of the fiber if simply unbent, whereas $\tau_{s}$ is the torsion of the fibre
centre line. Writing internal twist at node $i$ as
$(\theta_{i}-\theta_{i-1})/l_{0}$, and $\tau_{s,i}$ the torsion of the center
line at node $i$, we obtain the twist energy of the discrete fiber as:

$E^{(t)}=\frac{C~{}l_{0}}{2}\sum_{i=1}^{i=N-2}\bigl{(}\theta_{i}-\theta_{i-1}+l_{0}~{}\tau_{s,i}\bigr{)}^{2}$
(4)

The twist moment acting on segment joining nodes $i$ and $i+1$ is obtained by
differentiating (4) with respect to $\theta_{i}$ (see VI.2):

$\begin{split}{\bf
m}^{(t)}_{i}=\frac{C}{l_{0}}\bigl{[}&(\theta_{i+1}-\theta_{i}+l_{0}~{}\tau_{s,i+1})~{}{\bf
e}_{i+1}\\\ -&(\theta_{i}-\theta_{i-1}-l_{0}~{}\tau_{s,i})~{}{\bf
e}_{i-1}\bigr{]}\end{split}$ (5)

This elastic moment is split into one component $m^{(t)}_{i}={\bf
m}^{(t)}_{i}\cdot{\bf e}_{i}$ of the moment along the axis of the segment
$(i,i+1)$, and into forces acting on nodes (see VI.2).

###  II.3 Contact forces

Figure 2: Contact between two cylinders.

Contact between fibers may occur between segments of cylinders or spheres
belonging to identical or different fibers. The figure 2 shows the contact
between two sections of cylinders. The contact point ${\bf r}_{C}(t)$ is
located on the segment with ending points ${\bf r}_{H_{i}}^{(1)}$ and ${\bf
r}_{H_{i}}^{(2)}$ on axis cylinders which minimizes the distance between axis.
This segment is unique if the axis are not parallel. The determination of this
segment will be detailed in section III.3. Let $d(t)=\|{\bf
r}_{H_{i}}^{(2)}-{\bf r}_{H_{i}}^{(1)}\|$ this minimal distance, and ${\bf
n}(t)$ the normal unitary vector. We note $\delta(t)=(d_{1}+d_{2})/2-d(t)$ the
interpenetration of the two cylinders, and ${\bf r}_{C}(t)={{\bf
r}_{H_{i}}^{(1)}}+(r_{1}-\delta/2)~{}{\bf n}$ the contact point. We use the
Cundall-Strack model for the contact force [38]. The normal contact force
exerted by cylinder 1 on cylinder 2 is modeled as a spring-dashpot system:

${\bf
f}^{(c)}_{n}=-\bigl{[}k_{n}~{}\delta+\lambda_{n}~{}\dot{\delta}\bigr{]}{\bf
n}$ (6)

with $k_{n}$ the contact stiffness between the cylinders, and $\lambda_{n}$
the contact damping. This contact law is a simplified version of the elastic
contact force between 2 cylinders [39] which varies non-linearly with the
interpenetration $f^{(c)}_{n}\sim\delta^{3/2}$, and which depends on the angle
between cylinder axis. The tangential contact force is a Coulomb-Force:

${\bf f}^{(c)}_{t}=-Min\bigl{[}k_{t}~{}u_{t};\mu
k_{n}\delta\bigr{]}~{}\frac{{\bf u}_{t}}{u_{t}}$ (7)

where $k_{t}$ is the tangential stiffness, ${\bf u}_{t}$ the tangential
displacement, and $\mu$ the microscopic friction coefficient. The tangential
displacement ${\bf u}_{t}$ is initialized to $0$ when the contact is first
formed, and is evolves with time in the following way:

First, since the tangential displacement is expressed in the fixed frame,
${\bf u}_{t}$ is first rotated to take into account the rotation of the normal
vector. Lets $d\beta$ the angle between the normal at times $t-dt$ and $t$:
${\bf n}(t-dt)$ and ${\bf n}(t)$, and ${\bf k}=\bigl{[}{\bf n}(t-dt)\times{\bf
n}(t)\bigr{]}/\rvert{\bf n}(t-dt)\times{\bf n}(t)\lvert$ the axis rotation. We
name ${\bf u}^{rot}_{t}(t-dt)$ this rotated tangential displacement.

Then, the displacement is integrated as:

${\bf u}_{t}(t)={\bf u}^{rot}_{t}(t-dt)+({\bf v}^{(2)}-{\bf v}^{(1)})dt$ (8)

where ${\bf v}^{(1)}$ (and similar for ${\bf v}^{(2)}$) is the velocity of the
point of the cylinder $(1)$ coinciding with the contact point $C$. The
velocity is ${\bf v}^{(1)}=\dot{{\bf
r}}_{i}^{(1)}+{\bf\Omega}_{i}^{(1)}\times({\bf r}_{C}-{\bf r}_{i}^{(1)})$ with
${\bf\Omega}_{i}^{(1)}$ the rotational velocity of the segment $i$ of fiber
$1$. The rotation vector is separated into an axial and non-axial components
as: ${\bf\Omega}^{(1)}={\bf\Omega}^{(1)}_{\perp}+\dot{\theta}_{i}{\bf e}_{i}$,
where the non-axial component is
${\bf\Omega}_{i,\perp}^{(1)}=\frac{1}{l_{i}}{\bf e}_{i}\times(\dot{{\bf
r}}_{i+1}^{(1)}-\dot{{\bf r}}_{i}^{(1)})$.

Finally, the normal component $({\bf u}_{t}\cdot{\bf n})~{}{\bf n}$ is
removed. If $k_{t}u_{t}>\mu k_{n}\delta$, then the tangential displacement is
renormalized such that $u_{t}=\mu k_{n}\delta/k_{t}$.

The contact force ${\bf f}^{(c)}$ is then expressed as a system of forces
${\bf f}_{i}^{(c)}$ and ${\bf f}_{i+1}^{(c)}$ applied on nodes $i$ and $i+1$,
and a moment $m_{i}{\bf e}_{i}$ acting the cylinder connecting those nodes.
The conservation of the resultant and of the moment of contact force implies
that:

$\displaystyle{\bf f}_{i}^{(c)}+{\bf f}_{i+1}^{(c)}$ $\displaystyle={\bf
f}^{(c)}$ (9a) $\displaystyle({\bf r}_{i+1}-{\bf r}_{i})\times{\bf
f}_{i+1}+m_{i}{\bf e}_{i}$ $\displaystyle=({\bf r}_{C}-{\bf r}_{i})\times{\bf
f}^{(c)}$ (9b)

A possible choice for forces and moment is (see VI.3):

$\displaystyle m_{i}^{(c)}$ $\displaystyle=\bigl{[}({\bf r}_{C}-{\bf
r}_{i})\times{\bf f}^{(c)})\bigr{]}\cdot{\bf e}_{i}$ (10a) $\displaystyle{\bf
f}_{i}^{(c)}$ $\displaystyle=(1-s_{i}){\bf f}^{(c)}+\frac{R}{l_{i}}({\bf
f}^{(c)}\cdot{\bf e}_{i}){\bf n}$ (10b) $\displaystyle{\bf f}_{i+1}^{(c)}$
$\displaystyle=s_{i}{\bf f}^{(c)}-\frac{R}{l_{i}}({\bf f}^{(c)}\cdot{\bf
e}_{i}){\bf n}$ (10c)

It should be noticed that (9) does not set all the components of ${\bf
f}_{i}^{(c)}$ and ${\bf f}_{i+1}^{(c)}$, and that a supplementary conditions
expressed in appendix VI.3 must be added to obtain (10).

If the contact between two fibers involve one cylindrical segment of the fiber
and one sphere, or two spheres, the contact point is calculated accordingly to
the type of the surfaces in contact. The translation velocity of the sphere is
the velocity of the node. The rotation velocity of the sphere at node $i$ is
the rotation velocity of the cylinder joining nodes $i$ with $i+1$.

###  II.4 Miscellaneous forces.

In addition misc extra forces may be added. A global viscous damping force
${\bf f}_{i}^{(v)}=-\lambda_{v}~{}\dot{{\bf r}_{i}}$ may be added. It is
useful to damp transverse motion of fibers. Indeed, our mechanical model of
fiber does not include any dissipation for motion perpendicular to fiber axis
if there is no contacts. Volumetric forces such as gravity forces ${\bf
f}_{i}^{(g)}=m_{0}~{}{\bf g}$ with ${\bf g}$ the gravity field may be also
added. Other external forces may applied to fibers such pre-tension at ends of
fibers.

##  III Numerical implementation

###  III.1 Integration of equation of motions

The dynamical equations of motions writes as:

$\displaystyle M_{0}~{}\ddot{{\bf r}_{i}}$ $\displaystyle={\bf
f}_{i}^{(e)}+{\bf f}_{i}^{(b)}+{\bf f}_{i}^{(c)}+{\bf f}_{i}^{(v)}+{\bf
f}_{i}^{(g)}$ (11a) $\displaystyle J_{z}~{}\ddot{\theta_{i}}$
$\displaystyle=m_{i}^{(t)}+m_{i}^{(c)}$ (11b)

The second equation described the rotation of cylinder segment around its
axis. We did not consider in (11b) any elastic torque due to torsion of the
fiber, and the fiber is free to rotate around node ${\bf r}_{i}$. The
dynamical equations are integrated using a standard second-order Verlet
algorithm [40].

###  III.2 Physical parameters for simulations

####  III.2.1 Physical scales

We first define mass, length and stiffness scale for the simulation. The mass
scale $m_{0}$ is the point mass of nodes, and the length scale $l_{0}$ is the
equilibrium length of each segment, and the stiffness scale $k_{0}$ is the
elongation stiffness of spring. If all the fibers do not have identical
physical properties, those scales are chosen from the fibers of smallest
radius. For every physical quantities $x$, with a physical scale $x_{0}$, we
note the non-dimensional quantity as $x^{*}=x/x_{0}$.

The time scale is $t_{0}=(m_{0}/k_{0})^{1/2}$. For fibers of diameters $r=d/2$
made of an elastic material (Young modulus $E$, Poisson coefficient $\nu$) of
density $\rho$, we have $k_{0}=E\pi r^{2}/l_{0}$, $m_{0}=\rho\pi r^{2}l_{0}$,
and then $t_{0}=l_{0}~{}(\rho/E)^{1/2}$. The time scale $t_{0}$ is then the
time of propagation of compression waves through one segment of the fiber. The
force scale $f_{0}=k_{0}~{}l_{0}=E\pi r^{2}$ is the force that extend a
hypothetical perfectly elastic fiber by $100\%$.

####  III.2.2 Elastic forces and damping

When submitted to a traction force $f$, the relative expansion of the fibers
is $f/f_{0}=f^{*}$. It follows that if we want to stay in the limit of small
extension, we should keep $f^{*}\ll 1$. In practice the simulations are done
with $f^{*}\sim 10^{-5}-10^{-3}$. It should be noted that if $f^{*}$ is too
small, the propagation of transverse waves is very slow when no bending forces
are present. Indeed, the velocity of transverse wave $v_{t}$ in a string of
linear density $\rho_{l}$ under a tension $f$ is $v_{t}=(f/\rho_{l})^{1/2}$.
With $\rho_{l}=m_{0}/l_{0}$, we have $\rho_{l}^{*}=1$, and the non-dimensional
speed of transverse wave is
$v_{t}^{*}=(f^{*}/\rho_{l}^{*})^{1/2}=(f^{*})^{1/2}$ when no bending stiffness
are present.

The non-dimensional bending stiffness is $B^{*}=B/k_{0}l_{0}^{3}$. For an
elastic fiber as consider in III.2.1, we have $B=E\pi r^{4}/4$, and then
$B^{*}=(r^{*})^{2}/4$. Similarly, the non-dimensional torsional modulus is
$C^{*}=C/k_{0}l_{0}^{3}$. For an elastic fiber of radius $r$, we have $C=E\pi
r^{4}/2(1+\nu)$, and then $C^{*}=(r^{*})^{2}/2(1+\nu)$.

The longitudinal damping $\lambda$ is chosen to avoid compression waves that
travel continuously through the fibers, needing very long time to return to
equilibrium. We take $\lambda\sim(k_{0}m_{0})^{1/2}$, and then
$\lambda^{*}\sim 1$ for this.

####  III.2.3 Contact force

The value of the contact stiffness is fixed from a linearization of the
Hertzian contact between two elastic cylinders. If two cylinders of radius
$r$, with perpendicular axis are in contact, the problem is equivalent to the
the contact between a sphere of radius $r$ and a plane, and the normal force
is $f_{n}=(4/3)~{}E_{eff}~{}r^{1/2}~{}\delta^{3/2}$, with
$E_{eff}=E/(1-\nu^{2})$, $\nu$ being the Poisson ratio of the material. For
doing the linearization, we arbitrary set that the elastic energy of the
Hertzian contact $\sim E_{eff}~{}r^{1/2}~{}\delta^{5/2}$ is equal elastic
energy $k_{n}\delta^{2}/2$ of the spring for a normal force $f$ which is of
the order or the traction force that we applied on the fibers. Dropping
numerical factor of order $1$, we obtain $k_{n}=E^{2/3}~{}f^{1/3}~{}r^{1/3}$.
The non-dimensional stiffness may then be obtain as:

$k_{n}^{*}=\frac{(f^{*})^{1/3}}{r^{*}}$ (12)

where we again dropped constant term. $f^{*}$ is the typical non-dimensional
force (i.e. the non-dimensional traction applied to the fibers). This value of
$k_{n}^{*}$ is a reasonable choice for modeling contact, but evidently
different values may be set. In practice, since the tension is of order
$f^{*}\sim 10^{-5}-10^{-3}$, and typical radius are $r^{*}\sim 10^{-1}$, we
have $k_{n}^{*}\sim 1$. For sake of simplicity, the tangential stiffness is
taken as $k_{t}^{*}=k_{n}^{*}$.

Some damping of the normal force $\lambda_{n}$ may be introduced. We took
$\lambda_{n}^{*}\sim 1$ for rapid relaxation of oscillating motion of contact.

####  III.2.4 Time scale for simulation

The time step $dt$ for simulation is chosen such that the dynamic of length
relaxation and of contact establishment is correctly described. The length of
segment relaxes on a time scale $\sim(m_{0}/k_{0})^{1/2}=t_{0}$, whereas the
time scale for a contact to establish is
$\sim(m_{0}/k_{n})^{1/2}=t_{0}~{}(k_{0}/k_{n})^{1/2}$. The time step is chosen
as $dt=Min\bigl{[}t_{0};t_{0}~{}(k_{0}/k_{n})^{1/2}\bigr{]}/10$, leading to:

$dt^{*}=\frac{1}{10}~{}Min\bigl{[}1;(k_{n}^{*})^{-1/2}\bigr{]}$ (13)

such that both relaxat ions occur on at least $10$ time steps. In practice,
since $k_{n}^{*}\sim 1$, we take $dt^{*}=0.1$. For a given set of parameters,
it is checked that results are unchanged if time steps are divided by a factor
$2$.

###  III.3 Computation of contact points

The Discrete Element Method is mainly used in assemblies of spherical
particles. Due to the anisotropic shape of the segments, our algorithm for the
determination of the contact points has some particularities compared to
sphere-sphere contact that we discuss in this section.

####  III.3.1 Distance between fibers

Figure 3: (a) 2D view of two sets composed of one sphere and one cylinder. (b)
Motions of an external $C_{ext}$ and internal $C_{int}$ contact points at the
junction between two cylinders.

The distance between fibers is calculated in the following way. We first
consider a segment as a set composed of a sphere and a part of cylinder as
shown on figure 3(a). We first calculate the distance between the two parts of
cylinders following the method described in Appendix VI.4. If contact does not
occur along two cylinders, contact between spheres and cylinders are searched,
and finally between the two spheres. The hull of the fiber is therefore
composed of the external surface of the cylinders and of the spheres as shown
ib figure 3(b). The starting and the ending of fibers are finished by spheres.

####  III.3.2 Integration of displacement of contact point

The contact point is followed continuously during the motion of the fibers.
This may be done easily as long as the contact point between one segment and
one fiber is unique as in example the contact point $C_{ext}$ of 3(b). In this
case, the displacement of the contact point is continuously integrated along
the motion. In some case, two contact points may exist simultaneously as the
two points as in example the contact points $C_{int}^{(a)}$ and
$C_{int}^{(b)}$ of 3(b). When the contact at point $C_{int}^{(b)}$ occurs, its
tangential displacement is initially set to $0$ (as every new contact), and
this lower the tangential force. Since the fibers are weakly bend with $r\ll
l_{0}$, we expect that the number of such contacts are very small compared to
the total number of contact, producing very negligible errors. A possible
refinement may be to interpolate the two contact points as a single one,
allowing a continuous integration of displacement.

####  III.3.3 Neighbor search method

Figure 4: (a-b): Two possible choices for the size of the cells into which
collisions between fibers mat be searched: (a) size of the cell scales as the
length of segments; (b) size of the cell scales as the radius of segments. (c)
Simplification arising from the fact that the segments belonging to each
fibers are connected.

The search for contacts between discrete objects can significantly increase
the computation time of DEM algorithms. In our case, the algorithm for
measuring the distance between cylinders is slightly more complex than for
spheres, further increasing the computation time of collisions. Several
strategies are possible to significantly improve the computation time of the
collisions. They are based on the use of neighbor list (Verlet list) or on the
partition of the system in boxes (Linked Cell Method). We discuss here the
problem arising when using strongly anisotropic objects. In linked cell
method, the particles are assigned in cells, and the list of particles is each
cell is updated periodically. The collisions are searched only for particles
within the same or the neighboring cells. This strategy is very effective for
approximately monodisperse spheres. In the case of polydisperse spheres, the
size of the cell must be a multiple of the size of the largest particles, so
that the number of particles per box increases. As a consequence, the
computation time grows rapidly with the polydispersity as shown by Luding et
al [41]. The problem is very similar for strongly anistropic particles such
fibers, or segments of fibers. The figure 4(a) shows an assembly of fibers
with segments of size $l_{0}$. If collisions between segments are searched
within one or neighboring cells, the size of the cell should be $\sim
2~{}l_{0}$. For segments of section $\sim 4~{}r^{2}$, the number of segments
in each cell is $\sim 2~{}(l_{0}/r)^{2}$ for dense $3D$ system. Since
$l_{0}/r\gg 1$, sorting particles in cell of size $\sim l_{0}$ is not
efficient. A more convenient way to define cell may be considered. It
consists, as shown 4(b), of replacing segments by fictitious spheres of radius
$r$ inside each segment of length $l_{0}$, and to consider cells of size $\sim
4~{}r$. In this case for a system of $N_{f}$ fibers of $N$ segments each, the
total numbers of fictitious spheres is $\sim N_{f}N(l_{0}/2r)$. However those
two methods do not use the fact that different segments of one fiber are
linked together. Taking advantage of this knowledge may significantly speed up
the search of neighboring. The figure 4(c) shows two fibers, and we search
contact between segment $i$ of fiber $1$, with fiber $2$, by increasing $j$.
For a segment $j$, we calculate the distance $d(i,j)$. If this distance is
larger that $2r$ there is no contact, and we are sure that there is no contact
between the two fibers for $|j^{\prime}-j|\leq d(i,j)-2r$. So the next segment
where we need to search contact verifies $j^{\prime}>j+d(i,j)-2r$.

The optimal strategy to find contacts is expected dependent on the type of
fiber under studies. In case of fibers with numerous segments, taking
advantage of the constraint that the segment are linked as depicted in 4(c) is
presumably the better. At the opposite, in the case of an assembly of very
short fibers, such as a packing of one-segment needles, use of cell as 4(b)
should be preferred. The further study of such optimization is outside the
scope of this study.

##  IV illustration experiments.

The program has been tested on various simple geometries in order to check the
consistency with the theory, to verify the numerical stability of the
algorithm, and test the numerical precision. Those configurations were the
rolling or sliding of a cylinder on a inclined plane, the velocity of
transverse waves of a string, the static flexion of a fiber loaded at
extremity by a point force, the catenary shape of a massive string under
gravity. We present in the following four more complex situations. If not
otherwise specified, the simulation parameters are: time step $dt^{*}=0.1$,
internal damping $\lambda^{*}=2.8$, contact stiffness $k_{n}^{*}=k_{t}^{*}=1$,
contact damping $\lambda_{n}^{*}=1$, global viscous damping
$\lambda_{v}^{*}=0.001$, inertia momentum $J^{*}=m^{*}{r^{*}}^{2}/2$
(homogeneous cylinder).

###  IV.1 Elastic rods without contacts

The elastic rod model has been already tested in misc situations that do not
involve frictional contacts [24]. The test examples presented here are just
for checking the approximations used in II.2 and VI.

Figure 5: Force applied to a beam as a function of its deflection. (a)
Geometry. (b) Symbols: simulation results. Plain line: theoretical solution.
Dashed line: force in the small deflection limit $FL^{2}/B=3\delta/L$. (c)
Convergence: $\delta_{N}$ is the end deflection for a beam with $N$ nodes
($(N-1)$ cylinders), and $\delta_{th}$ the theoretical deflection. Applied
forces is $FL^{2}/B=10$. Symbols are relative error, and dashed line is a
$N^{-2}$ decay.

The first example is the deformation of a clamped elastic rod ($N=100$,
$B^{*}=0.1$) submitted to a point force applied at one end (see figure 5.a).
The clamping is imposed by fixing the first and second node of the rod. The
free rod length $L$ is then the number of cylinders $N-1$ minus one:
$L=N-2=98$. Results for different values of applied forces are shown on 5.b.
Those results may be compared to the deflection of an non-extensible rod. At
small deflections $\delta\ll L$, $FL^{2}/B\simeq 3\delta/L$. At large
deflections $\delta\sim L$ simulations agree well with the analytical solution
of Bisshopp et al. [42]. Simulations with beams made with different $N$ show
that the solution obtained with the discrete beam converges towards analytical
solution as $\sim N^{-2}$ (see fig.5.c). It may be noticed that since the
maximum force applied in those simulations are of order $F\sim 10~{}B/L^{2}$,
the maximum non-dimensional force $F^{*}\sim 10^{-4}<<1$, so that the beam
stretching is negligible.

Figure 6: Buckling of an elastic rod. (a) Geometry. (b) Torque applied at ends
at the torsional buckling threshold as a function of the rod length. Symbol:
simulation results with $P=0$. Line: $M_{t}^{(b)}=2\pi B/L$.

The second example is the buckling of a rod submitted to a compression and
applied torque at its ends (see fig.6.a). A numeric rod ($B^{*}=0.1$) is
submitted to a torque $M_{t}$ at its ends. The displacement of the ends
perpendicularly to the axis of the rod are blocked, and no compression forces
$P$ are applied. The torque is slowly increased until buckling of the beam
occurs. The buckling threshold is determined by measuring the displacement of
the ends along the axis of the rod. Those displacements are initially
negligible, and suddenly increases as buckling occurs. fig.6.b shows the
buckling torque $M_{t}^{(b)}$ as a function of the rod length. Stability
analysis of twisted rods leads to[43]: $M_{t}^{(b)}=2\pi B/L$. As shown on
fig.6.b, the numerical results are in correct agreement with this theoretical
law.

###  IV.2 Static without flexion : capstan

Figure 7: Tension in a rolled string around a cylinder: $T^{*}$ is the tension
in the string, and $\theta$ is the rolling angle.Circles are symbol, plain
line is an exponential fit. See for simulations parameters Inset: Schematic of
the experiment.

We simulate the tension along a string which is rolled up around a cylinder.
For this, we prepare a infinitely flexible spring ($B^{*}=0$, $N=200$,
$r^{*}=0.1$) which makes $5$ turns around a cylinder ($R^{*}=5$). The cylinder
had a huge mass and moment of inertia to prevent any motion. The friction
coefficient is $\mu=0.2$. We first apply an equal tension
$T_{1}^{*}=T_{2}^{*}=0.01$, with opposite directions, to the two ends of the
string. We let the system to reach equilibrium. Then, we slowly decreases
$T_{2}^{*}$ while keeping $T_{1}^{*}=0.01$. For a threshold value of
$T_{2}^{*}$, the sliding of the string occurs. We measure the tension in the
string using (1) at the onset of sliding. The fig.7 shows the decrease of the
tension $T^{*}$ along the string as a function of
$\theta=(s^{*}-s_{0}^{*})/(R^{*}+r^{*})$, with $s^{*}$ the abscissa along the
curve, and $s_{0}^{*}$ the abscissa of first contact contact point. The
solution of capstan problem with a finite thickness rod predicts that [44]:
$T^{*}/T_{1}^{*}=\exp(-\mu~{}\theta)$, which is the observed behavior on
figure 7. The measured decay is $\mu=0.198$ in agreement with the imposed
value $\mu=0.2$.

###  IV.3 Static with flexion : elastic knots

Figure 8: (a) Snapshot of a ($3_{1}$) knot. For seek of clarity, illustration
is made with $r^{*}=0.2$. (b) Tension as a function of
$\varepsilon=\sqrt{r/R}$ for frictionless and frictional strings. Symbols are
numerical data, and lines are theoretical expressions given by equation (14).

We consider the mechanical response of an elastic rod with an open knot. An
elastic fiber of length $L$, with a circular section of radius $r$, and
bending modulus $B$ is bent in an open trefoil knot ($3_{1}$). We then apply a
tension $T$ to the ends of the fibers. This experimental situation has been
addressed by Audoly et al. [16, 17]. When the tension is weak, the loop radius
$R$ is very large compared to $r$. In this limit, authors found analytical
solutions for the shape of the knot, either in the frictionless case, but also
for weak friction $\mu\ll 1$. This knot has been simulated very recently by
Choi et al. using an discrete rod model with an implicit solver for the
contact force [29].

We simulate numerically such knot by considering a flexible spring
($r^{*}=0.1$, $B^{*}=(r^{*})^{2}/4=2.5~{}10^{-3}$, $N=500$,
$\lambda_{v}^{*}=4.10^{-4}$) as shown on figure 8(a). We first knot the fiber
by setting $\mu=0$ and applying a tension $\pm~{}T^{*}{\bf e}_{z}$ at ends.
After this preparation stage, we set $\mu$ to its actual value, and we
increase or decrease $T^{*}$ depending if we tighten or lossen the knot. When
the knot begins to move, we measure the radius of curvature of the loop as
$R=<\|d{\bf t}_{i}/ds\|^{-1}>$, where $d{\bf t}_{i}/ds={\bf e}_{i+1}-{\bf
e}_{i}$ is the derivative of the tangent vector, and the average $<~{}>$ is
over all segments in the loop which are at a distance of at least one segment
from any contact point. Following [16, 17], we introduce
$\varepsilon=\sqrt{r/R}$. The figure 8(b) shows the tension $T^{*}$ as a
function of $\varepsilon$ for frictionless ($\mu=0$), and frictional
($\mu=0.1$) loosening and opening knots. The analytical solutions in the
limits $\varepsilon\ll 1$ and $\mu\ll 1$ are [16, 17]:

$\frac{Tr^{2}}{B}=\frac{\varepsilon^{4}}{2}\pm\mu\sigma\varepsilon^{3}$ (14)

where the sign $\pm$ depends if the knot is tightened $(+)$ or loosened $(-)$,
and $\sigma$ is a numerical constant which is $\sigma\simeq 0.492$ for trefoil
knot. As shown on figure 8(b), the numerical data agrees correctly with the
analytical one. In the frictionless case, we may observe deviations from the
scaling $T\sim\varepsilon^{4}$ when $\varepsilon\gtrsim 0.15$. Two possible
sources of deviations may be identified. First, the equation (14) is obtained
in the limit $\varepsilon\ll 1$, and deviations may arise from high order
$\varepsilon$ terms in equation (14). Second, for $\varepsilon\gtrsim 0.15$,
we have $R^{*}=R/l_{0}=r^{*}/\varepsilon^{2}\sim 4$, so that the
discretization of loop may then be an issue. The discret nature of the rod may
be clearly identified on numerical data form $\mu=0.1$ loosening, where some
steps in $\varepsilon$ are visible. For the frictional case, the model (14)
slightly underestimates the role of friction compared to numerical
simulations. It may be due to some departure from the hypothesis $\mu\ll 1$
which is used to obtain (14).

###  IV.4 Impact: falling chain.

Figure 9: (a) Experimental snapshots of the impact of a metallic chain on a
fixed perpendicular cylinder of radius $R=10~{}mm$. The perimeter of the
cylinder is underlined in red. Time $t=0$ is defined as the first contact
time. Chain length $L=190~{}mm$, mass $m=8.5~{}g$. (b) Simulated impact. See
text for physical parameters of the simulation.

We consider the dynamics of the impact of a metallic chain on a cylindrical
obstacle. We restrict this analysis to a qualitative analysis. A metallic
chain (length $L=190~{}mm$, mass $m=8.5~{}g$) is held at its extremities by
hands. The chain is released and its fall is recorded with a fast camera
operating at $200~{}fps$. The figure 9(a) show some snapshots of the impact.
The chain is simulated as a infinitely flexible spring $B^{*}=0$. We set the
length scale to $l_{0}=2~{}mm$, and $N=95$, so that $L=N~{}l_{0}$. The choice
of the time scales may be done in the following way. We want to simulate a
non-extensible chain, so we need that the non-dimensional typical force is
$<<1$. The gravity force is $F_{g}=Nm_{0}g$, with $m_{0}$ the mass scale of
one segment, and $g$ the gravity. The non-dimensional gravity force is then
$F_{g}^{*}=F_{g}/k_{0}l_{0}=Ng/l_{0}t_{0}^{-2}=Ng^{*}$. We take
$g^{*}=4.9~{}10^{-5}$ so that $N~{}g^{*}\simeq 5.~{}10^{-3}\ll 1$. This sets
the time scale $t_{0}=0.1~{}ms$. It should be noted that in the limit of a
non-extensible chain, the mass scale does not need to be specified. Other
parameters are $dt^{*}=0.1$, $R^{*}=5$, $r^{*}=0.1$, $\mu=0.1$,
$k_{n}^{*}=k_{t}^{*}=1$. Figure 9(b) shows the results of the simulations
which qualitatively agree with the experiments. We may remark that the
behaviour of the experimental chain is not symmetric in compression and in
extension (nearly infinite stiffness in extension, zero stiffness in
compression), whereas the numerical chain is symmetric (same stiffness in
compression and in extension). However, in impact experiment, the chain is
always in tension, and the lack of symmetry does not have importance.

###  IV.5 Multiple fibers: a yarn model.

Figure 10: (a) Assembly of of initially straight fibers. (b) Thread of fiber
after a torque is applied at ends. (c) Separation of the fibers due to applied
forces. (d) Force ratio necessary to separate the slivers as a function of the
twist angles. Red disks: $\mu=1$, blue triangles: $\mu=0.5$, green squares:
$\mu=0.2$. Lines are for guidelines. (e) Same data as (d) plotted as a
function of $\mathcal{H}=\mu\theta^{2}R/L$. Dotted line is
$0.75~{}\mu\theta^{2}R/L$. Simulation parameters are $N_{f}=20$, $N=30$,
$r^{*}=0.1$, $B^{*}=0$. For clarity, the sub-figures (a-c) are enlarged by a
factor 6 perpendicularly to $z$-axis.

In a recent study, Seguin et al. considered the situation of a staple yarn
made of twisted totally flexible fibers [7]. We present in this section some
numerical details about this simulation. The yarn is made of an assembly of
$N_{f}$ identical fibers of $N$ segments initially parallel to an axis $z$
(see figure 10(a). Their positions $(x_{i};y_{i})$ in the plane perpendicular
to the $z$ axis, with $1\leq i\leq N_{f}$ are the positions of a packing of
disks in 2D obtained from a separate simulation.

In a first phase of the simulation the fibers are twisted. The fiber are
submitted to a tension $T^{*}=10^{-4}$ along $z$, applied at both ends. A
torque $C^{*}{\bf e}_{z}$ is applied to both ends of the assembly of fibers.
For this, each fiber $i$ with $1\leq i\leq N_{f}$ is submitted at both ends
$j=0$ and $j=N$ to an external shear force:

${\bf\tau}_{i}(j)=\pm\frac{C^{*}}{\sum_{i}r_{i}^{2}(j)}\bigl{[}{\bf
e}_{z}\times{\bf r}_{i}(j)\bigr{]}$ (15)

where sign is $-$ for $j=0$, and $+$ for $j=N$ ends. The torque is gradually
increased until it reaches its target value and the shear forces are updated
at each time step. Under the action of this torque, the fibers twist and
becomes approximately helicoidal as shown on figure 10(b). During this
preparation, the friction coefficient is set to a low value $\mu=0.05$. This
is important in order to obtain a regular pitch along the thread. Indeed,
since the yarn is twisted by the application of torques at ends, the presence
of an important friction between fibers has the effect of concentrating the
twist near ends, with a central zone of low twist. This behavior is also
observed experimentally [7] if the twist is not homogenized along the yarn.
The duration of this preparation stage is $t^{*}=5.10^{5}$, and the total
twist $\theta$ is measured at the end of this phase.

In a second phase the fibers are separated. The friction is first set at its
target value. The fibers are randomly partitioned is two set -up- and -down-.
The tension of the up-fibers are multiplied by a factor $f>1$ at the up-
extremities: $T_{up}(j=0)=T^{*}$ and $T_{up}(j=N)=f~{}T^{*}$. Symmetrically,
$T_{down}(j=0)=f~{}T^{*}$ and $T_{down}(j=N)=T^{*}$. The factor is $f=1$ at
the beginning of the separating stage and is increased at a fixed rate
$(\Delta f/\Delta t^{*})=2.10^{-6}$. During this phase, the torque is kept
constant. The difference between the average positions of the -up and -down
fibers is measured. This difference stays constant, until a threshold value of
$f$ where the two slivers of fiber separates (see figure 10(c)). A mechanical
model of this problem developed in [7], show that the force necessary to
separate the two slivers is $\ln(1+f)\simeq 0.75~{}\mu\theta^{2}R/L$ which is
the behavior that is observed on figure 10(e).

##  V Conclusion

We have described a discrete element mechanics algorithm for the simulation of
flexible and frictional fibers. This algorithm is similar to the DEM type
algorithms widely used for the study of granular materials. The difference
arises from the type of surfaces in contact (cylinders and not spheres) and
from the elastic forces between the cylinders which are connected to form a
fiber. The algorithm has been tested on various configurations that can be
compared to experiments or to theoretical models.

The assumptions and approximations used to design this algorithm are quite
limited. The low bending assumption is not very compelling for many
applications, but could eventually be minimized by a finer discretization of
the fiber. The simplification of the Hertzian elastic contact law between the
cylindrical segments by a linear spring has probably a very small impact on
the modeling of real systems. An extension to non-linear contact laws should
not be a problem. Finally, the discretization of the fiber generates a
discontinuity of the displacement for some contact points at the passage
between successive segments of a fiber. A priori the number of such jumps is
negligible compared to the total number of contacts for thin and weakly bent
fibers, and this should not be an issue for simulations of real systems.

The main difference between this algorithm and those previously described to
simulate elastic fibers lies in the level of simplification of the mechanical
problem. Simulations of fibers with finite element algorithms are certainly of
high accuracy but can only simulate small systems. Implicit algorithms are
probably faster, but the indermination of forces in multi-contact cases is not
resolved by the dynamics of the system. The use of a DEM algorithm is a
compromise that allows to consider relatively complex assemblies of fibers and
that correctly handles the multiplicity of equilibrium solutions.

The potential applications of this algorithm are obviously multiple. The study
of complex knots between fibers of ropes, with or without bending energy is
possible. The mechanical response of fiber clusters in nests, cushions, or in
rigid needle stacks are also possible. For these studies, the contact search
should be optimized according to the aspect ratio of the fibers and the
geometry of the packing. The simulation of knitted or woven fabrics can also
be considered. For this, large systems can be simulated, but the introduction
of periodic boundary conditions should be more suitable. Finally, systems
mixing fibers and grains for the study of soils reinforced by fibers or roots
are also possible applications of this work.

###### Acknowledgements.

The author would like to thank Antoine Seguin and Sean McNamara for
discussions and careful reading of the manuscript, and Laurent Courbin for his
help in experiments on falling chains.

##  VI Appendices

###  VI.1 Flexion forces

The curvature $\kappa_{i}$ at a node $N-2\geq i\geq 1$ is first expressed as a
function of the positions of nodes ${\bf r}_{i-1}$, ${\bf r}_{i}$ and ${\bf
r}_{i+1}$. The radius $R_{i}=1/\kappa_{i}$ of the circle joining those three
points may be expressed as a function of the surface $S_{i}$ and the perimeter
$p_{i}$ of the triangle with vertices $({\bf r}_{i-1},{\bf r}_{i},{\bf
r}_{i+1})$ using the Heron formula. After elementary calculus, we obtain:

$\kappa_{i}^{2}=4~{}\frac{{\bf l}_{i-1}^{2}{\bf l}_{i}^{2}-({\bf
l}_{i-1}\cdot{\bf l}_{i})^{2}}{{\bf l}_{i-1}^{2}~{}{\bf l}_{i}^{2}~{}({\bf
l}_{i-1}+{\bf l}_{i})^{2}}$ (16)

where we noted ${\bf l}_{i}={\bf r}_{i}-{\bf r}_{i-1}$. The bending energy is

$E^{(b)}=\frac{B~{}l_{0}}{2}\sum_{i=1}^{i=N-2}\kappa_{i}^{2}$ (17)

The flexion force is then:

${\bf f}_{i}^{(b)}=\frac{B~{}l_{0}}{2}\frac{\partial}{\partial{\bf
r}_{i}}\bigl{[}\sum_{i^{\prime}=1}^{i^{\prime}=N-2}\kappa_{i^{\prime}}^{2}\bigr{]}$
(18)

First, we notice that for weakly bend fibers ${\bf l}_{i}\simeq{\bf l}_{i-1}$,
and for weakly extended fibers $l_{i}\simeq l_{0}$. Then, the denominator of
(16) is $\simeq 4~{}l_{0}^{6}$:

$\frac{\partial\kappa^{2}_{i^{\prime}}}{\partial{\bf
r}_{i}}\simeq\frac{1}{l_{0}^{6}}\frac{\partial}{\partial{\bf
r}_{i}}\bigl{[}{\bf l}_{i^{\prime}-1}^{2}{\bf l}_{i^{\prime}}^{2}-({\bf
l}_{i^{\prime}-1}\cdot{\bf l}_{i^{\prime}})^{2}\bigr{]}$ (19)

Using ${\bf l}_{i}={\bf r}_{i}-{\bf r}_{i-1}$, we obtain:

$\displaystyle\frac{\partial\kappa^{2}_{i-1}}{\partial{\bf r}_{i}}$
$\displaystyle\simeq\frac{1}{l_{0}^{4}}\bigl{[}2({\bf l}_{i-1}-{\bf
l}_{i-2})\bigr{]}$ (20a)
$\displaystyle\frac{\partial\kappa^{2}_{i}}{\partial{\bf r}_{i}}$
$\displaystyle\simeq\frac{1}{l_{0}^{4}}\bigl{[}-4({\bf l}_{i}-{\bf
l}_{i-1})\bigr{]}$ (20b)
$\displaystyle\frac{\partial\kappa^{2}_{i+1}}{\partial{\bf r}_{i}}$
$\displaystyle\simeq\frac{1}{l_{0}^{4}}\bigl{[}2({\bf l}_{i+1}-{\bf
l}_{i})\bigr{]}$ (20c)

and $(\partial\kappa^{2}_{j}/\partial{\bf r}_{i})=0$ if $|i-j|>1$. We obtain
finally:

$\displaystyle{\bf f}_{i}^{(b)}$ $\displaystyle=$
$\displaystyle\frac{B}{l_{0}^{3}}\bigl{[}-{\bf l}_{i-2}+3{\bf l}_{i-1}-3{\bf
l}_{i}+{\bf l}_{i+1}\bigr{]}$ (21) $\displaystyle=$
$\displaystyle\frac{B}{l_{0}^{3}}\bigl{[}{\bf r}_{i-2}-4{\bf r}_{i-1}+6{\bf
r}_{i}-4{\bf r}_{i+1}+{\bf r}_{i+2}\bigr{]}$ (22)

for $N-3\leq i\leq 2$. Expressions of the forces for $i<2$, and $i>N-3$ are
obtained by noticing that summation in (18) is for $i^{\prime}=1$ to
$i^{\prime}=N-2$.

$\displaystyle{\bf f}_{0}^{(b)}$
$\displaystyle=-\frac{B}{l_{0}^{3}}\bigl{[}{\bf r}_{0}-2{\bf r}_{1}+{\bf
r}_{2}\bigr{]}$ (23a) $\displaystyle{\bf f}_{1}^{(b)}$
$\displaystyle=-\frac{B}{l_{0}^{3}}\bigl{[}-2{\bf r}_{0}+5{\bf r}_{1}-4{\bf
r}_{2}+{\bf r}_{3}\bigr{]}$ (23b) $\displaystyle{\bf f}_{N-2}^{(b)}$
$\displaystyle=-\frac{B}{l_{0}^{3}}\bigl{[}{\bf r}_{N-4}-4{\bf r}_{N-3}+5{\bf
r}_{N-2}-2{\bf r}_{N-1}\bigr{]}$ (23c) $\displaystyle{\bf f}_{N-1}^{(b)}$
$\displaystyle=-\frac{B}{l_{0}^{3}}\bigl{[}{\bf r}_{N-3}-2{\bf r}_{N-2}+{\bf
r}_{N-1}\bigr{]}$ (23d)

###  VI.2 Twist moment and forces

The twist energy of the discrete rod is:

$E^{(t)}=\frac{C}{2~{}l_{0}}\sum_{i=1}^{i=N-2}(\theta_{i}-\theta_{i-1}+l_{0}~{}\tau_{s,i})^{2}$
(24)

where $\tau_{s,i}$ is the torsion of the center line at node $i$, and
$(\theta_{i}-\theta_{i-1})/l_{0}$ is the internal twist. The torsion
$\tau_{s}$ of the center line is obtained from Frenet-Serret equations as
$\tau_{s}=(d{\bf N}/ds)\cdot{\bf B}$, where $({\bf T},{\bf N},{\bf B})$ are
tangent, normal and bi-normal vector of the centre line of the fiber. They are
obtained by multiple differentiation of tangent vector ${\bf e}_{i}$, with
appropriate interpolations depending if derivatives are evaluated at nodes or
at cylinder.

${\bf
m}^{(t)}_{i-1;i}=\frac{C}{l_{0}}(\theta_{i}-\theta_{i-1}+l_{0}~{}\tau_{s,i}){\bf
e}_{i-1}$ (25)

Taking into account the torque acting from the segment $(i+1;i+2)$ on the
segment $(i;i+1)$, the total elastic twist torque acting on the segment
$(i;i+1)$ is:

$\begin{split}{\bf
m}^{(t)}_{i}=\frac{C}{l_{0}}\bigl{[}&(\theta_{i+1}-\theta_{i}+l_{0}~{}\tau_{s,i+1})~{}{\bf
e}_{i+1}\\\ -&(\theta_{i}-\theta_{i-1}-l_{0}~{}\tau_{s,i})~{}{\bf
e}_{i-1}\bigr{]}\end{split}$ (26)

This torque is split in two components. The axial (co-linear to ${\bf e}_{i}$)
component is:

$m^{(t)}_{i}={\bf m}^{(t)}_{i}\cdot{\bf e}_{i}$ (27)

whereas the remaining perpendicular component ${\bf
m}^{(t)}_{i}-m^{(t)}_{i}{\bf e}_{i}$ is written as a system of two points
forces ${\bf f}^{(t)}_{i}$ and ${\bf f}^{(t)}_{i+1}$ acting at points $i$ and
$i+1$ such that:

$\displaystyle{\bf f}^{(t)}_{i}+{\bf f}^{(t)}_{i+1}$ $\displaystyle=0$ (28a)
$\displaystyle({\bf r}_{i+1}-{\bf r}_{i})\times{\bf f}^{(t)}_{i+1}$
$\displaystyle={\bf m}^{(t)}_{i}-m^{(t)}_{i}{\bf e}_{i}$ (28b)
$\displaystyle{\bf f}^{(t)}_{i}\cdot{\bf e}_{i}$ $\displaystyle=0$ (28c)

(28a) ensures that the system of two points forces is a torque, (28b) assigns
the moment, and (28c) that those forces do not stretch the rod. Using $({\bf
r}_{i+1}-{\bf r}_{i})=l_{i}{\bf e}_{i}$, we finally obtain the two forces
acting on nodes:

${\bf f}^{(t)}_{i+1}=-{\bf f}^{(t)}_{i}=({\bf m}^{(t)}_{i}/l_{i})\times{\bf
e}_{i}$ (29)

###  VI.3 Contact forces distribution

Let’s a contact force ${\bf f}^{(c)}$ acting at point ${\bf r}_{C}$. We are
looking for two point forces ${\bf f}_{i}^{(c)}$ (respectively ${\bf
f}_{i+1}^{(c)}$) acting at point ${\bf r}_{i}$ (resp. ${\bf r}_{i+1}$) and a
moment $m_{i}{\bf e}_{i}$ such that:

$\displaystyle{\bf f}_{i}^{(c)}+{\bf f}_{i+1}^{(c)}$ $\displaystyle={\bf
f}^{(c)}$ (30a) $\displaystyle({\bf r}_{i+1}-{\bf r}_{i})\times{\bf
f}_{i+1}^{(c)}+m_{i}^{(c)}{\bf e}_{i}$ $\displaystyle=({\bf r}_{C}-{\bf
r}_{i})\times{\bf f}^{(c)}$ (30b)

Scalar product of (30b) with ${\bf e}_{i}$ gives:

$m_{i}^{(c)}=\bigl{[}({\bf r}_{C}-{\bf r}_{i})\times{\bf
f}^{(c)}\bigr{]}\cdot{\bf e}_{i}$ (31)

and cross product of (30b) with ${\bf e}_{i}$ gives:

$\displaystyle l_{i}{\bf f}_{i+1}^{(c)}-l_{i}\bigl{[}{\bf
f}_{i+1}^{(c)}\cdot{\bf e}_{i}\bigl{]}{\bf e}_{i}$ $\displaystyle={\bf
f}^{(c)}\bigl{[}{\bf e}_{i}\cdot({\bf r}_{C}-{\bf r}_{i})\bigr{]}$
$\displaystyle-({\bf r}_{C}-{\bf r}_{i})\bigl{[}{\bf e}_{i}\cdot{\bf
f}^{(c)}\bigl{]}$ (32)

Defining the parallel and perpendicular component of a force ${\bf f}$ with
respect to the cylinder axis as:

$\displaystyle{\bf f}^{\parallel}$ $\displaystyle=\bigr{[}{\bf f}\cdot{\bf
e}_{i}\bigr{]}{\bf e}_{i}$ (33a) $\displaystyle{\bf f}^{\perp}$
$\displaystyle={\bf f}-{\bf f}^{\parallel}$ (33b)

we obtain:

${\bf f}_{i+1}^{(c),\perp}=s_{i}{\bf
f}^{(c),\perp}-\frac{R}{l_{i}}f^{(c),\parallel}{\bf n}$ (34)

(34) determines only the components of ${\bf f}_{i+1}^{(c)}$ which are
perpendicular to the axis. The parallel component of ${\bf f}_{i+1}^{(c)}$ is
obtained in the following way. Consider the cylinder of length $l_{i}$ made of
an elastic material, and let’s $k$ the stiffness of the corresponding
compressing spring. This cylinder may be viewed as the reunion of one cylinder
of length $s_{i}l_{i}$ with stiffness $k/s_{i}$, and one cylinder of length
$(1-s_{i})l_{i}$ with stiffness $k/(1-s_{i})$. Let’s a force ${\bf
f}^{(c),\parallel}$ applied at the junction between cylinders. This force
moves the junction on a distance $\delta=\|{\bf
f}^{(c),\parallel}\|/[k/s_{i}+k/(1-s_{i})]$. This displacement deforms the
part of length $(1-s_{i})l_{i}$ and generates a force
$f_{i+1}^{(c),\parallel}=[k/s_{i}]\delta=s_{i}f^{(c),\parallel}$ on this
spring. Inserting this equation in (34), we finally obtain:

${\bf f}_{i+1}^{(c)}=s_{i}{\bf
f}^{(c)}-\frac{R}{l_{i}}\bigl{[}f^{(c)}\cdot{\bf e}_{i}\bigr{]}{\bf n}$ (35)

###  VI.4 Distance

Figure 11: Distance between two points located on two segments of line.

We consider two segments $1$ and $2$ whose axis are drawn on figure 11. On
each axis are located at abscissa $s=0$ a sphere of rayon $r$, and a segment
of cylinder of radius $r$ for $0\leq s\leq 1$. The distance between two points
at abscissa $s_{1}$ and $s_{2}$ is:

$d^{2}(s_{1},s_{2})=({\bf a}+s_{1}~{}{\bf l}_{1}+s_{2}~{}{\bf l}_{2})^{2}$
(36)

The distance is minimal for $s_{1}^{*}$ and $s_{2}^{*}$ which verify:

$\Bigl{(}\frac{\partial d^{2}(s_{1},s_{2})}{\partial
s_{1}}\Bigr{)}(s_{1}^{*},s_{2}^{*})=\Bigl{(}\frac{\partial
d^{2}(s_{1},s_{2})}{\partial s_{1}}\Bigr{)}(s_{1}^{*},s_{2}^{*})=0~{}$ (37)

Equation 37 is solved to obtain $(s_{1}^{*},s_{2}^{*})$, and the minimal
distance $d(s_{1}^{*},s_{2}^{*})$ is obtained. If $d(s_{1}^{*},s_{2}^{*})<2r$,
with $0\leq s_{1}^{*}\leq 1$ and $0\leq s_{2}^{*}\leq 1$, the contact is found
between the two cylinders.

It not, the contact is checked between the sphere located at $s_{1}=0$ and the
cylinder $2$. For this the minimal distance is obtained for $s_{2}^{*}$
verifying:

$\Bigl{(}\frac{\partial d^{2}(0,s_{2})}{\partial
s_{2}}\Bigr{)}(0,s_{2}^{*})=0~{}$ (38)

Equation 38 is solved to obtain $s_{2}^{*}$, and the minimal distance
$d(0,s_{2}^{*})$ is obtained. If $d(0,s_{2}^{*})<2r$, with $0\leq
s_{2}^{*}\leq 1$, the contact is found between the sphere (1) and the cylinder
(2).

The contact between sphere (2) and cylinder (1) is searched in a similar way.
If not, we check for a contact between the two spheres.

## References

  * [1] Brigitte Vigolo, Alain Pénicaud, Claude Coulon, Cédric Sauder, René Pailler, Catherine Journet, Patrick Bernier, and Philippe Poulin.  Macroscopic fibers and ribbons of oriented carbon nanotubes.  Science, 290(5495):1331–1334, 2000.
  * [2] Audrey Frenot and Ioannis S. Chronakis.  Polymer nanofibers assembled by electrospinning.  Current Opinion in Colloid & Interface Science, 8(1):64–75, 2003\.
  * [3] Staffan Toll.  Packing mechanics of fiber reinforcements.  Polymer Engineering & Science, 38(8):1337–1350, 1998.
  * [4] V. Negi and R. C. Picu.  Mechanical behavior of nonwoven non-crosslinked fibrous mats with adhesion and friction.  Soft Matter, 15:5951–5964, 2019.
  * [5] Ning Pan.  Exploring the significance of structural hierarchy in material systems—a review.  Applied Physics Reviews, 1(2):021302, 2014.
  * [6] Patrick B. Warren, Robin C. Ball, and Raymond E. Goldstein.  Why clothes don’t fall apart: Tension transmission in staple yarns.  Phys. Rev. Lett., 120:158001, Apr 2018.
  * [7] Antoine Seguin and Jérôme Crassous.  Twist-controlled force amplification and spinning tension transition in yarn.  Phys. Rev. Lett., 128:078002, Feb 2022.
  * [8] J. Bohr and K. Olsen.  The ancient art of laying rope.  EPL (Europhysics Letters), 93(6):60004, mar 2011.
  * [9] John WS Hearle, Percy Grosberg, and Stanley Backer.  Structural mechanics of fibers, yarns, and fabrics.  John Wiley & Sons Inc., 1969.
  * [10] Samuel Poincloux, Mokhtar Adda-Bedia, and Frédéric Lechenault.  Geometry and elasticity of a knitted fabric.  Phys. Rev. X, 8:021075, Jun 2018.
  * [11] Gautier Verhille, Sébastien Moulinet, Nicolas Vandenberghe, Mokhtar Adda-Bedia, and Patrice Le Gal.  Structure and mechanics of aegagropilae fiber network.  Proceedings of the National Academy of Sciences, 114(18):4607–4612, 2017.
  * [12] Andrade-Silva, Ignacio, Godefroy, Théo, Pouliquen, Olivier, and Marthelot, Joel.  Cohesion of bird nests.  EPJ Web Conf., 249:06014, 2021.
  * [13] N. Weiner, Y. Bhosale, M. Gazzola, and H. King.  Mechanics of randomly packed filaments—the “bird nest” as meta-material.  Journal of Applied Physics, 127(5):050902, 2020.
  * [14] Benjamin F. Bayman.  Theory of hitches.  American Journal of Physics, 45(2):185–190, 1977.
  * [15] M. K. Jawed, P. Dieleman, B. Audoly, and P. M. Reis.  Untangling the mechanics and topology in the frictional response of long overhand elastic knots.  Phys. Rev. Lett., 115:118302, Sep 2015.
  * [16] B. Audoly, N. Clauvelin, and S. Neukirch.  Elastic knots.  Phys. Rev. Lett., 99:164301, Oct 2007.
  * [17] N. Clauvelin, B. Audoly, and S. Neukirch.  Matched asymptotic expansions for twisted elastic knots: a self-contact problem with non-trivial contact topology.  Journal of the Mechanics and Physics of Solids, 57:1623—1656, 2009\.
  * [18] Paul Grandgeorge, Changyeob Baek, Harmeet Singh, Paul Johanns, Tomohiko G. Sano, Alastair Flynn, John H. Maddocks, and Pedro M. Reis.  Mechanics of two filaments in tight orthogonal contact.  Proceedings of the National Academy of Sciences, 118(15):e2021684118, 2021.
  * [19] Paul Johanns, Paul Grandgeorge, Changyeob Baek, Tomohiko G. Sano, John H. Maddocks, and Pedro M. Reis.  The shapes of physical trefoil knots.  Extreme Mechanics Letters, 43:101172, 2021.
  * [20] Changyeob Baek, Paul Johanns, Tomohiko G. Sano, Paul Grandgeorge, and Pedro M. Reis.  Finite Element Modeling of Tight Elastic Knots.  Journal of Applied Mechanics, 88(2), 11 2020.
  * [21] Henna Tangri, Yu Guo, and Jennifer S. Curtis.  Packing of cylindrical particles: Dem simulations and experimental measurements.  Powder Technology, 317:72–82, 2017.
  * [22] Paul Langston, Andrew Kennedy, and Hannah Constantin.  Discrete element modelling of flexible fibre packing.  Computational Materials Science, 96:108–116, 01 2015.
  * [23] Denis Dumont, Maurine Houze, Paul Rambach, Thomas Salez, Sylvain Patinet, and Pascal Damman.  Emergent strain stiffening in interlocked granular chains.  Phys. Rev. Lett., 120:088001, Feb 2018.
  * [24] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, and E. Grinspun.  Discrete elastic rods.  ACM Transactions on Graphics, 27(3):63:1–63:12, 2008.
  * [25] Damien Durville.  Simulation of the mechanical behaviour of woven fabrics at the scale of fibers.  International Journal of Material Forming, 3(2):1241–1251, Sep 2010\.
  * [26] Damien Durville.  Contact-friction modeling within elastic beam assemblies: an application to knot tightening.  Computational Mechanics, 49(6):687–707, Jun 2012.
  * [27] Florence Bertails-Descoubes, Florent Cadoux, Gilles Daviet, and Vincent Acary.  A Nonsmooth Newton Solver for Capturing Exact Coulomb Friction in Fiber Assemblies.  ACM Transactions on Graphics, 30(1):Article No. 6, January 2011\.
  * [28] Jean Jacques Moreau.  Indetermination due to dry friction in multibody dynamics.  In European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS 2004 proceedings, Jyväskylä, Finland, July 2004.
  * [29] Andrew Choi, Dezhong Tong, Mohammad K. Jawed, and Jungseock Joo.  Implicit Contact Model for Discrete Elastic Rods in Knot Tying.  Journal of Applied Mechanics, 88(5), 03 2021.  051010\.
  * [30] Mickaël Ly, Jean Jouve, Laurence Boissieux, and Florence Bertails-Descoubes.  Projective Dynamics with Dry Frictional Contact.  ACM Transactions on Graphics, 39(4):Article 57:1–8, 2020.
  * [31] Bruno Chareyre and Pascal Villard.  Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems.  Journal of Engineering Mechanics, 131(7):689–698, 2005.
  * [32] Franck Bourrier, François Kneib, Bruno Chareyre, and Thierry Fourcaud.  Discrete modeling of granular soils reinforcement by plant roots.  Ecological Engineering, 61:646–657, 2013.  Soil Bio- and Eco-Engineering: The Use of Vegetation to Improve Slope Stability.
  * [33] Anna Effeindzourou, Bruno Chareyre, Klaus Thoeni, Anna Giacomini, and François Kneib.  Modelling of deformable structures in the general framework of the discrete element method.  Geotextiles and Geomembranes, 44(2):143–156, 2016.
  * [34] D. Kunhappan, B. Harthong, B. Chareyre, G. Balarac, and P. J. J. Dumont.  Numerical modeling of high aspect ratio flexible fibers in inertial flows.  Physics of Fluids, 29(9):093302, 2017.
  * [35] A.E.H.Love.  A Treatise on the Mathematical Theory of Elasticity.  Cambridge University Press, third edition edition, 1920.
  * [36] Joel Langer and David A. Singer.  Lagrangian aspects of the kirchhoff elastic rod.  SIAM Review, 38(4):605–618, 1996.
  * [37] G. H. M. van der Heijden and J. M. T. Thompson.  Helical and localised buckling in twisted rods: A unified analysis of the symmetric case.  Nonlinear Dynamics, 21(1):71–99, Jan 2000.
  * [38] P. A. Cundall and O. D. L. Strack.  A discrete numerical model for granular assemblies.  Géotechnique, 29(1):47–65, 1979.
  * [39] Stephen Timoshenko and James N. Goodier.  Theory of Elasticity.  McGraw-Hill, third edition, 1970.
  * [40] Daan Frenkel and Berend Smit.  Chapter 4 - molecular dynamics simulations.  In Daan Frenkel and Berend Smit, editors, Understanding Molecular Simulation (Second Edition), pages 63–107. Academic Press, San Diego, second edition edition, 2002.
  * [41] B. Muth, M.-K. Müller, P. Eberhard, and Stefan Luding.  Collision detection and administration methods for many particles with different sizes.  In P. Cleary, editor, Discrete Element Methods, DEM 07, pages 1–18. Minerals Engineering Int., August 2007.  null ; Conference date: 27-08-2007 Through 29-08-2007.
  * [42] K. E. Bisshopp and D. C. Drucker.  Large deflection of cantilever beams.  Quart. Appl. Math., 3:272–275, 1945.
  * [43] James M. Gere Stephen P. Timoshenko.  Theory of Elastic Stability Elasticity.  Dover Publication, 2nd ed. edition edition, 2009.
  * [44] Jae Ho Jung, Ning Pan, and Tae Jin Kang.  Capstan equation including bending rigidity and non-linear frictional behavior.  Mechanism and Machine Theory, 43(6):661–675, 2008.