Jérôme Crassous

A Discrete Element Method model for frictional fibers

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 [START_REF] Vigolo | Macroscopic fibers and ribbons of oriented carbon nanotubes[END_REF] or polymer fibers [START_REF] Frenot | Polymer nanofibers assembled by electrospinning[END_REF] 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 [START_REF] Toll | Packing mechanics of fiber reinforcements[END_REF] such as cushions or non-woven fabrics [START_REF] Negi | Mechanical behavior of nonwoven non-crosslinked fibrous mats with adhesion and friction[END_REF]. Textile fibers can be twisted to produce yarns [START_REF] Ning Pan | Exploring the significance of structural hierarchy in material systems-a review[END_REF][START_REF] Warren | Why clothes don't fall apart: Tension transmission in staple yarns[END_REF][START_REF] Seguin | Twist-controlled force amplification and spinning tension transition in yarn[END_REF], which are then assembled into cords [START_REF] Bohr | The ancient art of laying rope[END_REF], woven [START_REF] John Ws Hearle | Structural mechanics of fibers, yarns, and fabrics[END_REF] or knitted fabrics [START_REF] John Ws Hearle | Structural mechanics of fibers, yarns, and fabrics[END_REF][START_REF] Poincloux | Geometry and elasticity of a knitted fabric[END_REF]. Cyclic mechanical stresses can form very compact natural structures [START_REF] Verhille | Structure and mechanics of aegagropilae fiber network[END_REF], and birds also assemble fibers to build their nests [START_REF] Andrade-Silva | Cohesion of bird nests[END_REF][START_REF] Weiner | Mechanics of randomly packed filaments-the "bird nest" as meta-material[END_REF]. The contacts between fibers play a fundamental role in describing the physics of knots, which is a subtle competition between tension and friction [START_REF] Benjamin | Theory of hitches[END_REF][START_REF] Jawed | Untangling the mechanics and topology in the frictional response of long overhand elastic knots[END_REF], as well as eventual bending of the fibers [START_REF] Audoly | Elastic knots[END_REF][START_REF] Clauvelin | Matched asymptotic expansions for twisted elastic knots: a self-contact problem with non-trivial contact topology[END_REF][START_REF] Grandgeorge | Mechanics of two filaments in tight orthogonal contact[END_REF][START_REF] Johanns | The shapes of physical trefoil knots[END_REF].

Several approaches have been proposed to numerically simulate these structures. One approach is to use finite element algorithms to discretize the fibers [START_REF] Baek | Finite Element Modeling of Tight Elastic Knots[END_REF]. This approach allows a complete solution of the elasticity equations in complex geometries such as nodes [START_REF] Grandgeorge | Mechanics of two filaments in tight orthogonal contact[END_REF], but is only possible for systems with small numbers of contacts. Another approach is to model the fibers as connected spheres [START_REF] Tangri | Packing of cylindrical particles: Dem simulations and experimental measurements[END_REF] or sphero-cylinders [START_REF] Langston | Discrete element modelling of flexible fibre pack-ing[END_REF] 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 [START_REF] Dumont | Emergent strain stiffening in interlocked granular chains[END_REF].

More realistic approaches are the simulations of fibers as discrete [START_REF] Bergou | Discrete elastic rods[END_REF] or continuous [START_REF] Durville | Simulation of the mechanical behaviour of woven fabrics at the scale of fibers[END_REF][START_REF] Durville | Contact-friction modeling within elastic beam assemblies: an application to knot tightening[END_REF] cylindrical elastic chains of circular cross-sections. The noninterpenetration 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 [START_REF] Durville | Simulation of the mechanical behaviour of woven fabrics at the scale of fibers[END_REF][START_REF] Durville | Contact-friction modeling within elastic beam assemblies: an application to knot tightening[END_REF][START_REF] Bertails-Descoubes | A Nonsmooth Newton Solver for Capturing Exact Coulomb Friction in Fiber Assemblies[END_REF]. In those algorithms, the fibers are moved in order * jerome.crassous@univ-rennes1.fr to find the positions of the surfaces that match the nonpenetration 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 [START_REF] Jacques | Indetermination due to dry friction in multibody dynamics[END_REF] simulations, and the solutions selected by iterative algorithms are not well controlled [START_REF] Jacques | Indetermination due to dry friction in multibody dynamics[END_REF], 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 [START_REF] Choi | Implicit Contact Model for Discrete Elastic Rods in Knot Tying[END_REF], or if qualitative simulations are needed as in computer graphic community [START_REF] Ly | Laurence Boissieux, and Florence Bertails-Descoubes[END_REF]. 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 [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF][START_REF] Bourrier | Discrete modeling of granular soils reinforcement by plant roots[END_REF], with geotextiles [START_REF] Effeindzourou | Modelling of de-formable structures in the general framework of the discrete element method[END_REF], and for the behavior of suspensions of frictional fibers in viscous flow [START_REF] Kunhappan | Numerical modeling of high aspect ratio flexible fibers in inertial flows[END_REF].

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. [START_REF] Bergou | Discrete elastic rods[END_REF]. 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. [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF]. 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

A. Description of the fiber Following [START_REF] Bergou | Discrete elastic rods[END_REF], we model a fiber as an ensemble of N connected points (see figure 1(a)). Let r i , with 0 ≤ i ≤ N -1 be the position of the point, and e i = (r i+1r i )/ r i+1 -r i , with 0 ≤ i ≤ N -2 the unit vector joining two successive points. We note l i = r i+1 -r i . The segment joining two successive points is the generatrix of a cylinder of circular basis of diameter d. In addition, each point 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 r i , and N -1 cylinder rotations θ 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.

B. 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 λ. The equilibrium length of the spring is l 0 , and the elongation force exerted by point i + 1 on the mass located at r i is:

f (e) i+1;i = k 0 (l i -l 0 ) + λ li e i (1)
Each point i is submitted to forces from points i -1 and i + 1 so that f

(e) i = f (e) i+1;i + f (e)
i-1;i , 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) = s (B/2) κ 2 ds with B the bending stiffness of the fiber, and κ the curvature. The bending energy of the discrete fiber is:

E (b) = B l 0 2 i=N -2 i=1 κ 2 i ( 2
)
where κ i is the curvature at node i, and the summation is extended to all nodes except ending ones. Writing the curvatures κ i as function of nodes positions r i , the flexion force f

(b) i = -(∂E (b) /∂r i
) acting on nodes i is (see Appendix A):

f (b) i = - B l 3 0 r i-2 -4r i-1 + 6r i -4r i+1 + r i+2 (3)
for (N -3) ≥ i ≥ 2. Expressions of the forces f calculation supposes that the fibers are weakly extended and bent (see Appendix A).

The internal elastic torque is obtained from the twisting energy [START_REF] Bergou | Discrete elastic rods[END_REF]: E (t) = s (C/2) τ 2 ds with C the torsion modulus of the fiber, and τ the twist of the fiber. The twist may be written as [START_REF] Love | A Treatise on the Mathematical Theory of Elasticity[END_REF][START_REF] Langer | Lagrangian aspects of the kirchhoff elastic rod[END_REF][START_REF] Van Der Heijden | Helical and localised buckling in twisted rods: A unified analysis of the symmetric case[END_REF]: τ = τ int + τ s . The internal twist τ int is the twist of the fiber if simply unbent, whereas τ s is the torsion of the fibre centre line. Writing internal twist at node i as (θ i -θ i-1 )/l 0 , and τ s,i the torsion of the center line at node i, we obtain the twist energy of the discrete fiber as:

E (t) = C l 0 2 i=N -2 i=1 θ i -θ i-1 + l 0 τ s,i 2 (4)
The twist moment acting on segment joining nodes i and i + 1 is obtained by differentiating (4) with respect to θ i (see VI B):

m (t) i = C l 0 (θ i+1 -θ i + l 0 τ s,i+1 ) e i+1 -(θ i -θ i-1 -l 0 τ s,i ) e i-1 (5)
This elastic moment is split into one component m

(t) i = m (t)
i • e i of the moment along the axis of the segment (i, i + 1), and into forces acting on nodes (see VI B).

C. Contact forces

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 r C (t) is located on the segment with ending points r (1)

Hi and r

(2) Hi 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 C. Let d(t) = r Hi + (r 1 -δ/2) n the contact point. We use the Cundall-Strack model for the contact force [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The normal contact force exerted by cylinder 1 on cylinder 2 is modeled as a spring-dashpot system:

f (c) n = -k n δ + λ n δ n (6)
with k n the contact stiffness between the cylinders, and λ n the contact damping. This contact law is a simplified version of the elastic contact force between 2 cylinders [START_REF] Timoshenko | Theory of Elasticity[END_REF] which varies non-linearly with the interpenetration f (c) n ∼ δ 3/2 , and which depends on the angle between cylinder axis. The tangential contact force is a Coulomb-Force:

f (c) t = -M in k t u t ; µk n δ u t u t ( 7
)
where k t is the tangential stiffness, u t the tangential displacement, and µ the microscopic friction coefficient. The tangential displacement 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, u t is first rotated to take into account the rotation of the normal vector. Lets dβ the angle between the normal at times t -dt and t: n(t -dt) and n(t), and k = n(t -dt) × n(t) /|n(t -dt) × n(t)| the axis rotation. We name u rot t (t -dt) this rotated tangential displacement.

Then, the displacement is integrated as:

u t (t) = u rot t (t -dt) + (v (2) -v (1) )dt (8)
where v (1) (and similar for v (2) ) is the velocity of the point of the cylinder (1) coinciding with the contact point C. The velocity is v

(1) = ṙ(1) i +Ω (1) i ×(r C -r (1)
i ) with Ω

(1) i the rotational velocity of the segment i of fiber 1. The rotation vector is separated into an axial and non-axial components as:

Ω (1) = Ω (1)
⊥ + θi e i , where the non-axial component is

Ω (1) i,⊥ = 1 li e i × (ṙ (1) i+1 - ṙ(1) i ). Finally, the normal component (u t •n) n is removed. If k t u t > µk n δ, then the tangential displacement is renor- malized such that u t = µk n δ/k t .
The contact force f (c) is then expressed as a system of forces f (c) i and f (c) i+1 applied on nodes i and i + 1, and a moment m i e i acting the cylinder connecting those nodes. The conservation of the resultant and of the moment of contact force implies that:

f (c) i + f (c) i+1 = f (c) (9a) (r i+1 -r i ) × f i+1 + m i e i = (r C -r i ) × f (c) (9b)
A possible choice for forces and moment is (see VI C):

m (c) i = (r C -r i ) × f (c) ) • e i (10a) f (c) i = (1 -s i )f (c) + R l i (f (c) • e i )n (10b) f (c) i+1 = s i f (c) - R l i (f (c) • e i )n (10c)
It should be noticed that [START_REF] John Ws Hearle | Structural mechanics of fibers, yarns, and fabrics[END_REF] does not set all the components of f

(c) i and f (c)
i+1 , and that a supplementary conditions expressed in appendix VI C must be added to obtain [START_REF] Poincloux | Geometry and elasticity of a knitted fabric[END_REF].

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.

D. Miscellaneous forces.

In addition misc extra forces may be added. A global viscous damping force f (v) i = -λ v ṙ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 f (g) i = m 0 g with g the gravity field may be also added. Other external forces may applied to fibers such pretension at ends of fibers.

III. NUMERICAL IMPLEMENTATION A. Integration of equation of motions

The dynamical equations of motions writes as:

M 0 ri = f (e) i + f (b) i + f (c) i + f (v) i + f (g) i (11a) J z θi = m (t) i + m (c) i (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 r i . The dynamical equations are integrated using a standard second-order Verlet algorithm [START_REF] Frenkel | Chapter 4 -molecular dynamics simulations[END_REF].

B. Physical parameters for simulations

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 ν) of density ρ, we have k 0 = Eπr 2 /l 0 , m 0 = ρπr 2 l 0 , and then t 0 = l 0 (ρ/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πr 2 is the force that extend a hypothetical perfectly elastic fiber by 100%.

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 * 1. In practice the simulations are done with f * ∼ 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 ρ l under a tension f is

v t = (f /ρ l ) 1/2 . With ρ l = m 0 /l 0 , we have ρ * l = 1, and the non-dimensional speed of trans- verse wave is v * t = (f * /ρ * l ) 1/2 = (f * ) 1/2
when no bending stiffness are present.

The non-dimensional bending stiffness is B * = B/k 0 l 3 0 . For an elastic fiber as consider in III B 1, we have B = Eπr 4 /4, and then B * = (r * ) 2 /4. Similarly, the nondimensional torsional modulus is C * = C/k 0 l 3 0 . For an elastic fiber of radius r, we have C = Eπr 4 /2(1 + ν), and then

C * = (r * ) 2 /2(1 + ν).
The longitudinal damping λ is chosen to avoid compression waves that travel continuously through the fibers, needing very long time to return to equilibrium. We take λ ∼ (k 0 m 0 ) 1/2 , and then λ * ∼ 1 for this.

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 ef f r 1/2 δ 3/2 , with E ef f = E/(1-ν 2 )
, ν being the Poisson ratio of the material. For doing the linearization, we arbitrary set that the elastic energy of the Hertzian contact ∼ E ef f r 1/2 δ 5/2 is equal elastic energy k n δ 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 = (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 * ∼ 10 -5 -10 -3 , and typical radius are r * ∼ 10 -1 , we have k * n ∼ 1. For sake of simplicity, the tangential stiffness is taken as k * t = k * n . Some damping of the normal force λ n may be introduced. We took λ * n ∼ 1 for rapid relaxation of oscillating motion of contact.

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 ∼ (m 0 /k 0 ) 1/2 = t 0 , whereas the time scale for a contact to establish is

∼ (m 0 /k n ) 1/2 = t 0 (k 0 /k n ) 1/2 . The time step is cho- sen as dt = M in t 0 ; t 0 (k 0 /k n ) 1/2 /10, leading to: dt * = 1 10 M in 1; (k * n ) -1/2 (13)
such that both relaxat ions occur on at least 10 time steps. In practice, since k * n ∼ 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.

C. 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.

Distance between fibers

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 D. 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.

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 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 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.

Neighbor search method

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 [START_REF] Muth | Collision detection and administration methods for many particles with different sizes[END_REF]. 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 ∼ 2 l 0 . For segments of section ∼ 4 r 2 , the number of segments in each cell is ∼ 2 (l 0 /r) 2 for dense 3D system. Since l 0 /r 1, sorting particles in cell of size ∼ 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 ∼ 4 r. In this case for a system of N f fibers of N segments each, the total numbers of fictitious spheres is ∼ 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 -j| ≤ d(i, j) -2r. So the next segment where we need to search contact verifies j > 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 λ * = 2.8, contact stiffness k * n = k * t = 1, contact damping λ * n = 1, global viscous damping λ * v = 0.001, inertia momentum J * = m * r * 2 /2 (homogeneous cylinder).

A. Elastic rods without contacts

The elastic rod model has been already tested in misc situations that do not involve frictional contacts [START_REF] Bergou | Discrete elastic rods[END_REF].

The test examples presented here are just for checking the approximations used in II B and VI.

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 small deflections δ L, F L 2 /B 3δ/L. At large deflections δ ∼ L simulations agree well with the analytical solution of Bisshopp et al. [START_REF] Bisshopp | Large deflection of cantilever beams[END_REF]. Simulations with beams made with different N show that the solution obtained with the discrete beam converges towards analytical solution as ∼ N -2 (see fig. 5.c). It may be noticed that since the maximum force applied in those simulations are of order F ∼ 10 B/L 2 , the maximum non-dimensional force F * ∼ 10 -4 << 1, so that the beam stretching is negligible.

The second example is the buckling of a rod submitted to a compression and applied torque at its ends (see fig.

B. Static without flexion : capstan

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 µ = 0.2. We first apply an

C. Static with flexion : elastic knots

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. [START_REF] Audoly | Elastic knots[END_REF][START_REF] Clauvelin | Matched asymptotic expansions for twisted elastic knots: a self-contact problem with non-trivial contact topology[END_REF]. 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 µ 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 [START_REF] Choi | Implicit Contact Model for Discrete Elastic Rods in Knot Tying[END_REF].

We simulate numerically such knot by considering a flexible spring (r * = 0.1, B * = (r * ) 2 /4 = 2.5 10 -3 , N = 500, λ * v = 4.10 -4 ) as shown on figure 8(a). We first knot the fiber by setting µ = 0 and applying a tension ± T * e z at ends. After this preparation stage, we set µ 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 =< dt i /ds -1 >, where dt i /ds = e i+1 -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 [START_REF] Audoly | Elastic knots[END_REF][START_REF] Clauvelin | Matched asymptotic expansions for twisted elastic knots: a self-contact problem with non-trivial contact topology[END_REF], we introduce ε = r/R. The figure 8(b) shows the tension T * as a function of ε for frictionless (µ = 0), and frictional (µ = 0.1) loosening and opening knots. The analytical solutions in the limits ε 1 and µ 1 are [START_REF] Audoly | Elastic knots[END_REF][START_REF] Clauvelin | Matched asymptotic expansions for twisted elastic knots: a self-contact problem with non-trivial contact topology[END_REF]:

T r 2 B = ε 4 2 ± µσε 3 (14)
where the sign ± depends if the knot is tightened (+) or loosened (-), and σ is a numerical constant which is σ 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 ∼ ε 4 when ε 0.15. Two possible sources of deviations may be identified. First, the equation ( 14) is obtained in the limit ε 1, and deviations may arise from high order ε terms in equation [START_REF] Benjamin | Theory of hitches[END_REF]. Second, for 0.15, we have R * = R/l 0 = r * /ε 2 ∼ 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 µ = 0.1 loosening, where some steps in ε 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 µ 1 which is used to obtain [START_REF] Benjamin | Theory of hitches[END_REF].

D. Impact: falling chain.

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 f ps. The figure 9(a) show some snapshots of the impact. The chain is simu-lated 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 = N m 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 = N g/l 0 t -2 0 = N g * . We take g * = 4.9 10 -5 so that N g * 5. 10 -3 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, µ = 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. In a recent study, Seguin et al. considered the situation of a staple yarn made of twisted totally flexible fibers [START_REF] Seguin | Twist-controlled force amplification and spinning tension transition in yarn[END_REF]. 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 ≤ i ≤ 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 * e z is applied to both ends of the assembly of fibers. For this, each fiber i with 1 ≤ i ≤ N f is submitted at both ends j = 0 and j = N to an external shear force:

τ i (j) = ± C * i r 2 i (j) e z × r i (j) ( 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 µ = 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 [START_REF] Seguin | Twist-controlled force amplification and spinning tension transition in yarn[END_REF] if the twist is not homogenized along the yarn. The duration of this preparation stage is t * = 5.10 5 , and the total twist θ 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 (∆f /∆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 [START_REF] Seguin | Twist-controlled force amplification and spinning tension transition in yarn[END_REF], show that the force necessary to separate the two slivers is ln(1 + f ) 0.75 µθ 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.

ACKNOWLEDGMENTS

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 A. Flexion forces

The curvature κ i at a node N -2 ≥ i ≥ 1 is first expressed as a function of the positions of nodes r i-1 , r i and r i+1 . The radius R i = 1/κ 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 (r i-1 , r i , r i+1 ) using the Heron formula. After elementary calculus, we obtain:

κ 2 i = 4 l 2 i-1 l 2 i -(l i-1 • l i ) 2 l 2 i-1 l 2 i (l i-1 + l i ) 2 (16)
where we noted l i = r i -r i-1 . The bending energy is

E (b) = B l 0 2 i=N -2 i=1 κ 2 i ( 17
)
The flexion force is then:

f (b) i = B l 0 2 ∂ ∂r i i =N -2 i =1 κ 2 i (18)
First, we notice that for weakly bend fibers l i l i-1 , and for weakly extended fibers l i l 0 . Then, the denominator of ( 16) is 4 l 6 0 :

∂κ 2 i ∂r i 1 l 6 0 ∂ ∂r i l 2 i -1 l 2 i -(l i -1 • l i ) 2 (19)
Using l i = r i -r i-1 , we obtain:

∂κ 2 i-1 ∂r i 1 l 4 0 2(l i-1 -l i-2 ) (20a) ∂κ 2 i ∂r i 1 l 4 0 -4(l i -l i-1 ) (20b) ∂κ 2 i+1 ∂r i 1 l 4 0 2(l i+1 -l i ) (20c)
and (∂κ 2 j /∂r i ) = 0 if |i -j| > 1. We obtain finally:

f (b) i = B l 3 0 -l i-2 + 3l i-1 -3l i + l i+1 (21)
= B l 3 0 r i-2 -4r i-1 + 6r i -4r i+1 + r i+2 (22)
for N -3 ≤ i ≤ 2. Expressions of the forces for i < 2, and i > N -3 are obtained by noticing that summation in ( 18) is for i

= 1 to i = N -2. f (b) 0 = - B l 3 0 r 0 -2r 1 + r 2 (23a) f (b) 1 = - B l 3 0 -2r 0 + 5r 1 -4r 2 + r 3 (23b) f (b) N -2 = - B l 3 0 r N -4 -4r N -3 + 5r N -2 -2r N -1 (23c) f (b) N -1 = - B l 3 0 r N -3 -2r N -2 + r N -1 (23d)

B. Twist moment and forces

The twist energy of the discrete rod is:

E (t) = C 2 l 0 i=N -2 i=1 (θ i -θ i-1 + l 0 τ s,i ) 2 (24)
where τ s,i is the torsion of the center line at node i, and (θ i -θ i-1 )/l 0 is the internal twist. The torsion τ s of the center line is obtained from Frenet-Serret equations as τ s = (dN/ds) • B, where (T, N, B) are tangent, normal and bi-normal vector of the centre line of the fiber. They are obtained by multiple differentiation of tangent vector e i , with appropriate interpolations depending if derivatives are evaluated at nodes or at cylinder.

m (t) i-1;i = C l 0 (θ i -θ i-1 + l 0 τ s,i )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:

m (t) i = C l 0 (θ i+1 -θ i + l 0 τ s,i+1 ) e i+1 -(θ i -θ i-1 -l 0 τ s,i ) e i-1 (26)
This torque is split in two components. The axial (colinear to e i ) component is:

m (t) i = m (t) i • e i ( 27
)
whereas the remaining perpendicular component m

(t) i - m (t)
i e i is written as a system of two points forces f (t) i and f (t) i+1 acting at points i and i + 1 such that:

f (t) i + f (t) i+1 = 0 (28a) (r i+1 -r i ) × f (t) i+1 = m (t) i -m (t) i e i (28b) f (t)
i • e i = 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 (r i+1 -r i ) = l i e i , we finally obtain the two forces acting on nodes:

f (t) i+1 = -f (t) i = (m (t) i /l i ) × e i (29)
C. Contact forces distribution

Let's a contact force f (c) acting at point r C . We are looking for two point forces f (c) i (respectively f (c) i+1 ) acting at point r i (resp. r i+1 ) and a moment m i e i such that:

f (c) i + f (c) i+1 = f (c) (30a) (r i+1 -r i ) × f (c) i+1 + m (c) i e i = (r C -r i ) × f (c) (30b)
Scalar product of (30b) with e i gives:

m (c) i = (r C -r i ) × f (c) • e i (31)
and cross product of (30b) with e i gives:

l i f (c) i+1 -l i f (c) i+1 • e i e i = f (c) e i • (r C -r i ) -(r C -r i ) e i • f (c) (32)
Defining the parallel and perpendicular component of a force f with respect to the cylinder axis as:

f = f • e i e i (33a)
f ⊥ = f -f (33b)
we obtain: 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 f (c), applied at the junction between cylinders. This force moves the junction on a distance δ = f (c), /[k/s i + k/(1 -s i )]. This displacement deforms the part of length (1 -s i )l i and generates a force f (c), i+1 = [k/s i ]δ = s i f (c), on this spring. Inserting this equation in [START_REF] Kunhappan | Numerical modeling of high aspect ratio flexible fibers in inertial flows[END_REF], we finally obtain:

f (c),⊥ i+1 = s i f (c),⊥ - R l i f (c), n (34)
f (c) i+1 = s i f (c) - R l i f (c) • e i n (35)

D. Distance

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 ≤ s ≤ 1. The distance between two points at abscissa s 1 and s 2 is: The distance is minimal for s * 1 and s * 2 which verify: 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: 2 ≤ 1, the contact is found between the sphere (1) and the cylinder [START_REF] Frenot | Polymer nanofibers assembled by electrospinning[END_REF].

d 2 (s 1 , s 2 ) = (a + s 1 l 1 + s 2 l 2 ) 2 (36)
∂d 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.

FIG. 1 .

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

  (b) i for i < 2 and i > (N -3) are given in Appendix A. The

FIG. 2 .

 2 FIG. 2. Contact between two cylinders.

( 2 )

 2 Hi -r(1) Hi this minimal distance, and n(t) the normal unitary vector. We note δ(t) = (d 1 + d 2 )/2 -d(t) the interpenetration of the two cylinders, and r C (t) = r[START_REF] Vigolo | Macroscopic fibers and ribbons of oriented carbon nanotubes[END_REF]

FIG. 3 .

 3 FIG. 3. (a) 2D view of two sets composed of one sphere and one cylinder. (b) Motions of an external Cext and internal Cint contact points at the junction between two cylinders.

FIG. 4 .

 4 FIG. 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.

FIG. 5 .

 5 FIG. 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 F L 2 /B = 3δ/L. (c) Convergence: δN is the end deflection for a beam with N nodes ((N -1) cylinders), and δ th the theoretical deflection. Applied forces is F L 2 /B = 10. Symbols are relative error, and dashed line is a N -2 decay.

  .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

FIG. 6 .

 6 FIG. 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 (b) t = 2πB/L.

  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 (b) t as a function of the rod length. Stability analysis of twisted rods leads to[43]: M (b) t = 2πB/L. As shown on fig.6.b, the numerical results are in correct agreement with this theoretical law.

FIG. 7 .

 7 FIG. 7. Tension in a rolled string around a cylinder: T * is the tension in the string, and θ is the rolling angle.Circles are symbol, plain line is an exponential fit. See for simulations parameters Inset: Schematic of the experiment.

FIG. 8 .

 8 FIG. 8. (a) Snapshot of a (31) knot. For seek of clarity, illustration is made with r * = 0.2. (b) Tension as a function of ε = r/R for frictionless and frictional strings. Symbols are numerical data, and lines are theoretical expressions given by equation (14).

FIG. 9 .

 9 FIG. 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.

  ε

E

  . Multiple fibers: a yarn model.

FIG. 10 .

 10 FIG. 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: µ = 1, blue triangles: µ = 0.5, green squares: µ = 0.2. Lines are for guidelines. (e) Same data as (d) plotted as a function of H = µθ 2 R/L. Dotted line is 0.75 µθ 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.

( 34 )

 34 determines only the components of f (c) i+1 which are perpendicular to the axis. The parallel component of f (c) i+1

FIG. 11 .

 11 FIG. 11. Distance between two points located on two segments of line.

  solved to obtain s * 2 , and the minimal distance d(0, s * 2 ) is obtained. If d(0, s * 2 ) < 2r, with 0 ≤ s *

  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 ≤ s * 1 ≤ 1 and 0 ≤ s * 2 ≤ 1, the contact is found between the two cylinders.

s 1 , s 2 ) ∂s 1 (s * 1 , s * 2 ) = ∂d 2 (s 1 , s 2 ) ∂s 1 (s * 1 , s * 2 ) = 0 (37)